dataanalysis2d.hh

来自「open lattice boltzmann project www.open」· HH 代码 · 共 686 行 · 第 1/2 页

/*  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 2D BlockStructures -- generic implementation.
 */ 

#ifndef DATA_ANALYSIS_2D_HH
#define DATA_ANALYSIS_2D_HH

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

namespace olb {

/////// Struct AnalysisFieldsImpl2D  ///////////////////////

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


/////// Struct AnalysisFields2D  ///////////////////////

template<typename T, template<typename U> class Lattice>
AnalysisFields2D<T,Lattice>::AnalysisFields2D (
        TensorField2D<T,2>& velField_,
        TensorField2D<T,2>& momentumField_,
        ScalarField2D<T>&   pressureField_,
        ScalarField2D<T>&   velNormField_,
        ScalarField2D<T>&   vortField_,
        TensorField2D<T,3>& strainRateField_,
        TensorField2D<T,3>& stressField_,
        ScalarField2D<T>&   divRhoUField_,
        ScalarField2D<T>&   poissonField_,
        TensorField2D<T,Lattice<T>::q >& populationField_ )
    :
      velField(velField_),
      momentumField(momentumField_),
      pressureField(pressureField_),
      velNormField(velNormField_),
      vortField(vortField_),
      strainRateField(strainRateField_),
      stressField(stressField_),
      divRhoUField(divRhoUField_),
      poissonField(poissonField_),
      populationField(populationField_)
{ }

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

/////// Class DataAnalysis2D  /////////////////////////////

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

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

template<typename T, template<typename U> class Lattice>
DataAnalysis2D<T,Lattice>::DataAnalysis2D(DataAnalysis2D<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>
DataAnalysis2D<T,Lattice>::~DataAnalysis2D() { }

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

template<typename T, template<typename U> class Lattice>
TensorFieldBase2D<T,2> const&
    DataAnalysis2D<T,Lattice>::getVelocity() const
{
    computeVelocityField();
    return fields.velField;
}

template<typename T, template<typename U> class Lattice>
TensorFieldBase2D<T,2> const&
    DataAnalysis2D<T,Lattice>::getMomentum() const
{
    computeMomentumField();
    return fields.momentumField;
}

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

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

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

template<typename T, template<typename U> class Lattice>
TensorFieldBase2D<T,3> const&
    DataAnalysis2D<T,Lattice>::getStrainRate() const
{
    computeStrainRateField();
    return fields.strainRateField;
}

template<typename T, template<typename U> class Lattice>
TensorFieldBase2D<T,3> const&
    DataAnalysis2D<T,Lattice>::getStrainRateFromStress() const
{
    computeStrainRateFieldFromStress();
    return fields.stressField;
}

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

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

template<typename T, template<typename U> class Lattice>
TensorFieldBase2D<T, Lattice<T>::q > const&
    DataAnalysis2D<T,Lattice>::getPopulations() const
{
    computePopulationField();
    return fields.populationField;
}

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

    T enstrophy = T();
    for (int iX=0; iX<nx; ++iX) {
        for (int iY=0; iY<ny; ++iY) {
           T dyux = (
                 fields.velField.get(iX  ,iY+1)[0]
               + fields.velField.get(iX+1,iY+1)[0]
               - fields.velField.get(iX  ,iY  )[0]
               - fields.velField.get(iX+1,iY  )[0] ) / (T)2;
           T dxuy = (
                 fields.velField.get(iX+1,iY  )[1]
               + fields.velField.get(iX+1,iY+1)[1]
               - fields.velField.get(iX  ,iY  )[1]
               - fields.velField.get(iX  ,iY+1)[1] ) / (T)2;
           T omega = dxuy - dyux;
           enstrophy += omega*omega;
        }
    }
    enstrophy /= (2*nx*ny);
    return enstrophy;
}

template<typename T, template<typename U> class Lattice>
T DataAnalysis2D<T,Lattice>::computeMeanEnstrophy2() const {
    computeVorticityField();
    int nx = fields.vortField.getNx();
    int ny = fields.vortField.getNy();

    T enstrophy = T();
    for (int iX=1; iX<nx-1; ++iX) {
        for (int iY=1; iY<ny-1; ++iY) {
           enstrophy += util::sqr(fields.vortField.get(iX,iY));
        }
    }

    for (int iX=1; iX<nx-1; ++iX) {
        enstrophy += 0.5* (
                util::sqr(fields.vortField.get(iX,0)) +
                util::sqr(fields.vortField.get(iX,ny-1)) );
    }
    for (int iY=1; iY<ny-1; ++iY) {
        enstrophy += 0.5* (
                util::sqr(fields.vortField.get(0,iY)) +
                util::sqr(fields.vortField.get(nx-1,iY)) );
    }
    enstrophy += 0.25 * util::sqr(fields.vortField.get(0,0));
    enstrophy += 0.25 * util::sqr(fields.vortField.get(0,ny-1));
    enstrophy += 0.25 * util::sqr(fields.vortField.get(nx-1,0));
    enstrophy += 0.25 * util::sqr(fields.vortField.get(nx-1,ny-1));

    enstrophy /= 2*(nx-1)*(ny-1);

    return enstrophy;
}

template<typename T, template<typename U> class Lattice>
void DataAnalysis2D<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) {
            block.get(iX,iY).computeU(fields.velField.get(iX,iY));
        }
    }
    flags.velFieldComputed = true;
}

template<typename T, template<typename U> class Lattice>
void DataAnalysis2D<T,Lattice>::computeMomentumField() const {
    if (flags.momentumFieldComputed) return;
    fields.momentumField.construct();
    for (int iX=0; iX<fields.momentumField.getNx(); ++iX) {
        for (int iY=0; iY<fields.momentumField.getNy(); ++iY) {
            T rho;
            block.get(iX,iY).computeRhoU(rho, fields.momentumField.get(iX,iY));
            for (int iD=0; iD<Lattice<T>::d; ++iD) {
                fields.momentumField.get(iX,iY)[iD] *= rho;
            }
        }
    }
    flags.momentumFieldComputed = true;
}

template<typename T, template<typename U> class Lattice>
void DataAnalysis2D<T,Lattice>::computePressureField() const {
    if (flags.pressureFieldComputed) return;
    fields.pressureField.construct();
    for (int iX=0; iX<fields.pressureField.getNx(); ++iX) {
        for (int iY=0; iY<fields.pressureField.getNy(); ++iY) {
            fields.pressureField.get(iX,iY) = block.get(iX,iY).computeRho();
            fields.pressureField.get(iX,iY) -= (T)1;
            fields.pressureField.get(iX,iY) /= Lattice<T>::invCs2;
        }
    }
    flags.pressureFieldComputed = true;
}

template<typename T, template<typename U> class Lattice>
void DataAnalysis2D<T,Lattice>::computeVelocityNormField() const {
    if (flags.velNormFieldComputed) return;
    fields.velNormField.construct();
    computeVelocityField();
    for (int iEl=0; iEl<fields.velNormField.getNx()*fields.velNormField.getNy(); ++iEl) {
        fields.velNormField[iEl] = sqrt(util::normSqr<T,2>(fields.velField[iEl]));
    }
    flags.velNormFieldComputed = true;
}

template<typename T, template<typename U> class Lattice>
T DataAnalysis2D<T,Lattice>::bulkVorticity(int iX, int iY) const {
    OLB_PRECONDITION(flags.velFieldComputed);
    OLB_PRECONDITION(iX>=1 && iX<=fields.vortField.getNx()-2);
    OLB_PRECONDITION(iY>=1 && iY<=fields.vortField.getNy()-2);

    T dyux = fd::centralGradient(fields.velField.get(iX,iY+1)[0],
                                 fields.velField.get(iX,iY-1)[0]);
    T dxuy = fd::centralGradient(fields.velField.get(iX+1,iY)[1],
                                 fields.velField.get(iX-1,iY)[1]);
    return dxuy - dyux;
}

template<typename T, template<typename U> class Lattice>
T DataAnalysis2D<T,Lattice>::bulkXderiv (