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

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());
    }
}

void averagingModel::setScalarSum
(
    volScalarField& field,
    double**& value,
    double**const& weight,
    double**const& mask
) const
{
    label cellI;
    scalar valueScal;
    scalar weightP;

    for(int index=0; index< particleCloud_.numberOfParticles(); index++)
    {
        if(mask[index][0])
        {
            for(int subCell=0;subCell<particleCloud_.voidFractionM().cellsPerParticle()[index][0];subCell++)
            {
                //Info << "subCell=" << subCell << endl;
                cellI = particleCloud_.cellIDs()[index][subCell];

                if (cellI >= 0)
                {
                    valueScal = value[index][0];
                    weightP = weight[index][subCell];
                    field[cellI] += valueScal*weightP;
                }
            }//forAllSubPoints
        }
    }

    // correct cell values to patches
    field.correctBoundaryConditions();
}

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 phaseChangeModel::correct
(
    const scalar dt,
    scalarField& availableMass,
    volScalarField& dMass,
    volScalarField& dEnergy
)
{
    if (!active())
    {
        return;
    }

    correctModel
    (
        dt,
        availableMass,
        dMass,
        dEnergy
    );

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

    availableMass -= dMass;
    dMass.correctBoundaryConditions();

    if (writeTime())
    {
        scalar phaseChangeMass = getModelProperty<scalar>("phaseChangeMass");
        phaseChangeMass += returnReduce(totalMassPC_, sumOp<scalar>());
        setModelProperty<scalar>("phaseChangeMass", phaseChangeMass);
        totalMassPC_ = 0.0;
    }
}

void Foam::MULES::explicitSolve
(
    const RhoType& rho,
    volScalarField& psi,
    const surfaceScalarField& phi,
    surfaceScalarField& phiPsi,
    const SpType& Sp,
    const SuType& Su,
    const scalar psiMax,
    const scalar psiMin
)
{
    const fvMesh& mesh = psi.mesh();

    psi.correctBoundaryConditions();

    if (fv::localEulerDdt::enabled(mesh))
    {
        const volScalarField& rDeltaT = fv::localEulerDdt::localRDeltaT(mesh);

        limit
        (
            rDeltaT,
            rho,
            psi,
            phi,
            phiPsi,
            Sp,
            Su,
            psiMax,
            psiMin,
            false
        );

        explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
    }
    else
    {
        const scalar rDeltaT = 1.0/mesh.time().deltaTValue();

        limit
        (
            rDeltaT,
            rho,
            psi,
            phi,
            phiPsi,
            Sp,
            Su,
            psiMax,
            psiMin,
            false
        );

        explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
    }
}

void Foam::MULES::explicitSolve
(
    const RhoType& rho,
    volScalarField& psi,
    const surfaceScalarField& phi,
    surfaceScalarField& phiPsi,
    const SpType& Sp,
    const SuType& Su,
    const scalar psiMax,
    const scalar psiMin
)
{
    psi.correctBoundaryConditions();
    limit(rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false);
    explicitSolve(rho, psi, phiPsi, Sp, Su);
}

void dense::setScalarAverage
(
    volScalarField& field,
    double**& value,
    double**& weight,
    volScalarField& weightField,
    double**const& mask
) const
{
    label cellI;
    scalar valueScal;
    scalar weightP;

    for(int index=0; index< particleCloud_.numberOfParticles(); index++)
    {
        if(mask[index][0])
        {
            for(int subCell=0;subCell<particleCloud_.voidFractionM().cellsPerParticle()[index][0];subCell++)
            {
                //Info << "subCell=" << subCell << endl;
                cellI = particleCloud_.cellIDs()[index][subCell];

                if (cellI >= 0)
                {
                    valueScal = value[index][0];
                    weightP = weight[index][0];

                    // first entry in this cell
                    if(weightField[cellI] == 0)
                    {
                        field[cellI] = valueScal;
                        weightField[cellI] = weightP;
                    }
                    else
                    {
                        field[cellI] = (field[cellI]*weightField[cellI]+valueScal*weightP)/(weightField[cellI]+weightP);
                        weightField[cellI] += weightP;
                    }
                }
            }
        }
    }

    // correct cell values to patches
    field.correctBoundaryConditions();
}

void Foam::MULES::explicitSolve
(
    const RhoType& rho,
    volScalarField& psi,
    const surfaceScalarField& phi,
    surfaceScalarField& phiPsi,
    const SpType& Sp,
    const SuType& Su,
    const scalar psiMax,
    const scalar psiMin
)
{
    const fvMesh& mesh = psi.mesh();
    const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
    psi.correctBoundaryConditions();
    limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false);
    explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
}

void Foam::MULES::explicitSolve
(
    const RdeltaTType& rDeltaT,
    const RhoType& rho,
    volScalarField& psi,
    const surfaceScalarField& phiPsi,
    const SpType& Sp,
    const SuType& Su
)
{
    Info<< "MULES: Solving for " << psi.name() << endl;

    const fvMesh& mesh = psi.mesh();

    scalarField& psiIf = psi;
    const scalarField& psi0 = psi.oldTime();

    psiIf = 0.0;
    fvc::surfaceIntegrate(psiIf, phiPsi);

    if (mesh.moving())
    {
        psiIf =
        (
            mesh.Vsc0()().field()*rho.oldTime().field()
           *psi0*rDeltaT/mesh.Vsc()().field()
          + Su.field()
          - psiIf
        )/(rho.field()*rDeltaT - Sp.field());
    }
    else
    {
        psiIf =
        (
            rho.oldTime().field()*psi0*rDeltaT
          + Su.field()
          - psiIf
        )/(rho.field()*rDeltaT - Sp.field());
    }

    psi.correctBoundaryConditions();
}

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();
}

void dilute::setScalarAverage
(
    volScalarField& field,
    double**& value,
    double**& weight,
    volScalarField& weightField,
    double**const& mask,
    double**const& weight2,       //allows the specification of a 2nd weight field
    bool      weightWithWeight2   //switch to activate 2nd weight field
) const
{
    label cellI;
    scalar valueScal;
    scalar weightP;

    if(weightWithWeight2) 
        FatalError << "dilute::setScalarAverage: attempt to weight with weight2, which is not implemented" << abort(FatalError);

    for(int index=0; index< particleCloud_.numberOfParticles(); index++)
    {
        for(int subCell=0;subCell<particleCloud_.cellsPerParticle()[index][0];subCell++)
        {
            //Info << "subCell=" << subCell << endl;
            cellI = particleCloud_.cellIDs()[index][subCell];

            if (cellI >= 0)
            {
                valueScal = value[index][0];
                weightP = weight[index][0];
                weightField[cellI] += weightP;

                field[cellI] = valueScal/weightP;
            }
        }
    }

    // correct cell values to patches
    field.correctBoundaryConditions();
}

bool Foam::cfdemCloud::evolve
(
    volScalarField& alpha,
    volVectorField& Us,
    volVectorField& U
)
{
    numberOfParticlesChanged_ = false;
    arraysReallocated_=false;
    bool doCouple=false;

    if(!ignore())
    {
        if (dataExchangeM().doCoupleNow())
        {
            Info << "\n Coupling..." << endl;
            dataExchangeM().couple(0);
            doCouple=true;

            // reset vol Fields
            clockM().start(16,"resetVolFields");
            if(verbose_)
            {
                Info << "couplingStep:" << dataExchangeM().couplingStep() 
                     << "\n- resetVolFields()" << endl;
            }
            averagingM().resetVectorAverage(averagingM().UsPrev(),averagingM().UsNext(),false);
            resetVoidFraction();
            averagingM().resetVectorAverage(forceM(0).impParticleForces(),forceM(0).impParticleForces(),true);
            averagingM().resetVectorAverage(forceM(0).expParticleForces(),forceM(0).expParticleForces(),true);
            averagingM().resetWeightFields();
            for (int i=0;i<momCoupleModels_.size(); i++)
                momCoupleM(i).resetMomSourceField();
            if(verbose_) Info << "resetVolFields done." << endl;
            clockM().stop("resetVolFields");

            if(verbose_) Info << "- getDEMdata()" << endl;
            clockM().start(17,"getDEMdata");
            getDEMdata();
            clockM().stop("getDEMdata");
            if(verbose_) Info << "- getDEMdata done." << endl;

            // search cellID of particles
            clockM().start(18,"findCell");
            if(verbose_) Info << "- findCell()" << endl;
            findCells();
            if(verbose_) Info << "findCell done." << endl;
            clockM().stop("findCell");

            // set void fraction field
            clockM().start(19,"setvoidFraction");
            if(verbose_) Info << "- setvoidFraction()" << endl;
            setVoidFraction();
            if(verbose_) Info << "setvoidFraction done." << endl;
            clockM().stop("setvoidFraction");

            // set average particles velocity field
            clockM().start(20,"setVectorAverage");
            setVectorAverages();


            //Smoothen "next" fields            
            smoothingM().dSmoothing();
            smoothingM().smoothen(voidFractionM().voidFractionNext());

            //only smoothen if we use implicit force coupling in cells void of particles
            //because we need unsmoothened Us field to detect cells for explicit 
            //force coupling
            if(!treatVoidCellsAsExplicitForce())
                smoothingM().smoothenReferenceField(averagingM().UsNext());
            
            clockM().stop("setVectorAverage");
        }
        
        //============================================
        //CHECK JUST TIME-INTERPOATE ALREADY SMOOTHENED VOIDFRACTIONNEXT AND UsNEXT FIELD 
        //      IMPLICIT FORCE CONTRIBUTION AND SOLVER USE EXACTLY THE SAME AVERAGED
        //      QUANTITIES AT THE GRID!
        Info << "\n timeStepFraction() = " << dataExchangeM().timeStepFraction() << endl;
        clockM().start(24,"interpolateEulerFields");

        // update voidFractionField
        setAlpha(alpha);
        if(dataExchangeM().couplingStep() < 2)
        {
            alpha.oldTime() = alpha; // supress volume src
            alpha.oldTime().correctBoundaryConditions();
        }
        alpha.correctBoundaryConditions();

        // calc ddt(voidfraction)
        calcDdtVoidfraction(alpha,Us);

        // update mean particle velocity Field
        Us = averagingM().UsInterp();
        Us.correctBoundaryConditions();

        clockM().stop("interpolateEulerFields");
        //============================================

//.........这里部分代码省略.........

void Foam::MULES::explicitSolve
(
    const RhoType& rho,
    volScalarField& psi,
    const surfaceScalarField& phi,
    surfaceScalarField& phiPsi,
    const SpType& Sp,
    const SuType& Su,
    const scalar psiMax,
    const scalar psiMin
)
{
    Info<< "MULES: Solving for " << psi.name() << endl;

    const fvMesh& mesh = psi.mesh();
    psi.correctBoundaryConditions();

    surfaceScalarField phiBD(upwind<scalar>(psi.mesh(), phi).flux(psi));

    surfaceScalarField& phiCorr = phiPsi;
    phiCorr -= phiBD;

    scalarField allLambda(mesh.nFaces(), 1.0);

    slicedSurfaceScalarField lambda
    (
        IOobject
        (
            "lambda",
            mesh.time().timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        mesh,
        dimless,
        allLambda,
        false   // Use slices for the couples
    );

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

    phiPsi = phiBD + lambda*phiCorr;

    scalarField& psiIf = psi;
    const scalarField& psi0 = psi.oldTime();
    const scalar deltaT = mesh.time().deltaTValue();

    psiIf = 0.0;
    fvc::surfaceIntegrate(psiIf, phiPsi);

    if (mesh.moving())
    {
        psiIf =
        (
            mesh.Vsc0()().field()*rho.oldTime().field()
           *psi0/(deltaT*mesh.Vsc()().field())
          + Su.field()
          - psiIf
        )/(rho.field()/deltaT - Sp.field());
    }
    else
    {
        psiIf =
        (
            rho.oldTime().field()*psi0/deltaT
          + Su.field()
          - psiIf
        )/(rho.field()/deltaT - Sp.field());
    }

    psi.correctBoundaryConditions();
}

//.........这里部分代码省略.........
            clockM().start(19,"setvoidFraction");
            if(verbose_) Info << "- setvoidFraction()" << endl;
            voidFractionM().setvoidFraction(NULL,voidfractions_,particleWeights_,particleVolumes_);
            if(verbose_) Info << "setvoidFraction done." << endl;
            clockM().stop("setvoidFraction");

            // set particles velocity field
            clockM().start(20,"setVectorAverage");
            setVectorAverages();
            clockM().stop("setVectorAverage");

            // set particles forces
            clockM().start(21,"setForce");
            if(verbose_) Info << "- setForce(forces_)" << endl;
            setForces();
            if(verbose_) Info << "setForce done." << endl;
            clockM().stop("setForce");

            // get next force field
            clockM().start(22,"setParticleForceField");
            if(verbose_) Info << "- setParticleForceField()" << endl;
            averagingM().setVectorSum
            (
                forceM(0).impParticleForces(),
                impForces_,
                particleWeights_,
                NULL //mask
            );
            averagingM().setVectorSum
            (
                forceM(0).expParticleForces(),
                expForces_,
                particleWeights_,
                NULL //mask
            );
            if(verbose_) Info << "- setParticleForceField done." << endl;
            clockM().stop("setParticleForceField");

            // write DEM data
            if(verbose_) Info << " -giveDEMdata()" << endl;
            clockM().start(23,"giveDEMdata");
            giveDEMdata();
            clockM().stop("giveDEMdata");
        }//end dataExchangeM().couple()
        Info << "\n timeStepFraction() = " << dataExchangeM().timeStepFraction() << endl;

        clockM().start(24,"interpolateEulerFields");
        // update smoothing model
        smoothingM().dSmoothing();

        //============================================
        // update voidFractionField V1
        alpha = voidFractionM().voidFractionInterp();
        smoothingM().smoothen(alpha);
        if(dataExchangeM().couplingStep() < 2)
        {
            alpha.oldTime() = alpha; // supress volume src
            alpha.oldTime().correctBoundaryConditions();
        }
        alpha.correctBoundaryConditions();

        // calc ddt(voidfraction)
        //calcDdtVoidfraction(voidFractionM().voidFractionNext());
        calcDdtVoidfraction(alpha);

        // update particle velocity Field
        Us = averagingM().UsInterp();
        //smoothingM().smoothenReferenceField(Us);
        Us.correctBoundaryConditions();
        /*//============================================
        // update voidFractionField
        volScalarField oldAlpha = alpha.oldTime(); //save old (smooth) alpha field
        alpha.oldTime().internalField() = voidFractionM().voidFractionInterp();
        smoothingM().smoothen(alpha);
        alpha.correctBoundaryConditions();
        alpha.oldTime() = oldAlpha; //set old (smooth) alpha field to allow correct computation of ddt

        // calc ddt(voidfraction)
        if (doCouple) calcDdtVoidfraction(alpha);
        //calcDdtVoidfraction(alpha); // alternative with scale=1! (does not see change in alpha?)

        // update particle velocity Field
        Us.oldTime().internalField() = averagingM().UsInterp();
        smoothingM().smoothenReferenceField(Us);
        Us.correctBoundaryConditions();
        //============================================*/
        clockM().stop("interpolateEulerFields");

        if(verbose_){
            #include "debugInfo.H"
        }

        clockM().start(25,"dumpDEMdata");
        // do particle IO
        IOM().dumpDEMdata();
        clockM().stop("dumpDEMdata");

    }//end ignore
    return doCouple;
}

void averagingModel::setDSauter
(
    volScalarField& dSauter,
    double**& weight,
    volScalarField& weightField,
    label myParticleType
) const
{
    label cellI;
    scalar valueScal;
    scalar weightP;
    scalar radius(-1);
    scalar radiusPow2(-1);
    scalar radiusPow3(-1);
    scalar volume(-1);

    scalar scale_ = particleCloud_.cg(); //scaling of parcel vs. primary particle diameter
    dSauter = 0.0 * dSauter; //set to zero, because we will use it to calc sum(wi*ri^3)
    volScalarField riPower2
    (
        IOobject
        (
            "dummy2",
            particleCloud_.mesh().time().timeName(),
            particleCloud_.mesh(),
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        particleCloud_.mesh(),
        dimensionedScalar("zero", dimensionSet(0, 0, 0, 0, 0),0)
    );

    for(int index=0; index< particleCloud_.numberOfParticles(); index++)
    {
        if(myParticleType!=0) //in case a particle type is specified, only consider particles of the right type
            if(myParticleType != particleCloud_.particleType(index)) continue; 

        radius         = particleCloud_.radii()[index][0] / scale_; //the primary particle diameter
        radiusPow2 = radius*radius;
        radiusPow3 = radiusPow2*radius;
        weightP      = weight[index][0];

        for(int subCell=0;subCell<particleCloud_.cellsPerParticle()[index][0];subCell++)
        {

            cellI = particleCloud_.cellIDs()[index][subCell];
            if (cellI >= 0)
            {
                // first entry in this cell
                if(weightField[cellI] == 0)
                {
                    dSauter[cellI]      = radiusPow3; //use dSauter to store sum(ri^3)
                    riPower2[cellI]    = radiusPow2;
                    weightField[cellI] = weightP;
                }
                else
                {
                    dSauter[cellI] = (dSauter[cellI]*weightField[cellI]+radiusPow3*weightP)
                                    /(weightField[cellI]+weightP);
                    riPower2[cellI] = (riPower2[cellI]*weightField[cellI]+radiusPow2*weightP)
                                    /(weightField[cellI]+weightP);
                    weightField[cellI] += weightP;
                }
            }
        }
    }

    // set value and correct cell values to patches
    dSauter=2.0*dSauter / (riPower2+1e-99);
    dSauter.correctBoundaryConditions();

    return;
}