From 2bc91ecff08a14677dd8ba61f247cff57c192985 Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Wed, 6 Dec 2023 12:39:34 +0000
Subject: [PATCH] phaseScalarTransport: Improved interface and documentation

This function now looks up an alphaPhi field by default and exits with
an error if the field cannot be found. In order to solve for alphaPhi a
new 'solveAlphaPhi' switch has to be set.

The documentation has been updated to reflect the fact that the VoF
solvers now store all the alphaPhi fluxes necessary for a in-phase
equation to be constructed.

The phaseScalarTransport function has been added to the damBreakLaminar
tutorial.
---
 .../phaseScalarTransport.C                    | 27 +++---
 .../phaseScalarTransport.H                    | 30 +++----
 .../damBreakLaminar/0/tracer.air.orig         | 39 ++++++++
 .../damBreakLaminar/0/tracer.water.orig       | 39 ++++++++
 .../damBreakLaminar/Allclean                  |  4 +-
 .../incompressibleVoF/damBreakLaminar/Allrun  |  4 -
 .../constant/momentumTransport                |  7 --
 .../damBreakLaminar/system/controlDict        | 62 +++++++++++++
 .../damBreakLaminar/system/fvSchemes          | 58 ++++++++++++
 .../damBreakLaminar/system/fvSolution         | 90 +++++++++++++++++++
 .../damBreakLaminar/system/setFieldsDict      | 53 +++++++++++
 11 files changed, 370 insertions(+), 43 deletions(-)
 create mode 100644 tutorials/incompressibleVoF/damBreakLaminar/0/tracer.air.orig
 create mode 100644 tutorials/incompressibleVoF/damBreakLaminar/0/tracer.water.orig
 create mode 100644 tutorials/incompressibleVoF/damBreakLaminar/system/controlDict
 create mode 100644 tutorials/incompressibleVoF/damBreakLaminar/system/fvSchemes
 create mode 100644 tutorials/incompressibleVoF/damBreakLaminar/system/fvSolution
 create mode 100644 tutorials/incompressibleVoF/damBreakLaminar/system/setFieldsDict

diff --git a/src/functionObjects/solvers/phaseScalarTransport/phaseScalarTransport.C b/src/functionObjects/solvers/phaseScalarTransport/phaseScalarTransport.C
index 83bd7be2dd..7d3eebc3e4 100644
--- a/src/functionObjects/solvers/phaseScalarTransport/phaseScalarTransport.C
+++ b/src/functionObjects/solvers/phaseScalarTransport/phaseScalarTransport.C
@@ -90,7 +90,7 @@ Foam::volScalarField& Foam::functionObjects::phaseScalarTransport::Phi()
             (
                 IOobject
                 (
-                    "Phi" + s_.name(),
+                    typedName(IOobject::groupName("Phi", phaseName_)),
                     time_.name(),
                     mesh_,
                     IOobject::READ_IF_PRESENT,
@@ -112,16 +112,11 @@ Foam::volScalarField& Foam::functionObjects::phaseScalarTransport::Phi()
 Foam::tmp<Foam::surfaceScalarField>
 Foam::functionObjects::phaseScalarTransport::alphaPhi()
 {
-    // If alphaPhi exists then return it
-    if (mesh_.foundObject<surfaceScalarField>(alphaPhiName_))
+    if (!solveAlphaPhi_)
     {
         return mesh_.lookupObject<surfaceScalarField>(alphaPhiName_);
     }
 
-    // Otherwise generate it ...
-    Info<< type() << ": " << surfaceScalarField::typeName << " "
-        << alphaPhiName_ << " was not found, so generating it" << endl;
-
     const volScalarField& alpha =
         mesh_.lookupObject<volScalarField>(alphaName_);
     const surfaceScalarField& phi =
@@ -288,8 +283,6 @@ Foam::functionObjects::phaseScalarTransport::phaseScalarTransport
     fvMeshFunctionObject(name, runTime, dict),
     fieldName_(dict.lookup("field")),
     phaseName_(IOobject::group(fieldName_)),
-    nCorr_(0),
-    residualAlpha_(rootSmall),
     s_
     (
         IOobject
@@ -330,18 +323,20 @@ bool Foam::functionObjects::phaseScalarTransport::read(const dictionary& dict)
 {
     fvMeshFunctionObject::read(dict);
 
+    solveAlphaPhi_ = dict.lookupOrDefault<bool>("solveAlphaPhi", false);
+
     alphaName_ =
         dict.lookupOrDefault<word>
         (
             "alpha",
             IOobject::groupName("alpha", phaseName_)
         );
+    const word defaultAlphaPhiName =
+        IOobject::groupName("alphaPhi", phaseName_);
     alphaPhiName_ =
-        dict.lookupOrDefault<word>
-        (
-            "alphaPhi",
-            IOobject::groupName("alphaPhi", phaseName_)
-        );
+        solveAlphaPhi_
+      ? typedName(defaultAlphaPhiName)
+      : dict.lookupOrDefault<word>("alphaPhi", defaultAlphaPhiName);
     phiName_ = dict.lookupOrDefault<word>("phi", "phi");
     rhoName_ =
         dict.lookupOrDefault<word>
@@ -356,8 +351,8 @@ bool Foam::functionObjects::phaseScalarTransport::read(const dictionary& dict)
     alphaD_ = dict.lookupOrDefault<scalar>("alphaD", 1);
     alphaDt_ = dict.lookupOrDefault<scalar>("alphaDt", 1);
 
-    dict.readIfPresent("nCorr", nCorr_);
-    dict.readIfPresent("residualAlpha", residualAlpha_);
+    nCorr_ = dict.lookupOrDefault<label>("nCorr", 0);
+    residualAlpha_ = dict.lookupOrDefault<scalar>("residualAlpha", rootSmall);
     writeAlphaField_ = dict.lookupOrDefault<bool>("writeAlphaField", true);
 
     return true;
diff --git a/src/functionObjects/solvers/phaseScalarTransport/phaseScalarTransport.H b/src/functionObjects/solvers/phaseScalarTransport/phaseScalarTransport.H
index 2a4ca1450d..e851e2bd3f 100644
--- a/src/functionObjects/solvers/phaseScalarTransport/phaseScalarTransport.H
+++ b/src/functionObjects/solvers/phaseScalarTransport/phaseScalarTransport.H
@@ -37,27 +37,25 @@ Description
     detailed below. Note that the phase-name will be determined by stripping
     the extension from the supplied field name.
 
-    If \c alphaPhi is specified and found in the database then it will be used
-    to transport the field. This is the likely mode of operation if this
-    function is used with Euler-Euler solvers, in which \c alphaPhi -like
-    fluxes are available. If \c alphaPhi is not found, then a pressure-like
-    equation will be solved in order to construct it so that it exactly matches
-    the time-derivative of the phase volume or mass. This is likely to be
-    necessary in volume-of-fluid solvers where \c alphaPhi is not part of the
-    solution procedure. The pressure field name will be required in this case.
+    If the solver does not provide an \c alphaPhi flux, or that flux is for
+    some reason unreliable, then the \c solveAlphaPhi switch can be used to
+    make this function solve a pressure-like equation from which \c alphaPhi is
+    recovered.
 
 Usage
     \table
-        Property        | Description                        | Req'd? | Default
-        alpha           | Name of the volume-fraction field  | no\\
+        Property        | Description                       | Req'd? | Default
+        alpha           | Name of the volume-fraction field | no\\
                         | alpha.<phase-name>
-        alphaPhi        | Name of the phase-flux field       | no\\
+        alphaPhi        | Name of the phase-flux field      | no\\
                         | alphaPhi.<phase-name>
-        p               | Name of the pressure field         | no     | p
+        solveAlphaPhi   | Solve for the alphaPhi flux, rather than looking it\\
+                          up                                | no     | false
+        p               | Name of the pressure field        | no     | p
         residualAlpha   | Small volume fraction used to stabilise the solution\\
-                        | no | rootSmall
+                                                            | no     | rootSmall
         writeAlphaField | Also write out alpha multiplied by the field\\
-                        | no | true
+                                                            |        no | true
     \endtable
 
     Example specification for incompressibleVoF:
@@ -68,7 +66,6 @@ Usage
         libs            ("libsolverFunctionObjects.so");
 
         field           s.water;
-        p               p_rgh;
     }
     \endverbatim
 
@@ -123,6 +120,9 @@ class phaseScalarTransport
         //- Name of the phase in which to solve
         const word phaseName_;
 
+        //- Solve for the alphaPhi flux, rather than looking it up (optional)
+        bool solveAlphaPhi_;
+
         //- Name of phase volume-fraction field (optional)
         word alphaName_;
 
diff --git a/tutorials/incompressibleVoF/damBreakLaminar/0/tracer.air.orig b/tutorials/incompressibleVoF/damBreakLaminar/0/tracer.air.orig
new file mode 100644
index 0000000000..6b05c88fc0
--- /dev/null
+++ b/tutorials/incompressibleVoF/damBreakLaminar/0/tracer.air.orig
@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      tracer.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    wall
+    {
+        type            zeroGradient;
+    }
+
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      $internalField;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressibleVoF/damBreakLaminar/0/tracer.water.orig b/tutorials/incompressibleVoF/damBreakLaminar/0/tracer.water.orig
new file mode 100644
index 0000000000..2244f7678c
--- /dev/null
+++ b/tutorials/incompressibleVoF/damBreakLaminar/0/tracer.water.orig
@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      tracer.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    wall
+    {
+        type            zeroGradient;
+    }
+
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      $internalField;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressibleVoF/damBreakLaminar/Allclean b/tutorials/incompressibleVoF/damBreakLaminar/Allclean
index 7b06dc70e3..9068aa2e07 100755
--- a/tutorials/incompressibleVoF/damBreakLaminar/Allclean
+++ b/tutorials/incompressibleVoF/damBreakLaminar/Allclean
@@ -6,10 +6,12 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 cleanVoFCase
 
-rm -rf 0 system
+find 0 -type f -not -name tracer.* -delete
 
 find constant -type f -not \( \
     -name fvModels -or -name momentumTransport \
     \) -delete
 
+rm -f system/blockMeshDict
+
 #------------------------------------------------------------------------------
diff --git a/tutorials/incompressibleVoF/damBreakLaminar/Allrun b/tutorials/incompressibleVoF/damBreakLaminar/Allrun
index 6bb259ecec..8064eb76ca 100755
--- a/tutorials/incompressibleVoF/damBreakLaminar/Allrun
+++ b/tutorials/incompressibleVoF/damBreakLaminar/Allrun
@@ -9,10 +9,6 @@ isTest "$@" && path=.. || path=$FOAM_TUTORIALS/incompressibleVoF
 cp -rn $path/damBreak/0 $path/damBreak/constant $path/damBreak/system .
 rm -f 0/alpha.water
 
-# Increase the number of limiter iterations
-runApplication foamDictionary system/fvSolution \
-    -entry solvers/"alpha.water.*"/nLimiterIter -set 5
-
 # Run
 runApplication blockMesh
 runApplication setFields
diff --git a/tutorials/incompressibleVoF/damBreakLaminar/constant/momentumTransport b/tutorials/incompressibleVoF/damBreakLaminar/constant/momentumTransport
index c7745a3880..666fcc4955 100644
--- a/tutorials/incompressibleVoF/damBreakLaminar/constant/momentumTransport
+++ b/tutorials/incompressibleVoF/damBreakLaminar/constant/momentumTransport
@@ -16,12 +16,5 @@ FoamFile
 
 simulationType  laminar;
 
-RAS
-{
-    model           kEpsilon;
-    turbulence      on;
-    printCoeffs     on;
-}
-
 
 // ************************************************************************* //
diff --git a/tutorials/incompressibleVoF/damBreakLaminar/system/controlDict b/tutorials/incompressibleVoF/damBreakLaminar/system/controlDict
new file mode 100644
index 0000000000..8e235f7f13
--- /dev/null
+++ b/tutorials/incompressibleVoF/damBreakLaminar/system/controlDict
@@ -0,0 +1,62 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     foamRun;
+
+solver          incompressibleVoF;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         1;
+
+deltaT          0.001;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.05;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxCo           1;
+maxAlphaCo      1;
+
+maxDeltaT       1;
+
+functions
+{
+    #includeFunc phaseScalarTransport(field=tracer.water)
+    #includeFunc phaseScalarTransport(field=tracer.air, residualAlpha=1e-6)
+}
+
+// ************************************************************************* //
diff --git a/tutorials/incompressibleVoF/damBreakLaminar/system/fvSchemes b/tutorials/incompressibleVoF/damBreakLaminar/system/fvSchemes
new file mode 100644
index 0000000000..d96420f79b
--- /dev/null
+++ b/tutorials/incompressibleVoF/damBreakLaminar/system/fvSchemes
@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    div(phi,alpha)  Gauss interfaceCompression vanLeer 1;
+
+    div(rhoPhi,U)   Gauss linearUpwind grad(U);
+
+    div(phi,k)      Gauss upwind;
+    div(phi,epsilon) Gauss upwind;
+    div(phi,omega)  Gauss upwind;
+
+    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
+
+    div(alphaPhi,tracer) Gauss vanLeer;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear uncorrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         uncorrected;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressibleVoF/damBreakLaminar/system/fvSolution b/tutorials/incompressibleVoF/damBreakLaminar/system/fvSolution
new file mode 100644
index 0000000000..ff81240933
--- /dev/null
+++ b/tutorials/incompressibleVoF/damBreakLaminar/system/fvSolution
@@ -0,0 +1,90 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    "alpha.water.*"
+    {
+        nAlphaCorr      2;
+        nAlphaSubCycles 1;
+
+        MULESCorr       yes;
+        nLimiterIter    5;
+
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-8;
+        relTol          0;
+    }
+
+    "pcorr.*"
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-5;
+        relTol          0;
+    }
+
+    p_rgh
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-07;
+        relTol          0.05;
+    }
+
+    p_rghFinal
+    {
+        $p_rgh;
+        relTol          0;
+    }
+
+    "(U|k|epsilon|omega).*"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-06;
+        relTol          0;
+        minIter         1;
+    }
+
+    "tracer.*"
+    {
+        solver          PBiCGStab;
+        preconditioner  DILU;
+        tolerance       1e-06;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    momentumPredictor   no;
+    nOuterCorrectors    1;
+    nCorrectors         3;
+    nNonOrthogonalCorrectors 0;
+}
+
+relaxationFactors
+{
+    equations
+    {
+        ".*"                1;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressibleVoF/damBreakLaminar/system/setFieldsDict b/tutorials/incompressibleVoF/damBreakLaminar/system/setFieldsDict
new file mode 100644
index 0000000000..936e153d7b
--- /dev/null
+++ b/tutorials/incompressibleVoF/damBreakLaminar/system/setFieldsDict
@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defaultFieldValues
+(
+    volScalarFieldValue alpha.water 0
+    volScalarFieldValue tracer.water 0
+    volScalarFieldValue tracer.air 0
+);
+
+regions
+(
+    boxToCell
+    {
+        box (0 0 -1) (0.1461 0.292 1);
+        fieldValues
+        (
+            volScalarFieldValue alpha.water 1
+        );
+    }
+    boxToCell
+    {
+        box (0 0 -1) (0.1461 0.146 1);
+        fieldValues
+        (
+            volScalarFieldValue tracer.water 1
+        );
+    }
+    boxToCell
+    {
+        box (0 0.292 -1) (0.1461 0.438 1);
+        fieldValues
+        (
+            volScalarFieldValue tracer.air 1
+        );
+    }
+);
+
+
+// ************************************************************************* //