From ba300c3c6f4c3985750e5c0c9d1bcd7a9ebd22c4 Mon Sep 17 00:00:00 2001
From: Vaggelis Papoutsis <vaggelisp@gmail.com>
Date: Tue, 1 Dec 2020 19:33:35 +0200
Subject: [PATCH] ENH: changed the treatment of fvOptions in primal solvers

fvOptions are no longer a member of incompressiblePrimalSolver but are
looked up from the registry in each iteration of each primal solver.
This means that the main system/fvOptions dictionary is read by ALL
instances of the primal solvers and the latter no longer have their
own fvOptions dict in optimisationDict. This is safe since each fvOption
is applied to a specific field and in case of many primal solvers, the
primal fields are named differently for each of them.

In addition, simple is now split in preLoop, loop and postLoop phase.
Furthermore, each SIMPLE iteration is broken down to
a preIter, mainIter and postIter phase, to allow for different behaviour
by derived classes.
---
 .../incompressiblePrimalSolver.C              | 13 ++-
 .../incompressiblePrimalSolver.H              |  4 -
 .../incompressible/simple/simple.C            | 88 ++++++++++++-------
 .../incompressible/simple/simple.H            | 19 +++-
 4 files changed, 80 insertions(+), 44 deletions(-)

diff --git a/src/optimisation/adjointOptimisation/adjoint/solvers/primalSolvers/incompressible/incompressiblePrimalSolver/incompressiblePrimalSolver.C b/src/optimisation/adjointOptimisation/adjoint/solvers/primalSolvers/incompressible/incompressiblePrimalSolver/incompressiblePrimalSolver.C
index f1f50530e36..3cffe23bfb3 100644
--- a/src/optimisation/adjointOptimisation/adjoint/solvers/primalSolvers/incompressible/incompressiblePrimalSolver/incompressiblePrimalSolver.C
+++ b/src/optimisation/adjointOptimisation/adjoint/solvers/primalSolvers/incompressible/incompressiblePrimalSolver/incompressiblePrimalSolver.C
@@ -30,6 +30,7 @@ License
 #include "incompressiblePrimalSolver.H"
 #include "adjustPhi.H"
 #include "adjointSolver.H"
+#include "fvOptions.H"
 #include "addToRunTimeSelectionTable.H"
 
 
@@ -67,8 +68,7 @@ Foam::incompressiblePrimalSolver::incompressiblePrimalSolver
     (
         dict.subOrEmptyDict("fieldReconstruction").
             getOrDefault<label>("iters", 10)
-    ),
-    fvOptions_(nullptr)
+    )
 {}
 
 
@@ -110,8 +110,6 @@ bool Foam::incompressiblePrimalSolver::readDict(const dictionary& dict)
 {
     if (primalSolver::readDict(dict))
     {
-        fvOptions_().read(dict.subOrEmptyDict("fvOptions"));
-
         return true;
     }
 
@@ -187,6 +185,7 @@ void Foam::incompressiblePrimalSolver::correctBoundaryConditions()
     volVectorField& U = vars.U();
     surfaceScalarField& phi = vars.phi();
     autoPtr<incompressible::turbulenceModel>& turbulence = vars.turbulence();
+    fv::options& fvOptions(fv::options::New(this->mesh_));
 
     scalar contError(GREAT), diff(GREAT);
     for (label iter = 0; iter < phiReconstructionIters_; ++iter)
@@ -200,11 +199,11 @@ void Foam::incompressiblePrimalSolver::correctBoundaryConditions()
             fvm::div(phi, U)
           + turbulence->divDevReff(U)
           ==
-            fvOptions_()(U)
+            fvOptions(U)
         );
         fvVectorMatrix& UEqn = tUEqn.ref();
         UEqn.relax();
-        fvOptions_().constrain(UEqn);
+        fvOptions.constrain(UEqn);
 
         // Pressure equation will give the Rhie-Chow correction
         //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -216,7 +215,7 @@ void Foam::incompressiblePrimalSolver::correctBoundaryConditions()
         surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
         adjustPhi(phiHbyA, U, p);
 
-        //fvOptions_().makeRelative(phiHbyA);
+        //fvOptions.makeRelative(phiHbyA);
 
         fvScalarMatrix pEqn
         (
diff --git a/src/optimisation/adjointOptimisation/adjoint/solvers/primalSolvers/incompressible/incompressiblePrimalSolver/incompressiblePrimalSolver.H b/src/optimisation/adjointOptimisation/adjoint/solvers/primalSolvers/incompressible/incompressiblePrimalSolver/incompressiblePrimalSolver.H
index 016d930d962..9def052a360 100644
--- a/src/optimisation/adjointOptimisation/adjoint/solvers/primalSolvers/incompressible/incompressiblePrimalSolver/incompressiblePrimalSolver.H
+++ b/src/optimisation/adjointOptimisation/adjoint/solvers/primalSolvers/incompressible/incompressiblePrimalSolver/incompressiblePrimalSolver.H
@@ -39,7 +39,6 @@ Description
 
 #include "primalSolver.H"
 #include "incompressibleVars.H"
-#include "fvOptionList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -77,9 +76,6 @@ protected:
         //- Max iterations  for reconstructing phi from U and p
         label phiReconstructionIters_;
 
-        //- optionList for source terms addition
-        autoPtr<fv::optionList> fvOptions_;
-
 
 public:
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/solvers/primalSolvers/incompressible/simple/simple.C b/src/optimisation/adjointOptimisation/adjoint/solvers/primalSolvers/incompressible/simple/simple.C
index 5919cb373de..90a33437094 100644
--- a/src/optimisation/adjointOptimisation/adjoint/solvers/primalSolvers/incompressible/simple/simple.C
+++ b/src/optimisation/adjointOptimisation/adjoint/solvers/primalSolvers/incompressible/simple/simple.C
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -33,6 +33,7 @@ License
 #include "constrainPressure.H"
 #include "adjustPhi.H"
 #include "Time.H"
+#include "fvOptions.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -101,10 +102,6 @@ Foam::simple::simple
     cumulativeContErr_(Zero),
     objectives_(0)
 {
-    fvOptions_.reset
-    (
-        new fv::optionList(mesh_, dict.subOrEmptyDict("fvOptions"))
-    );
     addExtraSchemes();
     setRefCell
     (
@@ -130,9 +127,20 @@ bool Foam::simple::readDict(const dictionary& dict)
 
 void Foam::simple::solveIter()
 {
-    const Time& time = mesh_.time();
-    Info<< "Time = " << time.timeName() << "\n" << endl;
+    preIter();
+    mainIter();
+    postIter();
+}
 
+
+void Foam::simple::preIter()
+{
+    Info<< "Time = " << mesh_.time().timeName() << "\n" << endl;
+}
+
+
+void Foam::simple::mainIter()
+{
     // Grab references
     volScalarField& p = incoVars_.pInst();
     volVectorField& U = incoVars_.UInst();
@@ -141,6 +149,7 @@ void Foam::simple::solveIter()
         incoVars_.turbulence();
     label&  pRefCell  = solverControl_().pRefCell();
     scalar& pRefValue = solverControl_().pRefValue();
+    fv::options& fvOptions(fv::options::New(this->mesh_));
 
     // Momentum predictor
     //~~~~~~~~~~~~~~~~~~~
@@ -153,21 +162,19 @@ void Foam::simple::solveIter()
       + MRF_.DDt(U)
       + turbulence->divDevReff(U)
       ==
-        fvOptions_()(U)
+        fvOptions(U)
     );
     fvVectorMatrix& UEqn = tUEqn.ref();
 
-    addOptimisationTypeSource(UEqn);
-
     UEqn.relax();
 
-    fvOptions_().constrain(UEqn);
+    fvOptions.constrain(UEqn);
 
     if (solverControl_().momentumPredictor())
     {
         Foam::solve(UEqn == -fvc::grad(p));
 
-        fvOptions_().correct(U);
+        fvOptions.correct(U);
     }
 
     // Pressure Eq
@@ -220,12 +227,16 @@ void Foam::simple::solveIter()
         // Momentum corrector
         U = HbyA - rAtU()*fvc::grad(p);
         U.correctBoundaryConditions();
-        fvOptions_().correct(U);
+        fvOptions.correct(U);
     }
 
     incoVars_.laminarTransport().correct();
     turbulence->correct();
+}
+
 
+void Foam::simple::postIter()
+{
     solverControl_().write();
 
     // Print objective values to screen and compute mean value
@@ -241,7 +252,7 @@ void Foam::simple::solveIter()
     incoVars_.computeMeanFields();
 
     // Print execution time
-    time.printExecutionTime(Info);
+    mesh_.time().printExecutionTime(Info);
 }
 
 
@@ -250,26 +261,12 @@ void Foam::simple::solve()
     // Iterate
     if (active_)
     {
-        // Get the objectives for this solver
-        if (objectives_.empty())
-        {
-            objectives_ = getObjectiveFunctions();
-        }
-
-        // Reset initial and mean fields before solving
-        restoreInitValues();
-        incoVars_.resetMeanFields();
-
-        // Validate turbulence model fields
-        incoVars_.turbulence()->validate();
-
+        preLoop();
         while (solverControl_().loop())
         {
             solveIter();
         }
-
-        // Safety
-        objectives_.clear();
+        postLoop();
     }
 }
 
@@ -286,6 +283,35 @@ void Foam::simple::restoreInitValues()
 }
 
 
+void Foam::simple::preLoop()
+{
+    // Get the objectives for this solver
+    if (objectives_.empty())
+    {
+        objectives_ = getObjectiveFunctions();
+    }
+
+    // Reset initial and mean fields before solving
+    restoreInitValues();
+    incoVars_.resetMeanFields();
+
+    // Validate turbulence model fields
+    incoVars_.turbulence()->validate();
+}
+
+
+void Foam::simple::postLoop()
+{
+    for (objective* obj : objectives_)
+    {
+        obj->writeInstantaneousSeparator();
+    }
+
+    // Safety
+    objectives_.clear();
+}
+
+
 bool Foam::simple::writeData(Ostream& os) const
 {
     os.writeEntry("averageIter", solverControl_().averageIter());
diff --git a/src/optimisation/adjointOptimisation/adjoint/solvers/primalSolvers/incompressible/simple/simple.H b/src/optimisation/adjointOptimisation/adjoint/solvers/primalSolvers/incompressible/simple/simple.H
index bc384375a16..335df8280fc 100644
--- a/src/optimisation/adjointOptimisation/adjoint/solvers/primalSolvers/incompressible/simple/simple.H
+++ b/src/optimisation/adjointOptimisation/adjoint/solvers/primalSolvers/incompressible/simple/simple.H
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -134,6 +134,15 @@ public:
             //- Execute one iteration of the solution algorithm
             virtual void solveIter();
 
+            //- Steps to be executed before each main SIMPLE iteration
+            virtual void preIter();
+
+            //- The main SIMPLE iter
+            virtual void mainIter();
+
+            //- Steps to be executed before each main SIMPLE iteration
+            virtual void postIter();
+
             //- Main control loop
             virtual void solve();
 
@@ -143,6 +152,12 @@ public:
             //- Restore initial field values if necessary
             virtual void restoreInitValues();
 
+            //- Functions to be called before loop
+            virtual void preLoop();
+
+            //- Functions to be called after loop
+            virtual void postLoop();
+
             //- Write average iteration
             virtual bool writeData(Ostream& os) const;
 };
-- 
GitLab