Foam::tmp<Foam::volScalarField>
Foam::diameterModels::IATEsources::dummy::R() const
{
    return volScalarField::New
    (
        "R",
        iate_.phase().U().mesh(),
        dimensionedScalar(dimless/dimTime, 0)
    );
}

void Foam::inhomogeneousMixture<ThermoType>::read(const dictionary& thermoDict)
{
    stoicRatio_ = dimensionedScalar
    (
        thermoDict.lookup("stoichiometricAirFuelMassRatio")
    );

    fuel_ = ThermoType(thermoDict.lookup("fuel"));
    oxidant_ = ThermoType(thermoDict.lookup("oxidant"));
    products_ = ThermoType(thermoDict.lookup("burntProducts"));
}

Foam::processorField::processorField
(
    const word& name,
    const objectRegistry& obr,
    const dictionary& dict,
    const bool loadFromFiles
)
:
    name_(name),
    obr_(obr),
    active_(true)
{
    // Check if the available mesh is an fvMesh otherise deactivate
    if (isA<fvMesh>(obr_))
    {
        read(dict);

        const fvMesh& mesh = refCast<const fvMesh>(obr_);

        volScalarField* procFieldPtr
        (
            new volScalarField
            (
                IOobject
                (
                    "processorID",
                    mesh.time().timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::NO_WRITE
                ),
                mesh,
                dimensionedScalar("0", dimless, 0.0)
            )
        );

        mesh.objectRegistry::store(procFieldPtr);
    }
    else
    {
        active_ = false;
        WarningIn
        (
            "processorField::processorField"
            "("
                "const word&, "
                "const objectRegistry&, "
                "const dictionary&, "
                "const bool"
            ")"
        )   << "No fvMesh available, deactivating " << name_
            << endl;
    }
}

Foam::tmp<Foam::volScalarField> Foam::JohnsonJacksonFrictionalStress::muf
(
    const volScalarField& alpha,
    const dimensionedScalar& alphaMax,
    const volScalarField& pf,
    const volSymmTensorField& D,
    const dimensionedScalar& phi
) const
{
    return dimensionedScalar("0.5", dimTime, 0.5)*pf*sin(phi);
}

void Foam::CellAverageParticleVelocity<CloudType>::preEvolve()
{
    if (UpPtr_.valid())
    {
        UpPtr_->internalField() = vector::zero;
    }
    else
    {
        const fvMesh& mesh = this->owner().mesh();

        UpPtr_.reset
        (
            new volVectorField
            (
                IOobject
                (
                    this->owner().name() + "Up",
                    mesh.time().timeName(),
                    mesh,
                    IOobject::MUST_READ,
                    IOobject::AUTO_WRITE
                ),
                mesh,
                dimensionedVector("zeroVector", dimVelocity, vector::zero)
            )
        );
    }

    if (pVolPtr_.valid())
    {
        pVolPtr_->internalField() = scalar(0);
    }
    else
    {
        const fvMesh& mesh = this->owner().mesh();

        pVolPtr_.reset
        (
            new volScalarField
            (
                IOobject
                (
                    this->owner().name() + "pVol",
                    mesh.time().timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::NO_WRITE
                ),
                mesh,
                dimensionedScalar("zeroVolume", dimVolume, scalar(0))
            )
        );
    }
}

Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::nu
(
    const volScalarField& alpha1,
    const dimensionedScalar& alphaMax,
    const volScalarField& pf,
    const volSymmTensorField& D
) const
{
    const scalar I2Dsmall = 1.0e-15;

    // Creating nu assuming it should be 0 on the boundary which may not be
    // true
    tmp<volScalarField> tnu
    (
        new volScalarField
        (
            IOobject
            (
                "Schaeffer:nu",
                alpha1.mesh().time().timeName(),
                alpha1.mesh(),
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            alpha1.mesh(),
            dimensionedScalar("nu", dimensionSet(0, 2, -1, 0, 0), 0.0)
        )
    );

    volScalarField& nuf = tnu();

    forAll (D, celli)
    {
        if (alpha1[celli] > alphaMax.value() - 5e-2)
        {
            nuf[celli] =
                0.5*pf[celli]*sin(phi_.value())
               /(
                    sqrt(1.0/6.0*(sqr(D[celli].xx() - D[celli].yy())
                  + sqr(D[celli].yy() - D[celli].zz())
                  + sqr(D[celli].zz() - D[celli].xx()))
                  + sqr(D[celli].xy()) + sqr(D[celli].xz())
                  + sqr(D[celli].yz())) + I2Dsmall
                );
        }
    }

    // Correct coupled BCs
    nuf.correctBoundaryConditions();

    return tnu;
}

void Foam::processorField::execute()
{
    if (active_)
    {
        const volScalarField& procField =
            obr_.lookupObject<volScalarField>("processorID");

        const_cast<volScalarField&>(procField) ==
            dimensionedScalar("procI", dimless, Pstream::myProcNo());
    }
}

//- Construct from objectRegistry arguments
Foam::engineTime::engineTime
(
    const word& name,
    const fileName& rootPath,
    const fileName& caseName,
    const fileName& systemName,
    const fileName& constantName,
    const fileName& dictName
)
:
    Time
    (
        name,
        rootPath,
        caseName,
        systemName,
        constantName
    ),
    dict_
    (
        IOobject
        (
            dictName,
            constant(),
            *this,
            IOobject::MUST_READ,
            IOobject::NO_WRITE,
            false
        )
    ),
    rpm_(dict_.lookup("rpm")),
    conRodLength_(dimensionedScalar("conRodLength", dimLength, 0)),
    bore_(dimensionedScalar("bore", dimLength, 0)),
    stroke_(dimensionedScalar("stroke", dimLength, 0)),
    clearance_(dimensionedScalar("clearance", dimLength, 0))
{
    // the geometric parameters are not strictly required for Time
    if (dict_.found("conRodLength"))
    {
        dict_.lookup("conRodLength") >> conRodLength_;
    }

volScalarField dynamicKEqn<BasicTurbulenceModel>::Ce() const
{
    const volSymmTensorField D(dev(symm(fvc::grad(this->U_))));

    volScalarField KK
    (
        0.5*(filter_(magSqr(this->U_)) - magSqr(filter_(this->U_)))
    );
    KK.max(dimensionedScalar("small", KK.dimensions(), SMALL));

    return Ce(D, KK);
}

Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::frictionalStressModels::JohnsonJackson::nu
(
    const phaseModel& phase,
    const dimensionedScalar& alphaMinFriction,
    const dimensionedScalar& alphaMax,
    const volScalarField& pf,
    const volSymmTensorField& D
) const
{
    return dimensionedScalar("0.5", dimTime, 0.5)*pf*sin(phi_);
}

tmp<volScalarField::Internal> kinematicSingleLayer::Srho
(
    const label i
) const
{
    return volScalarField::Internal::New
    (
        typeName + ":Srho(" + Foam::name(i) + ")",
        primaryMesh(),
        dimensionedScalar(dimMass/dimVolume/dimTime, 0)
    );
}

void kinematicSingleLayer::preEvolveRegion()
{
    if (debug)
    {
        InfoInFunction << endl;
    }

    surfaceFilmRegionModel::preEvolveRegion();

    transferPrimaryRegionThermoFields();

    correctThermoFields();

    transferPrimaryRegionSourceFields();

    // Reset transfer fields
    availableMass_ = mass();
    cloudMassTrans_ == dimensionedScalar(dimMass, 0);
    cloudDiameterTrans_ == dimensionedScalar(dimLength, 0);
    primaryMassTrans_ == dimensionedScalar(dimMass, 0);
}

standardRadiation::standardRadiation
(
    const surfaceFilmModel& owner,
    const dictionary& dict
)
    :
    filmRadiationModel(typeName, owner, dict),
    QrPrimary_
    (
       IOobject
       (
           "Qr", // same name as Qr on primary region to enable mapping
           owner.time().timeName(),
           owner.regionMesh(),
           IOobject::NO_READ,
           IOobject::NO_WRITE
       ),
       owner.regionMesh(),
       dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0),
       owner.mappedPushedFieldPatchTypes<scalar>()
    ),
    QrNet_
    (
       IOobject
       (
           "QrNet",
           owner.time().timeName(),
           owner.regionMesh(),
           IOobject::NO_READ,
           IOobject::NO_WRITE
       ),
       owner.regionMesh(),
       dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0),
       zeroGradientFvPatchScalarField::typeName
    ),
    delta_(owner.delta()),
    deltaMin_(readScalar(coeffs_.lookup("deltaMin"))),
    beta_(readScalar(coeffs_.lookup("beta"))),
    kappaBar_(readScalar(coeffs_.lookup("kappaBar")))
{}

void Foam::LESModels::IDDESDelta::calcDelta()
{
    const volScalarField& hmax = hmax_;
    const fvMesh& mesh = turbulenceModel_.mesh();

    // Wall-normal vectors
    const volVectorField& n = wallDist::New(mesh).n();

    tmp<volScalarField> tfaceToFacenMax
    (
        new volScalarField
        (
            IOobject
            (
                "faceToFaceMax",
                mesh.time().timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh,
            dimensionedScalar("zero", dimLength, 0.0)
        )
    );

    scalarField& faceToFacenMax = tfaceToFacenMax.ref().primitiveFieldRef();

    const cellList& cells = mesh.cells();
    const vectorField& faceCentres = mesh.faceCentres();

    forAll(cells, celli)
    {
        scalar maxDelta = 0.0;
        const labelList& cFaces = cells[celli];
        const vector nci = n[celli];

        forAll(cFaces, cFacei)
        {
            label facei = cFaces[cFacei];
            const point& fci = faceCentres[facei];

            forAll(cFaces, cFacej)
            {
                label facej = cFaces[cFacej];
                const point& fcj = faceCentres[facej];
                scalar ndfc = nci & (fcj - fci);

                if (ndfc > maxDelta)
                {
                    maxDelta = ndfc;
                }
            }

Foam::tmp<Foam::volScalarField> Foam::radiationModels::noRadiation::Rp() const
{
    return volScalarField::New
    (
        "Rp",
        mesh_,
        dimensionedScalar
        (
            constant::physicoChemical::sigma.dimensions()/dimLength,
            0
        )
    );
}

// * * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * * //
void Foam::ChengDiringer::update()
{
  Info<<"Updating Sigma Source terms"<<endl;

  //By the time, no quenching due to excesive strain is incorporated
  //(see the paper by Cheng & Diringer)

  ProdRateForSigma_ = rho_ * alphaSigma_ * turbulence_.epsilon() / 
    (turbulence_.k() + dimensionedScalar("tol", pow(dimVelocity,2), SMALL));

  DestrRateForSigma_ = rho_ * betaSigma_ * Su_ * Sigma_/(b_ + SMALL);

}

Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::
frictionalPressure
(
    const volScalarField& alpha1,
    const dimensionedScalar& alphaMinFriction,
    const dimensionedScalar& alphaMax
) const
{
    return
        dimensionedScalar("1e24", dimensionSet(1, -1, -2, 0, 0), 1e24)	//1e24*(alpha-alphaMinFriction)^10 if alpha1<Minfriction, then it is 0
       *pow(Foam::max(alpha1 - alphaMinFriction, scalar(0)), 10.0);
}

// Update the coefficients associated with the patch field
void smoluchowskiJumpTFvPatchScalarField::updateCoeffs()
{
    if (updated())
    {
        return;
    }

    const fvPatchScalarField& pmu =
        lookupPatchField<volScalarField, scalar>("mu");
    const fvPatchScalarField& prho =
        lookupPatchField<volScalarField, scalar>("rho");
    const fvPatchField<scalar>& ppsi =
        lookupPatchField<volScalarField, scalar>("psi");
    const fvPatchVectorField& pU =
        lookupPatchField<volVectorField, vector>("U");

    // Prandtl number reading consistent with rhoCentralFoam
    const dictionary& thermophysicalProperties =
        db().lookupObject<IOdictionary>("thermophysicalProperties");
    dimensionedScalar Pr = dimensionedScalar("Pr", dimless, 1.0);
    if (thermophysicalProperties.found("Pr"))
    {
        Pr = dimensionedScalar(thermophysicalProperties.lookup("Pr"));
    }

    Field<scalar> C2 = pmu/prho
        *sqrt(ppsi*mathematicalConstant::pi/2.0)
        *2.0*gamma_/Pr.value()/(gamma_ + 1.0)
        *(2.0 - accommodationCoeff_)/accommodationCoeff_;

    Field<scalar> aCoeff = prho.snGrad() - prho/C2;
    Field<scalar> KEbyRho = 0.5*magSqr(pU);

    valueFraction() = (1.0/(1.0 + patch().deltaCoeffs()*C2));
    refValue() = Twall_;
    refGrad() = 0.0;

    mixedFvPatchScalarField::updateCoeffs();
}

Foam::ConstantRateDevolatilisation<CloudType>::ConstantRateDevolatilisation
(
    const dictionary& dict,
    CloudType& owner
)
:
    DevolatilisationModel<CloudType>(dict, owner, typeName),
    A0_(dimensionedScalar(this->coeffDict().lookup("A0")).value()),
    volatileResidualCoeff_
    (
        readScalar(this->coeffDict().lookup("volatileResidualCoeff"))
    )
{}

singleStepCombustion<CombThermoType, ThermoType>::singleStepCombustion
(
    const word& modelType,
    const fvMesh& mesh,
    const word& combustionProperties,
    const word& phaseName
)
:
    CombThermoType(modelType, mesh, phaseName),
    singleMixturePtr_(nullptr),
    wFuel_
    (
        IOobject
        (
            IOobject::groupName("wFuel", phaseName),
            this->mesh().time().timeName(),
            this->mesh(),
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        this->mesh(),
        dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
    ),
    semiImplicit_(readBool(this->coeffs_.lookup("semiImplicit")))
{
    if (isA<singleStepReactingMixture<ThermoType>>(this->thermo()))
    {
        singleMixturePtr_ =
            &dynamic_cast<singleStepReactingMixture<ThermoType>&>
            (
                this->thermo()
            );
    }
    else
    {
        FatalErrorInFunction
            << "Inconsistent thermo package for " << this->type() << " model:\n"
            << "    " << this->thermo().type() << nl << nl
            << "Please select a thermo package based on "
            << "singleStepReactingMixture" << exit(FatalError);
    }

    if (semiImplicit_)
    {
        Info<< "Combustion mode: semi-implicit" << endl;
    }
    else
    {
        Info<< "Combustion mode: explicit" << endl;
    }
}