initial commit

OpenFOAM-v2112/applications/solvers/incompressible/simpleFoam/pEqn.H

+{
+    volScalargpuField rAU(1.0/UEqn.A());
+    volVectorgpuField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+    surfaceScalargpuField phiHbyA("phiHbyA", fvc::flux(HbyA));
+    MRF.makeRelative(phiHbyA);
+    adjustPhi(phiHbyA, U, p);
+
+    tmp<volScalargpuField> rAtU(rAU);
+
+    if (simple.consistent())
+    {
+        rAtU = 1.0/(1.0/rAU - UEqn.H1());
+        phiHbyA +=
+            fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*devicemesh.magSf();
+        HbyA -= (rAU - rAtU())*fvc::grad(p);
+    }
+
+    tUEqn.clear();
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, U, phiHbyA, rAtU(), MRF);
+
+    // Non-orthogonal pressure corrector loop
+    while (simple.correctNonOrthogonal())
+    {
+        gpufvScalarMatrix pEqn
+        (
+            fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
+        );
+
+        pEqn.setReference(pRefCell, pRefValue);
+
+        pEqn.solve();
+
+        if (simple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA - pEqn.flux();
+        }
+    }
+
+    #include "gpucontinuityErrs.H"
+
+    // Explicitly relax pressure for momentum corrector
+    p.relax();
+
+    // Momentum corrector
+    U = HbyA - rAtU()*fvc::grad(p);
+    U.correctBoundaryConditions();
+    fvOptions.correct(U);
+}

OpenFOAM-v2112/applications/solvers/incompressible/simpleFoam/simpleFoam.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    simpleFoam
+
+Group
+    grpIncompressibleSolvers
+
+Description
+    Steady-state solver for incompressible, turbulent flows.
+
+    \heading Solver details
+    The solver uses the SIMPLE algorithm to solve the continuity equation:
+
+        \f[
+            \div \vec{U} = 0
+        \f]
+
+    and momentum equation:
+
+        \f[
+            \div \left( \vec{U} \vec{U} \right) - \div \gvec{R}
+          = - \grad p + \vec{S}_U
+        \f]
+
+    Where:
+    \vartable
+        \vec{U} | Velocity
+        p       | Pressure
+        \vec{R} | Stress tensor
+        \vec{S}_U | Momentum source
+    \endvartable
+
+    \heading Required fields
+    \plaintable
+        U       | Velocity [m/s]
+        p       | Kinematic pressure, p/rho [m2/s2]
+        \<turbulence fields\> | As required by user selection
+    \endplaintable
+
+\*---------------------------------------------------------------------------*/
+
+#include "gpufvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "simpleControl.H"
+#include "gpufvOptions.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Steady-state solver for incompressible, turbulent flows."
+    );
+
+    #include "gpupostProcess.H"
+
+    #include "addCheckCaseOptions.H"
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "gpucreateMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+
+    turbulence->validate();
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (simple.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity SIMPLE corrector
+        {
+            #include "UEqn.H"
+            #include "pEqn.H"
+        }
+
+        laminarTransport.correct();
+        turbulence->correct();
+
+        runTime.write();
+
+        runTime.printExecutionTime(Info);
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

OpenFOAM-v2112/applications/solvers/multiphase/VoF/alphaCourantNo.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    alphaCourantNo
+
+Description
+    Calculates and outputs the mean and maximum Courant Numbers.
+
+\*---------------------------------------------------------------------------*/
+
+scalar maxAlphaCo
+(
+    runTime.controlDict().get<scalar>("maxAlphaCo")
+);
+
+scalar alphaCoNum = 0.0;
+scalar meanAlphaCoNum = 0.0;
+
+if (mesh.nInternalFaces())
+{
+    scalargpuField sumPhi
+    (
+        mixture.nearInterface()().primitiveField()
+       *fvc::surfaceSum(mag(phi))().primitiveField()
+    );
+
+    alphaCoNum = 0.5*gMax(sumPhi/devicemesh.V().field())*runTime.deltaTValue();
+
+    meanAlphaCoNum =
+        0.5*(gSum(sumPhi)/gSum(devicemesh.V().field()))*runTime.deltaTValue();
+}
+
+Info<< "Interface Courant Number mean: " << meanAlphaCoNum
+    << " max: " << alphaCoNum << endl;
+
+// ************************************************************************* //

OpenFOAM-v2112/applications/solvers/multiphase/VoF/alphaEqn.H

+{
+    word alphaScheme("div(phi,alpha)");
+    word alpharScheme("div(phirb,alpha)");
+
+    // Set the off-centering coefficient according to ddt scheme
+    scalar ocCoeff = 0;
+    {
+        tmp<fv::gpuddtScheme<scalar>> tddtAlpha
+        (
+            fv::gpuddtScheme<scalar>::New
+            (
+                devicemesh,
+                mesh.ddtScheme("ddt(alpha)")
+            )
+        );
+        const fv::gpuddtScheme<scalar>& ddtAlpha = tddtAlpha();
+
+        if
+        (
+            isType<fv::gpuEulerDdtScheme<scalar>>(ddtAlpha)
+         || isType<fv::gpulocalEulerDdtScheme<scalar>>(ddtAlpha)
+        )
+        {
+            ocCoeff = 0;
+        }
+        else if (isType<fv::gpuCrankNicolsonDdtScheme<scalar>>(ddtAlpha))
+        {
+            if (nAlphaSubCycles > 1)
+            {
+                FatalErrorInFunction
+                    << "Sub-cycling is not supported "
+                       "with the CrankNicolson ddt scheme"
+                    << exit(FatalError);
+            }
+
+            if
+            (
+                alphaRestart
+             || mesh.time().timeIndex() > mesh.time().startTimeIndex() + 1
+            )
+            {
+                ocCoeff =
+                    refCast<const fv::gpuCrankNicolsonDdtScheme<scalar>>(ddtAlpha)
+                   .ocCoeff();
+            }
+        }
+        else
+        {
+            FatalErrorInFunction
+                << "Only Euler and CrankNicolson ddt schemes are supported"
+                << exit(FatalError);
+        }
+    }
+
+    // Set the time blending factor, 1 for Euler
+    scalar cnCoeff = 1.0/(1.0 + ocCoeff);
+
+    // Standard face-flux compression coefficient
+    surfaceScalargpuField phic(mixture.cAlpha()*mag(phi/devicemesh.magSf()));
+
+    // Add the optional isotropic compression contribution
+    if (icAlpha > 0)
+    {
+        phic *= (1.0 - icAlpha);
+        phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
+    }
+
+    // Add the optional shear compression contribution
+    if (scAlpha > 0)
+    {
+        phic +=
+            scAlpha*mag(devicemesh.delta() & fvc::interpolate(symm(fvc::grad(U))));
+    }
+
+
+    surfaceScalargpuField::Boundary& phicBf =
+        phic.boundaryFieldRef();
+
+    // Do not compress interface at non-coupled boundary faces
+    // (inlets, outlets etc.)
+    forAll(phic.boundaryField(), patchi)
+    {
+        fvsPatchScalargpuField& phicp = phicBf[patchi];
+
+        if (!phicp.coupled())
+        {
+            phicp == 0;
+        }
+    }
+
+    tmp<surfaceScalargpuField> phiCN(phi);
+
+    // Calculate the Crank-Nicolson off-centred volumetric flux
+    if (ocCoeff > 0)
+    {
+        phiCN = cnCoeff*phi + (1.0 - cnCoeff)*phi.oldTime();
+    }
+
+    if (MULESCorr)
+    {
+        #include "alphaSuSp.H"
+
+        gpufvScalarMatrix alpha1Eqn
+        (
+            (
+                LTS
+              ? fv::gpulocalEulerDdtScheme<scalar>(devicemesh).fvmDdt(alpha1)
+              : fv::gpuEulerDdtScheme<scalar>(devicemesh).fvmDdt(alpha1)
+            )
+          + fv::gpugaussConvectionScheme<scalar>
+            (
+                devicemesh,
+                phiCN,
+                gpuupwind<scalar>(devicemesh, phiCN)
+            ).fvmDiv(phiCN, alpha1)
+       // - fvm::Sp(fvc::ddt(dimensionedScalar("1", dimless, 1), mesh)
+       //           + fvc::div(phiCN), alpha1)
+         ==
+            Su + fvm::Sp(Sp + divU, alpha1)
+        );
+
+        alpha1Eqn.solve();
+
+        Info<< "Phase-1 volume fraction = "
+            << alpha1.weightedAverage(devicemesh.Vsc()).value()
+            << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
+            << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
+            << endl;
+
+        tmp<surfaceScalargpuField> talphaPhi1UD(alpha1Eqn.flux());
+        alphaPhi10 = talphaPhi1UD();
+
+        if (alphaApplyPrevCorr && talphaPhi1Corr0.valid())
+        {
+            Info<< "Applying the previous iteration compression flux" << endl;
+            MULES::correct
+            (
+                geometricOneField(),
+                alpha1,
+                alphaPhi10,
+                talphaPhi1Corr0.ref(),
+                oneField(),
+                zeroField()
+            );
+
+            alphaPhi10 += talphaPhi1Corr0();
+        }
+
+        // Cache the upwind-flux
+        talphaPhi1Corr0 = talphaPhi1UD;
+
+        alpha2 = 1.0 - alpha1;
+
+        mixture.correct();
+    }
+
+
+    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
+    {
+        #include "alphaSuSp.H"
+
+        surfaceScalargpuField phir(phic*mixture.nHatf());
+
+        alphaPhiUn =
+            fvc::flux
+            (
+                phi,
+                alpha1,
+                alphaScheme
+            )
+          + fvc::flux
+            (
+               -fvc::flux(-phir, alpha2, alpharScheme),
+                alpha1,
+                alpharScheme
+            );
+
+        if (MULESCorr)
+        {
+            tmp<surfaceScalargpuField> talphaPhi1Corr(alphaPhiUn - alphaPhi10);
+            volScalargpuField alpha10("alpha10", alpha1);
+
+            MULES::correct
+            (
+                geometricOneField(),
+                alpha1,
+                alphaPhiUn,
+                talphaPhi1Corr.ref(),
+                Sp,
+                (-Sp*alpha1)(),
+                oneField(),
+                zeroField()
+            );
+
+            // Under-relax the correction for all but the 1st corrector
+            if (aCorr == 0)
+            {
+                alphaPhi10 += talphaPhi1Corr();
+            }
+            else
+            {
+                alpha1 = 0.5*alpha1 + 0.5*alpha10;
+                alphaPhi10 += 0.5*talphaPhi1Corr();
+            }
+        }
+        else
+        {
+            alphaPhi10 = alphaPhiUn;
+
+            MULES::explicitSolve
+            (
+                geometricOneField(),
+                alpha1,
+                phiCN,
+                alphaPhi10,
+                Sp,
+                (Su + divU*min(alpha1(), scalar(1)))(),
+                oneField(),
+                zeroField()
+            );
+        }
+
+        alpha2 = 1.0 - alpha1;
+
+        mixture.correct();
+    }
+
+    if (alphaApplyPrevCorr && MULESCorr)
+    {
+        talphaPhi1Corr0 = alphaPhi10 - talphaPhi1Corr0;
+        talphaPhi1Corr0.ref().rename("alphaPhi1Corr0");
+    }
+    else
+    {
+        talphaPhi1Corr0.clear();
+    }
+
+    #include "rhofs.H"
+
+    if
+    (
+        word(mesh.ddtScheme("ddt(rho,U)"))
+     == fv::gpuEulerDdtScheme<vector>::typeName
+     || word(mesh.ddtScheme("ddt(rho,U)"))
+     == fv::gpulocalEulerDdtScheme<vector>::typeName
+    )
+    {
+        rhoPhi = alphaPhi10*(rho1f - rho2f) + phiCN*rho2f;
+    }
+    else
+    {
+        if (ocCoeff > 0)
+        {
+            // Calculate the end-of-time-step alpha flux
+            alphaPhi10 =
+                (alphaPhi10 - (1.0 - cnCoeff)*alphaPhi10.oldTime())/cnCoeff;
+        }
+
+        // Calculate the end-of-time-step mass flux
+        rhoPhi = alphaPhi10*(rho1f - rho2f) + phi*rho2f;
+    }
+
+    Info<< "Phase-1 volume fraction = "
+        << alpha1.weightedAverage(devicemesh.Vsc()).value()
+        << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
+        << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
+        << endl;
+}

OpenFOAM-v2112/applications/solvers/multiphase/VoF/alphaEqnSubCycle.H

+if (nAlphaSubCycles > 1)
+{
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+    surfaceScalargpuField rhoPhiSum
+    (
+        IOobject
+        (
+            "rhoPhiSum",
+            runTime.timeName(),
+            mesh
+        ),
+        devicemesh,
+        dimensionedScalar(rhoPhi.dimensions(), Zero)
+    );
+
+    tmp<volScalargpuField> trSubDeltaT;
+
+    if (LTS)
+    {
+        trSubDeltaT =
+            fv::gpulocalEulerDdt::localRSubDeltaT(devicemesh, nAlphaSubCycles);
+    }
+
+    for
+    (
+        subCycle<volScalargpuField> alphaSubCycle(alpha1, nAlphaSubCycles);
+        !(++alphaSubCycle).end();
+    )
+    {
+        #include "alphaEqn.H"
+        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+    }
+
+    rhoPhi = rhoPhiSum;
+}
+else
+{
+    #include "alphaEqn.H"
+}
+
+rho == alpha1*rho1 + alpha2*rho2;

OpenFOAM-v2112/applications/solvers/multiphase/VoF/createAlphaFluxes.H

+IOobject alphaPhi10Header
+(
+    IOobject::groupName("alphaPhi0", alpha1.group()),
+    runTime.timeName(),
+    mesh,
+    IOobject::READ_IF_PRESENT,
+    IOobject::AUTO_WRITE
+);
+
+const bool alphaRestart =
+    alphaPhi10Header.typeHeaderOk<surfaceScalargpuField>(true);
+
+if (alphaRestart)
+{
+    Info << "Restarting alpha" << endl;
+}
+
+// MULES flux from previous time-step
+surfaceScalargpuField alphaPhi10
+(
+    alphaPhi10Header,
+    phi*fvc::interpolate(alpha1)
+);
+
+// MULES Correction
+tmp<surfaceScalargpuField> talphaPhi1Corr0;
+
+// MULES compressed flux is registered in case scalarTransport FO needs it.
+surfaceScalargpuField alphaPhiUn
+(
+    IOobject
+    (
+        "alphaPhiUn",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    devicemesh,
+    dimensionedScalar(phi.dimensions(), Zero)
+);

OpenFOAM-v2112/applications/solvers/multiphase/VoF/setDeltaT.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    setDeltaT
+
+Description
+    Reset the timestep to maintain a constant maximum courant Number.
+    Reduction of time-step is immediate, but increase is damped to avoid
+    unstable oscillations.
+
+\*---------------------------------------------------------------------------*/
+
+if (adjustTimeStep)
+{
+    scalar maxDeltaTFact =
+        min(maxCo/(CoNum + SMALL), maxAlphaCo/(alphaCoNum + SMALL));
+
+    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+
+    runTime.setDeltaT
+    (
+        min
+        (
+            deltaTFact*runTime.deltaTValue(),
+            maxDeltaT
+        )
+    );
+
+    Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
+}
+
+// ************************************************************************* //

OpenFOAM-v2112/applications/solvers/multiphase/VoF/setRDeltaT.H

+{
+    volScalargpuField& rDeltaT = trDeltaT.ref();
+
+    const dictionary& pimpleDict = pimple.dict();
+
+    scalar maxCo
+    (
+        pimpleDict.getOrDefault<scalar>("maxCo", 0.9)
+    );
+
+    scalar maxAlphaCo
+    (
+        pimpleDict.getOrDefault<scalar>("maxAlphaCo", 0.2)
+    );
+
+    scalar rDeltaTSmoothingCoeff
+    (
+        pimpleDict.getOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
+    );
+
+    label nAlphaSpreadIter
+    (
+        pimpleDict.getOrDefault<label>("nAlphaSpreadIter", 1)
+    );
+
+    scalar alphaSpreadDiff
+    (
+        pimpleDict.getOrDefault<scalar>("alphaSpreadDiff", 0.2)
+    );
+
+    scalar alphaSpreadMax
+    (
+        pimpleDict.getOrDefault<scalar>("alphaSpreadMax", 0.99)
+    );
+
+    scalar alphaSpreadMin
+    (
+        pimpleDict.getOrDefault<scalar>("alphaSpreadMin", 0.01)
+    );
+
+    label nAlphaSweepIter
+    (
+        pimpleDict.getOrDefault<label>("nAlphaSweepIter", 5)
+    );
+
+    scalar rDeltaTDampingCoeff
+    (
+        pimpleDict.getOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
+    );
+
+    scalar maxDeltaT
+    (
+        pimpleDict.getOrDefault<scalar>("maxDeltaT", GREAT)
+    );
+
+    volScalargpuField rDeltaT0("rDeltaT0", rDeltaT);
+
+    // Set the reciprocal time-step from the local Courant number
+    rDeltaT.ref() = max
+    (
+        1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
+        fvc::surfaceSum(mag(rhoPhi))()()
+       /((2*maxCo)*devicemesh.V()*rho())
+    );
+
+    if (maxAlphaCo < maxCo)
+    {
+        // Further limit the reciprocal time-step
+        // in the vicinity of the interface
+
+        volScalargpuField alpha1Bar(fvc::average(alpha1));
+
+        rDeltaT.ref() = max
+        (
+            rDeltaT(),
+            pos0(alpha1Bar() - alphaSpreadMin)
+           *pos0(alphaSpreadMax - alpha1Bar())
+           *fvc::surfaceSum(mag(phi))()()
+           /((2*maxAlphaCo)*devicemesh.V())
+        );
+    }
+
+    // Update tho boundary values of the reciprocal time-step
+    rDeltaT.correctBoundaryConditions();
+
+    Info<< "Flow time scale min/max = "
+        << gMin(1/rDeltaT.primitiveField())
+        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+
+    if (rDeltaTSmoothingCoeff < 1.0)
+    {
+        fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
+    }
+
+    if (nAlphaSpreadIter > 0)
+    {
+        fvc::spread
+        (
+            rDeltaT,
+            alpha1,
+            nAlphaSpreadIter,
+            alphaSpreadDiff,
+            alphaSpreadMax,
+            alphaSpreadMin
+        );
+    }
+
+    if (nAlphaSweepIter > 0)
+    {
+        fvc::sweep(rDeltaT, alpha1, nAlphaSweepIter, alphaSpreadDiff);
+    }
+
+    Info<< "Smoothed flow time scale min/max = "
+        << gMin(1/rDeltaT.primitiveField())
+        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+
+    // Limit rate of change of time scale
+    // - reduce as much as required
+    // - only increase at a fraction of old time scale
+    if
+    (
+        rDeltaTDampingCoeff < 1.0
+     && runTime.timeIndex() > runTime.startTimeIndex() + 1
+    )
+    {
+        rDeltaT = max
+        (
+            rDeltaT,
+            (scalar(1) - rDeltaTDampingCoeff)*rDeltaT0
+        );
+
+        Info<< "Damped flow time scale min/max = "
+            << gMin(1/rDeltaT.primitiveField())
+            << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+    }
+}

OpenFOAM-v2112/applications/solvers/multiphase/interFoam/Make/files

+interFoam.C
+
+EXE = $(FOAM_APPBIN)/interFoam

OpenFOAM-v2112/applications/solvers/multiphase/interFoam/Make/options

+EXE_INC = \
+    -I../VoF \
+    -I$(LIB_SRC)/phaseSystemModels-GPU/twoPhaseInter/incompressibleInterPhaseTransportModel/lnInclude \
+    -I$(LIB_SRC)/phaseSystemModels-GPU/twoPhaseInter/VoFphaseIncompressibleTurbulenceModels/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/gpuOpenFOAM/lnInclude \
+    -I$(LIB_SRC)/gpufiniteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/transportModels-GPU \
+    -I$(LIB_SRC)/transportModels-GPU/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels-GPU/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/transportModels-GPU/twoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels-GPU/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels-GPU/incompressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels-GPU/phaseIncompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels-GPU/immiscibleIncompressibleTwoPhaseMixture/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lgpuOpenFOAM \
+    -lgpufiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lsampling \
+    -ldynamicFvMesh \
+    -lincompressibleTransportModels-GPU \
+    -linterfaceProperties-GPU \
+    -limmiscibleIncompressibleTwoPhaseMixture-GPU \
+    -lturbulenceModels-GPU \
+    -lincompressibleTurbulenceModels-GPU \
+    -lVoFphaseTurbulentTransportModels-GPU \
+    -lincompressibleInterPhaseTransportModels-GPU

OpenFOAM-v2112/applications/solvers/multiphase/interFoam/UEqn.H

+    MRF.correctBoundaryVelocity(U);
+
+    gpufvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
+      + MRF.DDt(rho, U)
+      + turbulence->divDevRhoReff(rho, U)
+     ==
+        fvOptions(rho, U)
+    );
+
+    UEqn.relax();
+
+    fvOptions.constrain(UEqn);
+
+    if (pimple.momentumPredictor())
+    {
+        solve
+        (
+            UEqn
+         ==
+            fvc::reconstruct
+            (
+                (
+                    mixture.surfaceTensionForce()
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                ) * devicemesh.magSf()
+            )
+        );
+
+        fvOptions.correct(U);
+    }

OpenFOAM-v2112/applications/solvers/multiphase/interFoam/alphaSuSp.H

+zeroField Su;
+zeroField Sp;
+zeroField divU;

OpenFOAM-v2112/applications/solvers/multiphase/interFoam/correctPhi.H

+CorrectPhi
+(
+    U,
+    phi,
+    p_rgh,
+    surfaceScalargpuField("rAUf", fvc::interpolate(rAU())),
+    geometricZeroField(),
+    pimple
+);
+
+#include "gpucontinuityErrs.H"

OpenFOAM-v2112/applications/solvers/multiphase/interFoam/createFields.H

+#include "gpucreateRDeltaT.H"
+
+Info<< "Reading field p_rgh\n" << endl;
+volScalargpuField p_rgh
+(
+    IOobject
+    (
+        "p_rgh",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    devicemesh
+);
+
+Info<< "Reading field U\n" << endl;
+volVectorgpuField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    devicemesh
+);
+
+#include "gpucreatePhi.H"
+
+
+Info<< "Reading transportProperties\n" << endl;
+immiscibleIncompressibleTwoPhaseMixture mixture(U, phi);
+
+volScalargpuField& alpha1(mixture.alpha1());
+volScalargpuField& alpha2(mixture.alpha2());
+
+const dimensionedScalar& rho1 = mixture.rho1();
+const dimensionedScalar& rho2 = mixture.rho2();
+
+
+// Need to store rho for ddt(rho, U)
+volScalargpuField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT
+    ),
+    alpha1*rho1 + alpha2*rho2
+);
+rho.oldTime();
+
+
+// Mass flux
+surfaceScalargpuField rhoPhi
+(
+    IOobject
+    (
+        "rhoPhi",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    fvc::interpolate(rho)*phi
+);
+
+typedef incompressibleInterPhaseTransportModel
+    <
+        immiscibleIncompressibleTwoPhaseMixture
+    > transportModelType;
+
+autoPtr<transportModelType> turbulence
+(
+    new transportModelType
+    (
+        rho, U, phi, rhoPhi, mixture
+    )
+);
+
+#include "readGravitationalAcceleration.H"
+#include "readhRef.H"
+#include "gpugh.H"
+
+
+volScalargpuField p
+(
+    IOobject
+    (
+        "p",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    p_rgh + rho*gh
+);
+
+label pRefCell = 0;
+scalar pRefValue = 0.0;
+setRefCell
+(
+    p,
+    p_rgh,
+    pimple.dict(),
+    pRefCell,
+    pRefValue
+);
+
+if (p_rgh.needReference())
+{
+    p += dimensionedScalar
+    (
+        "p",
+        p.dimensions(),
+        pRefValue - getRefCellValue(p, pRefCell)
+    );
+    p_rgh = p - rho*gh;
+}
+
+mesh.setFluxRequired(p_rgh.name());
+mesh.setFluxRequired(alpha1.name());
+
+#include "gpucreateMRF.H"
+#include "gpucreateFvOptions.H"

OpenFOAM-v2112/applications/solvers/multiphase/interFoam/initCorrectPhi.H

+tmp<volScalargpuField> rAU;
+
+if (correctPhi)
+{
+    rAU = new volScalargpuField
+    (
+        IOobject
+        (
+            "rAU",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        devicemesh,
+        dimensionedScalar("rAU", dimTime/dimDensity, 1)
+    );
+
+    #include "correctPhi.H"
+}
+else
+{
+    CorrectPhi
+    (
+        U,
+        phi,
+        p_rgh,
+        dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
+        geometricZeroField(),
+        pimple
+    );
+
+    #include "gpucontinuityErrs.H"
+}

OpenFOAM-v2112/applications/solvers/multiphase/interFoam/interFoam.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    interFoam
+
+Group
+    grpMultiphaseSolvers
+
+Description
+    Solver for two incompressible, isothermal immiscible fluids using a VOF
+    (volume of fluid) phase-fraction based interface capturing approach,
+    with optional mesh motion and mesh topology changes including adaptive
+    re-meshing.
+
+\*---------------------------------------------------------------------------*/
+
+#include "gpufvCFD.H"
+#include "dynamicFvMesh.H"
+#include "gpuCMULES.H"
+#include "gpuEulerDdtScheme.H"
+#include "gpulocalEulerDdtScheme.H"
+#include "gpuCrankNicolsonDdtScheme.H"
+#include "subCycle.H"
+#include "immiscibleIncompressibleTwoPhaseMixture.H"
+#include "incompressibleInterPhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "pimpleControl.H"
+#include "gpufvOptions.H"
+#include "gpuCorrectPhi.H"
+#include "gpufvcSmooth.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Solver for two incompressible, isothermal immiscible fluids"
+        " using VOF phase-fraction based interface capturing.\n"
+        "With optional mesh motion and mesh topology changes including"
+        " adaptive re-meshing."
+    );
+
+    #include "gpupostProcess.H"
+
+    #include "addCheckCaseOptions.H"
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+	#include "gpucreateMesh.H"
+    #include "initContinuityErrs.H"
+    #include "createDyMControls.H"
+    #include "createFields.H"
+    #include "createAlphaFluxes.H"
+    #include "initCorrectPhi.H"
+    #include "gpucreateUfIfPresent.H"
+
+    if (!LTS)
+    {
+        #include "gpuCourantNo.H"
+        #include "setInitialDeltaT.H"
+    }
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readDyMControls.H"
+
+        if (LTS)
+        {
+            #include "setRDeltaT.H"
+        }
+        else
+        {
+            #include "gpuCourantNo.H"
+            #include "alphaCourantNo.H"
+            #include "setDeltaT.H"
+        }
+
+        ++runTime;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            if (pimple.firstIter() || moveMeshOuterCorrectors)
+            {
+                mesh.update();
+
+                if (mesh.changing())
+                {
+                	devicemesh.updateGpuGeom();
+                    // Do not apply previous time-step mesh compression flux
+                    // if the mesh topology changed
+                    if (mesh.topoChanging())
+                    {
+                        talphaPhi1Corr0.clear();
+                    }
+
+                    gh = (g & devicemesh.C()) - ghRef;
+                    ghf = (g & devicemesh.Cf()) - ghRef;
+
+                    MRF.update();
+
+                    if (correctPhi)
+                    {
+                        // Calculate absolute flux
+                        // from the mapped surface velocity
+                        phi = devicemesh.Sf() & Uf();
+
+                        #include "correctPhi.H"
+
+                        // Make the flux relative to the mesh motion
+                        fvc::makeRelative(phi, U);
+
+                        mixture.correct();
+                    }
+
+                    if (checkMeshCourantNo)
+                    {
+                        #include "gpumeshCourantNo.H"
+                    }
+                }
+            }
+
+            #include "alphaControls.H"
+            #include "alphaEqnSubCycle.H"
+
+            mixture.correct();
+
+            if (pimple.frozenFlow())
+            {
+                continue;
+            }
+
+            #include "UEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        runTime.write();
+
+        runTime.printExecutionTime(Info);
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

OpenFOAM-v2112/applications/solvers/multiphase/interFoam/pEqn.H

+{
+    if (correctPhi)
+    {
+        rAU.ref() = 1.0/UEqn.A();
+    }
+    else
+    {
+        rAU = 1.0/UEqn.A();
+    }
+
+    surfaceScalargpuField rAUf("rAUf", fvc::interpolate(rAU()));
+    volVectorgpuField HbyA(constrainHbyA(rAU()*UEqn.H(), U, p_rgh));
+    surfaceScalargpuField phiHbyA
+    (
+        "phiHbyA",
+        fvc::flux(HbyA)
+      + MRF.zeroFilter(fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf))
+    );
+    MRF.makeRelative(phiHbyA);
+
+    if (p_rgh.needReference())
+    {
+        fvc::makeRelative(phiHbyA, U);
+        adjustPhi(phiHbyA, U, p_rgh);
+        fvc::makeAbsolute(phiHbyA, U);
+    }
+
+    surfaceScalargpuField phig
+    (
+        (
+            mixture.surfaceTensionForce()
+          - ghf*fvc::snGrad(rho)
+        )*rAUf*devicemesh.magSf()
+    );
+
+    phiHbyA += phig;
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        gpufvScalarMatrix p_rghEqn
+        (
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+        );
+
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA - p_rghEqn.flux();
+
+            p_rgh.relax();
+
+            U = HbyA + rAU()*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
+            fvOptions.correct(U);
+        }
+    }
+
+    #include "gpucontinuityErrs.H"
+
+    // Correct Uf if the mesh is moving
+    fvc::correctUf(Uf, U, phi);
+
+    // Make the fluxes relative to the mesh motion
+    fvc::makeRelative(phi, U);
+
+    p == p_rgh + rho*gh;
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }
+
+    if (!correctPhi)
+    {
+        rAU.clear();
+    }
+}

OpenFOAM-v2112/applications/solvers/multiphase/interFoam/rhofs.H

+const dimensionedScalar& rho1f(rho1);
+const dimensionedScalar& rho2f(rho2);

OpenFOAM-v2112/applications/solvers/multiphase/multiphaseSolversDoc.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2015 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it under
+    the terms of the GNU General Public License as published by the Free
+    Software Foundation, either version 3 of the License, or (at your option)
+    any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT ANY
+    WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+    FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+    details.
+
+    You should have received a copy of the GNU General Public License along with
+    OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\defgroup grpMultiphaseSolvers Multiphase solvers
+@{
+    \ingroup grpSolvers
+    This group contains multiphase solvers.
+@}
+
+\*---------------------------------------------------------------------------*/

OpenFOAM-v2112/applications/solvers/rhoReactingCentralFoam/Allwclean

+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+set -x
+
+wclean libso BCs
+wclean
+wclean rhoCentralDyMFoam
+
+# ----------------------------------------------------------------- end-of-file