diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C index 922f9b8c59..6c9c92dad5 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C +++ b/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& 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* residualPtr = + lduMesh_.thisDb().lookupObjectRefPtr>(lookupName); + + if (residualPtr) + { + *residualPtr = residual; + } +} + + // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // Foam::Ostream& Foam::operator<<(Ostream& os, const lduMatrix& ldum) diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H index 38301514e9..e27782842a 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H +++ b/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& residual, + const word& fieldName, + const bool initial + ) const; template tmp> H(const Field&) const; diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C index dd23a802b2..58fc1dd1b1 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C +++ b/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; } diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C b/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C index f0e4a3feb2..5bc3a79dad 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C +++ b/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; } diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCGStab/PBiCGStab.C b/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCGStab/PBiCGStab.C index a189de6ecb..b7c1098afe 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCGStab/PBiCGStab.C +++ b/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; } diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C b/src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C index 283fd7c3df..6e70709670 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C +++ b/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 preconPtr = - lduMatrix::preconditioner::New - ( - *this, - controlDict_ - ); + lduMatrix::preconditioner::New + ( + *this, + controlDict_ + ); // --- Solver iteration do @@ -189,6 +191,8 @@ Foam::solverPerformance Foam::PCG::solve ); } + matrix().setResidualField(rA, fieldName_, false); + return solverPerf; } diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.C b/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.C index 31eb37add3..777fc5210e 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.C +++ b/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 - ( - (source - Apsi)(), - matrix().mesh().comm() - )/normFactor; + solverPerf.initialResidual() = + gSumMag(residual, matrix().mesh().comm())/normFactor; solverPerf.finalResidual() = solverPerf.initialResidual(); } @@ -170,9 +172,7 @@ Foam::solverPerformance Foam::smoothSolver::solve nSweeps_ ); - // Calculate the residual to check convergence - solverPerf.finalResidual() = gSumMag - ( + residual = 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; diff --git a/src/OpenFOAM/meshes/GeoMesh/GeoMesh.H b/src/OpenFOAM/meshes/GeoMesh/GeoMesh.H index a477876113..bc702273a0 100644 --- a/src/OpenFOAM/meshes/GeoMesh/GeoMesh.H +++ b/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 { diff --git a/src/OpenFOAM/meshes/lduMesh/lduMesh.H b/src/OpenFOAM/meshes/lduMesh/lduMesh.H index b871ee9566..04057b3902 100644 --- a/src/OpenFOAM/meshes/lduMesh/lduMesh.H +++ b/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; diff --git a/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H b/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H index 12a60773dc..d57c4afa46 100644 --- a/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H +++ b/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 { diff --git a/src/finiteVolume/fvMesh/fvMesh.H b/src/finiteVolume/fvMesh/fvMesh.H index a2bd7e403c..0504b59fc8 100644 --- a/src/finiteVolume/fvMesh/fvMesh.H +++ b/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 { diff --git a/src/functionObjects/utilities/residuals/residuals.C b/src/functionObjects/utilities/residuals/residuals.C index 0cdf2c6d97..0c20bb61c0 100644 --- a/src/functionObjects/utilities/residuals/residuals.C +++ b/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>(residualName)) + { + IOField* fieldPtr = + new IOField + ( + IOobject + ( + residualName, + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + Field(mesh_.nCells(), scalar(0)) + ); + + fieldPtr->store(); + } +} + + +void Foam::functionObjects::residuals::writeField(const word& fieldName) const +{ + const word residualName("initialResidual:" + fieldName); + + const IOField* residualPtr = + mesh_.lookupObjectPtr>(residualName); + + if (residualPtr) + { + volScalarField residual + ( + IOobject + ( + residualName, + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh_, + dimensionedScalar("0", dimless, 0), + zeroGradientFvPatchField::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; -} - - -bool Foam::functionObjects::residuals::write() -{ - if (Pstream::master()) + // Note: delaying initialisation until after first iteration so that + // we can find wildcard fields + if (!initialised_) { writeFileHeader(file()); - writeTime(file()); - - for (const word& fieldName : fieldSet_.selection()) + if (writeFields_) { - writeResidual(fieldName); - writeResidual(fieldName); - writeResidual(fieldName); - writeResidual(fieldName); - writeResidual(fieldName); + for (const word& fieldName : fieldSet_.selection()) + { + initialiseField(fieldName); + initialiseField(fieldName); + initialiseField(fieldName); + initialiseField(fieldName); + initialiseField(fieldName); + } } - file() << endl; + initialised_ = true; } return true; } +bool Foam::functionObjects::residuals::write() +{ + writeTime(file()); + + for (const word& fieldName : fieldSet_.selection()) + { + writeResidual(fieldName); + writeResidual(fieldName); + writeResidual(fieldName); + writeResidual(fieldName); + writeResidual(fieldName); + } + + file() << endl; + + return true; +} + + // ************************************************************************* // diff --git a/src/functionObjects/utilities/residuals/residuals.H b/src/functionObjects/utilities/residuals/residuals.H index 63b43e5fc9..d6c8ec9ed2 100644 --- a/src/functionObjects/utilities/residuals/residuals.H +++ b/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 void writeFileHeader(Ostream& os, const word& fileName) const; - //- Calculate the field min/max + //- Initialise a residual field + template + void initialiseField(const word& fieldName); + + //- Calculate the residual template void writeResidual(const word& fieldName); diff --git a/src/functionObjects/utilities/residuals/residualsTemplates.C b/src/functionObjects/utilities/residuals/residualsTemplates.C index 322b6bc52a..43c26a7f03 100644 --- a/src/functionObjects/utilities/residuals/residualsTemplates.C +++ b/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 +void Foam::functionObjects::residuals::initialiseField(const word& fieldName) +{ + typedef GeometricField volFieldType; + + if (foundObject(fieldName)) + { + const Foam::dictionary& solverDict = mesh_.solverPerformanceDict(); + + if (solverDict.found(fieldName)) + { + typename pTraits::labelType validComponents + ( + mesh_.validComponents() + ); + + for (direction cmpt=0; cmpt::nComponents; ++cmpt) + { + if (component(validComponents, cmpt) != -1) + { + const word resultName = + fieldName + word(pTraits::componentNames[cmpt]); + + createField(resultName); + } + } + } + } +} + + template void Foam::functionObjects::residuals::writeResidual(const word& fieldName) { - typedef GeometricField fieldType; + typedef GeometricField volFieldType; - if (foundObject(fieldName)) + if (foundObject(fieldName)) { const Foam::dictionary& solverDict = mesh_.solverPerformanceDict(); @@ -83,7 +115,7 @@ void Foam::functionObjects::residuals::writeResidual(const word& fieldName) mesh_.validComponents() ); - for (direction cmpt=0; cmpt::nComponents; cmpt++) + for (direction cmpt=0; cmpt::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::componentNames[cmpt]); setResult(resultName, r); + + writeField(resultName); } } }