ENH: refactoring and cleaning of optimisationType

Moved part common to all derived classes (e.g. update) to the base
class to avoid code duplication. Practically, only the protected
updateDesignVariables has to be overwritten in each derived class now.
steadyOptimisation was also affected in a minor way.

src/optimisation/adjointOptimisation/adjoint/optimisation/optimisationManager/singleRun/singleRun.H

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2020 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -107,8 +107,7 @@ public:
         virtual bool update();
 
         //- Update design variables.
-        //  Might employ a line search to find a correction satisfying the step
+        //  Does nothing
-        //  convergence criteria
         virtual void updateDesignVariables();
 };
 

src/optimisation/adjointOptimisation/adjoint/optimisation/optimisationManager/steadyOptimisation/steadyOptimisation.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2020 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -68,8 +68,13 @@ void Foam::steadyOptimisation::updateOptTypeSource()
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-void Foam::steadyOptimisation::lineSearchUpdate(scalarField& direction)
+void Foam::steadyOptimisation::lineSearchUpdate()
 {
+    // Compute direction of update
+    tmp<scalarField> tdirection = optType_->computeDirection();
+    scalarField& direction = tdirection.ref();
+
+    // Grab reference to line search
     autoPtr<lineSearch>& lineSrch = optType_->getLineSearch();
 
     // Store starting point
@@ -121,7 +126,7 @@ void Foam::steadyOptimisation::lineSearchUpdate(scalarField& direction)
         else
         {
             // If maximum number of iteration has been reached, continue
-            if (iter == lineSrch->maxIters()-1)
+            if (iter == lineSrch->maxIters() - 1)
             {
                 Info<< "Line search reached max. number of iterations.\n"
                     << "Proceeding to the next optimisation cycle" << endl;
@@ -141,15 +146,10 @@ void Foam::steadyOptimisation::lineSearchUpdate(scalarField& direction)
 }
 
 
-void Foam::steadyOptimisation::fixedStepUpdate(scalarField& direction)
+void Foam::steadyOptimisation::fixedStepUpdate()
 {
-    // Update based on fixed step
+    // Update design variables
-    optType_->update(direction);
+    optType_->update();
-
-    // If direction has been scaled (say by setting the initial eta), the
-    // old correction has to be updated
-    optType_->updateOldCorrection(direction);
-    optType_->write();
 
     // Solve primal equations
     solvePrimalEquations();
@@ -223,20 +223,15 @@ bool Foam::steadyOptimisation::update()
 
 void Foam::steadyOptimisation::updateDesignVariables()
 {
-    // Compute direction of update
-    tmp<scalarField> tdirection = optType_->computeDirection();
-    scalarField& direction = tdirection.ref();
-    autoPtr<lineSearch>& lineSrch = optType_->getLineSearch();
-
     // Update design variables using either a line-search scheme or
     // a fixed-step update
-    if (lineSrch.valid())
+    if (optType_->getLineSearch().valid())
     {
-        lineSearchUpdate(direction);
+        lineSearchUpdate();
     }
     else
     {
-        fixedStepUpdate(direction);
+        fixedStepUpdate();
     }
 
     // Reset adjoint sensitivities in all adjoint solver managers

src/optimisation/adjointOptimisation/adjoint/optimisation/optimisationManager/steadyOptimisation/steadyOptimisation.H

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2020 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -75,10 +75,10 @@ protected:
     // Protected Member Functions
 
         //- Update design variables using a line-search
-        void lineSearchUpdate(scalarField& direction);
+        void lineSearchUpdate();
 
         //- Update design variables using a fixed step
-        void fixedStepUpdate(scalarField& direction);
+        void fixedStepUpdate();
 
 
 public:

src/optimisation/adjointOptimisation/adjoint/optimisation/optimisationType/incompressible/optimisationType/optimisationTypeIncompressible.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2020 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -157,14 +157,79 @@ autoPtr<optimisationType> optimisationType::New
 
 // * * * * * * * * * * * * * * *  Member Functions   * * * * * * * * * * * * //
 
+void optimisationType::update()
+{
+    // Compute update of the design variables
+    tmp<scalarField> tcorrection(computeDirection());
+    scalarField& correction = tcorrection.ref();
+
+    // Update design variables given the correction
+    update(correction);
+
+    // If direction has been scaled (say by setting the initial eta), the
+    // old correction has to be updated
+    updateOldCorrection(correction);
+    write();
+}
+
+
+void optimisationType::update(scalarField& direction)
+{
+    // Compute eta if needed
+    computeEta(direction);
+
+    // Multiply with line search step, if necessary
+    scalarField correction(direction);
+    if (lineSearch_.valid())
+    {
+        correction *= lineSearch_->step();
+    }
+
+    // Update the design variables
+    updateDesignVariables(correction);
+}
+
+
 tmp<scalarField> optimisationType::computeDirection()
 {
-    // Sum contributions
+    // Sum contributions for sensitivities and objective/constraint values
     scalarField objectiveSens;
     PtrList<scalarField> constraintSens;
     scalar objectiveValue(Zero);
     scalarField constraintValues;
 
+    updateGradientsAndValues
+    (
+        objectiveSens,
+        constraintSens,
+        objectiveValue,
+        constraintValues
+    );
+
+    // Based on the sensitivities, return design variables correction
+    updateMethod_->setObjectiveDeriv(objectiveSens);
+    updateMethod_->setConstraintDeriv(constraintSens);
+    updateMethod_->setObjectiveValue(objectiveValue);
+    updateMethod_->setConstraintValues(constraintValues);
+    tmp<scalarField> tcorrection
+    (
+        new scalarField(objectiveSens.size(), Zero)
+    );
+    scalarField& correction = tcorrection.ref();
+    correction = updateMethod_->returnCorrection();
+
+    return tcorrection;
+}
+
+
+void optimisationType::updateGradientsAndValues
+(
+    scalarField& objectiveSens,
+    PtrList<scalarField>& constraintSens,
+    scalar& objectiveValue,
+    scalarField& constraintValues
+)
+{
     for (adjointSolverManager& adjSolvManager : adjointSolvManagers_)
     {
         const scalar opWeight = adjSolvManager.operatingPointWeight();
@@ -212,20 +277,6 @@ tmp<scalarField> optimisationType::computeDirection()
         }
         constraintValues += opWeight*cValues();
     }
-
-    // Based on the sensitivities, return design variables correction
-    updateMethod_->setObjectiveDeriv(objectiveSens);
-    updateMethod_->setConstraintDeriv(constraintSens);
-    updateMethod_->setObjectiveValue(objectiveValue);
-    updateMethod_->setConstraintValues(constraintValues);
-    tmp<scalarField> tcorrection
-    (
-        new scalarField(objectiveSens.size(), Zero)
-    );
-    scalarField& correction = tcorrection.ref();
-    correction = updateMethod_->returnCorrection();
-
-    return tcorrection;
 }
 
 

src/optimisation/adjointOptimisation/adjoint/optimisation/optimisationType/incompressible/optimisationType/optimisationTypeIncompressible.H

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2020 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -69,7 +69,11 @@ protected:
         autoPtr<volScalarField> sourcePtr_;
         autoPtr<lineSearch> lineSearch_;
 
-        virtual void computeEta(scalarField& correction)=0;
+        //- Update the design variables given their correction
+        virtual void updateDesignVariables(scalarField& correction) = 0;
+
+        //- Compute eta if not set in the first step
+        virtual void computeEta(scalarField& correction) = 0;
 
 
 private:
@@ -132,10 +136,10 @@ public:
         virtual ~optimisationType() = default;
 
         //- Update design variables
-        virtual void update() = 0;
+        virtual void update();
 
         //- Update design variables based on a given correction
-        virtual void update(scalarField& correction) = 0;
+        virtual void update(scalarField& correction);
 
         //- Store design variables, as the starting point for line search
         virtual void storeDesignVariables() = 0;
@@ -146,6 +150,15 @@ public:
         //- Compute update direction
         virtual tmp<scalarField> computeDirection();
 
+        //- Compute cumulative objective and constraint gradients
+        virtual void updateGradientsAndValues
+        (
+            scalarField& objectiveSens,
+            PtrList<scalarField>& constraintSens,
+            scalar& objectiveValue,
+            scalarField& constraintValues
+        );
+
         //- Compute the merit function of the optimisation problem.
         //  Could be different than the objective function in case of
         //  constraint optimisation

src/optimisation/adjointOptimisation/adjoint/optimisation/optimisationType/incompressible/shapeOptimisation/shapeOptimisationIncompressible.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2020 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -49,6 +49,39 @@ addToRunTimeSelectionTable
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
+void shapeOptimisation::updateDesignVariables(scalarField& correction)
+{
+    // Communicate the movement to optMeshMovement
+    optMeshMovement_->setCorrection(correction);
+
+    if (updateGeometry_)
+    {
+        // Update the mesh
+        optMeshMovement_->moveMesh();
+
+        if (writeEachMesh_)
+        {
+            Info<< "  Writing new mesh points " << endl;
+            pointIOField points
+            (
+                IOobject
+                (
+                   "points",
+                    mesh_.pointsInstance(),
+                    mesh_.meshSubDir,
+                    mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                mesh_.points()
+            );
+            points.write();
+        }
+    }
+}
+
+
 void shapeOptimisation::computeEta
 (
     scalarField& correction
@@ -146,144 +179,6 @@ shapeOptimisation::shapeOptimisation
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void shapeOptimisation::update()
-{
-    Info<< nl << "Moving Mesh..." << nl << endl;
-
-    // Sum contributions
-    scalarField objectiveSens(0);
-    PtrList<scalarField> constraintSens(0);
-    scalar objectiveValue(Zero);
-    scalarField constraintValues(0);
-    forAll(adjointSolvManagers_, amI)
-    {
-        adjointSolverManager& adjSolvManager(adjointSolvManagers_[amI]);
-        const scalar opWeight(adjSolvManager.operatingPointWeight());
-
-        // Allocate objective sens size if necessary
-        tmp<scalarField> tadjointSolverManagerSens =
-            adjSolvManager.aggregateSensitivities();
-
-        if (objectiveSens.empty())
-        {
-            objectiveSens.setSize(tadjointSolverManagerSens().size(), Zero);
-        }
-        objectiveSens += opWeight*tadjointSolverManagerSens();
-        objectiveValue += opWeight*adjSolvManager.objectiveValue();
-
-        // Allocate constraint sens size if necessary
-        PtrList<scalarField> adjointSolverManagerConstSens
-            = adjSolvManager.constraintSensitivities();
-        tmp<scalarField> cValues = adjSolvManager.constraintValues();
-        if (constraintSens.empty())
-        {
-            constraintSens.setSize(adjointSolverManagerConstSens.size());
-            forAll(constraintSens, cI)
-            {
-                constraintSens.set
-                (
-                    cI,
-                    new scalarField
-                    (
-                        adjointSolverManagerConstSens[cI].size(),
-                        Zero
-                    )
-                );
-                constraintValues.setSize(cValues().size());
-                constraintValues = Zero;
-            }
-        }
-
-        forAll(constraintSens, cI)
-        {
-            constraintSens[cI] += opWeight*adjointSolverManagerConstSens[cI];
-        }
-        constraintValues += opWeight*cValues();
-    }
-
-    // Based on the sensitivities, return design variables correction
-    updateMethod_->setObjectiveDeriv(objectiveSens);
-    updateMethod_->setConstraintDeriv(constraintSens);
-    updateMethod_->setObjectiveValue(objectiveValue);
-    updateMethod_->setConstraintValues(constraintValues);
-    scalarField& correction = updateMethod_->returnCorrection();
-
-    // Computed  eta if needed
-    computeEta(correction);
-    updateMethod_->writeCorrection();
-
-    // Communicate the movement to optMeshMovement
-    optMeshMovement_->setCorrection(correction);
-    if (updateGeometry_)
-    {
-        optMeshMovement_->moveMesh();
-
-        if (writeEachMesh_)
-        {
-            Info<< "  Writing new mesh points " << endl;
-            pointIOField points
-            (
-                IOobject
-                (
-                    "points",
-                     mesh_.pointsInstance(),
-                     mesh_.meshSubDir,
-                     mesh_,
-                     IOobject::NO_READ,
-                     IOobject::NO_WRITE,
-                     false
-                ),
-                mesh_.points()
-            );
-            points.write();
-        }
-    }
-}
-
-
-void shapeOptimisation::update(scalarField& direction)
-{
-    // Computed eta if needed
-    computeEta(direction);
-
-    // Multiply with line search step, if necessary
-    scalarField correction = direction;
-    if (lineSearch_.valid())
-    {
-        correction *= lineSearch_->step();
-    }
-
-    // Communicate the movement to optMeshMovement
-    optMeshMovement_->setCorrection(correction);
-
-    if (updateGeometry_)
-    {
-        // Update the mesh
-        optMeshMovement_->moveMesh();
-
-        if (writeEachMesh_)
-        {
-            Info<< "  Writing new mesh points " << endl;
-            pointIOField points
-            (
-                IOobject
-                (
-                   "points",
-                    mesh_.pointsInstance(),
-                    mesh_.meshSubDir,
-                    mesh_,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE,
-                    false
-                ),
-                mesh_.points()
-            );
-            points.write();
-        }
-    }
-}
-
-
 void shapeOptimisation::storeDesignVariables()
 {
     optMeshMovement_->storeDesignVariables();

src/optimisation/adjointOptimisation/adjoint/optimisation/optimisationType/incompressible/shapeOptimisation/shapeOptimisationIncompressible.H

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2007-2020 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -76,6 +76,10 @@ protected:
 
     // Protected Member Functions
 
+        //- Update the design variables given their correction
+        virtual void updateDesignVariables(scalarField& correction);
+
+        //- Compute eta if not set in the first step
         virtual void computeEta(scalarField& correction);
 
 
@@ -113,12 +117,6 @@ public:
 
     // Member Functions
 
-       //- Master function. Calls all the others
-       void update();
-
-       //- Update design variables based on a given correction
-       virtual void update(scalarField& correction);
-
         //- Store design variables, as the starting point for line search
         virtual void storeDesignVariables();