solvers::incompressibleFluid: New solver module for incompressible fluid flow

executed with foamRun for single region simulations of foamMultiRun for
multi-region simulations.  Replaces pimpleFoam, pisoFoam and simpleFoam and all
the corresponding tutorials have been updated and moved to
tutorials/modules/incompressibleFluid.

Class
    Foam::solvers::incompressibleFluid

Description
    Solver module for steady or transient turbulent flow of incompressible
    isothermal fluids with optional mesh motion and change.

    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
    pseudo-transient and steady simulations.

    Optional fvModels and fvConstraints are provided to enhance the simulation
    in many ways including adding various sources, constraining or limiting
    the solution.

    Reference:
    \verbatim
        Greenshields, C. J., & Weller, H. G. (2022).
        Notes on Computational Fluid Dynamics: General Principles.
        CFD Direct Ltd.: Reading, UK.
    \endverbatim

SourceFiles
    incompressibleFluid.C

See also
    Foam::solvers::fluidSolver
    Foam::solvers::isothermalFluid

applications/solvers/incompressible/SRFPimpleFoam/SRFPimpleFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/incompressible/SRFSimpleFoam/Make/options

@@ -1,5 +1,4 @@
 EXE_INC = \
-    -I.. \
     -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \
@@ -7,7 +6,6 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude
 
-
 EXE_LIBS = \
     -lmomentumTransportModels \
     -lincompressibleMomentumTransportModels \

applications/solvers/incompressible/SRFSimpleFoam/SRFSimpleFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/incompressible/pimpleFoam/Make/files

@@ -1,3 +0,0 @@
-pimpleFoam.C
-
-EXE = $(FOAM_APPBIN)/pimpleFoam

applications/solvers/incompressible/pimpleFoam/UEqn.H

@@ -1,24 +0,0 @@
-// Solve the Momentum equation
-
-MRF.correctBoundaryVelocity(U);
-
-tmp<fvVectorMatrix> tUEqn
-(
-    fvm::ddt(U) + fvm::div(phi, U)
-  + MRF.DDt(U)
-  + turbulence->divDevSigma(U)
- ==
-    fvModels.source(U)
-);
-fvVectorMatrix& UEqn = tUEqn.ref();
-
-UEqn.relax();
-
-fvConstraints.constrain(UEqn);
-
-if (pimple.momentumPredictor())
-{
-    solve(UEqn == -fvc::grad(p));
-
-    fvConstraints.constrain(U);
-}

applications/solvers/incompressible/pimpleFoam/correctPhi.H

@@ -1,21 +0,0 @@
-// Calculate absolute flux
-// from the mapped surface velocity
-phi = mesh.Sf() & Uf();
-
-correctUphiBCs(U, phi, true);
-
-CorrectPhi
-(
-    phi,
-    U,
-    p,
-    dimensionedScalar("rAUf", dimTime, 1),
-    geometricZeroField(),
-    pressureReference,
-    pimple
-);
-
-#include "continuityErrs.H"
-
-// Make the flux relative to the mesh motion
-fvc::makeRelative(phi, U);

applications/solvers/incompressible/pimpleFoam/createFields.H

@@ -1,47 +0,0 @@
-#include "createRDeltaT.H"
-
-Info<< "Reading field p\n" << endl;
-volScalarField p
-(
-    IOobject
-    (
-        "p",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-Info<< "Reading field U\n" << endl;
-volVectorField U
-(
-    IOobject
-    (
-        "U",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-#include "createPhi.H"
-
-pressureReference pressureReference(p, pimple.dict());
-
-mesh.schemes().setFluxRequired(p.name());
-
-
-autoPtr<viscosityModel> viscosity(viscosityModel::New(mesh));
-
-autoPtr<incompressible::momentumTransportModel> turbulence
-(
-    incompressible::momentumTransportModel::New(U, phi, viscosity)
-);
-
-#include "createMRF.H"
-#include "createFvModels.H"
-#include "createFvConstraints.H"

applications/solvers/incompressible/pimpleFoam/pEqn.H

@@ -1,72 +0,0 @@
-volScalarField rAU(1.0/UEqn.A());
-volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-surfaceScalarField phiHbyA
-(
-    "phiHbyA",
-    fvc::flux(HbyA)
-  + MRF.zeroFilter(fvc::interpolate(rAU)*fvc::ddtCorr(U, phi, Uf))
-);
-
-MRF.makeRelative(phiHbyA);
-
-if (p.needReference())
-{
-    fvc::makeRelative(phiHbyA, U);
-    adjustPhi(phiHbyA, U, p);
-    fvc::makeAbsolute(phiHbyA, U);
-}
-
-tmp<volScalarField> rAtU(rAU);
-
-if (pimple.consistent())
-{
-    rAtU = 1.0/max(1.0/rAU - UEqn.H1(), 0.1/rAU);
-    phiHbyA +=
-        fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
-    HbyA -= (rAU - rAtU())*fvc::grad(p);
-}
-
-if (pimple.nCorrPiso() <= 1)
-{
-    tUEqn.clear();
-}
-
-// Update the pressure BCs to ensure flux consistency
-constrainPressure(p, U, phiHbyA, rAtU(), MRF);
-
-// Non-orthogonal pressure corrector loop
-while (pimple.correctNonOrthogonal())
-{
-    fvScalarMatrix pEqn
-    (
-        fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
-    );
-
-    pEqn.setReference
-    (
-        pressureReference.refCell(),
-        pressureReference.refValue()
-    );
-
-    pEqn.solve();
-
-    if (pimple.finalNonOrthogonalIter())
-    {
-        phi = phiHbyA - pEqn.flux();
-    }
-}
-
-#include "continuityErrs.H"
-
-// Explicitly relax pressure for momentum corrector
-p.relax();
-
-U = HbyA - rAtU*fvc::grad(p);
-U.correctBoundaryConditions();
-fvConstraints.constrain(U);
-
-// Correct Uf if the mesh is moving
-fvc::correctUf(Uf, U, phi, MRF);
-
-// Make the fluxes relative to the mesh motion
-fvc::makeRelative(phi, U);

applications/solvers/incompressible/pimpleFoam/pimpleFoam.C

@@ -1,149 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-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
-    pimpleFoam
-
-Description
-    Transient solver for incompressible, turbulent flow of Newtonian fluids,
-    with optional mesh motion and mesh topology changes.
-
-    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "viscosityModel.H"
-#include "incompressibleMomentumTransportModels.H"
-#include "pimpleControl.H"
-#include "pressureReference.H"
-#include "CorrectPhi.H"
-#include "fvModels.H"
-#include "fvConstraints.H"
-#include "localEulerDdtScheme.H"
-#include "fvcSmooth.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #include "postProcess.H"
-
-    #include "setRootCaseLists.H"
-    #include "createTime.H"
-    #include "createMesh.H"
-    #include "initContinuityErrs.H"
-    #include "createDyMControls.H"
-    #include "createFields.H"
-    #include "createUfIfPresent.H"
-
-    turbulence->validate();
-
-    if (!LTS)
-    {
-        #include "CourantNo.H"
-        #include "setInitialDeltaT.H"
-    }
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (pimple.run(runTime))
-    {
-        #include "readDyMControls.H"
-
-        if (LTS)
-        {
-            #include "setRDeltaT.H"
-        }
-        else
-        {
-            #include "CourantNo.H"
-            #include "setDeltaT.H"
-        }
-
-        fvModels.preUpdateMesh();
-
-        // Update the mesh for topology change, mesh to mesh mapping
-        mesh.update();
-
-        runTime++;
-
-        Info<< "Time = " << runTime.userTimeName() << nl << endl;
-
-        // --- Pressure-velocity PIMPLE corrector loop
-        while (pimple.loop())
-        {
-            if (pimple.firstPimpleIter() || pimple.moveMeshOuterCorrectors())
-            {
-                // Move the mesh
-                mesh.move();
-
-                if (mesh.changing())
-                {
-                    MRF.update();
-
-                    if (correctPhi)
-                    {
-                        #include "correctPhi.H"
-                    }
-
-                    if (checkMeshCourantNo)
-                    {
-                        #include "meshCourantNo.H"
-                    }
-                }
-            }
-
-            fvModels.correct();
-
-            #include "UEqn.H"
-
-            // --- Pressure corrector loop
-            while (pimple.correct())
-            {
-                #include "pEqn.H"
-            }
-
-            if (pimple.turbCorr())
-            {
-                viscosity->correct();
-                turbulence->correct();
-            }
-        }
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/incompressible/pisoFoam/Make/files

@@ -1,3 +0,0 @@
-pisoFoam.C
-
-EXE = $(FOAM_APPBIN)/pisoFoam

applications/solvers/incompressible/pisoFoam/UEqn.H

@@ -1,23 +0,0 @@
-// Solve the Momentum equation
-
-MRF.correctBoundaryVelocity(U);
-
-fvVectorMatrix UEqn
-(
-    fvm::ddt(U) + fvm::div(phi, U)
-  + MRF.DDt(U)
-  + turbulence->divDevSigma(U)
- ==
-    fvModels.source(U)
-);
-
-UEqn.relax();
-
-fvConstraints.constrain(UEqn);
-
-if (piso.momentumPredictor())
-{
-    solve(UEqn == -fvc::grad(p));
-
-    fvConstraints.constrain(U);
-}

applications/solvers/incompressible/pisoFoam/pEqn.H

@@ -1,45 +0,0 @@
-volScalarField rAU(1.0/UEqn.A());
-volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-surfaceScalarField phiHbyA
-(
-    "phiHbyA",
-    fvc::flux(HbyA)
-  + MRF.zeroFilter(fvc::interpolate(rAU)*fvc::ddtCorr(U, phi))
-);
-
-MRF.makeRelative(phiHbyA);
-
-adjustPhi(phiHbyA, U, p);
-
-// Update the pressure BCs to ensure flux consistency
-constrainPressure(p, U, phiHbyA, rAU, MRF);
-
-// Non-orthogonal pressure corrector loop
-while (piso.correctNonOrthogonal())
-{
-    // Pressure corrector
-
-    fvScalarMatrix pEqn
-    (
-        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
-    );
-
-    pEqn.setReference
-    (
-        pressureReference.refCell(),
-        pressureReference.refValue()
-    );
-
-    pEqn.solve();
-
-    if (piso.finalNonOrthogonalIter())
-    {
-        phi = phiHbyA - pEqn.flux();
-    }
-}
-
-#include "continuityErrs.H"
-
-U = HbyA - rAU*fvc::grad(p);
-U.correctBoundaryConditions();
-fvConstraints.constrain(U);

applications/solvers/incompressible/pisoFoam/pisoFoam.C

@@ -1,99 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-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
-    pisoFoam
-
-Description
-    Transient solver for incompressible, turbulent flow, using the PISO
-    algorithm.
-
-    Sub-models include:
-    - turbulence modelling, i.e. laminar, RAS or LES
-    - run-time selectable MRF and finite volume options, e.g. explicit porosity
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "viscosityModel.H"
-#include "incompressibleMomentumTransportModels.H"
-#include "pisoControl.H"
-#include "pressureReference.H"
-#include "fvModels.H"
-#include "fvConstraints.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #include "postProcess.H"
-
-    #include "setRootCaseLists.H"
-    #include "createTime.H"
-    #include "createMesh.H"
-    #include "createControl.H"
-    #include "createFields.H"
-    #include "initContinuityErrs.H"
-
-    turbulence->validate();
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (runTime.loop())
-    {
-        Info<< "Time = " << runTime.userTimeName() << nl << endl;
-
-        #include "CourantNo.H"
-
-        // Pressure-velocity PISO corrector
-        {
-            fvModels.correct();
-
-            #include "UEqn.H"
-
-            // --- PISO loop
-            while (piso.correct())
-            {
-                #include "pEqn.H"
-            }
-        }
-
-        viscosity->correct();
-        turbulence->correct();
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/incompressible/porousSimpleFoam/Make/options

@@ -1,5 +1,4 @@
 EXE_INC = \
-    -I.. \
     -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \
@@ -7,7 +6,6 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude
 
-
 EXE_LIBS = \
     -lmomentumTransportModels \
     -lincompressibleMomentumTransportModels \

applications/solvers/incompressible/porousSimpleFoam/createFields.H

@@ -28,12 +28,10 @@ volVectorField U
 
 #include "createPhi.H"
 
-
+pressureReference pressureReference(p, simple.dict());
-pressureReference pressureReference(p, piso.dict());
 
 mesh.schemes().setFluxRequired(p.name());
 
-
 autoPtr<viscosityModel> viscosity(viscosityModel::New(mesh));
 
 autoPtr<incompressible::momentumTransportModel> turbulence

applications/solvers/incompressible/porousSimpleFoam/porousSimpleFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/modules/fluid/Allwmake

@@ -4,6 +4,8 @@ cd ${0%/*} || exit 1    # Run from this directory
 # Parse arguments for library compilation
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
 
+wmake $targetType fluidSolver
+wmake $targetType incompressibleFluid
 wmake $targetType isothermalFluid
 wmake $targetType fluid
 wmake $targetType multicomponentFluid

applications/solvers/modules/fluid/fluid/Make/options

@@ -1,5 +1,6 @@
 EXE_INC = \
-    -I../isothermalFluid/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/fluid/fluidSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/fluid/isothermalFluid/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
@@ -11,6 +12,7 @@ EXE_INC = \
     -I$(LIB_SRC)/sampling/lnInclude
 
 LIB_LIBS = \
+    -lfluidSolver \
     -lisothermalFluid \
     -lfluidThermophysicalModels \
     -lmomentumTransportModels \

applications/solvers/modules/fluid/fluid/fluid.H

@@ -27,11 +27,15 @@ Class
 Description
     Solver module for steady or transient turbulent flow of compressible fluids
     with heat-transfer for HVAC and similar applications, with optional
-    mesh motion and mesh topology changes.
+    mesh motion and change.
 
     Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
     pseudo-transient and steady simulations.
 
+    Optional fvModels and fvConstraints are provided to enhance the simulation
+    in many ways including adding various sources, Lagrangian particles,
+    radiation, surface film etc. and constraining or limiting the solution.
+
     Reference:
     \verbatim
         Greenshields, C. J., & Weller, H. G. (2022).

applications/solvers/modules/fluid/fluidSolver/Make/files

@@ -0,0 +1,3 @@
+fluidSolver.C
+
+LIB = $(FOAM_LIBBIN)/libfluidSolver

applications/solvers/modules/fluid/fluidSolver/Make/options

@@ -1,12 +1,11 @@
-EXE_INC = \
+EXE_INC = -ggdb3 \
     -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
-EXE_LIBS = \
+LIB_LIBS = \
     -lmomentumTransportModels \
     -lincompressibleMomentumTransportModels \
     -lphysicalProperties \

applications/solvers/modules/fluid/fluidSolver/fluidSolver.C

@@ -0,0 +1,220 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fluidSolver.H"
+#include "surfaceFields.H"
+#include "fvcDiv.H"
+#include "fvcSurfaceIntegrate.H"
+#include "fvcVolumeIntegrate.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solvers
+{
+    defineTypeNameAndDebug(fluidSolver, 0);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::solvers::fluidSolver::read()
+{
+    maxCo =
+        runTime.controlDict().lookupOrDefault<scalar>("maxCo", 1.0);
+
+    maxDeltaT_ =
+        runTime.controlDict().lookupOrDefault<scalar>("maxDeltaT", great);
+
+    correctPhi = pimple.dict().lookupOrDefault
+    (
+        "correctPhi",
+        mesh.dynamic()
+    );
+
+    checkMeshCourantNo = pimple.dict().lookupOrDefault
+    (
+        "checkMeshCourantNo",
+        false
+    );
+}
+
+
+void Foam::solvers::fluidSolver::meshCourantNo() const
+{
+    if (checkMeshCourantNo)
+    {
+        const scalarField sumPhi
+        (
+            fvc::surfaceSum(mag(mesh.phi()))().primitiveField()
+        );
+
+        const scalar meshCoNum
+        (
+            0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue()
+        );
+
+        const scalar meanMeshCoNum
+        (
+            0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue()
+        );
+
+        Info<< "Mesh Courant Number mean: " << meanMeshCoNum
+            << " max: " << meshCoNum << endl;
+    }
+}
+
+
+template<class RhoType>
+void Foam::solvers::fluidSolver::correctCoNum
+(
+    const RhoType& rho,
+    const surfaceScalarField& phi
+)
+{
+    const scalarField sumPhi
+    (
+        fvc::surfaceSum(mag(phi))().primitiveField()/rho.primitiveField()
+    );
+
+    CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+    const scalar meanCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+
+    Info<< "Courant Number mean: " << meanCoNum
+        << " max: " << CoNum << endl;
+}
+
+
+void Foam::solvers::fluidSolver::correctCoNum(const surfaceScalarField& phi)
+{
+    correctCoNum(geometricOneField(), phi);
+}
+
+
+void Foam::solvers::fluidSolver::correctCoNum
+(
+    const volScalarField& rho,
+    const surfaceScalarField& phi
+)
+{
+    correctCoNum<volScalarField>(rho, phi);
+}
+
+
+void Foam::solvers::fluidSolver::continuityErrors
+(
+    const surfaceScalarField& phi
+)
+{
+    const volScalarField contErr(fvc::div(phi));
+
+    const scalar sumLocalContErr =
+        runTime.deltaTValue()
+       *mag(contErr)().weightedAverage(mesh.V()).value();
+
+    const scalar globalContErr =
+        runTime.deltaTValue()
+       *contErr.weightedAverage(mesh.V()).value();
+
+    cumulativeContErr += globalContErr;
+
+    Info<< "time step continuity errors : sum local = " << sumLocalContErr
+        << ", global = " << globalContErr
+        << ", cumulative = " << cumulativeContErr
+        << endl;
+}
+
+
+void Foam::solvers::fluidSolver::continuityErrors
+(
+    const volScalarField& rho,
+    const volScalarField& thermoRho,
+    const surfaceScalarField& phi
+)
+{
+    if (mesh.schemes().steady())
+    {
+        continuityErrors(phi);
+    }
+    else
+    {
+        const dimensionedScalar totalMass = fvc::domainIntegrate(rho);
+
+        const scalar sumLocalContErr =
+            (fvc::domainIntegrate(mag(rho - thermoRho))/totalMass).value();
+
+        const scalar globalContErr =
+            (fvc::domainIntegrate(rho - thermoRho)/totalMass).value();
+
+        cumulativeContErr += globalContErr;
+
+        Info<< "time step continuity errors : sum local = " << sumLocalContErr
+            << ", global = " << globalContErr
+            << ", cumulative = " << cumulativeContErr
+            << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::solvers::fluidSolver::fluidSolver(fvMesh& mesh)
+:
+    solver(mesh),
+    cumulativeContErr(0),
+    CoNum(0)
+{
+    // Read the controls
+    read();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::solvers::fluidSolver::~fluidSolver()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::scalar Foam::solvers::fluidSolver::maxDeltaT() const
+{
+    if (CoNum > small)
+    {
+        const scalar deltaT = maxCo*runTime.deltaTValue()/CoNum;
+        return min(min(deltaT, fvModels().maxDeltaT()), maxDeltaT_);
+    }
+    else
+    {
+        return runTime.deltaTValue();
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/modules/fluid/fluidSolver/fluidSolver.H

@@ -0,0 +1,179 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+Class
+    Foam::solvers::fluidSolver
+
+Description
+    Base solver module for fluid solvers.
+
+    Provides Courant number time-step control and continuity checking.
+
+    Reference:
+    \verbatim
+        Greenshields, C. J., & Weller, H. G. (2022).
+        Notes on Computational Fluid Dynamics: General Principles.
+        CFD Direct Ltd.: Reading, UK.
+    \endverbatim
+
+SourceFiles
+    fluidSolver.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fluidSolver_H
+#define fluidSolver_H
+
+#include "solver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solvers
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class fluidSolver Declaration
+\*---------------------------------------------------------------------------*/
+
+class fluidSolver
+:
+    public solver
+{
+    // Control parameters
+
+        //- Maximum allowed Courant number
+        scalar maxCo;
+
+        //- Maximum time-step
+        scalar maxDeltaT_;
+
+        //- Switch to check the mesh Courant number after mesh change
+        bool checkMeshCourantNo;
+
+
+    // Continuity properties
+
+        //- Current cumulative continuity error
+        scalar cumulativeContErr;
+
+
+    // Cached Courant numbers for time-step control
+
+        //- Current maximum Courant number
+        scalar CoNum;
+
+        //- Current mean Courant number
+        scalar meanCoNum;
+
+
+    // Private Member Functions
+
+        //- Correct the cached Courant numbers
+        template<class RhoType>
+        inline void correctCoNum
+        (
+            const RhoType& rho,
+            const surfaceScalarField& phi
+        );
+
+
+protected:
+
+        //- Switch to correct the flux after mesh change
+        bool correctPhi;
+
+        //- Read controls
+        void read();
+
+        //- Check mesh Courant numbers for moving mesh cases
+        void meshCourantNo() const;
+
+        //- Correct the cached Courant numbers
+        void correctCoNum(const surfaceScalarField& phi);
+
+        //- Correct the cached Courant numbers
+        void correctCoNum
+        (
+            const volScalarField& rho,
+            const surfaceScalarField& phi
+        );
+
+        //- Calculate and print the continuity errors
+        void continuityErrors
+        (
+            const surfaceScalarField& phi
+        );
+
+        //- Calculate and print the continuity errors
+        void continuityErrors
+        (
+            const volScalarField& rho,
+            const volScalarField& thermoRho,
+            const surfaceScalarField& phi
+        );
+
+
+public:
+
+    //- Runtime type information
+    TypeName("incompresssibleFluid");
+
+
+    // Constructors
+
+        //- Construct from region mesh
+        fluidSolver(fvMesh& mesh);
+
+        //- Disallow default bitwise copy construction
+        fluidSolver(const fluidSolver&) = delete;
+
+
+    //- Destructor
+    virtual ~fluidSolver();
+
+
+    // Member Functions
+
+        //- Return the current maximum time-step for stable solution
+        virtual scalar maxDeltaT() const;
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const fluidSolver&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace solvers
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/modules/fluid/incompressibleFluid/Make/files

@@ -0,0 +1,7 @@
+setRDeltaT.C
+moveMesh.C
+momentumPredictor.C
+correctPressure.C
+incompressibleFluid.C
+
+LIB = $(FOAM_LIBBIN)/libincompressibleFluid

applications/solvers/modules/fluid/incompressibleFluid/Make/options

@@ -1,19 +1,18 @@
 EXE_INC = \
-    -I. \
+    -I$(FOAM_SOLVERS)/modules/fluid/fluidSolver/lnInclude \
-    -I.. \
     -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
-    -I$(LIB_SRC)/sampling/lnInclude
 
-EXE_LIBS = \
+LIB_LIBS = \
+    -lfluidSolver \
     -lmomentumTransportModels \
     -lincompressibleMomentumTransportModels \
     -lphysicalProperties \
     -lfiniteVolume \
-    -lmeshTools \
     -lfvModels \
     -lfvConstraints \
-    -lsampling
+    -lsampling \
+    -lmeshTools

applications/solvers/modules/fluid/incompressibleFluid/correctPressure.C

@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "incompressibleFluid.H"
+#include "constrainHbyA.H"
+#include "constrainPressure.H"
+#include "adjustPhi.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::solvers::incompressibleFluid::correctPressure()
+{
+    fvVectorMatrix& UEqn = tUEqn.ref();
+
+    volScalarField rAU(1.0/UEqn.A());
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::flux(HbyA)
+      + MRF.zeroFilter(fvc::interpolate(rAU)*fvc::ddtCorr(U, phi, Uf))
+    );
+
+    MRF.makeRelative(phiHbyA);
+
+    if (p.needReference())
+    {
+        fvc::makeRelative(phiHbyA, U);
+        adjustPhi(phiHbyA, U, p);
+        fvc::makeAbsolute(phiHbyA, U);
+    }
+
+    tmp<volScalarField> rAtU(rAU);
+
+    if (pimple.consistent())
+    {
+        rAtU = 1.0/max(1.0/rAU - UEqn.H1(), 0.1/rAU);
+        phiHbyA +=
+            fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
+        HbyA -= (rAU - rAtU())*fvc::grad(p);
+    }
+
+    if (pimple.nCorrPiso() <= 1)
+    {
+        tUEqn.clear();
+    }
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, U, phiHbyA, rAtU(), MRF);
+
+    // Non-orthogonal pressure corrector loop
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
+        );
+
+        pEqn.setReference
+        (
+            pressureReference.refCell(),
+            pressureReference.refValue()
+        );
+
+        pEqn.solve();
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA - pEqn.flux();
+        }
+    }
+
+    continuityErrors();
+
+    // Explicitly relax pressure for momentum corrector
+    p.relax();
+
+    U = HbyA - rAtU*fvc::grad(p);
+    U.correctBoundaryConditions();
+    fvConstraints().constrain(U);
+
+    // Correct Uf if the mesh is moving
+    fvc::correctUf(Uf, U, phi, MRF);
+
+    // Make the fluxes relative to the mesh motion
+    fvc::makeRelative(phi, U);
+}
+
+
+// ************************************************************************* //

applications/solvers/modules/fluid/incompressibleFluid/incompressibleFluid.C

@@ -0,0 +1,236 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "incompressibleFluid.H"
+#include "localEulerDdtScheme.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solvers
+{
+    defineTypeNameAndDebug(incompressibleFluid, 0);
+    addToRunTimeSelectionTable(solver, incompressibleFluid, fvMesh);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::solvers::incompressibleFluid::correctCoNum()
+{
+    fluidSolver::correctCoNum(phi);
+}
+
+
+void Foam::solvers::incompressibleFluid::continuityErrors()
+{
+    fluidSolver::continuityErrors(phi);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::solvers::incompressibleFluid::incompressibleFluid(fvMesh& mesh)
+:
+    fluidSolver(mesh),
+
+    p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    ),
+
+    pressureReference(p, pimple.dict()),
+
+    U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    ),
+
+    phi
+    (
+        IOobject
+        (
+            "phi",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        linearInterpolate(U) & mesh.Sf()
+    ),
+
+    viscosity(viscosityModel::New(mesh)),
+
+    turbulence
+    (
+        incompressible::momentumTransportModel::New
+        (
+            U,
+            phi,
+            viscosity
+        )
+    ),
+
+    MRF(mesh)
+{
+    mesh.schemes().setFluxRequired(p.name());
+
+    turbulence->validate();
+
+    if (mesh.dynamic())
+    {
+        Info<< "Constructing face momentum Uf" << endl;
+
+        Uf = new surfaceVectorField
+        (
+            IOobject
+            (
+                "Uf",
+                runTime.timeName(),
+                mesh,
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE
+            ),
+            fvc::interpolate(U)
+        );
+    }
+
+    if (transient())
+    {
+        correctCoNum();
+    }
+    else if (LTS)
+    {
+        Info<< "Using LTS" << endl;
+
+        trDeltaT = tmp<volScalarField>
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    fv::localEulerDdt::rDeltaTName,
+                    runTime.timeName(),
+                    mesh,
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::AUTO_WRITE
+                ),
+                mesh,
+                dimensionedScalar(dimless/dimTime, 1),
+                extrapolatedCalculatedFvPatchScalarField::typeName
+            )
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::solvers::incompressibleFluid::~incompressibleFluid()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::solvers::incompressibleFluid::preSolve()
+{
+    // Read the controls
+    read();
+
+    fvModels().preUpdateMesh();
+
+    // Update the mesh for topology change, mesh to mesh mapping
+    mesh.update();
+
+    if (transient())
+    {
+        correctCoNum();
+    }
+    else if (LTS)
+    {
+        setRDeltaT();
+    }
+}
+
+
+void Foam::solvers::incompressibleFluid::prePredictor()
+{
+    fvModels().correct();
+}
+
+
+void Foam::solvers::incompressibleFluid::thermophysicalPredictor()
+{}
+
+
+void Foam::solvers::incompressibleFluid::pressureCorrector()
+{
+    while (pimple.correct())
+    {
+        correctPressure();
+    }
+
+    tUEqn.clear();
+}
+
+
+void Foam::solvers::incompressibleFluid::momentumTransportCorrector()
+{
+    if (pimple.transportCorr())
+    {
+        viscosity->correct();
+        turbulence->correct();
+    }
+}
+
+
+void Foam::solvers::incompressibleFluid::thermophysicalTransportCorrector()
+{}
+
+
+void Foam::solvers::incompressibleFluid::postSolve()
+{}
+
+
+// ************************************************************************* //

applications/solvers/modules/fluid/incompressibleFluid/incompressibleFluid.H

@@ -0,0 +1,219 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+Class
+    Foam::solvers::incompressibleFluid
+
+Description
+    Solver module for steady or transient turbulent flow of incompressible
+    isothermal fluids with optional mesh motion and change.
+
+    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
+    pseudo-transient and steady simulations.
+
+    Optional fvModels and fvConstraints are provided to enhance the simulation
+    in many ways including adding various sources, constraining or limiting
+    the solution.
+
+    Reference:
+    \verbatim
+        Greenshields, C. J., & Weller, H. G. (2022).
+        Notes on Computational Fluid Dynamics: General Principles.
+        CFD Direct Ltd.: Reading, UK.
+    \endverbatim
+
+SourceFiles
+    incompressibleFluid.C
+
+See also
+    Foam::solvers::fluidSolver
+    Foam::solvers::isothermalFluid
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef incompressibleFluid_H
+#define incompressibleFluid_H
+
+#include "fluidSolver.H"
+#include "viscosityModel.H"
+#include "incompressibleMomentumTransportModels.H"
+#include "pressureReference.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solvers
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class incompressibleFluid Declaration
+\*---------------------------------------------------------------------------*/
+
+class incompressibleFluid
+:
+    public fluidSolver
+{
+
+protected:
+
+    // Pressure
+
+        //- Pressure field
+        volScalarField p;
+
+        //- Pressure reference
+        Foam::pressureReference pressureReference;
+
+
+    // Kinematic properties
+
+        //- Velocity field
+        volVectorField U;
+
+        //- Mass-flux field
+        surfaceScalarField phi;
+
+
+    // Momentum transport
+
+        //- Kinematic viscosity model
+        autoPtr<viscosityModel> viscosity;
+
+        //- Pointer to the momentum transport model
+        autoPtr<incompressible::momentumTransportModel> turbulence;
+
+
+    // Continuity properties
+
+        //- Current cumulative continuity error
+        scalar cumulativeContErr;
+
+
+    // Optional models
+
+        //- MRF zone list
+        IOMRFZoneList MRF;
+
+
+    // Cached temporary fields
+
+        //- Pointer to the surface momentum field
+        //  used to recreate the flux after mesh-change
+        autoPtr<surfaceVectorField> Uf;
+
+        //- Optional LTS reciprocal time-step field
+        tmp<volScalarField> trDeltaT;
+
+        //- Cached momentum matrix
+        //  shared between the momentum predictor and pressure corrector
+        tmp<fvVectorMatrix> tUEqn;
+
+
+private:
+
+    // Private Member Functions
+
+        //- Set rDeltaT for LTS
+        void setRDeltaT();
+
+        //- Correct the cached Courant numbers
+        void correctCoNum();
+
+        //- Calculate and print the continuity errors
+        void continuityErrors();
+
+        //- Construct the pressure equation
+        //  and correct the pressure and velocity
+        void correctPressure();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("incompressibleFluid");
+
+
+    // Constructors
+
+        //- Construct from region mesh
+        incompressibleFluid(fvMesh& mesh);
+
+        //- Disallow default bitwise copy construction
+        incompressibleFluid(const incompressibleFluid&) = delete;
+
+
+    //- Destructor
+    virtual ~incompressibleFluid();
+
+
+    // Member Functions
+
+        //- Called at the start of the time-step, before the PIMPLE loop
+        virtual void preSolve();
+
+        //- Called at the start of the PIMPLE loop to move the mesh
+        virtual bool moveMesh();
+
+        //- Called at the start of the PIMPLE loop
+        virtual void prePredictor();
+
+        //- Construct and optionally solve the momentum equation
+        virtual void momentumPredictor();
+
+        //- Construct and solve the energy equation,
+        //  convert to temperature
+        //  and update thermophysical and transport properties
+        virtual void thermophysicalPredictor();
+
+        //- Construct and solve the pressure equation in the PISO loop
+        virtual void pressureCorrector();
+
+        //- Correct the momentum transport modelling
+        //  Newtonian, non-Newtonian or turbulent
+        virtual void momentumTransportCorrector();
+
+        //- Correct the thermophysical transport modelling
+        virtual void thermophysicalTransportCorrector();
+
+        //- Called after the PIMPLE loop at the end of the time-step
+        virtual void postSolve();
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const incompressibleFluid&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace solvers
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/modules/fluid/incompressibleFluid/momentumPredictor.C

@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "incompressibleFluid.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::solvers::incompressibleFluid::momentumPredictor()
+{
+    MRF.correctBoundaryVelocity(U);
+
+    tUEqn =
+    (
+        fvm::ddt(U) + fvm::div(phi, U)
+      + MRF.DDt(U)
+      + turbulence->divDevSigma(U)
+     ==
+        fvModels().source(U)
+    );
+    fvVectorMatrix& UEqn = tUEqn.ref();
+
+    UEqn.relax();
+
+    fvConstraints().constrain(UEqn);
+
+    if (pimple.momentumPredictor())
+    {
+        solve(UEqn == -fvc::grad(p));
+
+        fvConstraints().constrain(U);
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/modules/fluid/incompressibleFluid/moveMesh.C

@@ -0,0 +1,75 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "incompressibleFluid.H"
+#include "CorrectPhi.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::solvers::incompressibleFluid::moveMesh()
+{
+    if (pimple.firstIter() || pimple.moveMeshOuterCorrectors())
+    {
+        // Move the mesh
+        mesh.move();
+
+        if (mesh.changing())
+        {
+            MRF.update();
+
+            if (correctPhi)
+            {
+                // Calculate absolute flux
+                // from the mapped surface velocity
+                phi = mesh.Sf() & Uf();
+
+                correctUphiBCs(U, phi, true);
+
+                CorrectPhi
+                (
+                    phi,
+                    U,
+                    p,
+                    dimensionedScalar("rAUf", dimTime, 1),
+                    geometricZeroField(),
+                    pressureReference,
+                    pimple
+                );
+
+                // Make the flux relative to the mesh motion
+                fvc::makeRelative(phi, U);
+            }
+
+            meshCourantNo();
+
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
+// ************************************************************************* //

applications/solvers/modules/fluid/incompressibleFluid/setRDeltaT.C

@@ -1,3 +1,34 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "incompressibleFluid.H"
+#include "fvcSmooth.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::solvers::incompressibleFluid::setRDeltaT()
 {
     volScalarField& rDeltaT = trDeltaT.ref();
 
@@ -77,3 +108,6 @@
             << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
     }
 }
+
+
+// ************************************************************************* //

applications/solvers/modules/fluid/isothermalFluid/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I$(FOAM_SOLVERS)/modules/fluid/fluidSolver/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
@@ -8,6 +9,7 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 LIB_LIBS = \
+    -lfluidSolver \
     -lfluidThermophysicalModels \
     -lmomentumTransportModels \
     -lcompressibleMomentumTransportModels \

applications/solvers/modules/fluid/isothermalFluid/isothermalFluid.C

@@ -42,79 +42,15 @@ namespace solvers
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::solvers::isothermalFluid::read()
-{
-    maxCo =
-        runTime.controlDict().lookupOrDefault<scalar>("maxCo", 1.0);
-
-    maxDeltaT_ =
-        runTime.controlDict().lookupOrDefault<scalar>("maxDeltaT", great);
-
-    correctPhi = pimple.dict().lookupOrDefault
-    (
-        "correctPhi",
-        mesh.dynamic()
-    );
-
-    checkMeshCourantNo = pimple.dict().lookupOrDefault
-    (
-        "checkMeshCourantNo",
-        false
-    );
-}
-
-
 void Foam::solvers::isothermalFluid::correctCoNum()
 {
-    const scalarField sumPhi
+    fluidSolver::correctCoNum(rho, phi);
-    (
-        fvc::surfaceSum(mag(phi))().primitiveField()/rho.primitiveField()
-    );
-
-    CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
-
-    const scalar meanCoNum =
-        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
-
-    Info<< "Courant Number mean: " << meanCoNum
-        << " max: " << CoNum << endl;
 }
 
 
 void Foam::solvers::isothermalFluid::continuityErrors()
 {
-    scalar sumLocalContErr = 0;
+    fluidSolver::continuityErrors(rho, thermo.rho(), phi);
-    scalar globalContErr = 0;
-
-    if (mesh.schemes().steady())
-    {
-        const volScalarField contErr(fvc::div(phi));
-
-        sumLocalContErr =
-            runTime.deltaTValue()
-           *mag(contErr)().weightedAverage(mesh.V()).value();
-
-        globalContErr =
-            runTime.deltaTValue()
-           *contErr.weightedAverage(mesh.V()).value();
-    }
-    else
-    {
-        const dimensionedScalar totalMass = fvc::domainIntegrate(rho);
-
-        sumLocalContErr =
-            (fvc::domainIntegrate(mag(rho - thermo.rho()))/totalMass).value();
-
-        globalContErr =
-            (fvc::domainIntegrate(rho - thermo.rho())/totalMass).value();
-    }
-
-    cumulativeContErr += globalContErr;
-
-    Info<< "time step continuity errors : sum local = " << sumLocalContErr
-        << ", global = " << globalContErr
-        << ", cumulative = " << cumulativeContErr
-        << endl;
 }
 
 
@@ -126,7 +62,7 @@ Foam::solvers::isothermalFluid::isothermalFluid
     autoPtr<fluidThermo> thermoPtr
 )
 :
-    solver(mesh),
+    fluidSolver(mesh),
 
     thermo_(thermoPtr),
     thermo(thermo_()),
@@ -211,11 +147,7 @@ Foam::solvers::isothermalFluid::isothermalFluid
 
     initialMass(fvc::domainIntegrate(rho)),
 
-    cumulativeContErr(0),
+    MRF(mesh)
-
-    MRF(mesh),
-
-    CoNum(0)
 {
     // Read the controls
     read();
@@ -301,20 +233,6 @@ Foam::solvers::isothermalFluid::~isothermalFluid()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-Foam::scalar Foam::solvers::isothermalFluid::maxDeltaT() const
-{
-    if (CoNum > small)
-    {
-        const scalar deltaT = maxCo*runTime.deltaTValue()/CoNum;
-        return min(min(deltaT, fvModels().maxDeltaT()), maxDeltaT_);
-    }
-    else
-    {
-        return runTime.deltaTValue();
-    }
-}
-
-
 void Foam::solvers::isothermalFluid::preSolve()
 {
     // Read the controls

applications/solvers/modules/fluid/isothermalFluid/isothermalFluid.H

@@ -26,11 +26,15 @@ Class
 
 Description
     Solver module for steady or transient turbulent flow of compressible
-    isothermal fluids with optional mesh motion and mesh topology changes.
+    isothermal fluids with optional mesh motion and change.
 
     Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
     pseudo-transient and steady simulations.
 
+    Optional fvModels and fvConstraints are provided to enhance the simulation
+    in many ways including adding various sources, Lagrangian
+    particles, surface film etc. and constraining or limiting the solution.
+
     Reference:
     \verbatim
         Greenshields, C. J., & Weller, H. G. (2022).
@@ -41,12 +45,16 @@ Description
 SourceFiles
     isothermalFluid.C
 
+See also
+    Foam::solvers::fluidSolver
+    Foam::solvers::incompressibleFluid
+
 \*---------------------------------------------------------------------------*/
 
 #ifndef isothermalFluid_H
 #define isothermalFluid_H
 
-#include "solver.H"
+#include "fluidSolver.H"
 #include "fluidThermo.H"
 #include "compressibleMomentumTransportModels.H"
 #include "buoyancy.H"
@@ -65,26 +73,11 @@ namespace solvers
 
 class isothermalFluid
 :
-    public solver
+    public fluidSolver
 {
 
 protected:
 
-    // Control parameters
-
-        //- Maximum allowed Courant number
-        scalar maxCo;
-
-        //- Maximum time-step
-        scalar maxDeltaT_;
-
-        //- Switch to correct the flux after mesh change
-        bool correctPhi;
-
-        //- Switch to check the mesh Courant number after mesh change
-        bool checkMeshCourantNo;
-
-
     // Thermophysical properties
 
         //- Pointer to the fluid thermophysical properties
@@ -145,9 +138,6 @@ protected:
         //- Initial mass in the region
         dimensionedScalar initialMass;
 
-        //- Current cumulative continuity error
-        scalar cumulativeContErr;
-
 
     // Optional models
 
@@ -155,15 +145,6 @@ protected:
         IOMRFZoneList MRF;
 
 
-    // Time-step controls
-
-        //- Current maximum Courant number
-        scalar CoNum;
-
-        //- Current mean Courant number
-        scalar meanCoNum;
-
-
     // Cached temporary fields
 
         //- Pointer to the surface momentum field
@@ -189,29 +170,14 @@ protected:
 
 private:
 
-    // Private Member Data
-
-        //- Pointer to the demand driven fvModels MeshObject
-        mutable Foam::fvModels* fvModelsPtr;
-
-        //- Pointer to the demand driven fvConstraints MeshObject
-        mutable Foam::fvConstraints* fvConstraintsPtr;
-
-
     // Private Member Functions
 
-        //- Read controls
-        void read();
-
         //- Set rDeltaT for LTS
         virtual void setRDeltaT();
 
         //- Correct the cached Courant numbers
         void correctCoNum();
 
-        //- Check mesh Courant numbers for moving mesh cases
-        void meshCourantNo() const;
-
         //- Calculate and print the continuity errors
         void continuityErrors();
 
@@ -251,9 +217,6 @@ public:
 
     // Member Functions
 
-        //- Return the current maximum time-step for stable solution
-        virtual scalar maxDeltaT() const;
-
         //- Called at the start of the time-step, before the PIMPLE loop
         virtual void preSolve();
 

applications/solvers/modules/fluid/isothermalFluid/moveMesh.C

@@ -28,31 +28,6 @@ License
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::solvers::isothermalFluid::meshCourantNo() const
-{
-    if (checkMeshCourantNo)
-    {
-        const scalarField sumPhi
-        (
-            fvc::surfaceSum(mag(mesh.phi()))().primitiveField()
-        );
-
-        const scalar meshCoNum
-        (
-            0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue()
-        );
-
-        const scalar meanMeshCoNum
-        (
-            0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue()
-        );
-
-        Info<< "Mesh Courant Number mean: " << meanMeshCoNum
-            << " max: " << meshCoNum << endl;
-    }
-}
-
-
 bool Foam::solvers::isothermalFluid::moveMesh()
 {
     if (pimple.firstIter() || pimple.moveMeshOuterCorrectors())

applications/solvers/modules/fluid/multicomponentFluid/Make/options

@@ -1,5 +1,6 @@
 EXE_INC = \
-    -I../isothermalFluid/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/fluid/fluidSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/fluid/isothermalFluid/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
     -I$(LIB_SRC)/ThermophysicalTransportModels/lnInclude \
@@ -17,6 +18,7 @@ EXE_INC = \
     -I$(LIB_SRC)/sampling/lnInclude
 
 LIB_LIBS = \
+    -lfluidSolver \
     -lisothermalFluid \
     -lfluidThermophysicalModels \
     -lspecie \

applications/solvers/modules/fluid/multicomponentFluid/multicomponentFluid.H

@@ -26,11 +26,16 @@ Class
 
 Description
     Solver module for steady or transient turbulent flow of compressible
-    reacting fluids with optional mesh motion and mesh topology changes.
+    multicomponent fluids with optional mesh motion and change.
 
     Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
     pseudo-transient and steady simulations.
 
+    Optional fvModels and fvConstraints are provided to enhance the simulation
+    in many ways including adding various sources, chemical reactions,
+    combustion, Lagrangian particles, radiation, surface film etc. and
+    constraining or limiting the solution.
+
     Reference:
     \verbatim
         Greenshields, C. J., & Weller, H. G. (2022).
@@ -41,6 +46,10 @@ Description
 SourceFiles
     multicomponentFluid.C
 
+See also
+    Foam::solvers::fluidSolver
+    Foam::solvers::isothermalFluid
+
 \*---------------------------------------------------------------------------*/
 
 #ifndef multicomponentFluid_H
@@ -70,12 +79,19 @@ class multicomponentFluid
 
 protected:
 
-    // Compositon
+    // Thermophysical properties
 
         fluidMulticomponentThermo& thermo;
 
+
+    // Compositon
+
         basicSpecieMixture& composition;
         PtrList<volScalarField>& Y;
+
+
+    // Reactions
+
         autoPtr<combustionModel> reaction;
         multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 

bin/buoyantFoam

@@ -27,14 +27,15 @@
 #
 # Description
 #     Script to inform the user that buoyantFoam has been superseded
-#     and replaced by the more general foamMultiRun solver.
+#     and replaced by the more general fluid solver module executed by the
+#     foamRun application.
 #
 #------------------------------------------------------------------------------
 
 cat <<EOF
 
 buoyantFoam has been superseded and replaced by the more general
-fluid solver module which will be executed using
+fluid solver module executed by the foamRun application:
 
     foamRun -solver fluid
 

bin/buoyantReactingFoam

@@ -27,14 +27,15 @@
 #
 # Description
 #     Script to inform the user that buoyantReactingFoam has been superseded
-#     and replaced by the more general foamMultiRun solver.
+#     and replaced by the more general multicomponentFluid solver module
+#     executed by the foamRun application.
 #
 #------------------------------------------------------------------------------
 
 cat <<EOF
 
 buoyantReactingFoam has been superseded and replaced by the more general
-fluid solver module which will be executed using
+multicomponentFluid solver module executed by the foamRun application:
 
     foamRun -solver multicomponentFluid
 

bin/pimpleFoam

@@ -0,0 +1,46 @@
+#!/bin/sh
+#------------------------------------------------------------------------------
+# =========                 |
+# \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+#  \\    /   O peration     | Website:  https://openfoam.org
+#   \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+#    \\/     M anipulation  |
+#------------------------------------------------------------------------------
+# 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/>.
+#
+# Script
+#     pimpleFoam
+#
+# Description
+#     Script to inform the user that pimpleFoam has been superseded
+#     and replaced by the more general incompressibleFluid solver module
+#     executed by the foamRun application.
+#
+#------------------------------------------------------------------------------
+
+cat <<EOF
+
+pimpleFoam has been superseded and replaced by the more general
+fluid solver module executed by the foamRun application:
+
+    foamRun -solver incompressibleFluid
+
+EOF
+
+foamRun -solver incompressibleFluid "$@"
+
+#------------------------------------------------------------------------------

bin/pisoFoam

@@ -0,0 +1,46 @@
+#!/bin/sh
+#------------------------------------------------------------------------------
+# =========                 |
+# \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+#  \\    /   O peration     | Website:  https://openfoam.org
+#   \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+#    \\/     M anipulation  |
+#------------------------------------------------------------------------------
+# 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/>.
+#
+# Script
+#     pisoFoam
+#
+# Description
+#     Script to inform the user that pisoFoam has been superseded
+#     and replaced by the more general incompressibleFluid solver module
+#     executed by the foamRun application.
+#
+#------------------------------------------------------------------------------
+
+cat <<EOF
+
+pisoFoam has been superseded and replaced by the more general
+fluid solver module executed by the foamRun application:
+
+    foamRun -solver incompressibleFluid
+
+EOF
+
+foamRun -solver incompressibleFluid "$@"
+
+#------------------------------------------------------------------------------

bin/reactingFoam

@@ -27,14 +27,15 @@
 #
 # Description
 #     Script to inform the user that reactingFoam has been superseded
-#     and replaced by the more general foamMultiRun solver.
+#     and replaced by the more general multicomponentFluid solver module
+#     executed by the foamRun application.
 #
 #------------------------------------------------------------------------------
 
 cat <<EOF
 
 reactingFoam has been superseded and replaced by the more general
-fluid solver module which will be executed using
+multicomponentFluid solver module executed by the foamRun application:
 
     foamRun -solver multicomponentFluid
 

bin/rhoPimpleFoam

@@ -27,14 +27,15 @@
 #
 # Description
 #     Script to inform the user that rhoPimpleFoam has been superseded
-#     and replaced by the more general foamMultiRun solver.
+#     and replaced by the more general fluid solver module executed by the
+#     foamRun application.
 #
 #------------------------------------------------------------------------------
 
 cat <<EOF
 
 rhoPimpleFoam has been superseded and replaced by the more general
-fluid solver module which will be executed using
+fluid solver module executed by the foamRun application:
 
     foamRun -solver fluid
 

bin/simpleFoam

@@ -0,0 +1,46 @@
+#!/bin/sh
+#------------------------------------------------------------------------------
+# =========                 |
+# \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+#  \\    /   O peration     | Website:  https://openfoam.org
+#   \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+#    \\/     M anipulation  |
+#------------------------------------------------------------------------------
+# 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/>.
+#
+# Script
+#     simpleFoam
+#
+# Description
+#     Script to inform the user that simpleFoam has been superseded
+#     and replaced by the more general incompressibleFluid solver module
+#     executed by the foamRun application.
+#
+#------------------------------------------------------------------------------
+
+cat <<EOF
+
+simpleFoam has been superseded and replaced by the more general
+fluid solver module executed by the foamRun application:
+
+    foamRun -solver incompressibleFluid
+
+EOF
+
+foamRun -solver incompressibleFluid "$@"
+
+#------------------------------------------------------------------------------

etc/config.sh/bash_completion

@@ -2614,56 +2614,6 @@ _PDRMesh_ ()
 }
 complete -o filenames -o nospace -F _PDRMesh_ PDRMesh
 
-_pimpleFoam_ ()
-{
-    local cur="${COMP_WORDS[COMP_CWORD]}"
-    local prev="${COMP_WORDS[COMP_CWORD-1]}"
-    local line=${COMP_LINE}
-    local used=$(echo "$line" | grep -oE "\-[a-zA-Z]+ ")
-
-    opts="-case -doc -fileHandler -help -hostRoots -libs -listFunctionObjects -listFvConstraints -listFvModels -listMomentumTransportModels -listScalarBCs -listSwitches -listVectorBCs -noFunctionObjects -parallel -postProcess -roots -srcDoc"
-    for o in $used ; do opts="${opts/$o/}" ; done
-    extra=""
-
-    [ "$COMP_CWORD" = 1 ] || \
-    case "$prev" in
-        -case)
-            opts="" ; extra="-d" ;;
-        -fileHandler)
-            opts="uncollated collated masterUncollated" ; extra="" ;;
-        -hostRoots|-libs|-roots)
-            opts="" ; extra="" ;;
-       *) ;;
-    esac
-    COMPREPLY=( $(compgen -W "${opts}" $extra -- ${cur}) )
-}
-complete -o filenames -o nospace -F _pimpleFoam_ pimpleFoam
-
-_pisoFoam_ ()
-{
-    local cur="${COMP_WORDS[COMP_CWORD]}"
-    local prev="${COMP_WORDS[COMP_CWORD-1]}"
-    local line=${COMP_LINE}
-    local used=$(echo "$line" | grep -oE "\-[a-zA-Z]+ ")
-
-    opts="-case -doc -fileHandler -help -hostRoots -libs -listFunctionObjects -listFvConstraints -listFvModels -listMomentumTransportModels -listScalarBCs -listSwitches -listVectorBCs -noFunctionObjects -parallel -postProcess -roots -srcDoc"
-    for o in $used ; do opts="${opts/$o/}" ; done
-    extra=""
-
-    [ "$COMP_CWORD" = 1 ] || \
-    case "$prev" in
-        -case)
-            opts="" ; extra="-d" ;;
-        -fileHandler)
-            opts="uncollated collated masterUncollated" ; extra="" ;;
-        -hostRoots|-libs|-roots)
-            opts="" ; extra="" ;;
-       *) ;;
-    esac
-    COMPREPLY=( $(compgen -W "${opts}" $extra -- ${cur}) )
-}
-complete -o filenames -o nospace -F _pisoFoam_ pisoFoam
-
 _plot3dToFoam_ ()
 {
     local cur="${COMP_WORDS[COMP_CWORD]}"
@@ -3322,7 +3272,7 @@ _simpleFoam_ ()
     local line=${COMP_LINE}
     local used=$(echo "$line" | grep -oE "\-[a-zA-Z]+ ")
 
-    opts="-case -doc -fileHandler -help -hostRoots -libs -listFunctionObjects -listFvConstraints -listFvModels -listMomentumTransportModels -listScalarBCs -listSwitches -listVectorBCs -noFunctionObjects -parallel -postProcess -roots -srcDoc"
+    opts="-case -doc -fileHandler -help -hostRoots -libs -noFunctionObjects -parallel -roots -solver -srcDoc"
     for o in $used ; do opts="${opts/$o/}" ; done
     extra=""
 
@@ -3332,7 +3282,7 @@ _simpleFoam_ ()
             opts="" ; extra="-d" ;;
         -fileHandler)
             opts="uncollated collated masterUncollated" ; extra="" ;;
-        -hostRoots|-libs|-roots)
+        -hostRoots|-libs|-roots|-solver)
             opts="" ; extra="" ;;
        *) ;;
     esac
@@ -5351,6 +5301,56 @@ _paraFoam_ ()
 }
 complete -o filenames -o nospace -F _paraFoam_ paraFoam
 
+_pimpleFoam_ ()
+{
+    local cur="${COMP_WORDS[COMP_CWORD]}"
+    local prev="${COMP_WORDS[COMP_CWORD-1]}"
+    local line=${COMP_LINE}
+    local used=$(echo "$line" | grep -oE "\-[a-zA-Z]+ ")
+
+    opts="-case -doc -fileHandler -help -hostRoots -libs -noFunctionObjects -parallel -roots -solver -srcDoc"
+    for o in $used ; do opts="${opts/$o/}" ; done
+    extra=""
+
+    [ "$COMP_CWORD" = 1 ] || \
+    case "$prev" in
+        -case)
+            opts="" ; extra="-d" ;;
+        -fileHandler)
+            opts="uncollated collated masterUncollated" ; extra="" ;;
+        -hostRoots|-libs|-roots|-solver)
+            opts="" ; extra="" ;;
+       *) ;;
+    esac
+    COMPREPLY=( $(compgen -W "${opts}" $extra -- ${cur}) )
+}
+complete -o filenames -o nospace -F _pimpleFoam_ pimpleFoam
+
+_pisoFoam_ ()
+{
+    local cur="${COMP_WORDS[COMP_CWORD]}"
+    local prev="${COMP_WORDS[COMP_CWORD-1]}"
+    local line=${COMP_LINE}
+    local used=$(echo "$line" | grep -oE "\-[a-zA-Z]+ ")
+
+    opts="-case -doc -fileHandler -help -hostRoots -libs -noFunctionObjects -parallel -roots -solver -srcDoc"
+    for o in $used ; do opts="${opts/$o/}" ; done
+    extra=""
+
+    [ "$COMP_CWORD" = 1 ] || \
+    case "$prev" in
+        -case)
+            opts="" ; extra="-d" ;;
+        -fileHandler)
+            opts="uncollated collated masterUncollated" ; extra="" ;;
+        -hostRoots|-libs|-roots|-solver)
+            opts="" ; extra="" ;;
+       *) ;;
+    esac
+    COMPREPLY=( $(compgen -W "${opts}" $extra -- ${cur}) )
+}
+complete -o filenames -o nospace -F _pisoFoam_ pisoFoam
+
 _postProcess_ ()
 {
     local cur="${COMP_WORDS[COMP_CWORD]}"
@@ -5455,6 +5455,31 @@ _rhoSimpleFoam_ ()
 }
 complete -o filenames -o nospace -F _rhoSimpleFoam_ rhoSimpleFoam
 
+_simpleFoam_ ()
+{
+    local cur="${COMP_WORDS[COMP_CWORD]}"
+    local prev="${COMP_WORDS[COMP_CWORD-1]}"
+    local line=${COMP_LINE}
+    local used=$(echo "$line" | grep -oE "\-[a-zA-Z]+ ")
+
+    opts="-case -doc -fileHandler -help -hostRoots -libs -noFunctionObjects -parallel -roots -solver -srcDoc"
+    for o in $used ; do opts="${opts/$o/}" ; done
+    extra=""
+
+    [ "$COMP_CWORD" = 1 ] || \
+    case "$prev" in
+        -case)
+            opts="" ; extra="-d" ;;
+        -fileHandler)
+            opts="uncollated collated masterUncollated" ; extra="" ;;
+        -hostRoots|-libs|-roots|-solver)
+            opts="" ; extra="" ;;
+       *) ;;
+    esac
+    COMPREPLY=( $(compgen -W "${opts}" $extra -- ${cur}) )
+}
+complete -o filenames -o nospace -F _simpleFoam_ simpleFoam
+
 _wclean_ ()
 {
     local cur="${COMP_WORDS[COMP_CWORD]}"

tutorials/modules/incompressibleFluid/T3A/system/controlDict

@@ -14,7 +14,9 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     simpleFoam;
+application     foamRun;
+
+solver          incompressibleFluid;
 
 startFrom       startTime;
 

tutorials/modules/incompressibleFluid/TJunction/system/controlDict

@@ -14,7 +14,9 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     pimpleFoam;
+application     foamRun;
+
+solver          incompressibleFluid;
 
 startFrom       startTime;