dataanalysis2d.hh

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

        int iX, int iY, int iD, TensorField2D<T,2> const& field) const
{
    OLB_PRECONDITION(iX>=1 && iX<=field.getNx()-2);
    OLB_PRECONDITION(iY>=1 && iY<=field.getNy()-2);

    T dxu = fd::centralGradient(field.get(iX+1,iY)[iD],
                                field.get(iX-1,iY)[iD]);
    return dxu;
}

template<typename T, template<typename U> class Lattice>
T DataAnalysis2D<T,Lattice>::bulkYderiv (
        int iX, int iY, int iD, TensorField2D<T,2> const& field) const
{
    OLB_PRECONDITION(iX>=1 && iX<=field.getNx()-2);
    OLB_PRECONDITION(iY>=1 && iY<=field.getNy()-2);

    T dyu = fd::centralGradient(field.get(iX,iY+1)[iD],
                                field.get(iX,iY-1)[iD]);
    return dyu;
}

template<typename T, template<typename U> class Lattice>
T DataAnalysis2D<T,Lattice>::bulkDeriv (
    int iX, int iY, int iAlpha, int iBeta,
    TensorField2D<T,2> const& field) const
{
    switch(iAlpha) {
        case 0:
            return bulkXderiv(iX,iY,iBeta,field);
        case 1:
            return bulkYderiv(iX,iY,iBeta,field);
        default:
            OLB_ASSERT( false, "iAlpha>1!");
            return T();
    }
}

template<typename T, template<typename U> class Lattice>
T DataAnalysis2D<T,Lattice>::bulkStrain (
    int iX, int iY, int iAlpha, int iBeta) const
{
    OLB_PRECONDITION( flags.momentumFieldComputed );
    return ( bulkDeriv(iX,iY,iAlpha,iBeta, fields.momentumField) +
             bulkDeriv(iX,iY,iBeta,iAlpha, fields.momentumField) ) / (T)2;

}

template<typename T, template<typename U> class Lattice>
T DataAnalysis2D<T,Lattice>::bulkDivRhoU(int iX, int iY) const {
    OLB_PRECONDITION( flags.momentumFieldComputed );
    return bulkDeriv(iX,iY, 0,0, fields.momentumField) +
           bulkDeriv(iX,iY, 1,1, fields.momentumField);
}

template<typename T, template<typename U> class Lattice>
T DataAnalysis2D<T,Lattice>::bulkPoisson(int iX, int iY) const {
    OLB_PRECONDITION( flags.velFieldComputed );
    T dxux = bulkDeriv(iX,iY,0,0, fields.velField);
    T dxuy = bulkDeriv(iX,iY,0,1, fields.velField);
    T dyux = bulkDeriv(iX,iY,1,0, fields.velField);
    T dyuy = bulkDeriv(iX,iY,1,1, fields.velField);
    return dxux*dxux + (T)2*dxuy*dyux + dyuy*dyuy;
}


template<typename T, template<typename U> class Lattice>
T DataAnalysis2D<T,Lattice>::boundaryXderiv (
        int iX, int iY, int iD, TensorField2D<T,2> const& field ) const
{
    OLB_PRECONDITION(iX>=0 && iX<=field.getNx()-1);
    OLB_PRECONDITION(iY>=0 && iY<=field.getNy()-1);

    T dxu;
    
    if (iX==0) {
        dxu = fd::boundaryGradient(field.get(iX,iY)[iD],
                                   field.get(iX+1,iY)[iD],
                                   field.get(iX+2,iY)[iD]);
    }
    else if (iX==fields.vortField.getNx()-1) {
        dxu = -fd::boundaryGradient(field.get(iX,iY)[iD],
                                    field.get(iX-1,iY)[iD],
                                    field.get(iX-2,iY)[iD]);
    }
    else {
        dxu = fd::centralGradient(field.get(iX+1,iY)[iD],
                                  field.get(iX-1,iY)[iD]);
    }

    return dxu;
}

template<typename T, template<typename U> class Lattice>
T DataAnalysis2D<T,Lattice>::boundaryYderiv (
        int iX, int iY, int iD, TensorField2D<T,2> const& field ) const
{
    OLB_PRECONDITION(iX>=0 && iX<=field.getNx()-1);
    OLB_PRECONDITION(iY>=0 && iY<=field.getNy()-1);

    T dyu;

    if (iY==0) {
        dyu = fd::boundaryGradient(field.get(iX,iY)[iD],
                                   field.get(iX,iY+1)[iD],
                                   field.get(iX,iY+2)[iD]);
    }
    else if (iY==fields.vortField.getNy()-1) {
        dyu = -fd::boundaryGradient(field.get(iX,iY)[iD],
                                    field.get(iX,iY-1)[iD],
                                    field.get(iX,iY-2)[iD]);
    }
    else {
        dyu = fd::centralGradient(field.get(iX,iY+1)[iD],
                                  field.get(iX,iY-1)[iD]);
    }
   
    return dyu;
}

template<typename T, template<typename U> class Lattice>
T DataAnalysis2D<T,Lattice>::boundaryDeriv (
    int iX, int iY, int iAlpha, int iBeta,
    TensorField2D<T,2> const& field ) const
{
    switch(iAlpha) {
        case 0:
            return boundaryXderiv(iX,iY,iBeta, field);
        case 1:
            return boundaryYderiv(iX,iY,iBeta, field);
        default:
            OLB_ASSERT( false, "iAlpha>1!");
            return T();
    }
}

template<typename T, template<typename U> class Lattice>
T DataAnalysis2D<T,Lattice>::boundaryStrain (
    int iX, int iY, int iAlpha, int iBeta) const
{
    OLB_PRECONDITION( flags.momentumFieldComputed );
    return ( boundaryDeriv(iX,iY,iAlpha,iBeta, fields.momentumField) +
             boundaryDeriv(iX,iY,iBeta,iAlpha, fields.momentumField) ) / (T)2;

}

template<typename T, template<typename U> class Lattice>
T DataAnalysis2D<T,Lattice>::boundaryDivRhoU(int iX, int iY) const {
    OLB_PRECONDITION( flags.momentumFieldComputed );
    return boundaryDeriv(iX,iY,0,0, fields.momentumField) +
           boundaryDeriv(iX,iY,1,1, fields.momentumField);
}

template<typename T, template<typename U> class Lattice>
T DataAnalysis2D<T,Lattice>::boundaryPoisson(int iX, int iY) const {
    OLB_PRECONDITION( flags.velFieldComputed );
    T dxux = boundaryDeriv(iX,iY,0,0, fields.velField);
    T dxuy = boundaryDeriv(iX,iY,0,1, fields.velField);
    T dyux = boundaryDeriv(iX,iY,1,0, fields.velField);
    T dyuy = boundaryDeriv(iX,iY,1,1, fields.velField);
    return dxux*dxux + (T)2*dxuy*dyux + dyuy*dyuy;
}

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

    T dyux = boundaryYderiv(iX,iY, 0, fields.velField);
    T dxuy = boundaryXderiv(iX,iY, 1, fields.velField);

    return dxuy - dyux;
}

template<typename T, template<typename U> class Lattice>
void DataAnalysis2D<T,Lattice>::computeVorticityField() const {
    if (flags.vortFieldComputed) return;
    fields.vortField.construct();
    computeVelocityField();

    int nx = fields.vortField.getNx();
    int ny = fields.vortField.getNy();

    for (int iX=1; iX<nx-1; ++iX) {
        for (int iY=1; iY<ny-1; ++iY) {
            fields.vortField.get(iX,iY) = bulkVorticity(iX,iY);
        }
    }

    for (int iX=1; iX<nx-1; ++iX) {
        fields.vortField.get(iX,0) = boundaryVorticity(iX,0);
        fields.vortField.get(iX,ny-1) = boundaryVorticity(iX,ny-1);
    }

    for (int iY=0; iY<ny; ++iY) {
        fields.vortField.get(0,iY) = boundaryVorticity(0,iY);
        fields.vortField.get(nx-1,iY) = boundaryVorticity(nx-1,iY);
    }

    flags.vortFieldComputed = true;
}

template<typename T, template<typename U> class Lattice>
void DataAnalysis2D<T,Lattice>::computeStrainRateField() const {
    if (flags.strainRateFieldComputed) return;
    fields.strainRateField.construct();
    computeMomentumField();

    int nx = fields.vortField.getNx();
    int ny = fields.vortField.getNy();

    int iPi = 0;
    for (int iAlpha=0; iAlpha<2; ++iAlpha) {
        for (int iBeta=iAlpha; iBeta<2; ++iBeta) {

            for (int iX=1; iX<nx-1; ++iX) {
                for (int iY=1; iY<ny-1; ++iY) {
                    fields.strainRateField.get(iX,iY)[iPi] =
                        bulkStrain(iX,iY, iAlpha,iBeta);
                }
            }

            for (int iX=1; iX<nx-1; ++iX) {
                fields.strainRateField.get(iX,0)[iPi] =
                    boundaryStrain(iX,0, iAlpha,iBeta);
                fields.strainRateField.get(iX,ny-1)[iPi] =
                    boundaryStrain(iX,ny-1, iAlpha,iBeta);
            }

            for (int iY=0; iY<ny; ++iY) {
                fields.strainRateField.get(0,iY)[iPi] =
                    boundaryStrain(0,iY, iAlpha,iBeta);
                fields.strainRateField.get(nx-1,iY)[iPi] =
                    boundaryStrain(nx-1,iY, iAlpha,iBeta);
            }

            ++iPi;
        }
    }

    flags.strainRateFieldComputed = true;
}

template<typename T, template<typename U> class Lattice>
void DataAnalysis2D<T,Lattice>::computeStrainRateFieldFromStress() const {
    if (flags.stressFieldComputed) return;
    fields.stressField.construct();

    for (int iX=0; iX<fields.velField.getNx(); ++iX) {
        for (int iY=0; iY<fields.velField.getNy(); ++iY) {
            block.get(iX,iY).computeStress(fields.stressField.get(iX,iY));
            T omega = block.get(iX,iY).getDynamics()->getOmega();
            for (int iPi=0; iPi<3; ++iPi) {
                fields.stressField.get(iX,iY)[iPi] *=
                    -omega / (T)2 * Lattice<T>::invCs2;
            }
        }
    }

    flags.stressFieldComputed = true;
}

template<typename T, template<typename U> class Lattice>
void DataAnalysis2D<T,Lattice>::computeDivRhoUField() const {
    if (flags.divRhoUFieldComputed) return;
    fields.divRhoUField.construct();
    computeMomentumField();

    int nx = fields.divRhoUField.getNx();
    int ny = fields.divRhoUField.getNy();

    for (int iX=1; iX<nx-1; ++iX) {
        for (int iY=1; iY<ny-1; ++iY) {
            fields.divRhoUField.get(iX,iY) = bulkDivRhoU(iX,iY);
        }
    }

    for (int iX=1; iX<nx-1; ++iX) {
        fields.divRhoUField.get(iX,0) = boundaryDivRhoU(iX,0);
        fields.divRhoUField.get(iX,ny-1) = boundaryDivRhoU(iX,ny-1);
    }

    for (int iY=0; iY<ny; ++iY) {
        fields.divRhoUField.get(0,iY) = boundaryDivRhoU(0,iY);
        fields.divRhoUField.get(nx-1,iY) = boundaryDivRhoU(nx-1,iY);
    }

    flags.divRhoUFieldComputed = true;
}


template<typename T, template<typename U> class Lattice>
void DataAnalysis2D<T,Lattice>::computePoissonTerm() const {
    if (flags.poissonFieldComputed) return;
    fields.poissonField.construct();
    computeVelocityField();

    int nx = fields.poissonField.getNx();
    int ny = fields.poissonField.getNy();

    for (int iX=1; iX<nx-1; ++iX) {
        for (int iY=1; iY<ny-1; ++iY) {
            fields.poissonField.get(iX,iY) = bulkPoisson(iX,iY);
        }
    }

    for (int iX=1; iX<nx-1; ++iX) {
        fields.poissonField.get(iX,0) = boundaryPoisson(iX,0);
        fields.poissonField.get(iX,ny-1) = boundaryPoisson(iX,ny-1);
    }

    for (int iY=0; iY<ny; ++iY) {
        fields.poissonField.get(0,iY) = boundaryPoisson(0,iY);
        fields.poissonField.get(nx-1,iY) = boundaryPoisson(nx-1,iY);
    }

    flags.poissonFieldComputed = true;
}

template<typename T, template<typename U> class Lattice>
void DataAnalysis2D<T,Lattice>::computePopulationField() const {
    if (flags.populationFieldComputed) return;
    fields.populationField.construct();

    int nx = fields.populationField.getNx();
    int ny = fields.populationField.getNy();
    for (int iX=0; iX<nx; ++iX) {
        for (int iY=0; iY<ny; ++iY) {
            for (int iPop=0; iPop<Lattice<T>::q; ++iPop) {
                fields.populationField.get(iX,iY)[iPop] = block.get(iX,iY)[iPop];
            }
        }
    }

    flags.populationFieldComputed = true;
}

}  // namespace olb


#endif