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 f1f50530e36e924c371e6c23092a43cdd92c8d83..3cffe23bfb36301cc6c7df340b8f6b35026a19f1 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 016d930d96220de062ed9f8694b87538fc492608..9def052a3607bc928ca2dcd322a1b443fa165bec 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 5919cb373def5e18546057c6533b447ba7c39243..90a3343709402298f4ea9bc2b788485d4ca3f82c 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 bc384375a16687be5e68c81e44bc54f582ea15d2..335df8280fc79f9045f4903c3a87f31973efd144 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;
 };