Handling consistency in solutionControls class

Note: additional machinery to enable easy top-level calls for
time/under-relaxation consistency in segregated solution algorithms.
Does not compile: still missing functions in ddtSchemes

applications/solvers/incompressible/icoFoam/icoFoam.C

@@ -60,6 +60,9 @@ int main(int argc, char *argv[])
 
 #       include "CourantNo.H"
 
+        // Time-derivative matrix
+        fvVectorMatrix ddtUEqn(fvm::ddt(U));
+
         // Convection-diffusion matrix
         fvVectorMatrix HUEqn
         (
@@ -69,7 +72,7 @@ int main(int argc, char *argv[])
 
         if (piso.momentumPredictor())
         {
-            solve(fvm::ddt(U) + HUEqn == -fvc::grad(p));
+            solve(ddtUEqn + HUEqn == -fvc::grad(p));
         }
 
         // Prepare clean 1/a_p without time derivative contribution
@@ -81,9 +84,9 @@ int main(int argc, char *argv[])
         {
             // Calculate U from convection-diffusion matrix
             U = rAU*HUEqn.H();
-            
+
-            // Consistently calculate flux and face velocity
+            // Consistently calculate flux
-            phi = piso.timeConsistentFlux(U, rAU);
+            piso.calculateTimeConsistentFlux(phi, U, rAU);
 
             adjustPhi(phi, U, p);
 
@@ -91,7 +94,14 @@ int main(int argc, char *argv[])
             {
                 fvScalarMatrix pEqn
                 (
-                    fvm::laplacian(rAU, p) == fvc::div(phi)
+                    fvm::laplacian
+                    (
+                        rAU/fvc::interpolate(piso.aCoeff()),
+                        p,
+                        "laplacian(rAU," + p.name() + ')'
+                    )
+                 ==
+                    fvc::div(phi)
                 );
 
                 pEqn.setReference(pRefCell, pRefValue);
@@ -109,8 +119,7 @@ int main(int argc, char *argv[])
 #           include "continuityErrs.H"
 
             // Consistently reconstruct velocity after pressure equation
-            U = piso.timeConsistentVelocity(U, rAU, p);
+            piso.reconstructVelocity(U, ddtUEqn, rAU, p, phi);
-            U.correctBoundaryConditions();
         }
 
         runTime.write();

src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C

@@ -25,6 +25,7 @@ License
 
 #include "solutionControl.H"
 #include "lduMatrix.H"
+#include "demandDrivenData.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -236,14 +237,244 @@ Foam::solutionControl::solutionControl(fvMesh& mesh, const word& algorithmName)
     transonic_(false),
     consistent_(false),
     corr_(0),
-    corrNonOrtho_(0)
+    corrNonOrtho_(0),
+    aCoeffPtr_(NULL),
+    faceUPtr_(NULL)
 {}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::solutionControl::~solutionControl()
-{}
+{
+    deleteDemandDrivenData(aCoeffPtr_);
+    deleteDemandDrivenData(faceUPtr_);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::solutionControl::calcTimeConsistentFlux
+(
+    surfaceScalarField& phi,
+    const volVectorField& U,
+    const volScalarField& rAU
+)
+{
+    // Store necessary fields for time consistency if they are not present
+    if (!aCoeffPtr_ && !faceUPtr_)
+    {
+        aCoeffPtr_ = new volScalarField
+        (
+            IOobject
+            (
+                "aCoeff",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimless, 0.0)
+        );
+
+        faceUPtr_ = new surfaceVectorField
+        (
+            IOobject
+            (
+                "faceU",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimVelocity, 0.0)
+        );
+    }
+    else if (!aCoeffPtr_ || !faceUPtr_)
+    {
+        FatalErrorIn
+        (
+            "tmp<surfaceScalarField> solutionControl::timeConsistentFlux"
+            "\n("
+            "\n    const volVectorField& U,"
+            "\n    const volScalarField& rAU"
+            "\n)"
+        )   << "Either aCoeffPtr_ or faceUPtr_ is allocated while the"
+            << " other is not. This must not happen."
+            << endl
+            << exit(FatalError);
+    }
+
+    // Now that we are sure that pointers are valid and the fields are there, we
+    // can do the calculation of the flux, while storing necessary data for
+    // future velocity reconstruction.
+    // Algorithm:
+    // 1. Update flux and aCoeff due to ddt discretisation
+    // 2. Update flux and aCoeff due to under-relaxation
+    // 3. Scale the flux with aCoeff, making sure that flux at fixed boundaries
+    //    remains is consistent
+
+    // Get fields that will be updated
+    volScalarField& aCoeff = *aCoeffPtr_;
+    surfaceVectorField& faceU = *faceUPtr_;
+
+    // Update face interpolated velocity field. Note: handling of oldTime faceU
+    // fields happens in ddt scheme when calling ddtConsistentPhiCorr
+    faceU = fvc::interpolate(U);
+
+    // Interpolate original rAU on the faces
+    const surfaceScalarField rAUf = fvc::interpolate(rAU);
+
+    // Calculate the ordinary part of the flux (H/A)
+    phi = (faceU & mesh_.Sf());
+
+    // STAGE 1: consistent ddt discretisation handling
+
+    // Add consistent flux contribution due to ddt discretisation
+    phi += fvc::ddtConsistentPhiCorr(U, rAUf, faceU);
+
+    // Reset aCoeff according to ddt discretisation
+    aCoeff = fvc::ddtRAUCorr(U, rAU); // Note: returns 1 + ddt(U).A()*rAU
+
+    // STAGE 2: consistent under-relaxation handling
+
+    // Get under-relaxation factor used in this iteration
+    const dimensionedScalar alphaU
+    (
+        "alphaU",
+        dimless,
+        mesh_.solutionDict().equationRelaxationFactor
+        (
+            U.select(this->finalIter())
+        )
+    );
+
+    // Helper variable
+    const dimensionedScalar urfCoeff = (1.0 - alphaU)/alphaU;
+
+    // Add consistent flux contribution due to possible under-relaxation
+    phi += urfCoeff*fvc::interpolate(aCoeff)*phi.prevIter();
+
+    // Add under-relaxation contribution to aCoeff
+    aCoeff += urfCoeff*aCoeff;
+
+    // Note: this is not 100% entirely consistent with the way FOAM handles
+    // under-relaxation. In fvMatrix::relax() member function, the diagonal
+    // dominance is first assured by limiting the diagonal and the matrix is
+    // then relaxed. However, this is a very minor inconsistency, especially
+    // since this term vanishes when non-linear iterations converge to a tight
+    // tolerance (achieving steady state in SIMPLE or converging PISO/PIMPLE in
+    // each time step). VV, 23/Dec/2016.
+
+    // STAGE 3: scale the flux and correct it at the boundaries
+
+    // Scale the flux
+    phi /= fvc::interpolate(aCoeff);
+
+    // Get necessary data
+    const volVectorField::GeometricBoundaryField& Ub = U.boundaryField();
+    const surfaceVectorField::GeometricBoundaryField& Sb =
+        mesh_.Sf().boundaryField();
+
+    surfaceScalarField::GeometricBoundaryField& phib = phi.boundaryField();
+
+    // Correct flux at boundaries depending on patch type. Note: it is possible
+    // that we missed some important boundary conditions that need special
+    // considerations when calculating the flux. Maybe we can reorganise this in
+    // future by relying on distinction between = and == operators for
+    // fvsPatchFields. However, that would require serious changes at the
+    // moment. VV, 23/Dec/2016.
+    forAll (phib, patchI)
+    {
+        if (Ub[patchI].fixesValue())
+        {
+            // This is fixed value patch, flux needs to be recalculated
+            // respecting the boundary condition
+            phib[patchI] == (Ub[patchI] & Sb[patchI]);
+        }
+        else if
+        (
+            isA<slipFvPatchVectorField>(Ub[patchI])
+         || isA<symmetryFvPatchVectorField>(Ub[patchI]])
+        )
+        {
+            // This is slip or symmetry, flux needs to be zero
+            phib[patchI] == 0.0;
+        }
+    }
+}
+
+
+void Foam::solutionControl::reconstructVelocity
+(
+    volVectorField& U,
+    const fvVectorMatrix& ddtUEqn,
+    const volScalarField& rAU,
+    const volScalarField& p,
+    const surfaceScalarField& phi
+)
+{
+    // Reconstruct the velocity using all the components from original equation
+    U = 1.0/(1.0/rAU + ddtUEqn.A())*
+        (
+            U/rAU + ddtUEqn.H() - fvc::grad(p)
+        );
+    U.correctBoundaryConditions();
+
+    // Update divergence free face velocity field, whose value will be used in
+    // the next time step
+    if (!faceUPtr_)
+    {
+        FatalErrorIn
+        (
+            "void solutionControl::reconstructVelocity"
+            "\n("
+            "\n    volVectorField& U,"
+            "\n    const volVectorField& ddtUEqn,"
+            "\n    const volScalarField& rAU,"
+            "\n    const volScalarField& p"
+            "\n) const"
+        )   << "faceUPtr_ not calculated. Make sure you have called"
+            << " calculateTimeConsistentFlux(...) before calling this function."
+            << endl;
+            << exit(FatalError);
+    }
+    surfaceVectorField& faceU = *faceUPtr_;
+
+    // First interpolate the reconstructed velocity on the faces
+    faceU = fvc::interpolate(U);
+
+    // Replace the normal component with conservative flux
+
+    const surfaceVectorField& Sf = mesh_.Sf();
+    const surfaceVectorField rSf = Sf/magSqr(Sf);
+
+    // Subtract interpolated normal component
+    faceU -= (Sf & faceU)*rSf;
+
+    // Now that the normal component is zero, add the normal component from
+    // conservative flux
+    faceU += phi*rSf;
+}
+
+
+const Foam::volScalarField& Foam::solutionControl::aCoeff() const
+{
+    if (!aCoeffPtr_)
+    {
+        FatalErrorIn
+        (
+            "const volScalarField& solutionControl::aCoeff() const"
+        )   << "aCoeffPtr_ not calculated. Make sure you have called"
+            << " calculateTimeConsistentFlux(...) before calling aCoeff()."
+            << endl;
+            << exit(FatalError);
+    }
+
+    return *aCoeffPtr_;
+}
 
 
 // ************************************************************************* //

src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H

@@ -25,7 +25,9 @@ Class
     Foam::solutionControl
 
 Description
-    Base class for solution control classes
+    Base class for solution control classes.
+    The class also provides additional member functions for calculation
+    of time and under-relaxation consistent flux and velocity.
 
 \*---------------------------------------------------------------------------*/
 
@@ -145,6 +147,19 @@ protected:
 
 private:
 
+    // Private data
+
+        // Fields necessary for time and under-relaxation consistency
+
+            //- A coeff (A^~) arising from consistency formulation
+            volScalarField* aCoeffPtr_;
+
+            //- Face velocity needed for consistent formulation
+            surfaceVectorField* faceUPtr_;
+
+
+    // Private member functions
+
         //- Disallow default bitwise copy construct
         solutionControl(const solutionControl&);
 
@@ -204,6 +219,38 @@ public:
             inline bool consistent() const;
 
 
+        // Time and under-relaxation consistency. Note: these functions could be
+        // made virtual and specific to SIMPLE, PISO and PIMPLE, but we will
+        // make them general in sense they account for arbitrary ddt scheme
+        // and arbitrary under-relaxation (e.g. if someone runs SIMPLE in
+        // transient mode or PISO with under-relaxation)
+
+            //- Calculate and return time consistent flux (before pressure
+            //  equation). At point of call, U = H/A where H and A come from
+            //  convection-diffusion equation only.
+            void calculateTimeConsistentFlux
+            (
+                surfaceScalarField& phi,
+                const volVectorField& U,
+                const volScalarField& rAU
+            );
+
+            //- Reconstruct velocity (after pressure equation). At point of
+            //  call, U = H/A where H and A come from convection-diffusion
+            //  equation only
+            void reconstructVelocity
+            (
+                volVectorField& U,
+                const fvVectorMatrix& ddtUEqn,
+                const volScalarField& rAU,
+                const volScalarField& p,
+                const surfaceScalarField& phi
+            );
+
+            //- Const access to aCoeff (needed for pressure equation)
+            const volScalarField& aCoeff() const;
+
+
         // Evolution
 
             //- Main control loop