diff --git a/src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/adjointSimple/adjointSimple.C b/src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/adjointSimple/adjointSimple.C
index 5b5eee62f6..8c850c5e07 100644
--- a/src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/adjointSimple/adjointSimple.C
+++ b/src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/adjointSimple/adjointSimple.C
@@ -181,9 +181,20 @@ bool Foam::adjointSimple::readDict(const dictionary& dict)
 
 void Foam::adjointSimple::solveIter()
 {
-    const Time& time = mesh_.time();
-    Info<< "Time = " << time.timeName() << "\n" << endl;
+    preIter();
+    mainIter();
+    postIter();
+}
 
+
+void Foam::adjointSimple::preIter()
+{
+    Info<< "Time = " << mesh_.time().timeName() << "\n" << endl;
+}
+
+
+void Foam::adjointSimple::mainIter()
+{
     // Grab primal references
     const surfaceScalarField& phi = primalVars_.phi();
     // Grab adjoint references
@@ -217,8 +228,8 @@ void Foam::adjointSimple::solveIter()
     // Add ATC term
     ATCModel_->addATC(UaEqn);
 
-    // Add source from optimisationType (e.g. topology)
-    addOptimisationTypeSource(UaEqn);
+    // Additional source terms (e.g. energy equation)
+    addMomentumSource(UaEqn);
 
     UaEqn.relax();
 
@@ -266,6 +277,8 @@ void Foam::adjointSimple::solveIter()
 
             paEqn.boundaryManipulate(pa.boundaryFieldRef());
 
+            addPressureSource(paEqn);
+
             paEqn.setReference(paRefCell, paRefValue);
 
             paEqn.solve();
@@ -292,19 +305,23 @@ void Foam::adjointSimple::solveIter()
 
     if (solverControl_().printMaxMags())
     {
-        dimensionedScalar maxUa = max(mag(Ua));
-        dimensionedScalar maxpa = max(mag(pa));
+        dimensionedScalar maxUa = gMax(mag(Ua)());
+        dimensionedScalar maxpa = gMax(mag(pa)());
         Info<< "Max mag of adjoint velocity = " << maxUa.value() << endl;
         Info<< "Max mag of adjoint pressure = " << maxpa.value() << endl;
     }
+}
 
+
+void Foam::adjointSimple::postIter()
+{
     solverControl_().write();
 
     // Average fields if necessary
     adjointVars_.computeMeanFields();
 
     // Print execution time
-    time.printExecutionTime(Info);
+    mesh_.time().printExecutionTime(Info);
 }
 
 
@@ -312,10 +329,7 @@ void Foam::adjointSimple::solve()
 {
     if (active_)
     {
-        // Reset mean fields before solving
-        adjointVars_.resetMeanFields();
-
-        // Iterate
+        preLoop();
         while (solverControl_().loop())
         {
             solveIter();
@@ -330,6 +344,13 @@ bool Foam::adjointSimple::loop()
 }
 
 
+void Foam::adjointSimple::preLoop()
+{
+    // Reset mean fields before solving
+    adjointVars_.resetMeanFields();
+}
+
+
 void Foam::adjointSimple::computeObjectiveSensitivities()
 {
     if (computeSensitivities_)
@@ -372,7 +393,7 @@ void Foam::adjointSimple::clearSensitivities()
 
 Foam::sensitivity& Foam::adjointSimple::getSensitivityBase()
 {
-    if (!adjointSensitivity_)
+    if (!adjointSensitivity_.valid())
     {
         FatalErrorInFunction
             << "Sensitivity object not allocated" << nl
@@ -386,6 +407,18 @@ Foam::sensitivity& Foam::adjointSimple::getSensitivityBase()
 }
 
 
+void Foam::adjointSimple::addMomentumSource(fvVectorMatrix& matrix)
+{
+    // Does nothing
+}
+
+
+void Foam::adjointSimple::addPressureSource(fvScalarMatrix& matrix)
+{
+    // Does nothing
+}
+
+
 void Foam::adjointSimple::updatePrimalBasedQuantities()
 {
     incompressibleAdjointSolver::updatePrimalBasedQuantities();
diff --git a/src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/adjointSimple/adjointSimple.H b/src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/adjointSimple/adjointSimple.H
index a5884247be..145e85b53e 100644
--- a/src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/adjointSimple/adjointSimple.H
+++ b/src/optimisation/adjointOptimisation/adjoint/solvers/adjointSolvers/incompressible/adjointSimple/adjointSimple.H
@@ -149,12 +149,24 @@ 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();
 
             //- Looper (advances iters, time step)
             virtual bool loop();
 
+            //- Functions to be called before loop
+            virtual void preLoop();
+
             //- Compute sensitivities of the underlaying objectives
             virtual void computeObjectiveSensitivities();
 
@@ -168,6 +180,15 @@ public:
             //- Return the base sensitivity object
             virtual sensitivity& getSensitivityBase();
 
+
+            // Source terms to be added to the adjoint equations
+
+                //- Source terms for the momentum equations
+                virtual void addMomentumSource(fvVectorMatrix& matrix);
+
+                //- Source terms for the continuity equation
+                virtual void addPressureSource(fvScalarMatrix& matrix);
+
             //- Update primal based quantities
             //- related to the objective functions.
             //  Also writes the objective function values to files.