dataanalysis3d.hh

open lattice boltzmann project www.openlb.org

/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2006, 2007 Jonas Latt and Bernd Stahl
 *  Address: Rue General Dufour 24,  1211 Geneva 4, Switzerland 
 *  E-mail: [e-mail]
 *
 *  This program is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU General Public License
 *  as published by the Free Software Foundation; either version 2
 *  of the License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public 
 *  License along with this program; if not, write to the Free 
 *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *  Boston, MA  02110-1301, USA.
*/

/** \file
 * Data analysis (formerly known as BlockStatistics) on 3D BlockStructures -- generic implementation.
 */ 

#ifndef DATA_ANALYSIS_3D_HH
#define DATA_ANALYSIS_3D_HH

#include <cmath>
#include "dataAnalysis3D.h"
#include "cell.h"
#include "finiteDifference.h"
#include "util.h"

namespace olb {

/////// Struct AnalysisFieldsImpl3D  ///////////////////////

template<typename T, template<typename U> class Lattice>
AnalysisFieldsImpl3D<T,Lattice>::AnalysisFieldsImpl3D(int nx, int ny, int nz)
    : velField(nx, ny, nz),
      momentumField(nx, ny, nz),
      pressureField(nx, ny, nz),
      velNormField(nx, ny, nz),
      vortField(nx, ny, nz),
      vortNormField(nx, ny, nz),
      strainRateField(nx, ny, nz),
      stressField(nx, ny, nz),
      divRhoUField(nx, ny, nz),
      poissonField(nx, ny, nz),
      populationField(nx, ny, nz)
{ }


/////// Struct AnalysisFields3D  ///////////////////////

template<typename T, template<typename U> class Lattice>
AnalysisFields3D<T,Lattice>::AnalysisFields3D (
        TensorField3D<T,3>& velField_,
        TensorField3D<T,3>& momentumField_,
        ScalarField3D<T>&   pressureField_,
        ScalarField3D<T>&   velNormField_,
        TensorField3D<T,3>& vortField_,
        ScalarField3D<T>&   vortNormField_,
        TensorField3D<T,6>& strainRateField_,
        TensorField3D<T,6>& stressField_,
        ScalarField3D<T>&   divRhoUField_,
        ScalarField3D<T>&   poissonField_,
        TensorField3D<T,Lattice<T>::q >& populationField_ )
    :
      velField(velField_),
      momentumField(momentumField_),
      pressureField(pressureField_),
      velNormField(velNormField_),
      vortField(vortField_),
      vortNormField(vortNormField_),
      strainRateField(strainRateField_),
      stressField(stressField_),
      divRhoUField(divRhoUField_),
      poissonField(poissonField_),
      populationField(populationField_)
{ }

template<typename T, template<typename U> class Lattice>
AnalysisFields3D<T,Lattice>::AnalysisFields3D(AnalysisFieldsImpl3D<T,Lattice>& impl)
    :
      velField(impl.velField),
      momentumField(impl.momentumField),
      pressureField(impl.pressureField),
      velNormField(impl.velNormField),
      vortField(impl.vortField),
      vortNormField(impl.vortNormField),
      strainRateField(impl.strainRateField),
      stressField(impl.stressField),
      divRhoUField(impl.divRhoUField),
      poissonField(impl.poissonField),
      populationField(impl.populationField)
{ }


/////// Class DataAnalysis3D  /////////////////////////////

template<typename T, template<typename U> class Lattice>
DataAnalysis3D<T,Lattice>::DataAnalysis3D (
    BlockStructure3D<T,Lattice> const& block_ )
    : block(block_),
      pointsToDefaultFields(true),
      defaultFields(block.getNx(), block.getNy(), block.getNz()),
      fields(defaultFields)
{
    flags.reset();
}

template<typename T, template<typename U> class Lattice>
DataAnalysis3D<T,Lattice>::DataAnalysis3D(BlockStructure3D<T,Lattice> const& block_,
                                          AnalysisFields3D<T,Lattice>& fields_ )
    : block(block_),
      pointsToDefaultFields(false),
      defaultFields(block.getNx(), block.getNy(), block.getNz()),
      fields(fields_ )
{
    flags.reset();
}

template<typename T, template<typename U> class Lattice>
DataAnalysis3D<T,Lattice>::DataAnalysis3D(DataAnalysis3D<T,Lattice> const& rhs)
    : block(rhs.block),
      pointsToDefaultFields(rhs.pointsToDefaultFields),
      defaultFields(rhs.defaultFields),
      fields( pointsToDefaultFields ? defaultFields : rhs.fields )
{ }


template<typename T, template<typename U> class Lattice>
DataAnalysis3D<T,Lattice>::~DataAnalysis3D() { }

template<typename T, template<typename U> class Lattice>
void DataAnalysis3D<T,Lattice>::reset() const {
    flags.reset();
}

template<typename T, template<typename U> class Lattice>
TensorFieldBase3D<T,3> const& DataAnalysis3D<T,Lattice>::getVelocity() const
{
    computeVelocityField();
    return fields.velField;
}

template<typename T, template<typename U> class Lattice>
TensorFieldBase3D<T,3> const&
    DataAnalysis3D<T,Lattice>::getMomentum() const
{
    computeMomentumField();
    return fields.momentumField;
}

template<typename T, template<typename U> class Lattice>
ScalarFieldBase3D<T> const&
    DataAnalysis3D<T,Lattice>::getPressure() const
{
    computePressureField();
    return fields.pressureField;
}

template<typename T, template<typename U> class Lattice>
TensorFieldBase3D<T,3> const&
    DataAnalysis3D<T,Lattice>::getVorticity() const
{
    computeVorticityField();
    return fields.vortField;
}

template<typename T, template<typename U> class Lattice>
ScalarFieldBase3D<T> const&
    DataAnalysis3D<T,Lattice>::getVelocityNorm() const
{
    computeVelocityNormField();
    return fields.velNormField;
}

template<typename T, template<typename U> class Lattice>
ScalarFieldBase3D<T> const&
    DataAnalysis3D<T,Lattice>::getVorticityNorm() const
{
    computeVorticityNormField();
    return fields.vortNormField;
}

template<typename T, template<typename U> class Lattice>
TensorFieldBase3D<T,6> const&
    DataAnalysis3D<T,Lattice>::getStrainRate() const
{
    computeStrainRateField();
    return fields.strainRateField;
}

template<typename T, template<typename U> class Lattice>
TensorFieldBase3D<T,6> const&
    DataAnalysis3D<T,Lattice>::getStrainRateFromStress() const
{
    computeStrainRateFieldFromStress();
    return fields.stressField;
}

template<typename T, template<typename U> class Lattice>
ScalarFieldBase3D<T> const&
    DataAnalysis3D<T,Lattice>::getDivRhoU() const
{
    computeDivRhoUField();
    return fields.divRhoUField;
}

template<typename T, template<typename U> class Lattice>
ScalarFieldBase3D<T> const&
    DataAnalysis3D<T,Lattice>::getPoissonTerm() const
{
    computePoissonTerm();
    return fields.poissonField;
}

template<typename T, template<typename U> class Lattice>
TensorFieldBase3D<T, Lattice<T>::q > const&
    DataAnalysis3D<T,Lattice>::getPopulations() const
{
    computePopulations();
    return fields.populationField;
}


template<typename T, template<typename U> class Lattice>
T DataAnalysis3D<T,Lattice>::computeMeanEnstrophy() const {
    computeVelocityField();
    int nx = fields.velField.getNx()-1;
    int ny = fields.velField.getNy()-1;
    int nz = fields.velField.getNz()-1;

    T enstrophy = T();
    for (int iX=0; iX<nx; ++iX) {
        for (int iY=0; iY<ny; ++iY) {
            for (int iZ=0; iZ<nz; ++iZ) {
                T dxuy = (
                      fields.velField.get(iX+1,iY+1,iZ  )[1]
                    + fields.velField.get(iX+1,iY  ,iZ  )[1]
                    + fields.velField.get(iX+1,iY+1,iZ+1)[1]
                    + fields.velField.get(iX+1,iY  ,iZ+1)[1]
                    - fields.velField.get(iX  ,iY+1,iZ  )[1]
                    - fields.velField.get(iX  ,iY  ,iZ  )[1]
                    - fields.velField.get(iX  ,iY+1,iZ+1)[1]
                    - fields.velField.get(iX  ,iY  ,iZ+1)[1] ) / (T)4;
                T dxuz = (
                      fields.velField.get(iX+1,iY+1,iZ  )[2]
                    + fields.velField.get(iX+1,iY  ,iZ  )[2]
                    + fields.velField.get(iX+1,iY+1,iZ+1)[2]
                    + fields.velField.get(iX+1,iY  ,iZ+1)[2]
                    - fields.velField.get(iX  ,iY+1,iZ  )[2]
                    - fields.velField.get(iX  ,iY  ,iZ  )[2]
                    - fields.velField.get(iX  ,iY+1,iZ+1)[2]
                    - fields.velField.get(iX  ,iY  ,iZ+1)[2] ) / (T)4;
                T dyux = (
                      fields.velField.get(iX  ,iY+1,iZ  )[0]
                    + fields.velField.get(iX+1,iY+1,iZ  )[0]
                    + fields.velField.get(iX  ,iY+1,iZ+1)[0]
                    + fields.velField.get(iX+1,iY+1,iZ+1)[0]
                    - fields.velField.get(iX  ,iY  ,iZ  )[0]
                    - fields.velField.get(iX+1,iY  ,iZ  )[0]
                    - fields.velField.get(iX  ,iY  ,iZ+1)[0]
                    - fields.velField.get(iX+1,iY  ,iZ+1)[0] ) / (T)4;
                T dyuz = (
                      fields.velField.get(iX  ,iY+1,iZ  )[2]
                    + fields.velField.get(iX+1,iY+1,iZ  )[2]
                    + fields.velField.get(iX  ,iY+1,iZ+1)[2]
                    + fields.velField.get(iX+1,iY+1,iZ+1)[2]
                    - fields.velField.get(iX  ,iY  ,iZ  )[2]
                    - fields.velField.get(iX+1,iY  ,iZ  )[2]
                    - fields.velField.get(iX  ,iY  ,iZ+1)[2]
                    - fields.velField.get(iX+1,iY  ,iZ+1)[2] ) / (T)4;
                T dzux = (
                      fields.velField.get(iX  ,iY  ,iZ+1)[0]
                    + fields.velField.get(iX+1,iY+1,iZ+1)[0]
                    + fields.velField.get(iX  ,iY+1,iZ+1)[0]
                    + fields.velField.get(iX+1,iY  ,iZ+1)[0]
                    - fields.velField.get(iX  ,iY  ,iZ  )[0]
                    - fields.velField.get(iX+1,iY  ,iZ  )[0]
                    - fields.velField.get(iX  ,iY+1,iZ  )[0]
                    - fields.velField.get(iX+1,iY+1,iZ  )[0] ) / (T)4;
                T dzuy = (
                      fields.velField.get(iX  ,iY  ,iZ+1)[1]
                    + fields.velField.get(iX+1,iY+1,iZ+1)[1]
                    + fields.velField.get(iX  ,iY+1,iZ+1)[1]
                    + fields.velField.get(iX+1,iY  ,iZ+1)[1]
                    - fields.velField.get(iX  ,iY  ,iZ  )[1]
                    - fields.velField.get(iX+1,iY  ,iZ  )[1]
                    - fields.velField.get(iX  ,iY+1,iZ  )[1]
                    - fields.velField.get(iX+1,iY+1,iZ  )[1] ) / (T)4;

                T omegaX = dyuz - dzuy;
                T omegaY = dzux - dxuz;
                T omegaZ = dxuy - dyux;

                enstrophy += omegaX*omegaX+omegaY*omegaY+omegaZ*omegaZ;
            }
        }
    }
    enstrophy /= (2*nx*ny*nz);
    return enstrophy;
}

template<typename T, template<typename U> class Lattice>
void DataAnalysis3D<T,Lattice>::computeVelocityField() const {
    if (flags.velFieldComputed) return;
    fields.velField.construct();
    for (int iX=0; iX<fields.velField.getNx(); ++iX) {
        for (int iY=0; iY<fields.velField.getNy(); ++iY) {
            for (int iZ=0; iZ<fields.velField.getNz(); ++iZ) {
                block.get(iX,iY,iZ).computeU(fields.velField.get(iX,iY,iZ));
            }