dataanalysis3d.hh

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

        }
    }
    flags.velFieldComputed = true;
}

template<typename T, template<typename U> class Lattice>
void DataAnalysis3D<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) {
            for (int iZ=0; iZ<fields.momentumField.getNz(); ++iZ) {
                T rho;
                block.get(iX,iY,iZ).computeRhoU (
                        rho, fields.momentumField.get(iX,iY,iZ) );
                for (int iD=0; iD<Lattice<T>::d; ++iD) {
                    fields.momentumField.get(iX,iY,iZ)[iD] *= rho;
                }
            }
        }
    }
    flags.momentumFieldComputed = true;
}

template<typename T, template<typename U> class Lattice>
void DataAnalysis3D<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) {
            for (int iZ=0; iZ<fields.pressureField.getNz(); ++iZ) {
                fields.pressureField.get(iX,iY,iZ) =
                    block.get(iX,iY,iZ).computeRho();
                fields.pressureField.get(iX,iY,iZ) -= (T)1;
                fields.pressureField.get(iX,iY,iZ) /= Lattice<T>::invCs2;
            }
        }
    }
    flags.pressureFieldComputed = true;
}

template<typename T, template<typename U> class Lattice>
void DataAnalysis3D<T,Lattice>::computeVelocityNormField() const {
    if (flags.velNormFieldComputed) return;
    fields.velNormField.construct();
    computeVelocityField();
    for (size_t iEl=0; iEl<fields.velNormField.getSize(); ++iEl) {
        fields.velNormField[iEl] = sqrt(util::normSqr<T,3>(fields.velField[iEl]));
    }
    flags.velNormFieldComputed = true;
}

template<typename T, template<typename U> class Lattice>
void DataAnalysis3D<T,Lattice>::computeVorticityNormField() const {
    if (flags.vortNormFieldComputed) return;
    fields.vortNormField.construct();
    computeVorticityField();
    for (size_t iEl=0; iEl<fields.vortNormField.getSize(); ++iEl) {
        fields.vortNormField[iEl] = sqrt(util::normSqr<T,3>(fields.vortField[iEl]));
    }
    flags.vortNormFieldComputed = true;
}

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

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

    return dyuz - dzuy;
}

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

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

    return dzux - dxuz;
}

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

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

    return dxuy - dyux;
}

template<typename T, template<typename U> class Lattice>
T DataAnalysis3D<T,Lattice>::bulkXderiv (
        int iX, int iY, int iZ, int iD,
        TensorField3D<T,3> const& field) const
{
    OLB_PRECONDITION(iX>=1 && iX<=field.getNx()-2);
    OLB_PRECONDITION(iY>=1 && iY<=field.getNy()-2);
    OLB_PRECONDITION(iZ>=1 && iZ<=field.getNz()-2);

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

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

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

template<typename T, template<typename U> class Lattice>
T DataAnalysis3D<T,Lattice>::bulkZderiv (
        int iX, int iY, int iZ, int iD,
        TensorField3D<T,3> const& field) const
{
    OLB_PRECONDITION(iX>=1 && iX<=field.getNx()-2);
    OLB_PRECONDITION(iY>=1 && iY<=field.getNy()-2);
    OLB_PRECONDITION(iZ>=1 && iZ<=field.getNz()-2);

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

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

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

}

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

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

    T dxux = bulkDeriv(iX,iY,iZ, 0,0, fields.velField);
    T dxuy = bulkDeriv(iX,iY,iZ, 0,1, fields.velField);
    T dxuz = bulkDeriv(iX,iY,iZ, 0,2, fields.velField);
    T dyux = bulkDeriv(iX,iY,iZ, 1,0, fields.velField);
    T dyuy = bulkDeriv(iX,iY,iZ, 1,1, fields.velField);
    T dyuz = bulkDeriv(iX,iY,iZ, 1,2, fields.velField);
    T dzux = bulkDeriv(iX,iY,iZ, 2,0, fields.velField);
    T dzuy = bulkDeriv(iX,iY,iZ, 2,1, fields.velField);
    T dzuz = bulkDeriv(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>::boundaryXderiv (
        int iX, int iY, int iZ, int iD,
        TensorField3D<T,3> const& field) const
{
    OLB_PRECONDITION(iX>=0 && iX<=field.getNx()-1);
    OLB_PRECONDITION(iY>=0 && iY<=field.getNy()-1);
    OLB_PRECONDITION(iZ>=0 && iZ<=field.getNz()-1);

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

    return dxu;
}

template<typename T, template<typename U> class Lattice>
T DataAnalysis3D<T,Lattice>::boundaryYderiv (
        int iX, int iY, int iZ, int iD,
        TensorField3D<T,3> const& field) const
{
    OLB_PRECONDITION(iX>=0 && iX<=field.getNx()-1);
    OLB_PRECONDITION(iY>=0 && iY<=field.getNy()-1);
    OLB_PRECONDITION(iZ>=0 && iZ<=field.getNz()-1);

    T dyu;

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

template<typename T, template<typename U> class Lattice>
T DataAnalysis3D<T,Lattice>::boundaryZderiv (
        int iX, int iY, int iZ, int iD,
        TensorField3D<T,3> const& field) const
{
    OLB_PRECONDITION(iX>=0 && iX<=field.getNx()-1);
    OLB_PRECONDITION(iY>=0 && iY<=field.getNy()-1);
    OLB_PRECONDITION(iZ>=0 && iZ<=field.getNz()-1);

    T dzu;

    if (iZ==0) {
        dzu = fd::boundaryGradient(field.get(iX,iY,iZ)[iD],
                                   field.get(iX,iY,iZ+1)[iD],
                                   field.get(iX,iY,iZ+2)[iD]);
    }
    else if (iZ==field.getNz()-1) {
        dzu = -fd::boundaryGradient(field.get(iX,iY,iZ)[iD],
                                    field.get(iX,iY,iZ-1)[iD],
                                    field.get(iX,iY,iZ-2)[iD]);
    }
    else {
        dzu = fd::centralGradient(field.get(iX,iY,iZ+1)[iD],
                                  field.get(iX,iY,iZ-1)[iD]);
    }
   
    return dzu;
}

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

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

}

template<typename T, template<typename U> class Lattice>
T DataAnalysis3D<T,Lattice>::boundaryDivRhoU (
        int iX, int iY, int iZ) const
{