dataanalysis3d.hh

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

    OLB_PRECONDITION( flags.momentumFieldComputed );
    return boundaryDeriv(iX,iY,iZ, 0,0, fields.momentumField) +
           boundaryDeriv(iX,iY,iZ, 1,1, fields.momentumField) +
           boundaryDeriv(iX,iY,iZ, 2,2, fields.momentumField);
}

template<typename T, template<typename U> class Lattice>
T DataAnalysis3D<T,Lattice>::boundaryPoisson (
        int iX, int iY, int iZ) const
{
    OLB_PRECONDITION( flags.velFieldComputed );

    T dxux = boundaryDeriv(iX,iY,iZ, 0,0, fields.velField);
    T dxuy = boundaryDeriv(iX,iY,iZ, 0,1, fields.velField);
    T dxuz = boundaryDeriv(iX,iY,iZ, 0,2, fields.velField);
    T dyux = boundaryDeriv(iX,iY,iZ, 1,0, fields.velField);
    T dyuy = boundaryDeriv(iX,iY,iZ, 1,1, fields.velField);
    T dyuz = boundaryDeriv(iX,iY,iZ, 1,2, fields.velField);
    T dzux = boundaryDeriv(iX,iY,iZ, 2,0, fields.velField);
    T dzuy = boundaryDeriv(iX,iY,iZ, 2,1, fields.velField);
    T dzuz = boundaryDeriv(iX,iY,iZ, 2,2, fields.velField);

    return dxux*dxux + dyuy*dyuy + dzuz*dzuz +
           (T)2*( dxuy*dyux + dxuz*dzux + dyuz*dzuy );
}

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

    T dyuz = boundaryYderiv(iX,iY,iZ, 2, fields.velField);
    T dzuy = boundaryZderiv(iX,iY,iZ, 1, fields.velField);

    return dyuz - dzuy;
}

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

    T dzux = boundaryZderiv(iX,iY,iZ, 0, fields.velField);
    T dxuz = boundaryXderiv(iX,iY,iZ, 2, fields.velField);

    return dzux - dxuz;
}

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

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

    return dxuy - dyux;
}

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

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

    for (int iX=1; iX<nx-1; ++iX) {
        for (int iY=1; iY<ny-1; ++iY) {
            for (int iZ=1; iZ<nz-1; ++iZ) {
                fields.vortField.get(iX,iY,iZ)[0] = bulkVorticityX(iX,iY,iZ);
                fields.vortField.get(iX,iY,iZ)[1] = bulkVorticityY(iX,iY,iZ);
                fields.vortField.get(iX,iY,iZ)[2] = bulkVorticityZ(iX,iY,iZ);
            }
        }
    }

    for (int iX=0; iX<nx-1; ++iX) {
        for (int iY=0; iY<ny-1; ++iY) {
            fields.vortField.get(iX,iY,0)[0] = boundaryVorticityX(iX,iY,0);
            fields.vortField.get(iX,iY,0)[1] = boundaryVorticityY(iX,iY,0);
            fields.vortField.get(iX,iY,0)[2] = boundaryVorticityZ(iX,iY,0);
            fields.vortField.get(iX,iY,nz-1)[0] = boundaryVorticityX(iX,iY,nz-1);
            fields.vortField.get(iX,iY,nz-1)[1] = boundaryVorticityY(iX,iY,nz-1);
            fields.vortField.get(iX,iY,nz-1)[2] = boundaryVorticityZ(iX,iY,nz-1);
        }
    }

    for (int iX=0; iX<nx-1; ++iX) {
        for (int iZ=0; iZ<nz-1; ++iZ) {
            fields.vortField.get(iX,0,iZ)[0] = boundaryVorticityX(iX,0,iZ);
            fields.vortField.get(iX,0,iZ)[1] = boundaryVorticityY(iX,0,iZ);
            fields.vortField.get(iX,0,iZ)[2] = boundaryVorticityZ(iX,0,iZ);
            fields.vortField.get(iX,ny-1,iZ)[0] = boundaryVorticityX(iX,ny-1,iZ);
            fields.vortField.get(iX,ny-1,iZ)[1] = boundaryVorticityY(iX,ny-1,iZ);
            fields.vortField.get(iX,ny-1,iZ)[2] = boundaryVorticityZ(iX,ny-1,iZ);
        }
    }

    for (int iY=0; iY<ny-1; ++iY) {
        for (int iZ=0; iZ<nz-1; ++iZ) {
            fields.vortField.get(0,iY,iZ)[0] = boundaryVorticityX(0,iY,iZ);
            fields.vortField.get(0,iY,iZ)[1] = boundaryVorticityY(0,iY,iZ);
            fields.vortField.get(0,iY,iZ)[2] = boundaryVorticityZ(0,iY,iZ);
            fields.vortField.get(nx-1,iY,iZ)[0] = boundaryVorticityX(nx-1,iY,iZ);
            fields.vortField.get(nx-1,iY,iZ)[1] = boundaryVorticityY(nx-1,iY,iZ);
            fields.vortField.get(nx-1,iY,iZ)[2] = boundaryVorticityZ(nx-1,iY,iZ);
        }
    }

    flags.vortFieldComputed = true;
}

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

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

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

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

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

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

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

            ++iPi;
        }
    }

    flags.strainRateFieldComputed = true;
}

template<typename T, template<typename U> class Lattice>
void DataAnalysis3D<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) {
            for (int iZ=0; iZ<fields.velField.getNz(); ++iZ) {
                block.get(iX,iY,iZ).computeStress(fields.stressField.get(iX,iY,iZ));
                T omega = block.get(iX,iY,iZ).getDynamics()->getOmega();
                for (int iPi=0; iPi<6; ++iPi) {
                    fields.stressField.get(iX,iY,iZ)[iPi] *= -omega / (T)2 * Lattice<T>::invCs2;
                }
            }
        }
    }

    flags.stressFieldComputed = true;
}

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

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

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

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

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

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

    flags.divRhoUFieldComputed = true;
}


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

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

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

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

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

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

    flags.poissonFieldComputed = true;
}

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

    int nx = fields.populationField.getNx();
    int ny = fields.populationField.getNy();
    int nz = fields.populationField.getNz();

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

    flags.populationFieldComputed = true;
}

}  // namespace olb


#endif