diff --git a/SU2_CFD/include/drivers/CDriverBase.hpp b/SU2_CFD/include/drivers/CDriverBase.hpp
index d720246efb74..c7123daff4b9 100644
--- a/SU2_CFD/include/drivers/CDriverBase.hpp
+++ b/SU2_CFD/include/drivers/CDriverBase.hpp
@@ -473,6 +473,25 @@ class CDriverBase {
         "MarkerSolutionTimeN of " + solver->GetSolverName(), false);
   }
 
+  /*!
+   * \brief Get a read/write view of the solution at time N-1 on all mesh nodes of a solver.
+   */
+  inline CPyWrapperMatrixView SolutionTimeN1(unsigned short iSolver) {
+    auto* solver = GetSolverAndCheckMarker(iSolver);
+    return CPyWrapperMatrixView(
+        solver->GetNodes()->GetSolution_time_n1(), "SolutionTimeN1 of " + solver->GetSolverName(), false);
+  }
+
+  /*!
+   * \brief Get a read/write view of the solution at time N-1 on the mesh nodes of a marker.
+   */
+  inline CPyWrapperMarkerMatrixView MarkerSolutionTimeN1(unsigned short iSolver, unsigned short iMarker) {
+    auto* solver = GetSolverAndCheckMarker(iSolver, iMarker);
+    return CPyWrapperMarkerMatrixView(
+        solver->GetNodes()->GetSolution_time_n1(), main_geometry->vertex[iMarker], main_geometry->GetnVertex(iMarker),
+        "MarkerSolutionTimeN1 of " + solver->GetSolverName(), false);
+  }
+
   /*!
    * \brief Get the flow solver primitive variable names with their associated indices.
    * These correspond to the column indices in the matrix returned by Primitives.
diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp
index aa75d07596db..a5beafd75e0f 100644
--- a/SU2_CFD/include/variables/CVariable.hpp
+++ b/SU2_CFD/include/variables/CVariable.hpp
@@ -500,6 +500,7 @@ class CVariable {
    * \return Pointer to the solution (at time n-1) vector.
    */
   inline su2double *GetSolution_time_n1(unsigned long iPoint) { return Solution_time_n1[iPoint]; }
+  inline MatrixType& GetSolution_time_n1() { return Solution_time_n1; }
 
   /*!
    * \brief Set the value of the old residual.
diff --git a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp
index f9cfb269a8ed..7731d8899bb0 100644
--- a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp
+++ b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp
@@ -129,6 +129,10 @@ CDiscAdjSinglezoneDriver::~CDiscAdjSinglezoneDriver() {
 
 void CDiscAdjSinglezoneDriver::Preprocess(unsigned long TimeIter) {
 
+  /*--- Set the current time iteration in the config and also in the driver
+   * because the python interface doesn't offer an explicit way of doing it. ---*/
+
+  this->TimeIter = TimeIter;
   config_container[ZONE_0]->SetTimeIter(TimeIter);
 
   /*--- Preprocess the adjoint iteration ---*/
diff --git a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp
index 3a3518d33ab1..964888a83fb0 100644
--- a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp
+++ b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp
@@ -110,8 +110,10 @@ void CSinglezoneDriver::StartSolver() {
 
 void CSinglezoneDriver::Preprocess(unsigned long TimeIter) {
 
-  /*--- Set the current time iteration in the config ---*/
+  /*--- Set the current time iteration in the config and also in the driver
+   * because the python interface doesn't offer an explicit way of doing it. ---*/
 
+  this->TimeIter = TimeIter;
   config_container[ZONE_0]->SetTimeIter(TimeIter);
 
   /*--- Store the current physical time in the config container, as
diff --git a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp
index 93e29ff3f30d..5b4519078651 100644
--- a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp
+++ b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp
@@ -41,74 +41,68 @@ void CDiscAdjFluidIteration::Preprocess(COutput* output, CIntegration**** integr
   const bool grid_IsMoving = config[iZone]->GetGrid_Movement();
   const bool species = config[iZone]->GetKind_Species_Model() != SPECIES_MODEL::NONE;
   const bool heat = config[iZone]->GetWeakly_Coupled_Heat();
+  const bool radiation = config[iZone]->AddRadiation();
 
   auto solvers0 = solver[iZone][iInst][MESH_0];
   auto geometries = geometry[iZone][iInst];
 
+  unsigned long nSolvers = 1;
+  std::array<int, 5> solversToProcess{{FLOW_SOL}};
+  if (turbulent) solversToProcess[nSolvers++] = TURB_SOL;
+  if (species) solversToProcess[nSolvers++] = SPECIES_SOL;
+  if (heat) solversToProcess[nSolvers++] = HEAT_SOL;
+  if (radiation) solversToProcess[nSolvers++] = RAD_SOL;
+
   /*--- For the unsteady adjoint, load direct solutions from restart files. ---*/
 
   if (config[iZone]->GetTime_Marching() != TIME_MARCHING::STEADY) {
-    const int Direct_Iter = static_cast<int>(config[iZone]->GetUnst_AdjointIter()) - static_cast<int>(TimeIter) - 2 + dual_time;
+    const int Direct_Iter = static_cast<int>(config[iZone]->GetUnst_AdjointIter()) -
+                            static_cast<int>(TimeIter) - 2 + dual_time;
 
-    /*--- For dual-time stepping we want to load the already converged solution at timestep n ---*/
+    /*--- For dual-time stepping we want to load the already converged solution at previous timesteps.
+     * In general we only load one file and shift the previously loaded solutions, on the first we
+     * load one or two more (depending on dual time order). ---*/
 
-    if (TimeIter == 0) {
-      if (dual_time_2nd) {
-        /*--- Load solution at timestep n-2 ---*/
-        LoadUnsteady_Solution(geometry, solver, config, iZone, iInst, Direct_Iter - 2);
+    if (dual_time_2nd) {
+      LoadUnsteady_Solution(geometry, solver, config, iZone, iInst, Direct_Iter - 2);
+    } else if (dual_time_1st) {
+      LoadUnsteady_Solution(geometry, solver, config, iZone, iInst, Direct_Iter - 1);
+    }
 
-        /*--- Push solution back to correct array ---*/
+    if (TimeIter == 0) {
+      /*--- Push solution back one level. ---*/
 
+      if (dual_time) {
         for (auto iMesh = 0u; iMesh <= config[iZone]->GetnMGLevels(); iMesh++) {
-          auto solvers = solver[iZone][iInst][iMesh];
-
-          solvers[FLOW_SOL]->GetNodes()->Set_Solution_time_n();
-          solvers[FLOW_SOL]->GetNodes()->Set_Solution_time_n1();
-          if (turbulent) {
-            solvers[TURB_SOL]->GetNodes()->Set_Solution_time_n();
-            solvers[TURB_SOL]->GetNodes()->Set_Solution_time_n1();
-          }
-          if (species) {
-            solvers[SPECIES_SOL]->GetNodes()->Set_Solution_time_n();
-            solvers[SPECIES_SOL]->GetNodes()->Set_Solution_time_n1();
-          }
-          if (heat) {
-            solvers[HEAT_SOL]->GetNodes()->Set_Solution_time_n();
-            solvers[HEAT_SOL]->GetNodes()->Set_Solution_time_n1();
+          for (auto iSol = 0ul; iSol < nSolvers; ++iSol) {
+            solver[iZone][iInst][iMesh][solversToProcess[iSol]]->GetNodes()->Set_Solution_time_n();
           }
           if (grid_IsMoving) {
             geometries[iMesh]->nodes->SetCoord_n();
-            geometries[iMesh]->nodes->SetCoord_n1();
           }
           if (config[iZone]->GetDynamic_Grid()) {
             geometries[iMesh]->nodes->SetVolume_n();
-            geometries[iMesh]->nodes->SetVolume_nM1();
           }
         }
       }
-      if (dual_time) {
-        /*--- Load solution at timestep n-1 ---*/
-        LoadUnsteady_Solution(geometry, solver, config, iZone, iInst, Direct_Iter - 1);
 
-        /*--- Push solution back to correct array ---*/
+      /*--- If required load another time step. Push the previous time step to n-1 and the
+       loaded time step to n. ---*/
 
-        for (auto iMesh = 0u; iMesh <= config[iZone]->GetnMGLevels(); iMesh++) {
-          auto solvers = solver[iZone][iInst][iMesh];
+      if (dual_time_2nd) {
+        LoadUnsteady_Solution(geometry, solver, config, iZone, iInst, Direct_Iter - 1);
 
-          solvers[FLOW_SOL]->GetNodes()->Set_Solution_time_n();
-          if (turbulent) {
-            solvers[TURB_SOL]->GetNodes()->Set_Solution_time_n();
-          }
-          if (species) {
-            solvers[SPECIES_SOL]->GetNodes()->Set_Solution_time_n();
-          }
-          if (heat) {
-            solvers[HEAT_SOL]->GetNodes()->Set_Solution_time_n();
+        for (auto iMesh = 0u; iMesh <= config[iZone]->GetnMGLevels(); iMesh++) {
+          for (auto iSol = 0ul; iSol < nSolvers; ++iSol) {
+            solver[iZone][iInst][iMesh][solversToProcess[iSol]]->GetNodes()->Set_Solution_time_n1();
+            solver[iZone][iInst][iMesh][solversToProcess[iSol]]->GetNodes()->Set_Solution_time_n();
           }
           if (grid_IsMoving) {
+            geometries[iMesh]->nodes->SetCoord_n1();
             geometries[iMesh]->nodes->SetCoord_n();
           }
           if (config[iZone]->GetDynamic_Grid()) {
+            geometries[iMesh]->nodes->SetVolume_nM1();
             geometries[iMesh]->nodes->SetVolume_n();
           }
         }
@@ -123,15 +117,13 @@ void CDiscAdjFluidIteration::Preprocess(COutput* output, CIntegration**** integr
       }
 
     }
-    else if ((TimeIter > 0) && dual_time) {
-      /*---
-      Here the primal solutions (only working variables) are loaded and put in the correct order
-      into containers. For ALE the mesh coordinates have to be put into the
-      correct containers as well, i.e. follow the same logic for the solution.
-      Afterwards the GridVelocity is computed based on the Coordinates.
-      ---*/
-
-      /*--- Temporarily store the loaded volumes into old containers ---*/
+    if ((TimeIter > 0) && dual_time) {
+      /*--- Here the primal solutions (only working variables) are loaded and put in the correct order
+       into containers. For ALE the mesh coordinates have to be put into the correct containers as well,
+       i.e. follow the same logic for the solution. Afterwards the GridVelocity is computed based on
+       the Coordinates. ---*/
+
+      /*--- Temporarily store the loaded volumes into old containers. ---*/
       if (config[iZone]->GetDynamic_Grid()) {
         for (auto iMesh=0; iMesh<=config[iZone]->GetnMGLevels();iMesh++) {
           geometries[iMesh]->nodes->SetVolume_Old();
@@ -140,155 +132,80 @@ void CDiscAdjFluidIteration::Preprocess(COutput* output, CIntegration**** integr
         }
       }
 
-      /*-- Load mesh solver ---*/
+      /*--- Load mesh solver. ---*/
       if (config[iZone]->GetDeform_Mesh()) {
         solvers0[MESH_SOL]->LoadRestart(geometries, solver[iZone][iInst], config[iZone], Direct_Iter, true);
       }
 
-      /*--- Load solution timestep n-1 | n-2 for DualTimestepping 1st | 2nd order ---*/
-      if (dual_time_1st) {
-        LoadUnsteady_Solution(geometry, solver, config, iZone, iInst, Direct_Iter - 1);
-      } else {
-        LoadUnsteady_Solution(geometry, solver, config, iZone, iInst, Direct_Iter - 2);
-
-        /*--- Set volumes into correct containers ---*/
-        if (config[iZone]->GetDynamic_Grid()) {
-          for (auto iMesh=0; iMesh<=config[iZone]->GetnMGLevels();iMesh++) {
-            /*--- If negative iteration number, set default ---*/
-            if (Direct_Iter - 2 < 0) {
-              for(auto iPoint=0ul; iPoint<geometries[iMesh]->GetnPoint();iPoint++) {
-                geometries[iMesh]->nodes->SetVolume(iPoint,0.0);
-              }
+      /*--- Set volumes into correct containers ---*/
+      if (config[iZone]->GetDynamic_Grid()) {
+        for (auto iMesh = 0; iMesh <= config[iZone]->GetnMGLevels(); iMesh++) {
+          /*--- If negative iteration number, set default. ---*/
+          if (Direct_Iter - 1 - dual_time_2nd < 0) {
+            for(auto iPoint = 0ul; iPoint < geometries[iMesh]->GetnPoint(); iPoint++) {
+              geometries[iMesh]->nodes->SetVolume(iPoint, 0.0);
             }
+          }
 
-            /*--- Set currently loaded volume to Volume_nM1 ---*/
-            geometries[iMesh]->nodes->SetVolume_n();
-            geometries[iMesh]->nodes->SetVolume_nM1();
+          /*--- Set currently loaded volume to Volume_nM1 ---*/
+          geometries[iMesh]->nodes->SetVolume_n();
+          geometries[iMesh]->nodes->SetVolume_nM1();
 
-            /*--- Set Volume_n and Volume from old containers ---*/
-            geometries[iMesh]->nodes->SetVolume_n_from_OldnM1();
-            geometries[iMesh]->nodes->SetVolume_from_Oldn();
-          }
+          /*--- Set Volume_n and Volume from old containers ---*/
+          geometries[iMesh]->nodes->SetVolume_n_from_OldnM1();
+          geometries[iMesh]->nodes->SetVolume_from_Oldn();
         }
       }
 
-      /*--- Temporarily store the loaded solution in the Solution_Old array ---*/
+      /*--- Temporarily store the loaded solution and coordinates in the "Old" arrays. ---*/
 
       for (auto iMesh = 0u; iMesh <= config[iZone]->GetnMGLevels(); iMesh++) {
-        auto solvers = solver[iZone][iInst][iMesh];
-
-        solvers[FLOW_SOL]->Set_OldSolution();
-        if (turbulent) {
-          solvers[TURB_SOL]->Set_OldSolution();
-        }
-        if (species) {
-          solvers[SPECIES_SOL]->Set_OldSolution();
-        }
-        if (heat) {
-          solvers[HEAT_SOL]->Set_OldSolution();
+        for (auto iSol = 0ul; iSol < nSolvers; ++iSol) {
+          solver[iZone][iInst][iMesh][solversToProcess[iSol]]->Set_OldSolution();
         }
         if (grid_IsMoving) {
           geometries[iMesh]->nodes->SetCoord_Old();
         }
       }
 
-      /*--- Set Solution at timestep n to solution at n-1 ---*/
+      /*--- Move timestep n to current solution. ---*/
 
       for (auto iMesh = 0u; iMesh <= config[iZone]->GetnMGLevels(); iMesh++) {
-        auto solvers = solver[iZone][iInst][iMesh];
-
-        for (auto iPoint = 0ul; iPoint < geometries[iMesh]->GetnPoint(); iPoint++) {
-          solvers[FLOW_SOL]->GetNodes()->SetSolution(iPoint, solvers[FLOW_SOL]->GetNodes()->GetSolution_time_n(iPoint));
-
-          if (grid_IsMoving) {
-            geometries[iMesh]->nodes->SetCoord(iPoint, geometries[iMesh]->nodes->GetCoord_n(iPoint));
-          }
-          if (turbulent) {
-            solvers[TURB_SOL]->GetNodes()->SetSolution(iPoint, solvers[TURB_SOL]->GetNodes()->GetSolution_time_n(iPoint));
-          }
-          if (species) {
-            solvers[SPECIES_SOL]->GetNodes()->SetSolution(iPoint, solvers[SPECIES_SOL]->GetNodes()->GetSolution_time_n(iPoint));
-          }
-          if (heat) {
-            solvers[HEAT_SOL]->GetNodes()->SetSolution(iPoint, solvers[HEAT_SOL]->GetNodes()->GetSolution_time_n(iPoint));
+        for (auto iSol = 0ul; iSol < nSolvers; ++iSol) {
+          auto* s = solver[iZone][iInst][iMesh][solversToProcess[iSol]];
+          for (auto iPoint = 0ul; iPoint < geometries[iMesh]->GetnPoint(); iPoint++) {
+            s->GetNodes()->SetSolution(iPoint, s->GetNodes()->GetSolution_time_n(iPoint));
           }
         }
-      }
-      if (dual_time_1st) {
-        /*--- Set Solution at timestep n-1 to the previously loaded solution ---*/
-        for (auto iMesh = 0u; iMesh <= config[iZone]->GetnMGLevels(); iMesh++) {
-          auto solvers = solver[iZone][iInst][iMesh];
-
+        if (grid_IsMoving) {
           for (auto iPoint = 0ul; iPoint < geometries[iMesh]->GetnPoint(); iPoint++) {
-            solvers[FLOW_SOL]->GetNodes()->Set_Solution_time_n(
-                iPoint, solvers[FLOW_SOL]->GetNodes()->GetSolution_Old(iPoint));
-
-            if (grid_IsMoving) {
-              geometries[iMesh]->nodes->SetCoord_n(iPoint, geometries[iMesh]->nodes->GetCoord_Old(iPoint));
-            }
-            if (turbulent) {
-              solvers[TURB_SOL]->GetNodes()->Set_Solution_time_n(
-                  iPoint, solvers[TURB_SOL]->GetNodes()->GetSolution_Old(iPoint));
-            }
-            if (species) {
-              solvers[SPECIES_SOL]->GetNodes()->Set_Solution_time_n(
-                  iPoint, solvers[SPECIES_SOL]->GetNodes()->GetSolution_Old(iPoint));
-            }
-            if (heat) {
-              solvers[HEAT_SOL]->GetNodes()->Set_Solution_time_n(
-                  iPoint, solvers[HEAT_SOL]->GetNodes()->GetSolution_Old(iPoint));
-            }
+            geometries[iMesh]->nodes->SetCoord(iPoint, geometries[iMesh]->nodes->GetCoord_n(iPoint));
           }
         }
       }
-      if (dual_time_2nd) {
-        /*--- Set Solution at timestep n-1 to solution at n-2 ---*/
-        for (auto iMesh = 0u; iMesh <= config[iZone]->GetnMGLevels(); iMesh++) {
-          auto solvers = solver[iZone][iInst][iMesh];
 
-          for (auto iPoint = 0ul; iPoint < geometries[iMesh]->GetnPoint(); iPoint++) {
-            solvers[FLOW_SOL]->GetNodes()->Set_Solution_time_n(
-                iPoint, solvers[FLOW_SOL]->GetNodes()->GetSolution_time_n1(iPoint));
+      /*--- Finally, place the loaded solution in the correct place (n or n-1 depending on order). ---*/
 
-            if (grid_IsMoving) {
-              geometries[iMesh]->nodes->SetCoord_n(iPoint, geometries[iMesh]->nodes->GetCoord_n1(iPoint));
-            }
-            if (turbulent) {
-              solvers[TURB_SOL]->GetNodes()->Set_Solution_time_n(
-                  iPoint, solvers[TURB_SOL]->GetNodes()->GetSolution_time_n1(iPoint));
-            }
-            if (species) {
-              solvers[SPECIES_SOL]->GetNodes()->Set_Solution_time_n(
-                  iPoint, solvers[SPECIES_SOL]->GetNodes()->GetSolution_time_n1(iPoint));
-            }
-            if (heat) {
-              solvers[HEAT_SOL]->GetNodes()->Set_Solution_time_n(
-                  iPoint, solvers[HEAT_SOL]->GetNodes()->GetSolution_time_n1(iPoint));
+      for (auto iMesh = 0u; iMesh <= config[iZone]->GetnMGLevels(); iMesh++) {
+        for (auto iSol = 0ul; iSol < nSolvers; ++iSol) {
+          auto* s = solver[iZone][iInst][iMesh][solversToProcess[iSol]];
+          for (auto iPoint = 0ul; iPoint < geometries[iMesh]->GetnPoint(); iPoint++) {
+            if (dual_time_2nd) {
+              /*--- If required also move timestep n-1 to timestep n. ---*/
+              s->GetNodes()->Set_Solution_time_n(iPoint, s->GetNodes()->GetSolution_time_n1(iPoint));
+              s->GetNodes()->Set_Solution_time_n1(iPoint, s->GetNodes()->GetSolution_Old(iPoint));
+            } else {
+              s->GetNodes()->Set_Solution_time_n(iPoint, s->GetNodes()->GetSolution_Old(iPoint));
             }
           }
         }
-        /*--- Set Solution at timestep n-2 to the previously loaded solution ---*/
-        for (auto iMesh = 0u; iMesh <= config[iZone]->GetnMGLevels(); iMesh++) {
-          auto solvers = solver[iZone][iInst][iMesh];
-
+        if (grid_IsMoving) {
           for (auto iPoint = 0ul; iPoint < geometries[iMesh]->GetnPoint(); iPoint++) {
-            solvers[FLOW_SOL]->GetNodes()->Set_Solution_time_n1(
-                iPoint, solvers[FLOW_SOL]->GetNodes()->GetSolution_Old(iPoint));
-
-            if (grid_IsMoving) {
+            if (dual_time_2nd) {
+              geometries[iMesh]->nodes->SetCoord_n(iPoint, geometries[iMesh]->nodes->GetCoord_n1(iPoint));
               geometries[iMesh]->nodes->SetCoord_n1(iPoint, geometries[iMesh]->nodes->GetCoord_Old(iPoint));
-            }
-            if (turbulent) {
-              solvers[TURB_SOL]->GetNodes()->Set_Solution_time_n1(
-                  iPoint, solvers[TURB_SOL]->GetNodes()->GetSolution_Old(iPoint));
-            }
-            if (species) {
-              solvers[SPECIES_SOL]->GetNodes()->Set_Solution_time_n1(
-                  iPoint, solvers[SPECIES_SOL]->GetNodes()->GetSolution_Old(iPoint));
-            }
-            if (heat) {
-              solvers[HEAT_SOL]->GetNodes()->Set_Solution_time_n1(
-                  iPoint, solvers[HEAT_SOL]->GetNodes()->GetSolution_Old(iPoint));
+            } else {
+              geometries[iMesh]->nodes->SetCoord_n(iPoint, geometries[iMesh]->nodes->GetCoord_Old(iPoint));
             }
           }
         }
@@ -297,9 +214,10 @@ void CDiscAdjFluidIteration::Preprocess(COutput* output, CIntegration**** integr
     }  // else if TimeIter > 0
 
     /*--- Compute & set Grid Velocity via finite differences of the Coordinates. ---*/
-    if (grid_IsMoving)
+    if (grid_IsMoving) {
       for (auto iMesh = 0u; iMesh <= config[iZone]->GetnMGLevels(); iMesh++)
         geometries[iMesh]->SetGridVelocity(config[iZone]);
+    }
 
   }  // if unsteady
 
@@ -308,36 +226,28 @@ void CDiscAdjFluidIteration::Preprocess(COutput* output, CIntegration**** integr
   /*--- Store flow solution also in the adjoint solver in order to be able to reset it later ---*/
 
   if (TimeIter == 0 || dual_time) {
-    for (auto iMesh = 0u; iMesh <= config[iZone]->GetnMGLevels(); iMesh++) {
-      auto solvers = solver[iZone][iInst][iMesh];
+    auto SetSolutionDirect = [&](CSolver** solvers, int adj, int primal, unsigned long nPoint) {
       SU2_OMP_FOR_STAT(1024)
-      for (auto iPoint = 0ul; iPoint < geometries[iMesh]->GetnPoint(); iPoint++)
-        solvers[ADJFLOW_SOL]->GetNodes()->SetSolution_Direct(iPoint, solvers[FLOW_SOL]->GetNodes()->GetSolution(iPoint));
+      for (auto iPoint = 0ul; iPoint < nPoint; iPoint++)
+        solvers[adj]->GetNodes()->SetSolution_Direct(iPoint, solvers[primal]->GetNodes()->GetSolution(iPoint));
       END_SU2_OMP_FOR
+    };
+
+    for (auto iMesh = 0u; iMesh <= config[iZone]->GetnMGLevels(); iMesh++) {
+      auto solvers = solver[iZone][iInst][iMesh];
+      SetSolutionDirect(solvers, ADJFLOW_SOL, FLOW_SOL, geometries[iMesh]->GetnPoint());
     }
     if (turbulent && !config[iZone]->GetFrozen_Visc_Disc()) {
-      SU2_OMP_FOR_STAT(1024)
-      for (auto iPoint = 0ul; iPoint < geometries[MESH_0]->GetnPoint(); iPoint++)
-        solvers0[ADJTURB_SOL]->GetNodes()->SetSolution_Direct(iPoint, solvers0[TURB_SOL]->GetNodes()->GetSolution(iPoint));
-      END_SU2_OMP_FOR
+      SetSolutionDirect(solvers0, ADJTURB_SOL, TURB_SOL, geometries[MESH_0]->GetnPoint());
     }
     if (species) {
-      SU2_OMP_FOR_STAT(1024)
-      for (auto iPoint = 0ul; iPoint < geometries[MESH_0]->GetnPoint(); iPoint++)
-        solvers0[ADJSPECIES_SOL]->GetNodes()->SetSolution_Direct(iPoint, solvers0[SPECIES_SOL]->GetNodes()->GetSolution(iPoint));
-      END_SU2_OMP_FOR
+      SetSolutionDirect(solvers0, ADJSPECIES_SOL, SPECIES_SOL, geometries[MESH_0]->GetnPoint());
     }
     if (heat) {
-      SU2_OMP_FOR_STAT(1024)
-      for (auto iPoint = 0ul; iPoint < geometries[MESH_0]->GetnPoint(); iPoint++)
-        solvers0[ADJHEAT_SOL]->GetNodes()->SetSolution_Direct(iPoint, solvers0[HEAT_SOL]->GetNodes()->GetSolution(iPoint));
-      END_SU2_OMP_FOR
+      SetSolutionDirect(solvers0, ADJHEAT_SOL, HEAT_SOL, geometries[MESH_0]->GetnPoint());
     }
-    if (config[iZone]->AddRadiation()) {
-      SU2_OMP_FOR_STAT(1024)
-      for (auto iPoint = 0ul; iPoint < geometries[MESH_0]->GetnPoint(); iPoint++)
-        solvers0[ADJRAD_SOL]->GetNodes()->SetSolution_Direct(iPoint, solvers0[RAD_SOL]->GetNodes()->GetSolution(iPoint));
-      END_SU2_OMP_FOR
+    if (radiation) {
+      SetSolutionDirect(solvers0, ADJRAD_SOL, RAD_SOL, geometries[MESH_0]->GetnPoint());
     }
   }
 
@@ -356,7 +266,7 @@ void CDiscAdjFluidIteration::Preprocess(COutput* output, CIntegration**** integr
     solvers0[ADJHEAT_SOL]->Preprocessing(geometries[MESH_0], solvers0, config[iZone],
                                          MESH_0, 0, RUNTIME_ADJHEAT_SYS, false);
   }
-  if (config[iZone]->AddRadiation()) {
+  if (radiation) {
     solvers0[ADJRAD_SOL]->Preprocessing(geometries[MESH_0], solvers0, config[iZone],
                                         MESH_0, 0, RUNTIME_ADJRAD_SYS, false);
   }
@@ -369,22 +279,23 @@ void CDiscAdjFluidIteration::LoadUnsteady_Solution(CGeometry**** geometry, CSolv
                                                    unsigned short iZone, unsigned short iInst, int DirectIter) {
 
   auto solvers = solver[iZone][iInst];
+  auto geometries = geometry[iZone][iInst];
   const bool species = config[iZone]->GetKind_Species_Model() != SPECIES_MODEL::NONE;
 
   if (DirectIter >= 0) {
     if (rank == MASTER_NODE)
       cout << " Loading flow solution from direct iteration " << DirectIter << " for zone " << iZone << "." << endl;
 
-    solvers[MESH_0][FLOW_SOL]->LoadRestart(geometry[iZone][iInst], solvers, config[iZone], DirectIter, true);
+    solvers[MESH_0][FLOW_SOL]->LoadRestart(geometries, solvers, config[iZone], DirectIter, true);
 
     if (turbulent) {
-      solvers[MESH_0][TURB_SOL]->LoadRestart(geometry[iZone][iInst], solvers, config[iZone], DirectIter, false);
+      solvers[MESH_0][TURB_SOL]->LoadRestart(geometries, solvers, config[iZone], DirectIter, false);
     }
     if (species) {
-      solvers[MESH_0][SPECIES_SOL]->LoadRestart(geometry[iZone][iInst], solvers, config[iZone], DirectIter, false);
+      solvers[MESH_0][SPECIES_SOL]->LoadRestart(geometries, solvers, config[iZone], DirectIter, false);
     }
     if (config[iZone]->GetWeakly_Coupled_Heat()) {
-      solvers[MESH_0][HEAT_SOL]->LoadRestart(geometry[iZone][iInst], solvers, config[iZone], DirectIter, false);
+      solvers[MESH_0][HEAT_SOL]->LoadRestart(geometries, solvers, config[iZone], DirectIter, false);
     }
   } else {
     /*--- If there is no solution file we set the freestream condition ---*/
@@ -393,19 +304,19 @@ void CDiscAdjFluidIteration::LoadUnsteady_Solution(CGeometry**** geometry, CSolv
 
     for (auto iMesh = 0u; iMesh <= config[iZone]->GetnMGLevels(); iMesh++) {
       solvers[iMesh][FLOW_SOL]->SetFreeStream_Solution(config[iZone]);
-      solvers[iMesh][FLOW_SOL]->Preprocessing(geometry[iZone][iInst][iMesh], solvers[iMesh], config[iZone], iMesh,
+      solvers[iMesh][FLOW_SOL]->Preprocessing(geometries[iMesh], solvers[iMesh], config[iZone], iMesh,
                                               DirectIter, RUNTIME_FLOW_SYS, false);
       if (turbulent) {
         solvers[iMesh][TURB_SOL]->SetFreeStream_Solution(config[iZone]);
-        solvers[iMesh][TURB_SOL]->Postprocessing(geometry[iZone][iInst][iMesh], solvers[iMesh], config[iZone], iMesh);
+        solvers[iMesh][TURB_SOL]->Postprocessing(geometries[iMesh], solvers[iMesh], config[iZone], iMesh);
       }
       if (species) {
         solvers[iMesh][SPECIES_SOL]->SetFreeStream_Solution(config[iZone]);
-        solvers[iMesh][SPECIES_SOL]->Postprocessing(geometry[iZone][iInst][iMesh], solvers[iMesh], config[iZone], iMesh);
+        solvers[iMesh][SPECIES_SOL]->Postprocessing(geometries[iMesh], solvers[iMesh], config[iZone], iMesh);
       }
       if (config[iZone]->GetWeakly_Coupled_Heat()) {
         solvers[iMesh][HEAT_SOL]->SetFreeStream_Solution(config[iZone]);
-        solvers[iMesh][HEAT_SOL]->Postprocessing(geometry[iZone][iInst][iMesh], solvers[iMesh], config[iZone], iMesh);
+        solvers[iMesh][HEAT_SOL]->Postprocessing(geometries[iMesh], solvers[iMesh], config[iZone], iMesh);
       }
     }
   }
@@ -413,29 +324,31 @@ void CDiscAdjFluidIteration::LoadUnsteady_Solution(CGeometry**** geometry, CSolv
 
 void CDiscAdjFluidIteration::IterateDiscAdj(CGeometry**** geometry, CSolver***** solver, CConfig** config,
                                             unsigned short iZone, unsigned short iInst, bool CrossTerm) {
+  auto solvers0 = solver[iZone][iInst][MESH_0];
+  auto geometry0 = geometry[iZone][iInst][MESH_0];
 
-  SU2_OMP_PARALLEL_(if(solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->GetHasHybridParallel())) {
+  SU2_OMP_PARALLEL_(if(solvers0[ADJFLOW_SOL]->GetHasHybridParallel())) {
 
   /*--- Extract the adjoints of the conservative input variables and store them for the next iteration ---*/
 
   if (config[iZone]->GetFluidProblem()) {
-    solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone], CrossTerm);
+    solvers0[ADJFLOW_SOL]->ExtractAdjoint_Solution(geometry0, config[iZone], CrossTerm);
 
-    solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->ExtractAdjoint_Variables(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJFLOW_SOL]->ExtractAdjoint_Variables(geometry0, config[iZone]);
   }
   if (turbulent && !config[iZone]->GetFrozen_Visc_Disc()) {
-    solver[iZone][iInst][MESH_0][ADJTURB_SOL]->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone], CrossTerm);
+    solvers0[ADJTURB_SOL]->ExtractAdjoint_Solution(geometry0, config[iZone], CrossTerm);
   }
   if (config[iZone]->GetKind_Species_Model() != SPECIES_MODEL::NONE) {
-    solver[iZone][iInst][MESH_0][ADJSPECIES_SOL]->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone], CrossTerm);
+    solvers0[ADJSPECIES_SOL]->ExtractAdjoint_Solution(geometry0, config[iZone], CrossTerm);
   }
   if (config[iZone]->GetWeakly_Coupled_Heat()) {
-    solver[iZone][iInst][MESH_0][ADJHEAT_SOL]->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone], CrossTerm);
+    solvers0[ADJHEAT_SOL]->ExtractAdjoint_Solution(geometry0, config[iZone], CrossTerm);
   }
   if (config[iZone]->AddRadiation()) {
-    solver[iZone][iInst][MESH_0][ADJRAD_SOL]->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone], CrossTerm);
+    solvers0[ADJRAD_SOL]->ExtractAdjoint_Solution(geometry0, config[iZone], CrossTerm);
 
-    solver[iZone][iInst][MESH_0][ADJRAD_SOL]->ExtractAdjoint_Variables(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJRAD_SOL]->ExtractAdjoint_Variables(geometry0, config[iZone]);
   }
 
   }
@@ -444,33 +357,35 @@ void CDiscAdjFluidIteration::IterateDiscAdj(CGeometry**** geometry, CSolver*****
 
 void CDiscAdjFluidIteration::InitializeAdjoint(CSolver***** solver, CGeometry**** geometry, CConfig** config,
                                                unsigned short iZone, unsigned short iInst) {
+  auto solvers0 = solver[iZone][iInst][MESH_0];
+  auto geometry0 = geometry[iZone][iInst][MESH_0];
 
-  SU2_OMP_PARALLEL_(if(solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->GetHasHybridParallel())) {
+  SU2_OMP_PARALLEL_(if(solvers0[ADJFLOW_SOL]->GetHasHybridParallel())) {
 
   /*--- Initialize the adjoints the conservative variables ---*/
 
   if (config[iZone]->GetFluidProblem()) {
-    solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->SetAdjoint_Output(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJFLOW_SOL]->SetAdjoint_Output(geometry0, config[iZone]);
   }
 
   if (turbulent && !config[iZone]->GetFrozen_Visc_Disc()) {
-    solver[iZone][iInst][MESH_0][ADJTURB_SOL]->SetAdjoint_Output(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJTURB_SOL]->SetAdjoint_Output(geometry0, config[iZone]);
   }
 
   if (config[iZone]->GetKind_Species_Model() != SPECIES_MODEL::NONE) {
-    solver[iZone][iInst][MESH_0][ADJSPECIES_SOL]->SetAdjoint_Output(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJSPECIES_SOL]->SetAdjoint_Output(geometry0, config[iZone]);
   }
 
   if (config[iZone]->GetWeakly_Coupled_Heat()) {
-    solver[iZone][iInst][MESH_0][ADJHEAT_SOL]->SetAdjoint_Output(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJHEAT_SOL]->SetAdjoint_Output(geometry0, config[iZone]);
   }
 
   if (config[iZone]->AddRadiation()) {
-    solver[iZone][iInst][MESH_0][ADJRAD_SOL]->SetAdjoint_Output(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJRAD_SOL]->SetAdjoint_Output(geometry0, config[iZone]);
   }
 
   if (config[iZone]->GetFluidProblem() && !config[iZone]->GetMultizone_Problem()) {
-    solver[iZone][iInst][MESH_0][FLOW_SOL]->SetVertexTractionsAdjoint(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[FLOW_SOL]->SetVertexTractionsAdjoint(geometry0, config[iZone]);
   }
 
   }
@@ -479,46 +394,48 @@ void CDiscAdjFluidIteration::InitializeAdjoint(CSolver***** solver, CGeometry***
 
 void CDiscAdjFluidIteration::RegisterInput(CSolver***** solver, CGeometry**** geometry, CConfig** config,
                                            unsigned short iZone, unsigned short iInst, RECORDING kind_recording) {
+  auto solvers0 = solver[iZone][iInst][MESH_0];
+  auto geometry0 = geometry[iZone][iInst][MESH_0];
 
-  SU2_OMP_PARALLEL_(if(solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->GetHasHybridParallel())) {
+  SU2_OMP_PARALLEL_(if(solvers0[ADJFLOW_SOL]->GetHasHybridParallel())) {
 
   if (kind_recording == RECORDING::SOLUTION_VARIABLES || kind_recording == RECORDING::SOLUTION_AND_MESH) {
     /*--- Register flow and turbulent variables as input ---*/
 
     if (config[iZone]->GetFluidProblem()) {
-      solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]);
+      solvers0[ADJFLOW_SOL]->RegisterSolution(geometry0, config[iZone]);
 
-      solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->RegisterVariables(geometry[iZone][iInst][MESH_0], config[iZone]);
+      solvers0[ADJFLOW_SOL]->RegisterVariables(geometry0, config[iZone]);
     }
 
     if (turbulent && !config[iZone]->GetFrozen_Visc_Disc()) {
-      solver[iZone][iInst][MESH_0][ADJTURB_SOL]->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]);
+      solvers0[ADJTURB_SOL]->RegisterSolution(geometry0, config[iZone]);
     }
     if (config[iZone]->GetKind_Species_Model() != SPECIES_MODEL::NONE) {
-      solver[iZone][iInst][MESH_0][ADJSPECIES_SOL]->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]);
+      solvers0[ADJSPECIES_SOL]->RegisterSolution(geometry0, config[iZone]);
     }
     if (config[iZone]->GetWeakly_Coupled_Heat()) {
-      solver[iZone][iInst][MESH_0][ADJHEAT_SOL]->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]);
+      solvers0[ADJHEAT_SOL]->RegisterSolution(geometry0, config[iZone]);
     }
     if (config[iZone]->AddRadiation()) {
-      solver[iZone][iInst][MESH_0][ADJRAD_SOL]->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]);
+      solvers0[ADJRAD_SOL]->RegisterSolution(geometry0, config[iZone]);
 
-      solver[iZone][iInst][MESH_0][ADJRAD_SOL]->RegisterVariables(geometry[iZone][iInst][MESH_0], config[iZone]);
+      solvers0[ADJRAD_SOL]->RegisterVariables(geometry0, config[iZone]);
     }
   }
 
   if (kind_recording == RECORDING::MESH_COORDS || kind_recording == RECORDING::SOLUTION_AND_MESH) {
     /*--- Register node coordinates as input ---*/
 
-    geometry[iZone][iInst][MESH_0]->RegisterCoordinates();
+    geometry0->RegisterCoordinates();
   }
 
   if (kind_recording == RECORDING::MESH_DEFORM) {
     /*--- Undeformed mesh coordinates ---*/
-    solver[iZone][iInst][MESH_0][ADJMESH_SOL]->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJMESH_SOL]->RegisterSolution(geometry0, config[iZone]);
 
     /*--- Boundary displacements ---*/
-    solver[iZone][iInst][MESH_0][ADJMESH_SOL]->RegisterVariables(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJMESH_SOL]->RegisterVariables(geometry0, config[iZone]);
   }
   }
   END_SU2_OMP_PARALLEL
@@ -527,6 +444,8 @@ void CDiscAdjFluidIteration::RegisterInput(CSolver***** solver, CGeometry**** ge
 void CDiscAdjFluidIteration::SetDependencies(CSolver***** solver, CGeometry**** geometry, CNumerics****** numerics,
                                              CConfig** config, unsigned short iZone, unsigned short iInst,
                                              RECORDING kind_recording) {
+  auto solvers0 = solver[iZone][iInst][MESH_0];
+  auto geometry0 = geometry[iZone][iInst][MESH_0];
 
   if ((kind_recording == RECORDING::MESH_COORDS) ||
       (kind_recording == RECORDING::CLEAR_INDICES) ||
@@ -538,71 +457,68 @@ void CDiscAdjFluidIteration::SetDependencies(CSolver***** solver, CGeometry****
     CGeometry::ComputeWallDistance(config, geometry);
   }
 
-  SU2_OMP_PARALLEL_(if(solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->GetHasHybridParallel())) {
+  SU2_OMP_PARALLEL_(if(solvers0[ADJFLOW_SOL]->GetHasHybridParallel())) {
 
   /*--- Compute coupling between flow, turbulent and species equations ---*/
-  solver[iZone][iInst][MESH_0][FLOW_SOL]->Preprocessing(geometry[iZone][iInst][MESH_0], solver[iZone][iInst][MESH_0],
-                                                        config[iZone], MESH_0, NO_RK_ITER, RUNTIME_FLOW_SYS, true);
-  solver[iZone][iInst][MESH_0][FLOW_SOL]->InitiateComms(geometry[iZone][iInst][MESH_0], config[iZone], SOLUTION);
-  solver[iZone][iInst][MESH_0][FLOW_SOL]->CompleteComms(geometry[iZone][iInst][MESH_0], config[iZone], SOLUTION);
+  solvers0[FLOW_SOL]->Preprocessing(geometry0, solvers0, config[iZone], MESH_0, NO_RK_ITER, RUNTIME_FLOW_SYS, true);
+  solvers0[FLOW_SOL]->InitiateComms(geometry0, config[iZone], SOLUTION);
+  solvers0[FLOW_SOL]->CompleteComms(geometry0, config[iZone], SOLUTION);
 
   if (turbulent && !config[iZone]->GetFrozen_Visc_Disc()) {
-    solver[iZone][iInst][MESH_0][TURB_SOL]->Postprocessing(geometry[iZone][iInst][MESH_0], solver[iZone][iInst][MESH_0],
+    solvers0[TURB_SOL]->Postprocessing(geometry0, solvers0,
                                                            config[iZone], MESH_0);
-    solver[iZone][iInst][MESH_0][TURB_SOL]->InitiateComms(geometry[iZone][iInst][MESH_0], config[iZone], SOLUTION);
-    solver[iZone][iInst][MESH_0][TURB_SOL]->CompleteComms(geometry[iZone][iInst][MESH_0], config[iZone], SOLUTION);
+    solvers0[TURB_SOL]->InitiateComms(geometry0, config[iZone], SOLUTION);
+    solvers0[TURB_SOL]->CompleteComms(geometry0, config[iZone], SOLUTION);
   }
   if (config[iZone]->GetKind_Species_Model() != SPECIES_MODEL::NONE) {
-    solver[iZone][iInst][MESH_0][SPECIES_SOL]->Preprocessing(geometry[iZone][iInst][MESH_0], solver[iZone][iInst][MESH_0],
-                                                             config[iZone], MESH_0, NO_RK_ITER, RUNTIME_FLOW_SYS, true);
-    solver[iZone][iInst][MESH_0][SPECIES_SOL]->InitiateComms(geometry[iZone][iInst][MESH_0], config[iZone], SOLUTION);
-    solver[iZone][iInst][MESH_0][SPECIES_SOL]->CompleteComms(geometry[iZone][iInst][MESH_0], config[iZone], SOLUTION);
+    solvers0[SPECIES_SOL]->Preprocessing(geometry0, solvers0, config[iZone], MESH_0, NO_RK_ITER, RUNTIME_FLOW_SYS, true);
+    solvers0[SPECIES_SOL]->InitiateComms(geometry0, config[iZone], SOLUTION);
+    solvers0[SPECIES_SOL]->CompleteComms(geometry0, config[iZone], SOLUTION);
   }
 
   }
   END_SU2_OMP_PARALLEL
 
   if (config[iZone]->GetWeakly_Coupled_Heat()) {
-    solver[iZone][iInst][MESH_0][HEAT_SOL]->Set_Heatflux_Areas(geometry[iZone][iInst][MESH_0], config[iZone]);
-    solver[iZone][iInst][MESH_0][HEAT_SOL]->Preprocessing(geometry[iZone][iInst][MESH_0], solver[iZone][iInst][MESH_0],
-                                                          config[iZone], MESH_0, NO_RK_ITER, RUNTIME_HEAT_SYS, true);
-    solver[iZone][iInst][MESH_0][HEAT_SOL]->Postprocessing(geometry[iZone][iInst][MESH_0], solver[iZone][iInst][MESH_0],
-                                                           config[iZone], MESH_0);
-    solver[iZone][iInst][MESH_0][HEAT_SOL]->InitiateComms(geometry[iZone][iInst][MESH_0], config[iZone], SOLUTION);
-    solver[iZone][iInst][MESH_0][HEAT_SOL]->CompleteComms(geometry[iZone][iInst][MESH_0], config[iZone], SOLUTION);
+    solvers0[HEAT_SOL]->Set_Heatflux_Areas(geometry0, config[iZone]);
+    solvers0[HEAT_SOL]->Preprocessing(geometry0, solvers0, config[iZone], MESH_0, NO_RK_ITER, RUNTIME_HEAT_SYS, true);
+    solvers0[HEAT_SOL]->Postprocessing(geometry0, solvers0, config[iZone], MESH_0);
+    solvers0[HEAT_SOL]->InitiateComms(geometry0, config[iZone], SOLUTION);
+    solvers0[HEAT_SOL]->CompleteComms(geometry0, config[iZone], SOLUTION);
   }
   if (config[iZone]->AddRadiation()) {
-    solver[iZone][iInst][MESH_0][RAD_SOL]->Postprocessing(geometry[iZone][iInst][MESH_0], solver[iZone][iInst][MESH_0],
-                                                          config[iZone], MESH_0);
-    solver[iZone][iInst][MESH_0][RAD_SOL]->InitiateComms(geometry[iZone][iInst][MESH_0], config[iZone], SOLUTION);
-    solver[iZone][iInst][MESH_0][RAD_SOL]->CompleteComms(geometry[iZone][iInst][MESH_0], config[iZone], SOLUTION);
+    solvers0[RAD_SOL]->Postprocessing(geometry0, solvers0, config[iZone], MESH_0);
+    solvers0[RAD_SOL]->InitiateComms(geometry0, config[iZone], SOLUTION);
+    solvers0[RAD_SOL]->CompleteComms(geometry0, config[iZone], SOLUTION);
   }
 }
 
 void CDiscAdjFluidIteration::RegisterOutput(CSolver***** solver, CGeometry**** geometry, CConfig** config,
                                             unsigned short iZone, unsigned short iInst) {
+  auto solvers0 = solver[iZone][iInst][MESH_0];
+  auto geometry0 = geometry[iZone][iInst][MESH_0];
 
-  SU2_OMP_PARALLEL_(if(solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->GetHasHybridParallel())) {
+  SU2_OMP_PARALLEL_(if(solvers0[ADJFLOW_SOL]->GetHasHybridParallel())) {
 
   /*--- Register conservative variables as output of the iteration ---*/
 
   if (config[iZone]->GetFluidProblem()) {
-    solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJFLOW_SOL]->RegisterOutput(geometry0, config[iZone]);
   }
   if (turbulent && !config[iZone]->GetFrozen_Visc_Disc()) {
-    solver[iZone][iInst][MESH_0][ADJTURB_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJTURB_SOL]->RegisterOutput(geometry0, config[iZone]);
   }
   if (config[iZone]->GetKind_Species_Model() != SPECIES_MODEL::NONE) {
-    solver[iZone][iInst][MESH_0][ADJSPECIES_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJSPECIES_SOL]->RegisterOutput(geometry0, config[iZone]);
   }
   if (config[iZone]->GetWeakly_Coupled_Heat()) {
-    solver[iZone][iInst][MESH_0][ADJHEAT_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJHEAT_SOL]->RegisterOutput(geometry0, config[iZone]);
   }
   if (config[iZone]->AddRadiation()) {
-    solver[iZone][iInst][MESH_0][ADJRAD_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJRAD_SOL]->RegisterOutput(geometry0, config[iZone]);
   }
   if (config[iZone]->GetFluidProblem() && !config[iZone]->GetMultizone_Problem()) {
-    solver[iZone][iInst][MESH_0][FLOW_SOL]->RegisterVertexTractions(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[FLOW_SOL]->RegisterVertexTractions(geometry0, config[iZone]);
   }
 
   }
diff --git a/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp
index a12cb43bdc07..a7935c92b19d 100644
--- a/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp
+++ b/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp
@@ -43,79 +43,64 @@ void CDiscAdjHeatIteration::Preprocess(COutput* output, CIntegration**** integra
   /*--- For the unsteady adjoint, load direct solutions from restart files. ---*/
 
   if (config[val_iZone]->GetTime_Marching() != TIME_MARCHING::STEADY) {
-    const int Direct_Iter = static_cast<int>(config[val_iZone]->GetUnst_AdjointIter()) - static_cast<int>(TimeIter) - 2 + dual_time;
+    const int Direct_Iter = static_cast<int>(config[val_iZone]->GetUnst_AdjointIter()) -
+                            static_cast<int>(TimeIter) - 2 + dual_time;
 
-    /*--- For dual-time stepping we want to load the already converged solution at timestep n ---*/
+    /*--- For dual-time stepping we want to load the already converged solution at previous timesteps.
+     * In general we only load one file and shift the previously loaded solutions, on the first we
+     * load one or two more (depending on dual time order). ---*/
 
-    if (TimeIter == 0) {
-      if (dual_time_2nd) {
-        /*--- Load solution at timestep n-2 ---*/
-
-        LoadUnsteady_Solution(geometry, solver, config, val_iZone, val_iInst, Direct_Iter - 2);
-
-        /*--- Push solution back to correct array ---*/
+    if (dual_time_2nd) {
+      LoadUnsteady_Solution(geometry, solver, config, val_iZone, val_iInst, Direct_Iter - 2);
+    } else if (dual_time_1st) {
+      LoadUnsteady_Solution(geometry, solver, config, val_iZone, val_iInst, Direct_Iter - 1);
+    }
 
+    if (TimeIter == 0) {
+      /*--- Push solution back one level. ---*/
+      if (dual_time) {
         for (auto iMesh = 0u; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++) {
           solvers[iMesh][HEAT_SOL]->GetNodes()->Set_Solution_time_n();
-          solvers[iMesh][HEAT_SOL]->GetNodes()->Set_Solution_time_n1();
         }
       }
-      if (dual_time) {
-        /*--- Load solution at timestep n-1 ---*/
 
+      /*--- If required load another time step. Push the previous time step to n-1 and the
+       loaded time step to n. ---*/
+      if (dual_time_2nd) {
         LoadUnsteady_Solution(geometry, solver, config, val_iZone, val_iInst, Direct_Iter - 1);
 
-        /*--- Push solution back to correct array ---*/
-
         for (auto iMesh = 0u; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++) {
+          solvers[iMesh][HEAT_SOL]->GetNodes()->Set_Solution_time_n1();
           solvers[iMesh][HEAT_SOL]->GetNodes()->Set_Solution_time_n();
         }
       }
 
-      /*--- Load solution timestep n ---*/
-
+      /*--- Load current solution. ---*/
       LoadUnsteady_Solution(geometry, solver, config, val_iZone, val_iInst, Direct_Iter);
     }
     if ((TimeIter > 0) && dual_time) {
-      /*--- Load solution timestep n - 2 ---*/
-
-      LoadUnsteady_Solution(geometry, solver, config, val_iZone, val_iInst, Direct_Iter - 2);
-
       /*--- Temporarily store the loaded solution in the Solution_Old array ---*/
-
       for (auto iMesh = 0u; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++)
         solvers[iMesh][HEAT_SOL]->Set_OldSolution();
 
-      /*--- Set Solution at timestep n to solution at n-1 ---*/
-
+      /*--- Move timestep n to current solution. ---*/
       for (auto iMesh = 0u; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++) {
         for (auto iPoint = 0ul; iPoint < geometry[val_iZone][val_iInst][iMesh]->GetnPoint(); iPoint++) {
           solvers[iMesh][HEAT_SOL]->GetNodes()->SetSolution(
               iPoint, solvers[iMesh][HEAT_SOL]->GetNodes()->GetSolution_time_n(iPoint));
         }
       }
-      if (dual_time_1st) {
-        /*--- Set Solution at timestep n-1 to the previously loaded solution ---*/
-        for (auto iMesh = 0u; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++) {
-          for (auto iPoint = 0ul; iPoint < geometry[val_iZone][val_iInst][iMesh]->GetnPoint(); iPoint++) {
-            solvers[iMesh][HEAT_SOL]->GetNodes()->Set_Solution_time_n(
-                iPoint, solvers[iMesh][HEAT_SOL]->GetNodes()->GetSolution_time_n1(iPoint));
-          }
-        }
-      }
-      if (dual_time_2nd) {
-        /*--- Set Solution at timestep n-1 to solution at n-2 ---*/
-        for (auto iMesh = 0u; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++) {
-          for (auto iPoint = 0ul; iPoint < geometry[val_iZone][val_iInst][iMesh]->GetnPoint(); iPoint++) {
-            solvers[iMesh][HEAT_SOL]->GetNodes()->Set_Solution_time_n(
-                iPoint, solvers[iMesh][HEAT_SOL]->GetNodes()->GetSolution_time_n1(iPoint));
-          }
-        }
-        /*--- Set Solution at timestep n-2 to the previously loaded solution ---*/
-        for (auto iMesh = 0u; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++) {
-          for (auto iPoint = 0ul; iPoint < geometry[val_iZone][val_iInst][iMesh]->GetnPoint(); iPoint++) {
-            solvers[iMesh][HEAT_SOL]->GetNodes()->Set_Solution_time_n1(
-                iPoint, solvers[iMesh][HEAT_SOL]->GetNodes()->GetSolution_Old(iPoint));
+
+      /*--- Finally, place the loaded solution in the correct place (n or n-1 depending on order). ---*/
+      for (auto iMesh = 0u; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++) {
+        auto* heatSol = solvers[iMesh][HEAT_SOL];
+        for (auto iPoint = 0ul; iPoint < geometry[val_iZone][val_iInst][iMesh]->GetnPoint(); iPoint++) {
+          if (dual_time_2nd) {
+            /*--- If required also move timestep n-1 to timestep n. ---*/
+            heatSol->GetNodes()->Set_Solution_time_n(iPoint, heatSol->GetNodes()->GetSolution_time_n1(iPoint));
+            heatSol->GetNodes()->Set_Solution_time_n1(iPoint, heatSol->GetNodes()->GetSolution_Old(iPoint));
+          } else {
+            heatSol->GetNodes()->Set_Solution_time_n(iPoint, heatSol->GetNodes()->GetSolution_Old(iPoint));
           }
         }
       }
@@ -140,13 +125,14 @@ void CDiscAdjHeatIteration::Preprocess(COutput* output, CIntegration**** integra
 void CDiscAdjHeatIteration::LoadUnsteady_Solution(CGeometry**** geometry, CSolver***** solver, CConfig** config,
                                                   unsigned short val_iZone, unsigned short val_iInst,
                                                   int val_DirectIter) {
+  auto solvers = solver[val_iZone][val_iInst];
+  auto geometries = geometry[val_iZone][val_iInst];
 
   if (val_DirectIter >= 0) {
     if (rank == MASTER_NODE)
       cout << " Loading heat solution from direct iteration " << val_DirectIter << " for zone " << val_iZone << "." << endl;
 
-    solver[val_iZone][val_iInst][MESH_0][HEAT_SOL]->LoadRestart(
-        geometry[val_iZone][val_iInst], solver[val_iZone][val_iInst], config[val_iZone], val_DirectIter, false);
+    solvers[MESH_0][HEAT_SOL]->LoadRestart(geometries, solvers, config[val_iZone], val_DirectIter, false);
   }
   else {
     /*--- If there is no solution file we set the freestream condition ---*/
@@ -154,9 +140,8 @@ void CDiscAdjHeatIteration::LoadUnsteady_Solution(CGeometry**** geometry, CSolve
       cout << " Setting freestream conditions at direct iteration " << val_DirectIter << " for zone " << val_iZone << "." << endl;
 
     for (auto iMesh = 0u; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++) {
-      solver[val_iZone][val_iInst][iMesh][HEAT_SOL]->SetFreeStream_Solution(config[val_iZone]);
-      solver[val_iZone][val_iInst][iMesh][HEAT_SOL]->Postprocessing(
-          geometry[val_iZone][val_iInst][iMesh], solver[val_iZone][val_iInst][iMesh], config[val_iZone], iMesh);
+      solvers[iMesh][HEAT_SOL]->SetFreeStream_Solution(config[val_iZone]);
+      solvers[iMesh][HEAT_SOL]->Postprocessing(geometries[iMesh], solvers[iMesh], config[val_iZone], iMesh);
     }
   }
 }
@@ -177,26 +162,28 @@ void CDiscAdjHeatIteration::InitializeAdjoint(CSolver***** solver, CGeometry****
 
 void CDiscAdjHeatIteration::RegisterInput(CSolver***** solver, CGeometry**** geometry, CConfig** config,
                                           unsigned short iZone, unsigned short iInst, RECORDING kind_recording) {
+  auto solvers0 = solver[iZone][iInst][MESH_0];
+  auto geometry0 = geometry[iZone][iInst][MESH_0];
 
   if (kind_recording == RECORDING::SOLUTION_VARIABLES || kind_recording == RECORDING::SOLUTION_AND_MESH) {
     /*--- Register flow and turbulent variables as input ---*/
 
-    solver[iZone][iInst][MESH_0][ADJHEAT_SOL]->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJHEAT_SOL]->RegisterSolution(geometry0, config[iZone]);
 
-    solver[iZone][iInst][MESH_0][ADJHEAT_SOL]->RegisterVariables(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJHEAT_SOL]->RegisterVariables(geometry0, config[iZone]);
   }
   else if (kind_recording == RECORDING::MESH_COORDS) {
     /*--- Register node coordinates as input ---*/
 
-    geometry[iZone][iInst][MESH_0]->RegisterCoordinates();
+    geometry0->RegisterCoordinates();
   }
   else if (kind_recording == RECORDING::MESH_DEFORM) {
     /*--- Register the variables of the mesh deformation ---*/
     /*--- Undeformed mesh coordinates ---*/
-    solver[iZone][iInst][MESH_0][ADJMESH_SOL]->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJMESH_SOL]->RegisterSolution(geometry0, config[iZone]);
 
     /*--- Boundary displacements ---*/
-    solver[iZone][iInst][MESH_0][ADJMESH_SOL]->RegisterVariables(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJMESH_SOL]->RegisterVariables(geometry0, config[iZone]);
   }
 }
 
diff --git a/SU2_CFD/src/output/CAdjElasticityOutput.cpp b/SU2_CFD/src/output/CAdjElasticityOutput.cpp
index e30f2895f654..d27b96c990c7 100644
--- a/SU2_CFD/src/output/CAdjElasticityOutput.cpp
+++ b/SU2_CFD/src/output/CAdjElasticityOutput.cpp
@@ -244,19 +244,19 @@ void CAdjElasticityOutput::SetVolumeOutputFields(CConfig *config){
 
   /*--- Sensitivities with respect to initial conditions. ---*/
 
-  AddVolumeOutput("SENS_DISP-X", "SensInitialDisp_x", "SENSITIVITY_T0", "sensitivity to the initial x displacement");
-  AddVolumeOutput("SENS_DISP-Y", "SensInitialDisp_y", "SENSITIVITY_T0", "sensitivity to the initial y displacement");
+  AddVolumeOutput("SENS_DISP-X", "SensitivityDispN_x", "SENSITIVITY_N", "sensitivity to the previous x displacement");
+  AddVolumeOutput("SENS_DISP-Y", "SensitivityDispN_y", "SENSITIVITY_N", "sensitivity to the previous y displacement");
   if (nDim == 3)
-    AddVolumeOutput("SENS_DISP-Z", "SensInitialDisp_z", "SENSITIVITY_T0", "sensitivity to the initial z displacement");
+    AddVolumeOutput("SENS_DISP-Z", "SensitivityDispN_z", "SENSITIVITY_N", "sensitivity to the previous z displacement");
 
-  AddVolumeOutput("SENS_VEL-X", "SensInitialVel_x", "SENSITIVITY_T0", "sensitivity to the initial x velocity");
-  AddVolumeOutput("SENS_VEL-Y", "SensInitialVel_y", "SENSITIVITY_T0", "sensitivity to the initial y velocity");
+  AddVolumeOutput("SENS_VEL-X", "SensitivityVelN_x", "SENSITIVITY_N", "sensitivity to the previous x velocity");
+  AddVolumeOutput("SENS_VEL-Y", "SensitivityVelN_y", "SENSITIVITY_N", "sensitivity to the previous y velocity");
   if (nDim == 3)
-    AddVolumeOutput("SENS_VEL-Z", "SensInitialVel_z", "SENSITIVITY_T0", "sensitivity to the initial z velocity");
+    AddVolumeOutput("SENS_VEL-Z", "SensitivityVelN_z", "SENSITIVITY_N", "sensitivity to the previous z velocity");
 
-  AddVolumeOutput("SENS_ACCEL-X", "SensInitialAccel_x", "SENSITIVITY_T0", "sensitivity to the initial x acceleration");
-  AddVolumeOutput("SENS_ACCEL-Y", "SensInitialAccel_y", "SENSITIVITY_T0", "sensitivity to the initial y acceleration");
+  AddVolumeOutput("SENS_ACCEL-X", "SensitivityAccelN_x", "SENSITIVITY_N", "sensitivity to the previous x acceleration");
+  AddVolumeOutput("SENS_ACCEL-Y", "SensitivityAccelN_y", "SENSITIVITY_N", "sensitivity to the previous y acceleration");
   if (nDim == 3)
-    AddVolumeOutput("SENS_ACCEL-Z", "SensInitialAccel_z", "SENSITIVITY_T0", "sensitivity to the initial z acceleration");
+    AddVolumeOutput("SENS_ACCEL-Z", "SensitivityAccelN_z", "SENSITIVITY_N", "sensitivity to the previous z acceleration");
 
 }
diff --git a/SU2_CFD/src/output/CAdjHeatOutput.cpp b/SU2_CFD/src/output/CAdjHeatOutput.cpp
index 65d5679b8cad..889275f436b4 100644
--- a/SU2_CFD/src/output/CAdjHeatOutput.cpp
+++ b/SU2_CFD/src/output/CAdjHeatOutput.cpp
@@ -176,6 +176,14 @@ void CAdjHeatOutput::SetVolumeOutputFields(CConfig *config){
   AddVolumeOutput("SENSITIVITY", "Surface_Sensitivity", "SENSITIVITY", "sensitivity in normal direction");
   /// END_GROUP
 
+  if (!config->GetTime_Domain()) return;
+
+  /*--- Sensitivities with respect to initial conditions. ---*/
+
+  AddVolumeOutput("SENS_TEMP_N", "SensitivityTempN", "SENSITIVITY_N", "sensitivity to the previous temperature");
+  if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) {
+    AddVolumeOutput("SENS_TEMP_N1", "SensitivityTempN1", "SENSITIVITY_N", "sensitivity to the previous-1 temperature");
+  }
 }
 
 void CAdjHeatOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint){
@@ -199,6 +207,13 @@ void CAdjHeatOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve
   if (nDim == 3)
     SetVolumeOutputValue("SENSITIVITY-Z", iPoint, Node_AdjHeat->GetSensitivity(iPoint, 2));
 
+  if (!config->GetTime_Domain()) return;
+
+  SetVolumeOutputValue("SENS_TEMP_N", iPoint, Node_AdjHeat->GetSolution_time_n(iPoint, 0));
+  if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) {
+    SetVolumeOutputValue("SENS_TEMP_N1", iPoint, Node_AdjHeat->GetSolution_time_n1(iPoint, 0));
+  }
+
 }
 
 void CAdjHeatOutput::LoadSurfaceData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint, unsigned short iMarker, unsigned long iVertex){
diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py
index 61c79202fb80..8efd8f7a6d72 100644
--- a/TestCases/parallel_regression_AD.py
+++ b/TestCases/parallel_regression_AD.py
@@ -421,6 +421,19 @@ def main():
     test_list.append(pywrapper_Unst_FEA_AD)
     pass_list.append(pywrapper_Unst_FEA_AD.run_test())
 
+    # Heat solver unsteady AD
+    pywrapper_Unst_Heat_AD = TestCase('pywrapper_Unst_Heat_AD')
+    pywrapper_Unst_Heat_AD.cfg_dir = "py_wrapper/custom_heat_flux"
+    pywrapper_Unst_Heat_AD.cfg_file = "run_ad.py"
+    pywrapper_Unst_Heat_AD.test_iter = 100
+    pywrapper_Unst_Heat_AD.test_vals = [0.776365, 0.776430, 1.000003]
+    pywrapper_Unst_Heat_AD.command = TestCase.Command("mpirun -n 2", "python", "run_ad.py")
+    pywrapper_Unst_Heat_AD.timeout = 1600
+    pywrapper_Unst_Heat_AD.tol = 0.00001
+    pywrapper_Unst_Heat_AD.new_output = False
+    test_list.append(pywrapper_Unst_Heat_AD)
+    pass_list.append(pywrapper_Unst_Heat_AD.run_test())
+
     # Flow AD Mesh Displacement Sensitivity
     pywrapper_CFD_AD_MeshDisp               = TestCase('pywrapper_CFD_AD_MeshDisp')
     pywrapper_CFD_AD_MeshDisp.cfg_dir       = "py_wrapper/disc_adj_flow/mesh_disp_sens"
diff --git a/TestCases/py_wrapper/custom_heat_flux/run_ad.py b/TestCases/py_wrapper/custom_heat_flux/run_ad.py
new file mode 100644
index 000000000000..15d2ed8d845c
--- /dev/null
+++ b/TestCases/py_wrapper/custom_heat_flux/run_ad.py
@@ -0,0 +1,246 @@
+#!/usr/bin/env python
+
+## \file run.py
+#  \brief Unsteady adjoint heat transfer case with custom heat flux.
+#  \version 7.5.1 "Blackbird"
+#
+# SU2 Project Website: https://su2code.github.io
+#
+# The SU2 Project is maintained by the SU2 Foundation
+# (http://su2foundation.org)
+#
+# Copyright 2012-2023, SU2 Contributors (cf. AUTHORS.md)
+#
+# SU2 is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License, or (at your option) any later version.
+#
+# SU2 is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+
+import pysu2ad
+from mpi4py import MPI
+
+common_settings = """
+SOLVER= HEAT_EQUATION
+INC_NONDIM= DIMENSIONAL
+
+TIME_DOMAIN= YES
+TIME_STEP= 0.05
+TIME_MARCHING= DUAL_TIME_STEPPING-2ND_ORDER
+
+FREESTREAM_TEMPERATURE= 345
+MATERIAL_DENSITY= 2700
+SPECIFIC_HEAT_CP = 870
+% Very high value to make the case more interesting.
+THERMAL_CONDUCTIVITY_CONSTANT= 1000
+
+MARKER_HEATFLUX= ( x_minus, 0, x_plus, 0, y_minus, 0, y_plus, 0 )
+MARKER_PYTHON_CUSTOM= ( x_minus, x_plus, y_minus, y_plus )
+MARKER_MONITORING= ( x_minus, x_plus, y_minus, y_plus )
+
+NUM_METHOD_GRAD= GREEN_GAUSS
+CFL_NUMBER= 1e8
+
+LINEAR_SOLVER= CONJUGATE_GRADIENT
+DISCADJ_LIN_SOLVER= CONJUGATE_GRADIENT
+LINEAR_SOLVER_PREC= ILU
+DISCADJ_LIN_PREC= ILU
+LINEAR_SOLVER_ERROR= 1E-5
+LINEAR_SOLVER_ITER= 100
+
+MESH_FORMAT= RECTANGLE
+MESH_BOX_SIZE= ( 33, 33, 0 )
+MESH_BOX_LENGTH= ( __SIZE__, __SIZE__, 0 )
+
+OUTPUT_FILES= RESTART, PARAVIEW
+OUTPUT_WRT_FREQ= 1
+OBJECTIVE_FUNCTION= AVG_TEMPERATURE
+
+INNER_ITER= 20
+CONV_RESIDUAL_MINVAL= -4
+CONV_STARTITER= 2
+MAX_TIME= 2.0
+TIME_ITER= 41
+"""
+
+primal_settings = """
+MATH_PROBLEM= DIRECT
+SCREEN_OUTPUT= TIME_ITER, CUR_TIME, INNER_ITER, RMS_RES, LINSOL, AVG_TEMPERATURE, TAVG_AVG_TEMPERATURE
+HISTORY_OUTPUT= ITER, RMS_RES, HEAT, TAVG_HEAT
+"""
+
+adjoint_settings = """
+MATH_PROBLEM= DISCRETE_ADJOINT
+SCREEN_OUTPUT= TIME_ITER, CUR_TIME, INNER_ITER, RMS_RES, LINSOL
+UNST_ADJOINT_ITER= 41
+% Only interested in the final value of the objective.
+ITER_AVERAGE_OBJ= 1
+SOLUTION_FILENAME= restart.dat
+VOLUME_OUTPUT= SOLUTION, SENSITIVITY, SENSITIVITY_N
+"""
+
+
+def ApplyHeatFlux(time, driver, marker_ids):
+  """Applies a heat flux on all boundaries for a period of time and then removes it."""
+  for marker_id in marker_ids:
+    if marker_id < 0:
+      continue
+    hf = (-1, 0)[time > 0.2]
+    for i_vertex in range(driver.GetNumberMarkerNodes(marker_id)):
+      driver.SetMarkerCustomNormalHeatFlux(marker_id, i_vertex, hf)
+
+
+def RunPrimal(size):
+  """
+  Run the heat solver with a custom heat (function of time) flux on all boundaries.
+  Returns the final average boundary temperature.
+  """
+  comm = MPI.COMM_WORLD
+
+  with open('config_unsteady.cfg', 'w') as f:
+    f.write(common_settings.replace('__SIZE__', str(size)) + primal_settings)
+
+  # Initialize the primal driver of SU2, this includes solver preprocessing.
+  try:
+    driver = pysu2ad.CSinglezoneDriver('config_unsteady.cfg', 1, comm)
+  except TypeError as exception:
+    print('A TypeError occured in pysu2ad.CSinglezoneDriver : ', exception)
+    raise
+
+  # Get the ID of the markers where the heat flux is applied.
+  all_marker_ids = driver.GetMarkerIndices()
+  marker_names = ['x_minus', 'x_plus', 'y_minus', 'y_plus']
+  marker_ids = []
+  for name in marker_names:
+    marker_ids.append(all_marker_ids[name] if name in all_marker_ids else -1)
+
+  # Run the time loop in python to vary the heat flux.
+  dt = driver.GetUnsteadyTimeStep()
+
+  for time_iter in range(driver.GetNumberTimeIter()):
+    # Custom heat flux.
+    ApplyHeatFlux(time_iter * dt, driver, marker_ids)
+
+    driver.Preprocess(time_iter)
+
+    # Run one time iteration.
+    driver.Run()
+    driver.Postprocess()
+    driver.Update()
+
+    # Monitor the solver and output solution to file if required.
+    driver.Monitor(time_iter)
+    driver.Output(time_iter)
+
+  # Get the final average temperature.
+  avg_temperature = driver.GetOutputValue('AVG_TEMPERATURE')
+
+  # Finalize the solver and exit cleanly.
+  driver.Finalize()
+
+  return avg_temperature
+
+
+def RunAdjoint(size):
+  """
+  Runs the adjoint heat solver and returns the sensitivity of the objective function to
+  size of the domain and to the initial temperature.
+  """
+  comm = MPI.COMM_WORLD
+
+  with open('config_unsteady_ad.cfg', 'w') as f:
+    f.write(common_settings.replace('__SIZE__', str(size)) + adjoint_settings)
+
+  # Initialize the adjoint driver of SU2, this includes solver preprocessing.
+  try:
+    driver = pysu2ad.CDiscAdjSinglezoneDriver('config_unsteady_ad.cfg', 1, comm)
+  except TypeError as exception:
+    print('A TypeError occured in pysu2ad.CDiscAdjSinglezoneDriver : ', exception)
+    raise
+
+  # Get the ID of the markers where the heat flux is applied.
+  all_marker_ids = driver.GetMarkerIndices()
+  marker_names = ['x_minus', 'x_plus', 'y_minus', 'y_plus']
+  marker_ids = []
+  for name in marker_names:
+    marker_ids.append(all_marker_ids[name] if name in all_marker_ids else -1)
+
+  # Run the time loop in python to vary the heat flux.
+  dt = driver.GetUnsteadyTimeStep()
+
+  # Run the time loop in python to extract sensitivities at each step.
+  for time_iter in range(driver.GetNumberTimeIter()):
+    # Note that time runs in reverse for the adjoint solver.
+    ApplyHeatFlux((driver.GetNumberTimeIter() - time_iter - 1) * dt, driver, marker_ids)
+
+    # Preprocess adjoint iteration (AD recording).
+    driver.Preprocess(time_iter)
+
+    # Run one time iteration.
+    driver.Run()
+    driver.Postprocess()
+    driver.Update()
+
+    # Monitor the solver and output solution to file if required.
+    driver.Monitor(time_iter)
+    driver.Output(time_iter)
+
+  size_sens = 0.0
+  temp_sens = 0.0
+  sensitivity = driver.Sensitivity(driver.GetSolverIndices()['ADJ.HEAT'])
+  sens_n = driver.SolutionTimeN(driver.GetSolverIndices()['ADJ.HEAT'])
+  sens_n1 = driver.SolutionTimeN1(driver.GetSolverIndices()['ADJ.HEAT'])
+  coords = driver.Coordinates()
+
+  for i_node in range(driver.GetNumberNodes() - driver.GetNumberHaloNodes()):
+    x, y = coords.Get(i_node)
+    dx_ds = x / size
+    dy_ds = y / size
+    size_sens += dx_ds * sensitivity(i_node, 0) + dy_ds * sensitivity(i_node, 1)
+    # For 2nd order BDF the initial conditions apply to time n and n-1,
+    # hence we combine the two sensitivities. n-1 is multiplied by 2 because
+    # the initial conditions are used as the n-1 solution twice (in the 1st and
+    # second iterations). However, note that this is an approximation, the
+    # correct value would be sens_n_t0 + sens_n1_t0 + sens_n1_t1 (we approximate
+    # sens_n1 at t1 as sens_n1 at t0).
+    temp_sens += sens_n(i_node, 0) + 2 * sens_n1(i_node, 0)
+
+  # Finalize the solver and exit cleanly.
+  driver.Finalize()
+
+  return comm.allreduce(size_sens), comm.allreduce(temp_sens)
+
+
+def main():
+  comm = MPI.COMM_WORLD
+  rank = comm.Get_rank()
+
+  obj_pert_size = RunPrimal(0.100005)
+  obj = RunPrimal(0.1)
+  sens_size_fd = (obj_pert_size - obj) / 0.000005
+
+  sens_size, sens_temp = RunAdjoint(0.1)
+
+  if rank == 0:
+    print("     Finite Differences\tDiscrete Adjoint")
+    print(f"Size {sens_size_fd}\t{sens_size}")
+
+  assert abs(sens_size / sens_size_fd - 1) < 1e-4, "Error in geometric derivatives."
+  # We expect the final average temperature to be directly proportional to the initial
+  # temperature since the applied heat flux is not a function of temperature.
+  assert abs(sens_temp - 1) < 1e-5, "Error in initial condition derivatives."
+
+  # Print results for the regression script to check.
+  if rank == 0:
+    print("\n------------------------------ Begin Solver -----------------------------")
+    print(100, 100, sens_size_fd / 100, sens_size / 100, sens_temp)
+
+if __name__ == '__main__':
+  main()
diff --git a/TestCases/py_wrapper/custom_load_fea/run_ad.py b/TestCases/py_wrapper/custom_load_fea/run_ad.py
index dd64982b51b4..8a44b1e033f1 100644
--- a/TestCases/py_wrapper/custom_load_fea/run_ad.py
+++ b/TestCases/py_wrapper/custom_load_fea/run_ad.py
@@ -84,7 +84,7 @@
 UNST_ADJOINT_ITER= 21
 ITER_AVERAGE_OBJ= 0
 SOLUTION_FILENAME= restart.dat
-VOLUME_OUTPUT= SENSITIVITY_T0
+VOLUME_OUTPUT= SOLUTION, SENSITIVITY, SENSITIVITY_N
 """