ENH: residuals function object - extended to write residual fields

Residual fields can be written using the new 'writeFields' entry, e.g.

    functions
    {
        residual
        {
            type            residuals;
            libs            ("libutilityFunctionObjects.so");
            fields          (".*");
            writeControl    writeTime;
            writeFields     true;
        }
    }

Fields currently correspond to the initial residual for the last solver
iteration.

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C

@@ -26,6 +26,9 @@ License
 #include "lduMatrix.H"
 #include "IOstreams.H"
 #include "Switch.H"
+#include "objectRegistry.H"
+#include "IOField.H"
+#include "Time.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -312,6 +315,38 @@ const Foam::scalarField& Foam::lduMatrix::upper() const
 }
 
 
+void Foam::lduMatrix::setResidualField
+(
+    const Field<scalar>& residual,
+    const word& fieldName,
+    const bool initial
+) const
+{
+    if (!lduMesh_.hasDb())
+    {
+        return;
+    }
+
+    word lookupName;
+    if (initial)
+    {
+        lookupName = word("initialResidual:" + fieldName);
+    }
+    else
+    {
+        lookupName = word("residual:" + fieldName);
+    }
+
+    IOField<scalar>* residualPtr =
+        lduMesh_.thisDb().lookupObjectRefPtr<IOField<scalar>>(lookupName);
+
+    if (residualPtr)
+    {
+        *residualPtr = residual;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
 Foam::Ostream& Foam::operator<<(Ostream& os, const lduMatrix& ldum)

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H

@@ -124,6 +124,7 @@ public:
 
             profilingTrigger profiling_;
 
+
         // Protected Member Functions
 
             //- Read the control parameters from the controlDict_
@@ -258,7 +259,7 @@ public:
             ) const = 0;
 
             //- Return the matrix norm used to normalise the residual for the
-            //  stopping criterion
+            //- stopping criterion
             scalar normFactor
             (
                 const scalarField& psi,
@@ -485,7 +486,7 @@ public:
         // Member functions
 
             //- Read and reset the preconditioner parameters
-            //  from the given stream
+            //- from the given stream
             virtual void read(const dictionary&)
             {}
 
@@ -498,7 +499,7 @@ public:
             ) const = 0;
 
             //- Return wT the transpose-matrix preconditioned form of
-            //  residual rT.
+            //- residual rT.
             //  This is only required for preconditioning asymmetric matrices.
             virtual void preconditionT
             (
@@ -531,7 +532,7 @@ public:
         lduMatrix(lduMatrix&, bool reuse);
 
         //- Construct given an LDU addressed mesh and an Istream
-        //  from which the coefficients are read
+        //- from which the coefficients are read
         lduMatrix(const lduMesh&, Istream&);
 
 
@@ -669,7 +670,7 @@ public:
 
 
             //- Initialise the update of interfaced interfaces
-            //  for matrix operations
+            //- for matrix operations
             void initMatrixInterfaces
             (
                 const bool add,
@@ -691,6 +692,14 @@ public:
                 const direction cmpt
             ) const;
 
+            //- Set the residual field using an IOField on the object registry
+            //- if it exists
+            void setResidualField
+            (
+                const Field<scalar>& residual,
+                const word& fieldName,
+                const bool initial
+            ) const;
 
             template<class Type>
             tmp<Field<Type>> H(const Field<Type>&) const;

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C

@@ -59,6 +59,8 @@ Foam::solverPerformance Foam::GAMGSolver::solve
     // Calculate initial finest-grid residual field
     scalarField finestResidual(source - Apsi);
 
+    matrix().setResidualField(finestResidual, fieldName_, true);
+
     // Calculate normalised residual for convergence test
     solverPerf.initialResidual() = gSumMag
     (
@@ -143,6 +145,8 @@ Foam::solverPerformance Foam::GAMGSolver::solve
         );
     }
 
+    matrix().setResidualField(finestResidual, fieldName_, false);
+
     return solverPerf;
 }
 

src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C

@@ -93,6 +93,8 @@ Foam::solverPerformance Foam::PBiCG::solve
     scalarField rA(source - wA);
     scalar* __restrict__ rAPtr = rA.begin();
 
+    matrix().setResidualField(rA, fieldName_, true);
+
     // --- Calculate normalisation factor
     const scalar normFactor = this->normFactor(psi, source, wA, pA);
 
@@ -218,6 +220,8 @@ Foam::solverPerformance Foam::PBiCG::solve
             << exit(FatalError);
     }
 
+    matrix().setResidualField(rA, fieldName_, false);
+
     return solverPerf;
 }
 

src/OpenFOAM/matrices/lduMatrix/solvers/PBiCGStab/PBiCGStab.C

@@ -96,6 +96,8 @@ Foam::solverPerformance Foam::PBiCGStab::solve
     scalarField rA(source - yA);
     scalar* __restrict__ rAPtr = rA.begin();
 
+    matrix().setResidualField(rA, fieldName_, true);
+
     // --- Calculate normalisation factor
     const scalar normFactor = this->normFactor(psi, source, yA, pA);
 
@@ -248,6 +250,8 @@ Foam::solverPerformance Foam::PBiCGStab::solve
         );
     }
 
+    matrix().setResidualField(rA, fieldName_, false);
+
     return solverPerf;
 }
 

src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C

@@ -96,6 +96,8 @@ Foam::solverPerformance Foam::PCG::solve
     scalarField rA(source - wA);
     scalar* __restrict__ rAPtr = rA.begin();
 
+    matrix().setResidualField(rA, fieldName_, true);
+
     // --- Calculate normalisation factor
     scalar normFactor = this->normFactor(psi, source, wA, pA);
 
@@ -119,11 +121,11 @@ Foam::solverPerformance Foam::PCG::solve
     {
         // --- Select and construct the preconditioner
         autoPtr<lduMatrix::preconditioner> preconPtr =
-        lduMatrix::preconditioner::New
+            lduMatrix::preconditioner::New
-        (
+            (
-            *this,
+                *this,
-            controlDict_
+                controlDict_
-        );
+            );
 
         // --- Solver iteration
         do
@@ -189,6 +191,8 @@ Foam::solverPerformance Foam::PCG::solve
         );
     }
 
+    matrix().setResidualField(rA, fieldName_, false);
+
     return solverPerf;
 }
 

src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.C

@@ -113,6 +113,7 @@ Foam::solverPerformance Foam::smoothSolver::solve
     else
     {
         scalar normFactor = 0;
+        scalarField residual;
 
         {
             scalarField Apsi(psi.size());
@@ -124,12 +125,13 @@ Foam::solverPerformance Foam::smoothSolver::solve
             // Calculate normalisation factor
             normFactor = this->normFactor(psi, source, Apsi, temp);
 
+            residual = source - Apsi;
+
+            matrix().setResidualField(residual, fieldName_, true);
+
             // Calculate residual magnitude
-            solverPerf.initialResidual() = gSumMag
+            solverPerf.initialResidual() =
-            (
+                gSumMag(residual, matrix().mesh().comm())/normFactor;
-                (source - Apsi)(),
-                matrix().mesh().comm()
-            )/normFactor;
             solverPerf.finalResidual() = solverPerf.initialResidual();
         }
 
@@ -170,9 +172,7 @@ Foam::solverPerformance Foam::smoothSolver::solve
                     nSweeps_
                 );
 
-                // Calculate the residual to check convergence
+                residual =
-                solverPerf.finalResidual() = gSumMag
-                (
                     matrix_.residual
                     (
                         psi,
@@ -180,9 +180,11 @@ Foam::solverPerformance Foam::smoothSolver::solve
                         interfaceBouCoeffs_,
                         interfaces_,
                         cmpt
-                    )(),
+                    );
-                    matrix().mesh().comm()
+
-                )/normFactor;
+                // Calculate the residual to check convergence
+                solverPerf.finalResidual() =
+                    gSumMag(residual, matrix().mesh().comm())/normFactor;
             } while
             (
                 (
@@ -192,6 +194,8 @@ Foam::solverPerformance Foam::smoothSolver::solve
              || solverPerf.nIterations() < minIter_
             );
         }
+
+        matrix().setResidualField(residual, fieldName_, false);
     }
 
     return solverPerf;

src/OpenFOAM/meshes/GeoMesh/GeoMesh.H

@@ -74,6 +74,12 @@ public:
 
     // Member Functions
 
+        //- Return true if thisDb() is a valid DB - here = false
+        bool hasDb() const
+        {
+            return true;
+        }
+
         //- Return the object registry
         const objectRegistry& thisDb() const
         {

src/OpenFOAM/meshes/lduMesh/lduMesh.H

@@ -77,6 +77,9 @@ public:
 
         // Access
 
+            //- Return true if thisDb() is a valid DB
+            virtual bool hasDb() const = 0;
+
             //- Return the object registry
             virtual const objectRegistry& thisDb() const;
 

src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H

@@ -73,6 +73,7 @@ class lduPrimitiveMesh
         //- Communicator to use for any parallel communication
         const label comm_;
 
+
     // Private Member Functions
 
         //- Get size of all meshes
@@ -92,6 +93,7 @@ public:
         // Declare name of the class and its debug switch
         ClassName("lduPrimitiveMesh");
 
+
     // Constructors
 
         //- Construct from components but without interfaces. Add interfaces
@@ -170,6 +172,12 @@ public:
 
         // Access
 
+            //- Return true if thisDb() is a valid DB
+            virtual bool hasDb() const
+            {
+                return false;
+            }
+
             //- Return ldu addressing
             virtual const lduAddressing& lduAddr() const
             {

src/finiteVolume/fvMesh/fvMesh.H

@@ -239,6 +239,12 @@ public:
                 return polyMesh::time();
             }
 
+            //- Return true if thisDb() is a valid DB
+            virtual bool hasDb() const
+            {
+                return true;
+            }
+
             //- Return the object registry - resolve conflict polyMesh/lduMesh
             virtual const objectRegistry& thisDb() const
             {

src/functionObjects/utilities/residuals/residuals.C

@@ -79,6 +79,69 @@ void Foam::functionObjects::residuals::writeFileHeader(Ostream& os)
 }
 
 
+void Foam::functionObjects::residuals::createField(const word& fieldName)
+{
+    if (!writeFields_)
+    {
+        return;
+    }
+
+    const word residualName("initialResidual:" + fieldName);
+
+    if (!mesh_.foundObject<IOField<scalar>>(residualName))
+    {
+        IOField<scalar>* fieldPtr =
+            new IOField<scalar>
+            (
+                IOobject
+                (
+                    residualName,
+                    mesh_.time().timeName(),
+                    mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                Field<scalar>(mesh_.nCells(), scalar(0))
+            );
+
+        fieldPtr->store();
+    }
+}
+
+
+void Foam::functionObjects::residuals::writeField(const word& fieldName) const
+{
+    const word residualName("initialResidual:" + fieldName);
+
+    const IOField<scalar>* residualPtr =
+        mesh_.lookupObjectPtr<IOField<scalar>>(residualName);
+
+    if (residualPtr)
+    {
+        volScalarField residual
+        (
+            IOobject
+            (
+                residualName,
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("0", dimless, 0),
+            zeroGradientFvPatchField<scalar>::typeName
+        );
+
+        residual.primitiveFieldRef() = *residualPtr;
+        residual.correctBoundaryConditions();
+
+        residual.write();
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::functionObjects::residuals::residuals
@@ -90,7 +153,9 @@ Foam::functionObjects::residuals::residuals
 :
     fvMeshFunctionObject(name, runTime, dict),
     writeFile(obr_, name, typeName, dict),
-    fieldSet_(mesh_)
+    fieldSet_(mesh_),
+    writeFields_(false),
+    initialised_(false)
 {
     read(dict);
 }
@@ -109,6 +174,9 @@ bool Foam::functionObjects::residuals::read(const dictionary& dict)
     if (fvMeshFunctionObject::read(dict))
     {
         fieldSet_.read(dict);
+
+        writeFields_ = dict.lookupOrDefault("writeFields", false);
+
         return true;
     }
 
@@ -118,32 +186,48 @@ bool Foam::functionObjects::residuals::read(const dictionary& dict)
 
 bool Foam::functionObjects::residuals::execute()
 {
-    return true;
+    // Note: delaying initialisation until after first iteration so that
-}
+    // we can find wildcard fields
-
+    if (!initialised_)
-
-bool Foam::functionObjects::residuals::write()
-{
-    if (Pstream::master())
     {
         writeFileHeader(file());
 
-        writeTime(file());
+        if (writeFields_)
-
-        for (const word& fieldName : fieldSet_.selection())
         {
-            writeResidual<scalar>(fieldName);
+            for (const word& fieldName : fieldSet_.selection())
-            writeResidual<vector>(fieldName);
+            {
-            writeResidual<sphericalTensor>(fieldName);
+                initialiseField<scalar>(fieldName);
-            writeResidual<symmTensor>(fieldName);
+                initialiseField<vector>(fieldName);
-            writeResidual<tensor>(fieldName);
+                initialiseField<sphericalTensor>(fieldName);
+                initialiseField<symmTensor>(fieldName);
+                initialiseField<tensor>(fieldName);
+            }
         }
 
-        file() << endl;
+        initialised_ = true;
     }
 
     return true;
 }
 
 
+bool Foam::functionObjects::residuals::write()
+{
+    writeTime(file());
+
+    for (const word& fieldName : fieldSet_.selection())
+    {
+        writeResidual<scalar>(fieldName);
+        writeResidual<vector>(fieldName);
+        writeResidual<sphericalTensor>(fieldName);
+        writeResidual<symmTensor>(fieldName);
+        writeResidual<tensor>(fieldName);
+    }
+
+    file() << endl;
+
+    return true;
+}
+
+
 // ************************************************************************* //

src/functionObjects/utilities/residuals/residuals.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -87,17 +87,33 @@ protected:
         //- Fields to write residuals
         solverFieldSelection fieldSet_;
 
+        //- Flag to write the residual as a vol field
+        bool writeFields_;
+
+        //- Initialisation flag
+        bool initialised_;
+
 
     // Protected Member Functions
 
         //- Output file header information
         void writeFileHeader(Ostream& os);
 
+        //- Create and store a residual field on the mesh database
+        void createField(const word& fieldName);
+
+        //- Write a residual field
+        void writeField(const word& fieldName) const;
+
         //- Output file header information per primitive type value
         template<class Type>
         void writeFileHeader(Ostream& os, const word& fileName) const;
 
-        //- Calculate the field min/max
+        //- Initialise a residual field
+        template<class Type>
+        void initialiseField(const word& fieldName);
+
+        //- Calculate the residual
         template<class Type>
         void writeResidual(const word& fieldName);
 

src/functionObjects/utilities/residuals/residualsTemplates.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,6 +26,7 @@ License
 #include "residuals.H"
 #include "volFields.H"
 #include "ListOps.H"
+#include "zeroGradientFvPatchField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -60,12 +61,43 @@ void Foam::functionObjects::residuals::writeFileHeader
 }
 
 
+template<class Type>
+void Foam::functionObjects::residuals::initialiseField(const word& fieldName)
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> volFieldType;
+
+    if (foundObject<volFieldType>(fieldName))
+    {
+        const Foam::dictionary& solverDict = mesh_.solverPerformanceDict();
+
+        if (solverDict.found(fieldName))
+        {
+            typename pTraits<Type>::labelType validComponents
+            (
+                mesh_.validComponents<Type>()
+            );
+
+            for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; ++cmpt)
+            {
+                if (component(validComponents, cmpt) != -1)
+                {
+                    const word resultName =
+                        fieldName + word(pTraits<Type>::componentNames[cmpt]);
+
+                    createField(resultName);
+                }
+            }
+        }
+    }
+}
+
+
 template<class Type>
 void Foam::functionObjects::residuals::writeResidual(const word& fieldName)
 {
-    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+    typedef GeometricField<Type, fvPatchField, volMesh> volFieldType;
 
-    if (foundObject<fieldType>(fieldName))
+    if (foundObject<volFieldType>(fieldName))
     {
         const Foam::dictionary& solverDict = mesh_.solverPerformanceDict();
 
@@ -83,7 +115,7 @@ void Foam::functionObjects::residuals::writeResidual(const word& fieldName)
                 mesh_.validComponents<Type>()
             );
 
-            for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
+            for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; ++cmpt)
             {
                 if (component(validComponents, cmpt) != -1)
                 {
@@ -94,6 +126,8 @@ void Foam::functionObjects::residuals::writeResidual(const word& fieldName)
                     const word resultName =
                         fieldName + word(pTraits<Type>::componentNames[cmpt]);
                     setResult(resultName, r);
+
+                    writeField(resultName);
                 }
             }
         }