ENH: changes towards a machine-accurate continuation

Each solver now writes its sensitivity derivatives to its dictionary,
enabling also a binary format. If present, the sensitivities are then
re-read from the dictionary, avoiding thus possible loss of information
due to re-computation.

As a side-effect, sensitivities are computed after the completion of
each adjoint solver, instead of being computed after all adjoint solvers
have been completed.

applications/solvers/incompressible/adjointOptimisationFoam/adjointOptimisationFoam.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-2021 PCOpt/NTUA
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -71,9 +71,6 @@ int main(int argc, char *argv[])
 
         // Solve all adjoint equations
         om.solveAdjointEquations();
-
-        // Compute all sensitivities
-        om.computeSensitivities();
     }
 
     Info<< "End\n" << endl;

src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/adjointSolver/adjointSolver.C

@@ -151,6 +151,16 @@ Foam::objectiveManager& Foam::adjointSolver::getObjectiveManager()
 }
 
 
+void Foam::adjointSolver::postLoop()
+{
+    computeObjectiveSensitivities();
+    // The solver dictionary has been already written after the termination
+    // of the adjoint loop. Force re-writing it to include the sensitivities
+    // as well
+    regIOobject::write(true);
+}
+
+
 bool Foam::adjointSolver::isConstraint()
 {
     return isConstraint_;
@@ -169,4 +179,14 @@ void Foam::adjointSolver::updatePrimalBasedQuantities()
 }
 
 
+bool Foam::adjointSolver::writeData(Ostream& os) const
+{
+    if (sensitivities_.valid())
+    {
+        sensitivities_().writeEntry("sensitivities", os);
+    }
+    return true;
+}
+
+
 // ************************************************************************* //

src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/adjointSolver/adjointSolver.H

@@ -81,7 +81,7 @@ protected:
         autoPtr<objectiveManager> objectiveManagerPtr_;
 
         //- Sensitivities field
-        autoPtr<scalarField> sensitivities_;
+        tmp<scalarField> sensitivities_;
 
         //- Are sensitivities computed
         bool computeSensitivities_;
@@ -172,6 +172,10 @@ public:
 
         // Evolution
 
+            //- Functions to be called after loop
+            //  Computes adjoint sensitivities
+            virtual void postLoop();
+
             //- Compute sensitivities of the underlaying objectives
             virtual void computeObjectiveSensitivities() = 0;
 
@@ -191,6 +195,9 @@ public:
             //- in adjoint turbulence models
             //  Does nothing in the base
             virtual void updatePrimalBasedQuantities();
+
+            //- Write the sensitivity derivatives
+            virtual bool writeData(Ostream& os) const;
 };
 
 

src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/adjointSimple/adjointSimple.C

@@ -151,6 +151,19 @@ Foam::adjointSimple::adjointSimple
                 objectiveManagerPtr_()
             ).ptr()
         );
+        // Read stored sensitivities, if they exist
+        // Need to know the size of the sensitivity field, retrieved after the
+        // allocation of the corresponding object
+        if (dictionary::found("sensitivities"))
+        {
+            sensitivities_ =
+                tmp<scalarField>::New
+                (
+                    "sensitivities",
+                    *this,
+                    adjointSensitivity_().getSensitivities().size()
+                );
+        }
     }
 }
 
@@ -336,6 +349,7 @@ void Foam::adjointSimple::solve()
         {
             solveIter();
         }
+        postLoop();
     }
 }
 
@@ -363,7 +377,7 @@ void Foam::adjointSimple::computeObjectiveSensitivities()
         {
             sensitivities_.reset(new scalarField(sens.size(), Zero));
         }
-        *sensitivities_ = sens;
+        sensitivities_.ref() = sens;
     }
     else
     {
@@ -433,8 +447,7 @@ void Foam::adjointSimple::updatePrimalBasedQuantities()
 bool Foam::adjointSimple::writeData(Ostream& os) const
 {
     os.writeEntry("averageIter", solverControl_().averageIter());
-
-    return true;
+    return adjointSolver::writeData(os);
 }