diff --git a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/O_mesh_cylinder.geo b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/O_mesh_cylinder.geo
new file mode 100644
index 000000000000..b048536e3420
--- /dev/null
+++ b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/O_mesh_cylinder.geo
@@ -0,0 +1,147 @@
+// ----------------------------------------------------------------------------------- //
+// Kattmann 20.08.2018, O mesh for CHT vortex shedding behind cylinder
+// The O mesh around the cylinder consists out of two half cylinders.
+// The inner pin is hollow. 
+// ----------------------------------------------------------------------------------- //
+
+// Create fluid and solid mesh seperately and merge together e.g. by hand.
+
+// Which domain part should be handled
+Which_Mesh_Part= 1;// 0=all, 1=Fluid, 2=Solid
+// Evoke Meshing Algorithm?
+Do_Meshing= 1; // 0=false, 1=true
+// Write Mesh files in .su2 format
+Write_mesh= 1; // 0=false, 1=true
+
+//Geometric inputs
+cylinder_diameter = 1;
+cylinder_radius = cylinder_diameter/2;
+mesh_radius = 20 * cylinder_diameter;
+inner_pin_d = 0.5;
+inner_pin_r = inner_pin_d/2;
+
+// ----------------------------------------------------------------------------------- //
+//Mesh inputs
+gridsize = 0.1;
+Ncylinder = 40;
+Nradial = 50;
+Rradial = 1.15;
+
+NPinRadial = 10;
+RPinRadial = 0.91;
+
+// Each zone is self-sufficient (i.e. has all of its own Points/Lines etc.)
+// ----------------------------------------------------------------------------------- //
+// Fluid zone
+If (Which_Mesh_Part == 0 || Which_Mesh_Part == 1)
+
+    // Geometry definition
+    // Points
+    Point(1) = {-mesh_radius, 0, 0, gridsize};
+    Point(2) = {-cylinder_radius, 0, 0, gridsize};
+    Point(3) = {cylinder_radius, 0, 0, gridsize};
+    Point(4) = {mesh_radius, 0, 0, gridsize};
+    Point(5) = {0, 0, 0, gridsize};
+
+    //helping point to know height of first layer
+    //Point(6) = {-cylinder_radius - 0.002, 0, 0, gridsize};
+
+    // Lines
+    Line(1) = {1, 2}; // to the left
+    Line(2) = {3, 4}; // to the right
+
+    Circle(3) = {2, 5, 3}; // lower inner
+    Circle(4) = {1, 5, 4}; // lower outer
+    Circle(5) = {3, 5, 2}; // upper inner
+    Circle(6) = {4, 5, 1}; // upper outer
+
+    // Lineloops and surfaces
+    Line Loop(1) = {1, 3, 2, -4}; Plane Surface(1) = {1}; // lower half cylinder
+    Line Loop(2) = {1, -5, 2, 6}; Plane Surface(2) = {2}; // upper half cylinder
+
+    // ----------------------------------------------------------------------------------- //
+    // Mesh definition
+    // make structured mesh with transfinite Lines
+
+    // lower
+    Transfinite Line{3, 4} = Ncylinder;
+    Transfinite Line{-1, 2} = Nradial Using Progression Rradial;
+
+    // upper
+    Transfinite Line{-5, -6} = Ncylinder;
+    Transfinite Line{-1, 2} = Nradial Using Progression Rradial;
+
+    // Physical Groups
+    Physical Line("cylinder_fluid") = {3, 5};
+    Physical Line("farfield") = {4, 6};
+    Physical Surface("surface_mesh") = {1, 2};
+
+EndIf
+
+// ----------------------------------------------------------------------------------- //
+// Pin zone
+If (Which_Mesh_Part == 0 || Which_Mesh_Part == 2)
+
+     // Geometry definition
+    // Points
+    Point(11) = {-cylinder_radius, 0, 0, gridsize};
+    Point(12) = {-inner_pin_r, 0, 0, gridsize};
+    Point(13) = {inner_pin_r, 0, 0, gridsize};
+    Point(14) = {cylinder_radius, 0, 0, gridsize};
+    Point(15) = {0, 0, 0, gridsize};
+
+    // Lines
+    Line(11) = {11, 12}; // to the left
+    Line(12) = {13, 14}; // to the right
+
+    Circle(13) = {12, 15, 13}; // lower inner
+    Circle(14) = {11, 15, 14}; // lower outer
+    Circle(15) = {13, 15, 12}; // upper inner
+    Circle(16) = {14, 15, 11}; // upper outer
+
+    // Lineloops and surfaces
+    Line Loop(11) = {11, 13, 12, -14}; Plane Surface(11) = {11}; // lower half cylinder
+    Line Loop(12) = {11, -15, 12, 16}; Plane Surface(12) = {12}; // upper half cylinder
+
+    // ----------------------------------------------------------------------------------- //
+    // Mesh definition
+    // make structured mesh with transfinite Lines
+
+    // lower
+    Transfinite Line{13, 14} = Ncylinder;
+    Transfinite Line{-11, 12} = NPinRadial Using Progression RPinRadial;
+
+    // upper
+    Transfinite Line{-15, -16} = Ncylinder;
+    Transfinite Line{-11, 12} = NPinRadial Using Progression RPinRadial;
+
+    // Physical Groups
+    Physical Line("inner_pin") = {13, 15};
+    Physical Line("cylinder_solid") = {14, 16};
+    Physical Surface("surface_mesh") = {11, 12};
+
+EndIf
+
+// ----------------------------------------------------------------------------------- //
+Transfinite Surface "*";
+Recombine Surface "*";
+
+If (Do_Meshing == 1)
+    Mesh 1; Mesh 2;
+EndIf
+
+// ----------------------------------------------------------------------------------- //
+// Write .su2 meshfile
+If (Write_mesh == 1)
+
+    Mesh.Format = 42; // .su2 mesh format, 
+    If (Which_Mesh_Part == 1)
+        Save "fluid.su2";
+    ElseIf (Which_Mesh_Part == 2)
+        Save "solid.su2";
+    Else
+        Printf("Invalid Which_Mesh_Part variable.");
+        Abort;
+    EndIf
+
+EndIf
diff --git a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/README.md b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/README.md
new file mode 100644
index 000000000000..7023e9819430
--- /dev/null
+++ b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/README.md
@@ -0,0 +1,30 @@
+# Unsteady CHT Adjoint Testcase
+
+## Short Description
+This is a 2D cylinder in freestream testcase. The flow is incompressible and laminar at Re=200.
+A uniform vortex shedding forms behind the cylinder and each shedding cycle is resolved by 54 timesteps.
+The pin is heated on the inner surface.
+
+## Mesh
+The mesh is for testing purposes only and contains ~4000 elements for the flow and ~1000 for the heat domain.
+A gmsh .geo file is added such that the mesh can be recreated and modified.
+
+## Recreating full primal
+The primal for a full cycle can be restarted with the `solution_*_00000.dat` and `solution_*_00001.dat`.
+The primal solution is necessary for the Discrete Adjoint sweep and for the gradient of the full
+shedding cycle the full primal is necessary. The necessary changes to `chtMaster.cfg` are documented
+in the config itself.
+
+## Discrete Adjoint
+In the regression testcase of SU2 only 2 reverse steps are taken.
+For that, the solution files 52-55 for the adjoint are added.
+The objective Function is the average temperature on the inner pin surface, averaged over the full time.
+
+## Gradient validation via Finite Differences using FADO
+In order to validate the Discrete Adjoint gradient a Finite Differences python script `gradient_validation.py`
+using [FADO](www.github.com/su2code/FADO) is added.
+Note that the script can be used with the files as they are. Necessary adaptions are made by FADO itself.
+The script deforms the mesh and runs the primal for each of the 18 Design Variables.
+Afterwards the baseline mesh is evaluated and then the Discrete Adjoint.
+Use `postprocess.py` to print the absolute difference and relative difference in percent to screen.
+The relative differences in percent are <0.15% for all Design Variables (2021-05-14).
diff --git a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/chtMaster.cfg b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/chtMaster.cfg
new file mode 100644
index 000000000000..6821416134a5
--- /dev/null
+++ b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/chtMaster.cfg
@@ -0,0 +1,151 @@
+% 2021-03-11 TobiKattmann
+%
+SOLVER= MULTIPHYSICS
+%
+% Set RESTART_SOL=YES for primal runs including the FD sweep
+RESTART_SOL= NO
+RESTART_ITER= 2
+READ_BINARY_RESTART= YES
+SOLUTION_FILENAME= solution
+RESTART_FILENAME= solution
+%
+CONFIG_LIST = ( fluid.cfg, solid.cfg )
+%
+MARKER_ZONE_INTERFACE= ( cylinder_fluid, cylinder_solid )
+MARKER_CHT_INTERFACE=  ( cylinder_fluid, cylinder_solid )
+%
+% ------------------------- UNSTEADY SIMULATION -------------------------------%
+%
+TIME_DOMAIN= YES
+TIME_MARCHING= DUAL_TIME_STEPPING-2ND_ORDER
+%
+TIME_STEP= 500
+%
+MAX_TIME= 1e9
+% For a primal restart change TIME_ITER=56 for the correct number of steps. 54 is for the adjoint run.
+TIME_ITER= 54
+% For the primal (and therefore FD sweep) OUTER_ITER=200 is suitable.
+% For an accurate adjont run set OUTER_ITER=500. 100 is for the regression test.
+OUTER_ITER= 100
+%INNER_ITER= 1
+%
+UNST_ADJOINT_ITER= 56
+%
+ITER_AVERAGE_OBJ= 54
+%
+% ------------------------- INPUT/OUTPUT FILE INFORMATION --------------------------%
+%
+MESH_FILENAME= MeshCHT.su2
+%
+SCREEN_OUTPUT= (TIME_ITER, OUTER_ITER, BGS_ADJ_PRESSURE[0], BGS_ADJ_VELOCITY-X[0], BGS_ADJ_VELOCITY-Y[0], BGS_ADJ_TEMPERATURE[0], BGS_ADJ_TEMPERATURE[1], SENS_TEMP[0], SENS_GEO[0], SENS_GEO[1] )
+% Suitable output for primal simulations
+%SCREEN_OUTPUT= (TIME_ITER, OUTER_ITER, WALL_TIME, BGS_PRESSURE[0], BGS_TEMPERATURE[0], BGS_TEMPERATURE[1], DRAG[0], AVG_TEMPERATURE[1] )
+SCREEN_WRT_FREQ_OUTER= 50
+%
+HISTORY_OUTPUT= ( ITER, BGS_RES[0], RMS_RES[0], BGS_RES[1], RMS_RES[1],\
+ FLOW_COEFF[0], HEAT[0], AERO_COEFF[0], HEAT[1],\
+ LINSOL[0], LINSOL[1])
+%
+OUTPUT_FILES= (RESTART, PARAVIEW)
+OUTPUT_WRT_FREQ= 1
+VOLUME_FILENAME= flow
+WRT_PERFORMANCE= YES
+%
+SOLUTION_ADJ_FILENAME= solution_adj
+RESTART_ADJ_FILENAME= solution_adj
+VOLUME_ADJ_FILENAME= flow_adj
+%
+TABULAR_FORMAT= CSV
+GRAD_OBJFUNC_FILENAME= of_grad.csv
+OUTPUT_PRECISION=16
+%
+% -------------------- FREE-FORM DEFORMATION PARAMETERS -----------------------%
+%
+FFD_TOLERANCE= 1E-10
+FFD_ITERATIONS= 500
+%
+% FFD box definition: 2D case (FFD_BoxTag, X1, Y1, 0.0, X2, Y2, 0.0, X3, Y3, 0.0, X4, Y4, 0.0,
+%                              0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
+% Counterclockwise definition of FFD cornerpoints
+FFD_DEFINITION= (BOX,\
+ -0.6,-0.6,0.0,\
+  0.6,-0.6,0.0,\
+  0.6, 0.6,0.0,\
+ -0.6, 0.6,0.0,\
+ 0.0,0.0,0.0,  0.0,0.0,0.0  0.0,0.0,0.0,  0.0,0.0,0.0 )
+%
+% FFD box degree: 2D case (x_degree, y_degree, 0)
+FFD_DEGREE= (8, 1, 0)
+%
+% Surface grid continuity at the intersection with the faces of the FFD boxes.
+% To keep a particular level of surface continuity, SU2 automatically freezes the right
+% number of control point planes (NO_DERIVATIVE, 1ST_DERIVATIVE, 2ND_DERIVATIVE, USER_INPUT)
+FFD_CONTINUITY= NO_DERIVATIVE
+%
+% ----------------------- DESIGN VARIABLE PARAMETERS --------------------------%
+%
+%DV_KIND= FFD_SETTING
+% First 9 are upper, second 9 are lower DV's
+DV_KIND= FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D, FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D,FFD_CONTROL_POINT_2D
+%
+% Marker of the surface in which we are going apply the shape deformation
+DV_MARKER= ( cylinder_fluid, cylinder_solid )
+%
+% Parameters of the shape deformation
+% - FFD_SETTING ( 1.0 )
+% - FFD_CONTROL_POINT_2D ( FFD_BoxTag, i_Ind, j_Ind, x_Disp, y_Disp )
+%DV_PARAM= ( 1.0 )
+DV_PARAM= \
+( BOX, 0, 1, 0.0, 1.0);\
+( BOX, 1, 1, 0.0, 1.0);\
+( BOX, 2, 1, 0.0, 1.0);\
+( BOX, 3, 1, 0.0, 1.0);\
+( BOX, 4, 1, 0.0, 1.0);\
+( BOX, 5, 1, 0.0, 1.0);\
+( BOX, 6, 1, 0.0, 1.0);\
+( BOX, 7, 1, 0.0, 1.0);\
+( BOX, 8, 1, 0.0, 1.0);\
+( BOX, 0, 0, 0.0, 1.0);\
+( BOX, 1, 0, 0.0, 1.0);\
+( BOX, 2, 0, 0.0, 1.0);\
+( BOX, 3, 0, 0.0, 1.0);\
+( BOX, 4, 0, 0.0, 1.0);\
+( BOX, 5, 0, 0.0, 1.0);\
+( BOX, 6, 0, 0.0, 1.0);\
+( BOX, 7, 0, 0.0, 1.0);\
+( BOX, 8, 0, 0.0, 1.0)
+%
+% Value of the shape deformation
+DV_VALUE= 1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0, 1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0
+%
+% ------------------------ GRID DEFORMATION PARAMETERS ------------------------%
+%
+DEFORM_LINEAR_SOLVER= FGMRES
+DEFORM_LINEAR_SOLVER_PREC= ILU
+DEFORM_LINEAR_SOLVER_ERROR= 1E-14
+DEFORM_NONLINEAR_ITER= 1
+DEFORM_LINEAR_SOLVER_ITER= 1000
+%
+DEFORM_CONSOLE_OUTPUT= YES
+DEFORM_STIFFNESS_TYPE= WALL_DISTANCE
+%
+DEFINITION_DV= \
+( 19, 1.0 | cylinder_fluid, cylinder_solid | BOX, 0, 1, 0.0, 1.0 ); \
+( 19, 1.0 | cylinder_fluid, cylinder_solid | BOX, 1, 1, 0.0, 1.0 ); \
+( 19, 1.0 | cylinder_fluid, cylinder_solid | BOX, 2, 1, 0.0, 1.0 ); \
+( 19, 1.0 | cylinder_fluid, cylinder_solid | BOX, 3, 1, 0.0, 1.0 ); \
+( 19, 1.0 | cylinder_fluid, cylinder_solid | BOX, 4, 1, 0.0, 1.0 ); \
+( 19, 1.0 | cylinder_fluid, cylinder_solid | BOX, 5, 1, 0.0, 1.0 ); \
+( 19, 1.0 | cylinder_fluid, cylinder_solid | BOX, 6, 1, 0.0, 1.0 ); \
+( 19, 1.0 | cylinder_fluid, cylinder_solid | BOX, 7, 1, 0.0, 1.0 ); \
+( 19, 1.0 | cylinder_fluid, cylinder_solid | BOX, 8, 1, 0.0, 1.0 ); \
+( 19, 1.0 | cylinder_fluid, cylinder_solid | BOX, 0, 0, 0.0, 1.0 ); \
+( 19, 1.0 | cylinder_fluid, cylinder_solid | BOX, 1, 0, 0.0, 1.0 ); \
+( 19, 1.0 | cylinder_fluid, cylinder_solid | BOX, 2, 0, 0.0, 1.0 ); \
+( 19, 1.0 | cylinder_fluid, cylinder_solid | BOX, 3, 0, 0.0, 1.0 ); \
+( 19, 1.0 | cylinder_fluid, cylinder_solid | BOX, 4, 0, 0.0, 1.0 ); \
+( 19, 1.0 | cylinder_fluid, cylinder_solid | BOX, 5, 0, 0.0, 1.0 ); \
+( 19, 1.0 | cylinder_fluid, cylinder_solid | BOX, 6, 0, 0.0, 1.0 ); \
+( 19, 1.0 | cylinder_fluid, cylinder_solid | BOX, 7, 0, 0.0, 1.0 ); \
+( 19, 1.0 | cylinder_fluid, cylinder_solid | BOX, 8, 0, 0.0, 1.0 )
+
diff --git a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/fluid.cfg b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/fluid.cfg
new file mode 100644
index 000000000000..b4fe5b28efb8
--- /dev/null
+++ b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/fluid.cfg
@@ -0,0 +1,80 @@
+% TobiKattmann 2021-03-11
+%
+SOLVER= INC_NAVIER_STOKES
+%
+OBJECTIVE_FUNCTION= AVG_TEMPERATURE
+%OBJECTIVE_FUNCTION= DRAG
+OBJECTIVE_WEIGHT= 0.0
+%
+% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------%
+%
+INC_NONDIM= DIMENSIONAL
+% ~Water at room temperature
+INC_DENSITY_MODEL= CONSTANT
+INC_DENSITY_INIT= 1.0e3
+%
+% Re = rho*v*d/mu_dyn = 1e3*v*1/1e-3 = 200 => v = 2.0e-4
+INC_VELOCITY_INIT= ( 2.0e-4, 0.0, 0.0 )
+%
+INC_ENERGY_EQUATION= YES
+INC_TEMPERATURE_INIT= 288.15
+SPECIFIC_HEAT_CP= 4180
+%
+VISCOSITY_MODEL= CONSTANT_VISCOSITY
+MU_CONSTANT= 1.0e-3
+%
+% Pr = c_p * mu_dyn / lambda
+% lambda_water = 0.6 => Pr = 4180 * 1e-3 / 0.6 = 6.96
+CONDUCTIVITY_MODEL= CONSTANT_PRANDTL
+PRANDTL_LAM = 6.96
+%
+% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
+%
+MARKER_FAR= ( farfield )
+%MARKER_HEATFLUX= ( cylinder_fluid, 1e2 )
+%
+% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
+%
+NUM_METHOD_GRAD= GREEN_GAUSS
+%
+CFL_NUMBER= 1e3
+CFL_ADAPT= NO
+CFL_ADAPT_PARAM= ( 0.8, 1.4, 10.0, 1e10, 0.001 )
+%
+RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 )
+%
+% ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
+%
+LINEAR_SOLVER= FGMRES
+LINEAR_SOLVER_PREC= ILU
+%LINEAR_SOLVER_ERROR= 1E-4
+LINEAR_SOLVER_ERROR= 1E-8
+LINEAR_SOLVER_ITER= 20
+%
+CONV_NUM_METHOD_FLOW= FDS
+MUSCL_FLOW= YES
+SLOPE_LIMITER_FLOW= NONE
+TIME_DISCRE_FLOW= EULER_IMPLICIT
+%
+% --------------------------- CONVERGENCE PARAMETERS --------------------------%
+%
+CONV_RESIDUAL_MINVAL= -20
+CONV_STARTITER= 2000
+CONV_CAUCHY_ELEMS= 100
+CONV_CAUCHY_EPS= 1E-6
+%
+% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
+%
+%MESH_FILENAME= fluid.su2
+%
+%SCREEN_OUTPUT= (TIME_ITER, WALL_TIME, INNER_ITER, RMS_PRESSURE, RMS_VELOCITY-X, DRAG, LIFT, AVG_TEMPERATURE, TOTAL_HEATFLUX )
+%SCREEN_WRT_FREQ_INNER= 100
+%
+%HISTORY_OUTPUT= (ITER, RMS_RES, AERO_COEFF, HEAT, FLOW_COEFF, LINSOL)
+%HISTORY_WRT_FREQ_INNER= 100
+MARKER_MONITORING= ( cylinder_fluid )
+MARKER_ANALYZE= ( cylinder_fluid )
+MARKER_ANALYZE_AVERAGE= AREA
+%
+%OUTPUT_FILES= ( RESTART, PARAVIEW )
+
diff --git a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/gradient_validation.py b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/gradient_validation.py
new file mode 100644
index 000000000000..003f742ce090
--- /dev/null
+++ b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/gradient_validation.py
@@ -0,0 +1,139 @@
+# FADO script: Finite Differences of unsteady CHT and adjoint run
+
+from FADO import *
+
+# Design variables ----------------------------------------------------- #
+
+nDV = 18
+ffd = InputVariable(0.0,PreStringHandler("DV_VALUE="),nDV)
+
+# Parameters ----------------------------------------------------------- #
+
+# The master config `chtMaster.cfg` serves as an SU2 adjoint regression test.
+# For a correct gradient validation we need to exchange some options
+
+time_iter_primal = Parameter(["TIME_ITER= 56"],\
+                LabelReplacer("TIME_ITER= 54"))
+outer_iter_primal = Parameter(["OUTER_ITER= 200"],\
+                 LabelReplacer("OUTER_ITER= 100"))
+restart_sol_primal = Parameter(["RESTART_SOL= YES"],\
+                  LabelReplacer("RESTART_SOL= NO"))
+
+outer_iter_adjoint = Parameter(["OUTER_ITER= 500"],\
+                  LabelReplacer("OUTER_ITER= 100"))
+
+# Evaluations ---------------------------------------------------------- #
+
+# Note that correct SU2 version needs to be in PATH
+
+def_command = "SU2_DEF chtMaster.cfg"
+cfd_command = "mpirun -n 12 SU2_CFD chtMaster.cfg"
+
+cfd_ad_command = "mpirun -n 12 SU2_CFD_AD chtMaster.cfg"
+dot_ad_command = "mpirun -n 12 SU2_DOT_AD chtMaster.cfg"
+
+max_tries = 1
+
+# mesh deformation
+deform = ExternalRun("DEFORM",def_command,True) # True means sym links are used for addData
+deform.setMaxTries(max_tries)
+deform.addConfig("chtMaster.cfg")
+deform.addData("fluid.cfg") # zonal cfg's can be symlinked as they are unchanged
+deform.addData("solid.cfg")
+deform.addData("MeshCHT.su2")
+deform.addExpected("mesh_out.su2")
+
+# direct run
+direct = ExternalRun("DIRECT",cfd_command,True)
+direct.setMaxTries(max_tries)
+direct.addConfig("chtMaster.cfg")
+direct.addData("fluid.cfg")
+direct.addData("solid.cfg")
+direct.addData("DEFORM/mesh_out.su2",destination="MeshCHT.su2")
+direct.addData("solution_0_00000.dat")
+direct.addData("solution_0_00001.dat")
+direct.addData("solution_1_00000.dat")
+direct.addData("solution_1_00001.dat")
+direct.addExpected("solution_0_00055.dat")
+direct.addExpected("solution_1_00055.dat")
+direct.addParameter(time_iter_primal)
+direct.addParameter(outer_iter_primal)
+direct.addParameter(restart_sol_primal)
+
+# adjoint run
+adjoint = ExternalRun("ADJOINT",cfd_ad_command,True)
+adjoint.setMaxTries(max_tries)
+adjoint.addConfig("chtMaster.cfg")
+adjoint.addData("fluid.cfg") # zonal cfg's can be symlinked
+adjoint.addData("solid.cfg")
+adjoint.addData("DEFORM/mesh_out.su2", destination="MeshCHT.su2")
+# add all primal solution files
+for timeIter in range(56): #
+    if timeIter < 10:
+        timeIter = "0" + str(timeIter)
+    adjoint.addData("DIRECT/solution_0_000" + str(timeIter) + ".dat")
+    adjoint.addData("DIRECT/solution_1_000" + str(timeIter) + ".dat")
+#end
+# replace OUTER_ITER= by 500 and TIME_ITER= 56
+adjoint.addParameter(outer_iter_adjoint)
+adjoint.addExpected("solution_adj_avtp_0_00053.dat")
+adjoint.addExpected("solution_adj_avtp_1_00053.dat")
+
+# gradient projection
+dot = adjont = ExternalRun("DOT",dot_ad_command,True)
+dot.setMaxTries(max_tries)
+dot.addConfig("chtMaster.cfg")
+dot.addData("fluid.cfg") # zonal cfg's can be symlinked
+dot.addData("solid.cfg")
+dot.addData("DEFORM/mesh_out.su2", destination="MeshCHT.su2")
+# add all adjoint solution files
+for timeIter in range(54):
+    if timeIter < 10:
+        timeIter = "0" + str(timeIter)
+    dot.addData("ADJOINT/solution_adj_avtp_0_000" + str(timeIter) + ".dat")
+    dot.addData("ADJOINT/solution_adj_avtp_1_000" + str(timeIter) + ".dat")
+#end
+dot.addExpected("of_grad.csv")
+
+# Functions ------------------------------------------------------------ #
+
+tavgT = Function("tavgT", "DIRECT/chtMaster.csv",LabeledTableReader("\"tavg[AvgTemp[1]]\""))
+tavgT.addInputVariable(ffd,"DOT/of_grad.csv",TableReader(None,0,(1,0))) # all rows, col 0, don't read the header
+tavgT.addValueEvalStep(deform)
+tavgT.addValueEvalStep(direct)
+tavgT.addGradientEvalStep(adjoint)
+tavgT.addGradientEvalStep(dot)
+
+# Driver --------------------------------------------------------------- #
+
+driver = ExteriorPenaltyDriver(0.005)
+driver.addObjective("min", tavgT)
+
+driver.setWorkingDirectory("DOE")
+driver.preprocessVariables()
+driver.setStorageMode(True,"DESIGN_")
+
+his = open("doe.his","w",1)
+driver.setHistorian(his)
+
+# Simulation Runs ------------------------------------------------------ #
+
+# Primal simulation for each deformed DV
+for iLoop in range(0, nDV, 1):
+    x = driver.getInitial()
+    x[iLoop] = 1e-4 # DV_VALUE
+    driver.fun(x)
+#end
+
+# Undeformed/initial primal last in order to have the correct solution in
+# the WorkindDirectory for the following adjoint
+x = driver.getInitial()
+driver.fun(x) # baseline evaluation
+
+# Compute discrete adjoint gradient
+driver.grad(x)
+
+his.close()
+
+# For results run `python postprocess.py` to get screen output
+# of the differences between primal and adjoint simulation.
diff --git a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/postprocess.py b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/postprocess.py
new file mode 100644
index 000000000000..8e59310c1a94
--- /dev/null
+++ b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/postprocess.py
@@ -0,0 +1,26 @@
+# Compute and print absolute difference between Discrete Adjoint
+# and Finite Difference gradient. Prints also percentage difference.
+#
+# Run this script after `python gradient_validation.py` successfully finished
+
+import pandas as pd
+
+# load files
+DAgrad = pd.read_csv("DOE/DOT/of_grad.csv")
+FDvals = pd.read_csv("doe.his")
+
+# additional values
+FDstep = 1e-4
+FDstring = '  tavgT'
+DAstring = 'AVG_TEMPERATURE gradient '
+
+# create FD gradient
+FDgrad = (FDvals[FDstring].iloc[:18] - FDvals[FDstring].iloc[18]) / FDstep
+
+# absolute difference
+absoluteDiff = DAgrad[DAstring] - FDgrad
+print("DAgrad - FDgrad\n", absoluteDiff)
+
+# relative difference in percent
+relDiffPercent = abs(DAgrad[DAstring] - FDgrad)/DAgrad[DAstring] * 100
+print("(DAgrad - FDgrad) / DAgrad * 100\n", relDiffPercent)
diff --git a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/solid.cfg b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/solid.cfg
new file mode 100644
index 000000000000..8af4dae4a7d7
--- /dev/null
+++ b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/solid.cfg
@@ -0,0 +1,59 @@
+% 2021-03-11 TobiKattmann
+%
+SOLVER= HEAT_EQUATION
+%
+OBJECTIVE_FUNCTION= AVG_TEMPERATURE
+%OBJECTIVE_FUNCTION= DRAG
+OBJECTIVE_WEIGHT= 1.0
+%
+% ---------------- (SOLIDS) CONDUCTION CONDITION DEFINITION -------------------%
+%
+% ~Aluminium at room temperature
+INC_NONDIM= DIMENSIONAL
+%
+SOLID_TEMPERATURE_INIT= 340.0
+SOLID_DENSITY= 2719
+SPECIFIC_HEAT_CP = 0.8710
+SOLID_THERMAL_CONDUCTIVITY= 200
+%
+% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
+%
+MARKER_HEATFLUX= ( inner_pin, 1e3 )
+%
+% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
+%
+NUM_METHOD_GRAD= GREEN_GAUSS
+%
+CFL_NUMBER= 1e4
+CFL_ADAPT= NO
+CFL_ADAPT_PARAM= ( 0.1, 2.0, 10.0, 1e10, 0.001 )
+%
+BETA_FACTOR= 50
+MAX_DELTA_TIME= 1.0
+%
+% ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
+%
+LINEAR_SOLVER= FGMRES
+LINEAR_SOLVER_PREC= ILU
+LINEAR_SOLVER_ERROR= 1E-15
+LINEAR_SOLVER_ITER= 20
+%
+% --------------------------- CONVERGENCE PARAMETERS --------------------------%
+%
+CONV_RESIDUAL_MINVAL= -20
+CONV_STARTITER= 10000
+%
+% -------------------- HEAT NUMERICAL METHOD DEFINITION -----------------------%
+%
+CONV_NUM_METHOD_HEAT = SPACE_CENTERED
+TIME_DISCRE_HEAT= EULER_IMPLICIT
+%
+% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
+%
+%MESH_FILENAME= solid.su2
+%
+%HISTORY_OUTPUT= (ITER, RMS_RES, HEAT, LINSOL)
+HISTORY_OUTPUT= (HEAT)
+MARKER_MONITORING= ( inner_pin )
+%MARKER_ANALYZE= ( inner_pin )
+
diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py
index 0759ab1913db..bb4813591a10 100644
--- a/TestCases/parallel_regression_AD.py
+++ b/TestCases/parallel_regression_AD.py
@@ -335,6 +335,19 @@ def main():
     da_sp_pinArray_cht_2d_dp_hf.multizone = True
     test_list.append(da_sp_pinArray_cht_2d_dp_hf)
 
+    # 2D unsteady CHT vortex shedding at RE=200. TAVG_Temperature OF
+    da_unsteadyCHT_cylinder           = TestCase('da_unsteadyCHT_cylinder')
+    da_unsteadyCHT_cylinder.cfg_dir   = "coupled_cht/disc_adj_unsteadyCHT_cylinder"
+    da_unsteadyCHT_cylinder.cfg_file  = "chtMaster.cfg"
+    da_unsteadyCHT_cylinder.test_iter = 2
+    da_unsteadyCHT_cylinder.test_vals = [-3.521358, -4.312658, -4.271025, -9.846075, -7.967741, 0.0000e+00, 3.6840e+00, 2.9483e-01]
+    da_unsteadyCHT_cylinder.su2_exec  = "mpirun -n 2 SU2_CFD_AD"
+    da_unsteadyCHT_cylinder.timeout   = 1600
+    da_unsteadyCHT_cylinder.tol       = 0.00001
+    da_unsteadyCHT_cylinder.unsteady  = True
+    da_unsteadyCHT_cylinder.multizone = True
+    test_list.append(da_unsteadyCHT_cylinder)
+
     ######################################
     ### RUN TESTS                      ###
     ######################################