C++ volScalarField::internalField方法代码示例

void Foam::boundMinMax
(
    volScalarField& vsf,
    const dimensionedScalar& vsf0,
    const dimensionedScalar& vsf1
)
{
    scalar minVsf = min(vsf).value();
    scalar maxVsf = max(vsf).value();

    if (minVsf < vsf0.value() || maxVsf > vsf1.value())
    {
        Info<< "bounding " << vsf.name()
            << ", min: " << gMin(vsf.internalField())
            << " max: " << gMax(vsf.internalField())
            << " average: " << gAverage(vsf.internalField())
            << endl;
    }

    if (minVsf < vsf0.value())
    {
        vsf.internalField() = max
        (
            max
            (
                vsf.internalField(),
                fvc::average(max(vsf, vsf0))().internalField()
                *pos(vsf0.value() - vsf.internalField())
            ),
            vsf0.value()
        );

        vsf.correctBoundaryConditions();
        vsf.boundaryField() = max(vsf.boundaryField(), vsf0.value());
    }

    if (maxVsf > vsf1.value())
    {
        vsf.internalField() = min
        (
            min
            (
                vsf.internalField(),
                fvc::average(min(vsf, vsf1))().internalField()
                *neg(vsf1.value() - vsf.internalField())
                // This is needed when all values are above max
                // HJ, 18/Apr/2009
              + pos(vsf1.value() - vsf.internalField())*vsf1.value()
            ),
            vsf1.value()
        );

        vsf.correctBoundaryConditions();
        vsf.boundaryField() = min(vsf.boundaryField(), vsf1.value());
    }
}

tmp<GeometricField<Type, fvPatchField, volMesh> >
EulerLocalDdtScheme<Type>::fvcDdt
(
    const volScalarField& rho,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const objectRegistry& registry = this->mesh();

    // get access to the scalar beta[i]
    const scalarField& beta =
        registry.lookupObject<scalarField>(deltaTName_);

    volScalarField rDeltaT =
        1.0/(beta[0]*registry.lookupObject<volScalarField>(deltaTauName_));

    IOobject ddtIOobject
    (
        "ddt("+rho.name()+','+vf.name()+')',
        mesh().time().timeName(),
        mesh()
    );

    if (mesh().moving())
    {
        return tmp<GeometricField<Type, fvPatchField, volMesh> >
        (
            new GeometricField<Type, fvPatchField, volMesh>
            (
                ddtIOobject,
                mesh(),
                rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
                rDeltaT.internalField()*
                (
                    rho.internalField()*vf.internalField()
                  - rho.oldTime().internalField()
                   *vf.oldTime().internalField()*mesh().V0()/mesh().V()
                ),
                rDeltaT.boundaryField()*
                (
                    rho.boundaryField()*vf.boundaryField()
                  - rho.oldTime().boundaryField()
                   *vf.oldTime().boundaryField()
                )
            )
        );
    }
    else
    {
        return tmp<GeometricField<Type, fvPatchField, volMesh> >
        (
            new GeometricField<Type, fvPatchField, volMesh>
            (
                ddtIOobject,
                rDeltaT*(rho*vf - rho.oldTime()*vf.oldTime())
            )
        );
    }
}

Foam::tmp<Foam::volVectorField>
Foam::fv::faceLimitedGrad<Foam::scalar>::calcGrad
(
    const volScalarField& vsf,
    const word& name
) const
{
    const fvMesh& mesh = vsf.mesh();

    tmp<volVectorField> tGrad = basicGradScheme_().calcGrad(vsf, name);

    if (k_ < SMALL)
    {
        return tGrad;
    }

    volVectorField& g = tGrad();

    const labelUList& owner = mesh.owner();
    const labelUList& neighbour = mesh.neighbour();

    const volVectorField& C = mesh.C();
    const surfaceVectorField& Cf = mesh.Cf();

    // create limiter
    scalarField limiter(vsf.internalField().size(), 1.0);

    scalar rk = (1.0/k_ - 1.0);

    forAll(owner, facei)
    {
        label own = owner[facei];
        label nei = neighbour[facei];

        scalar vsfOwn = vsf[own];
        scalar vsfNei = vsf[nei];

        scalar maxFace = max(vsfOwn, vsfNei);
        scalar minFace = min(vsfOwn, vsfNei);
        scalar maxMinFace = rk*(maxFace - minFace);
        maxFace += maxMinFace;
        minFace -= maxMinFace;

        // owner side
        limitFace
        (
            limiter[own],
            maxFace - vsfOwn, minFace - vsfOwn,
            (Cf[facei] - C[own]) & g[own]
        );

        // neighbour side
        limitFace
        (
            limiter[nei],
            maxFace - vsfNei, minFace - vsfNei,
            (Cf[facei] - C[nei]) & g[nei]
        );
    }

void makeGraph
(
    const scalarField& x,
    const volScalarField& vsf,
    const word& name,
    const word& graphFormat
)
{
    makeGraph(x, vsf.internalField(), name, vsf.path(), graphFormat);
}

void Foam::bound(volScalarField& vsf, const dimensionedScalar& vsf0)
{
    scalar minVsf = min(vsf).value();

    if (minVsf < vsf0.value())
    {
        Info<< "bounding " << vsf.name()
            << ", min: " << gMin(vsf.internalField())
            << " max: " << gMax(vsf.internalField())
            << " average: " << gAverage(vsf.internalField())
            << endl;

        vsf.internalField() = max
        (
            max
            (
                vsf.internalField(),
                fvc::average(max(vsf, vsf0))().internalField()
                // Bug fix: was assuming bound on zero.  HJ, 25/Nov/2008
                *pos(vsf0.value() - vsf.internalField())
            ),
            vsf0.value()
        );

        vsf.correctBoundaryConditions();
        vsf.boundaryField() = max(vsf.boundaryField(), vsf0.value());
    }
}

void Foam::fv::limitTemperature::correct(volScalarField& he)
{
    const basicThermo& thermo =
        mesh_.lookupObject<basicThermo>(basicThermo::dictName);

    scalarField Tmin(cells_.size(), Tmin_);
    scalarField Tmax(cells_.size(), Tmax_);

    scalarField heMin(thermo.he(thermo.p(), Tmin, cells_));
    scalarField heMax(thermo.he(thermo.p(), Tmax, cells_));

    scalarField& hec = he.internalField();

    forAll(cells_, i)
    {
        label cellI = cells_[i];
        hec[cellI]= max(min(hec[cellI], heMax[i]), heMin[i]);
    }

tmp<fvMatrix<Type> >
EulerLocalDdtScheme<Type>::fvmDdt
(
    const volScalarField& rho,
    GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const objectRegistry& registry = this->mesh();

    // get access to the scalar beta[i]
    const scalarField& beta =
        registry.lookupObject<scalarField>(deltaTName_);

    tmp<fvMatrix<Type> > tfvm
    (
        new fvMatrix<Type>
        (
            vf,
            rho.dimensions()*vf.dimensions()*dimVol/dimTime
        )
    );
    fvMatrix<Type>& fvm = tfvm();

    scalarField rDeltaT =
        1.0/(beta[0]*registry.lookupObject<volScalarField>(deltaTauName_).internalField());

    fvm.diag() = rDeltaT*rho.internalField()*mesh().V();

    if (mesh().moving())
    {
        fvm.source() = rDeltaT
            *rho.oldTime().internalField()
            *vf.oldTime().internalField()*mesh().V0();
    }
    else
    {
        fvm.source() = rDeltaT
            *rho.oldTime().internalField()
            *vf.oldTime().internalField()*mesh().V();
    }

    return tfvm;
}

void makeGraph
(
    const scalarField& x,
    const volScalarField& vsf,
    const word& name,
    const word& graphFormat
)
{
    fileName path(vsf.rootPath()/vsf.caseName()/"graphs"/vsf.instance());
    mkDir(path);

    makeGraph
    (
        x,
        vsf.internalField(),
        name,
        path,
        graphFormat
    );
}

Foam::tmp<Foam::volVectorField>
Foam::fv::cellLimitedGrad<Foam::scalar>::calcGrad
(
    const volScalarField& vsf,
    const word& name
) const
{
    const fvMesh& mesh = vsf.mesh();

    tmp<volVectorField> tGrad = basicGradScheme_().calcGrad(vsf, name);

    if (k_ < SMALL)
    {
        return tGrad;
    }

    volVectorField& g = tGrad();

    const labelUList& owner = mesh.owner();
    const labelUList& neighbour = mesh.neighbour();

    const volVectorField& C = mesh.C();
    const surfaceVectorField& Cf = mesh.Cf();

    scalarField maxVsf(vsf.internalField());
    scalarField minVsf(vsf.internalField());

    forAll(owner, facei)
    {
        label own = owner[facei];
        label nei = neighbour[facei];

        scalar vsfOwn = vsf[own];
        scalar vsfNei = vsf[nei];

        maxVsf[own] = max(maxVsf[own], vsfNei);
        minVsf[own] = min(minVsf[own], vsfNei);

        maxVsf[nei] = max(maxVsf[nei], vsfOwn);
        minVsf[nei] = min(minVsf[nei], vsfOwn);
    }

void phaseChangeModel::correct
(
    const scalar dt,
    scalarField& availableMass,
    volScalarField& dMass,
    volScalarField& dEnergy
)
{
    correctModel
    (
        dt,
        availableMass,
        dMass,
        dEnergy
    );

    latestMassPC_ = sum(dMass.internalField());
    totalMassPC_ += latestMassPC_;

    availableMass -= dMass;
    dMass.correctBoundaryConditions();
}

tmp<volVectorField> cellMDLimitedGrad<scalar>::grad
(
    const volScalarField& vsf
) const
{
    const fvMesh& mesh = vsf.mesh();

    tmp<volVectorField> tGrad = basicGradScheme_().grad(vsf);

    if (k_ < SMALL)
    {
        return tGrad;
    }

    volVectorField& g = tGrad();

    const unallocLabelList& owner = mesh.owner();
    const unallocLabelList& neighbour = mesh.neighbour();

    const volVectorField& C = mesh.C();
    const surfaceVectorField& Cf = mesh.Cf();

    scalarField maxVsf(vsf.internalField());
    scalarField minVsf(vsf.internalField());

    forAll(owner, facei)
    {
        label own = owner[facei];
        label nei = neighbour[facei];

        scalar vsfOwn = vsf[own];
        scalar vsfNei = vsf[nei];

        maxVsf[own] = max(maxVsf[own], vsfNei);
        minVsf[own] = min(minVsf[own], vsfNei);

        maxVsf[nei] = max(maxVsf[nei], vsfOwn);
        minVsf[nei] = min(minVsf[nei], vsfOwn);
    }

//.........这里部分代码省略.........
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        mesh,
        dimless,
        allLambda,
        false   // Use slices for the couples
    );

    linear<scalar> CDs(mesh);
    upwind<scalar> UDs(mesh, phi);
    //fv::uncorrectedSnGrad<scalar> snGrads(mesh);

    fvScalarMatrix psiConvectionDiffusion
    (
        fvm::ddt(rho, psi)
      + fv::gaussConvectionScheme<scalar>(mesh, phi, UDs).fvmDiv(phi, psi)
        //- fv::gaussLaplacianScheme<scalar, scalar>(mesh, CDs, snGrads)
        //.fvmLaplacian(Dpsif, psi)
      - fvm::Sp(Sp, psi)
      - Su
    );

    surfaceScalarField phiBD(psiConvectionDiffusion.flux());

    surfaceScalarField& phiCorr = phiPsi;
    phiCorr -= phiBD;

    for (label i=0; i<maxIter; i++)
    {
        if (i != 0 && i < 4)
        {
            allLambda = allCoLambda;
        }

        limiter
        (
            allLambda,
            rho,
            psi,
            phiBD,
            phiCorr,
            Sp,
            Su,
            psiMax,
            psiMin,
            nLimiterIter
        );

        solve
        (
            psiConvectionDiffusion + fvc::div(lambda*phiCorr),
            MULEScontrols
        );

        scalar maxPsiM1 = gMax(psi.internalField()) - 1.0;
        scalar minPsi = gMin(psi.internalField());

        scalar unboundedness = max(max(maxPsiM1, 0.0), -min(minPsi, 0.0));

        if (unboundedness < maxUnboundedness)
        {
            break;
        }
        else
        {
            Info<< "MULES: max(" << psi.name() << " - 1) = " << maxPsiM1
                << " min(" << psi.name() << ") = " << minPsi << endl;

            phiBD = psiConvectionDiffusion.flux();

            /*
            word gammaScheme("div(phi,gamma)");
            word gammarScheme("div(phirb,gamma)");

            const surfaceScalarField& phir =
                mesh.lookupObject<surfaceScalarField>("phir");

            phiCorr =
                fvc::flux
                (
                    phi,
                    psi,
                    gammaScheme
                )
              + fvc::flux
                (
                    -fvc::flux(-phir, scalar(1) - psi, gammarScheme),
                    psi,
                    gammarScheme
                )
                - phiBD;
            */
        }
    }

    phiPsi = psiConvectionDiffusion.flux() + lambda*phiCorr;
}

void Foam::fv::interRegionExplicitPorositySource::addSup
(
    const volScalarField& rho,
    fvMatrix<vector>& eqn,
    const label fieldI
)
{
    initialise();

    const fvMesh& nbrMesh = mesh_.time().lookupObject<fvMesh>(nbrRegionName_);

    const volVectorField& U = eqn.psi();

    volVectorField UNbr
    (
        IOobject
        (
            name_ + ":UNbr",
            nbrMesh.time().timeName(),
            nbrMesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        nbrMesh,
        dimensionedVector("zero", U.dimensions(), vector::zero)
    );

    // map local velocity onto neighbour region
    meshInterp().mapSrcToTgt
    (
        U.internalField(),
        plusEqOp<vector>(),
        UNbr.internalField()
    );

    fvMatrix<vector> nbrEqn(UNbr, eqn.dimensions());

    volScalarField rhoNbr
    (
        IOobject
        (
            "rho:UNbr",
            nbrMesh.time().timeName(),
            nbrMesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        nbrMesh,
        dimensionedScalar("zero", dimDensity, 0.0)
    );

    volScalarField muNbr
    (
        IOobject
        (
            "mu:UNbr",
            nbrMesh.time().timeName(),
            nbrMesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        nbrMesh,
        dimensionedScalar("zero", dimViscosity, 0.0)
    );

    const volScalarField& mu =
        mesh_.lookupObject<volScalarField>(muName_);

    // map local rho onto neighbour region
    meshInterp().mapSrcToTgt
    (
        rho.internalField(),
        plusEqOp<scalar>(),
        rhoNbr.internalField()
    );

    // map local mu onto neighbour region
    meshInterp().mapSrcToTgt
    (
        mu.internalField(),
        plusEqOp<scalar>(),
        muNbr.internalField()
    );

    porosityPtr_->addResistance(nbrEqn, rhoNbr, muNbr);

    // convert source from neighbour to local region
    fvMatrix<vector> porosityEqn(U, eqn.dimensions());
    scalarField& Udiag = porosityEqn.diag();
    vectorField& Usource = porosityEqn.source();

    Udiag.setSize(eqn.diag().size(), 0.0);
    Usource.setSize(eqn.source().size(), vector::zero);

    meshInterp().mapTgtToSrc(nbrEqn.diag(), plusEqOp<scalar>(), Udiag);
    meshInterp().mapTgtToSrc(nbrEqn.source(), plusEqOp<vector>(), Usource);

    eqn -= porosityEqn;
}

// solve vegetation model
void modifiedVegetationModel::solve(volVectorField& U, volScalarField& T, volScalarField& w)
{
    // solve radiation within vegetation
    radiation();

    const double p_ = 101325;

    // Magnitude of velocity
    volScalarField magU("magU", mag(U));

    // Bounding velocity
    bound(magU, UMin_);

    // solve aerodynamic, stomatal resistance
    volScalarField new_Tl("new_Tl", Tl_);

    // info
    Info << "    max leaf temp tl=" << max(T.internalField())
         << "k, iteration i=0" << endl;

    scalar maxError, maxRelError;
    int i;

    // solve leaf temperature, iteratively.
    int maxIter = 500;
    for (i=1; i<=maxIter; i++)
    {
        // Solve aerodynamc, stomatal resistance
        resistance(magU, T, w, new_Tl);

        forAll(LAD_, cellI)
        {
            if (LAD_[cellI] > 10*SMALL)
            {
                // Initial leaf temperature
                if (i==1)
                    Tl_[cellI];// = T[cellI];//*0. + 300.;//T[cellI];

                // Calculate saturated density, specific humidity
                rhosat_[cellI] = calc_rhosat(Tl_[cellI]);
                evsat_[cellI] = calc_evsat(Tl_[cellI]);
                wsat_[cellI] = 0.621945*(evsat_[cellI]/(p_-evsat_[cellI])); // ASHRAE 1, eq.23

                // Calculate transpiration rate]);
                E_[cellI] = nEvapSides_.value()*LAD_[cellI]*rhoa_.value()*(wsat_[cellI]-w[cellI])/(ra_[cellI]+rs_[cellI]);
                //E_[cellI] = 0.0; // No evapotranspiration

                // Calculate latent heat flux
                Ql_[cellI] = lambda_.value()*E_[cellI];

                // Calculate new leaf temperature
                new_Tl[cellI] = T[cellI] + (Rn_[cellI] - Ql_[cellI])*(ra_[cellI]/(2.0*rhoa_.value()*cpa_.value()*LAD_[cellI]));

            }
        }

        // info
        Info << "    max leaf temp tl=" << gMax(new_Tl.internalField())
             << " K, iteration i="   << i << endl;

        // Check rel. L-infinity error
        maxError = gMax(mag(new_Tl.internalField()-Tl_.internalField()));
        maxRelError = maxError/gMax(mag(new_Tl.internalField()));

        // update leaf temp.
        forAll(Tl_, cellI)
            Tl_[cellI] = 0.5*Tl_[cellI]+0.5*new_Tl[cellI];

         // convergence check
         if (maxRelError < 1e-8)
             break;
    }
    Tl_.correctBoundaryConditions();

    // Iteration info
    Info << "Vegetation model:  Solving for Tl, Final residual = " << maxError
         << ", Final relative residual = " << maxRelError
         << ", No Iterations " << i << endl;

    Info << "temperature parameters: max Tl = " << gMax(Tl_)
         << ", min T = " << gMin(T) << ", max T = " << gMax(T) << endl;

    Info << "resistances: max rs = " << gMax(rs_)
         << ", max ra = " << gMax(ra_) << endl;

    // Final: Solve aerodynamc, stomatal resistance
    resistance(magU, T, w, Tl_);

    // Final: Update sensible and latent heat flux
    forAll(LAD_, cellI)
    {
        if (LAD_[cellI] > 10*SMALL)
        {
            // Calculate saturated density, specific humidity
            rhosat_[cellI] = calc_rhosat(Tl_[cellI]);
            evsat_[cellI] = calc_evsat(Tl_[cellI]);
            wsat_[cellI] = 0.621945*(evsat_[cellI]/(p_-evsat_[cellI])); // ASHRAE 1, eq.23

            // Calculate transpiration rate
//.........这里部分代码省略.........

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void scalarGeneralExchange::manipulateScalarField(volScalarField& explicitEulerSource,
                                                  volScalarField& implicitEulerSource,
                                                  int speciesID) const
{

    // reset Scalar field
    explicitEulerSource.internalField() = 0.0;
    implicitEulerSource.internalField() = 0.0;

    if(speciesID>=0 && particleSpeciesValue_[speciesID]<0.0)    //skip if species is not active
        return;

    //Set the names of the exchange fields
    word    fieldName;
    word    partDatName;
    word    partFluxName;
    word    partTransCoeffName;
    word    partFluidName;
    scalar  transportParameter;

    // realloc the arrays to pull particle data


    if(speciesID<0) //this is the temperature - always pull from LIGGGHTS
    {
        fieldName          = tempFieldName_;
        partDatName        = partTempName_;
        partFluxName       = partHeatFluxName_;
        partTransCoeffName = partHeatTransCoeffName_;
        partFluidName      = partHeatFluidName_;
        transportParameter = lambda_;

        allocateMyArrays(0.0);
        particleCloud_.dataExchangeM().getData(partDatName,"scalar-atom", partDat_);
    }
    else
    {
        fieldName          = eulerianFieldNames_[speciesID];
        partDatName        = partSpeciesNames_[speciesID];
        partFluxName       = partSpeciesFluxNames_[speciesID]; 
        partTransCoeffName = partSpeciesTransCoeffNames_[speciesID]; 
        partFluidName      = partSpeciesFluidNames_[speciesID]; 
        transportParameter = DMolecular_[speciesID];

        allocateMyArrays(0.0);
        if(particleSpeciesValue_[speciesID]>ALARGECONCENTRATION)   
            particleCloud_.dataExchangeM().getData(partDatName,"scalar-atom", partDat_);
    }

    if (scaleDia_ > 1)
        Info << typeName << " using scale = " << scaleDia_ << endl;
    else if (particleCloud_.cg() > 1)
    {
        scaleDia_=particleCloud_.cg();
        Info << typeName << " using scale from liggghts cg = " << scaleDia_ << endl;
    }

    //==============================
    // get references
    const volScalarField& voidfraction_(particleCloud_.mesh().lookupObject<volScalarField> (voidfractionFieldName_));    // ref to voidfraction field
    const volVectorField& U_(particleCloud_.mesh().lookupObject<volVectorField> (velFieldName_));
    const volScalarField& fluidScalarField_(particleCloud_.mesh().lookupObject<volScalarField> (fieldName));            // ref to scalar field
    const volScalarField& nufField = forceSubM(0).nuField();
    //==============================

    // calc La based heat flux
    vector position(0,0,0);
    scalar voidfraction(1);
    vector Ufluid(0,0,0);
    scalar fluidValue(0);
    label  cellI=0;
    vector Us(0,0,0);
    vector Ur(0,0,0);
    scalar dscaled(0);
    scalar nuf(0);
    scalar magUr(0);
    scalar As(0);
    scalar Rep(0);
    scalar Pr(0);

    scalar sDth(scaleDia_*scaleDia_*scaleDia_);

    interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
    interpolationCellPoint<vector> UInterpolator_(U_);
    interpolationCellPoint<scalar> fluidScalarFieldInterpolator_(fluidScalarField_);

    #include "setupProbeModel.H"

    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
    {
            cellI = particleCloud_.cellIDs()[index][0];
            if(cellI >= 0)
            {
                if(forceSubM(0).interpolation())
                {
	                position = particleCloud_.position(index);
                    voidfraction = voidfractionInterpolator_.interpolate(position,cellI);
                    Ufluid = UInterpolator_.interpolate(position,cellI);
                    fluidValue = fluidScalarFieldInterpolator_.interpolate(position,cellI);
//.........这里部分代码省略.........