From 3548413ed7ab7b312b32f3def9894e88095ad43b Mon Sep 17 00:00:00 2001
From: TobiKattmann <Tobias.Kattmann@de.bosch.com>
Date: Thu, 13 May 2021 01:45:11 +0200
Subject: [PATCH 1/4] Add files for unsteady CHT adjoint and primal regression
 test.

---
 .../O_mesh_cylinder.geo                       | 147 +++++++++++++++++
 .../disc_adj_unsteadyCHT_cylinder/README.md   |  28 ++++
 .../chtMaster.cfg                             | 151 ++++++++++++++++++
 .../disc_adj_unsteadyCHT_cylinder/findiff.py  |  78 +++++++++
 .../disc_adj_unsteadyCHT_cylinder/fluid.cfg   |  80 ++++++++++
 .../postprocess.py                            |  28 ++++
 .../disc_adj_unsteadyCHT_cylinder/solid.cfg   |  59 +++++++
 TestCases/parallel_regression_AD.py           |  13 ++
 8 files changed, 584 insertions(+)
 create mode 100644 TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/O_mesh_cylinder.geo
 create mode 100644 TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/README.md
 create mode 100644 TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/chtMaster.cfg
 create mode 100644 TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/findiff.py
 create mode 100644 TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/fluid.cfg
 create mode 100644 TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/postprocess.py
 create mode 100644 TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/solid.cfg

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..35403c1e17cb
--- /dev/null
+++ b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/README.md
@@ -0,0 +1,28 @@
+# 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.
+
+## Finite Differences using FADO
+In order to validate the Discrete Adjoint gradient a Finite Differences python script `findiff.py`
+using [FADO](www.github.com/su2code/FADO) is added.
+The script starts a baseline simulation and an additional for each of the 18 Design Variables.
+Using the `postprocess.py` to print the absolute difference and relative difference in percent to screen.
+Make sure to adapt the `chtMaster.cfg` as for the primal.
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..d4be542e2903
--- /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] )
+% 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/findiff.py b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/findiff.py
new file mode 100644
index 000000000000..5fc711393041
--- /dev/null
+++ b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/findiff.py
@@ -0,0 +1,78 @@
+# FADO script: Finite Differences of unsteady CHT
+
+from FADO import *
+import subprocess
+
+# Design variables ----------------------------------------------------- #
+
+nDV = 18
+ffd = InputVariable(0.0,PreStringHandler("DV_VALUE="),nDV)
+
+# Parameters ----------------------------------------------------------- #
+
+# switch input mesh to perform deformation
+mesh_in = Parameter(["MESH_FILENAME= mesh_out.su2"],\
+       LabelReplacer("MESH_FILENAME= MeshCHT.su2"))
+
+# Evaluations ---------------------------------------------------------- #
+
+def_command = "SU2_DEF chtMaster.cfg" # Note that correct SU2 version needs to be in PATH
+cfd_command = "mpirun -n 12 SU2_CFD 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")
+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(mesh_in)
+
+# Functions ------------------------------------------------------------ #
+
+tavgT = Function("tavgT", "DIRECT/chtMaster.csv",LabeledTableReader("\"tavg[AvgTemp[1]]\""))
+tavgT.addInputVariable(ffd,"",None)
+tavgT.addValueEvalStep(deform)
+tavgT.addValueEvalStep(direct)
+
+# 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)
+
+# Finite Differences
+x = driver.getInitial()
+driver.fun(x) # baseline evaluation
+
+for iLoop in range(0, nDV, 1):
+    x = driver.getInitial()
+    x[iLoop] = 1e-4 # DV_VALUE
+    driver.fun(x)
+#end
+
+his.close()
+
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/postprocess.py b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/postprocess.py
new file mode 100644
index 000000000000..d03bea21c983
--- /dev/null
+++ b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/postprocess.py
@@ -0,0 +1,28 @@
+# Compute and print absolute difference between Discrete Adjoint
+# and Finite Difference gradient. Prints also percentage difference.
+
+import pandas as pd
+
+# load files
+DAgrad = pd.read_csv("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[1:] - FDvals[FDstring].iloc[0]) / FDstep
+# The above operation creates a pd.series that starts at 1 because the first FDval is the baseline.
+# This is reset to match the DAgrad
+FDgrad = FDgrad.reset_index(inplace=False)
+
+# absolute difference
+absoluteDiff = DAgrad[DAstring] - FDgrad[FDstring]
+print("DAgrad - FDgrad\n", absoluteDiff)
+
+# relative difference in percent
+relDiffPercent = abs(DAgrad[DAstring] - FDgrad[FDstring])/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 8688cc2d455f..000b16f0fb90 100644
--- a/TestCases/parallel_regression_AD.py
+++ b/TestCases/parallel_regression_AD.py
@@ -348,6 +348,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]
+    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                      ###
     ######################################

From 9ae02f0b0c4434b2a2c05464bb53e236c86c14cb Mon Sep 17 00:00:00 2001
From: TobiKattmann <31306376+TobiKattmann@users.noreply.github.com>
Date: Thu, 13 May 2021 11:25:50 +0200
Subject: [PATCH 2/4] remove unnecessary python import. LGTM warning

---
 TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/findiff.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/findiff.py b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/findiff.py
index 5fc711393041..c1672bfc9de8 100644
--- a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/findiff.py
+++ b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/findiff.py
@@ -1,7 +1,6 @@
 # FADO script: Finite Differences of unsteady CHT
 
 from FADO import *
-import subprocess
 
 # Design variables ----------------------------------------------------- #
 
@@ -75,4 +74,3 @@
 #end
 
 his.close()
-

From 49f2d50c6390c8141100ed5eee1990f9112ebb11 Mon Sep 17 00:00:00 2001
From: TobiKattmann <Tobias.Kattmann@de.bosch.com>
Date: Fri, 14 May 2021 16:10:59 +0200
Subject: [PATCH 3/4] Addtional reg test screen output. Adapt unst cht adjoint
 FADO script to run adjoint as well.

---
 .../disc_adj_unsteadyCHT_cylinder/README.md   |  14 +-
 .../chtMaster.cfg                             |   2 +-
 .../disc_adj_unsteadyCHT_cylinder/findiff.py  |  76 ----------
 .../gradient_validation.py                    | 139 ++++++++++++++++++
 .../postprocess.py                            |  14 +-
 TestCases/parallel_regression_AD.py           |   2 +-
 6 files changed, 155 insertions(+), 92 deletions(-)
 delete mode 100644 TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/findiff.py
 create mode 100644 TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/gradient_validation.py

diff --git a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/README.md b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/README.md
index 35403c1e17cb..7440d8d4896e 100644
--- a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/README.md
+++ b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/README.md
@@ -13,16 +13,18 @@ A gmsh .geo file is added such that the mesh can be recreated and modified.
 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
+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.
 
-## Finite Differences using FADO
-In order to validate the Discrete Adjoint gradient a Finite Differences python script `findiff.py`
+## 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.
-The script starts a baseline simulation and an additional for each of the 18 Design Variables.
-Using the `postprocess.py` to print the absolute difference and relative difference in percent to screen.
-Make sure to adapt the `chtMaster.cfg` as for the primal.
+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.
+
diff --git a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/chtMaster.cfg b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/chtMaster.cfg
index d4be542e2903..6821416134a5 100644
--- a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/chtMaster.cfg
+++ b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/chtMaster.cfg
@@ -37,7 +37,7 @@ ITER_AVERAGE_OBJ= 54
 %
 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] )
+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
diff --git a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/findiff.py b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/findiff.py
deleted file mode 100644
index c1672bfc9de8..000000000000
--- a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/findiff.py
+++ /dev/null
@@ -1,76 +0,0 @@
-# FADO script: Finite Differences of unsteady CHT
-
-from FADO import *
-
-# Design variables ----------------------------------------------------- #
-
-nDV = 18
-ffd = InputVariable(0.0,PreStringHandler("DV_VALUE="),nDV)
-
-# Parameters ----------------------------------------------------------- #
-
-# switch input mesh to perform deformation
-mesh_in = Parameter(["MESH_FILENAME= mesh_out.su2"],\
-       LabelReplacer("MESH_FILENAME= MeshCHT.su2"))
-
-# Evaluations ---------------------------------------------------------- #
-
-def_command = "SU2_DEF chtMaster.cfg" # Note that correct SU2 version needs to be in PATH
-cfd_command = "mpirun -n 12 SU2_CFD 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")
-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(mesh_in)
-
-# Functions ------------------------------------------------------------ #
-
-tavgT = Function("tavgT", "DIRECT/chtMaster.csv",LabeledTableReader("\"tavg[AvgTemp[1]]\""))
-tavgT.addInputVariable(ffd,"",None)
-tavgT.addValueEvalStep(deform)
-tavgT.addValueEvalStep(direct)
-
-# 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)
-
-# Finite Differences
-x = driver.getInitial()
-driver.fun(x) # baseline evaluation
-
-for iLoop in range(0, nDV, 1):
-    x = driver.getInitial()
-    x[iLoop] = 1e-4 # DV_VALUE
-    driver.fun(x)
-#end
-
-his.close()
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
index d03bea21c983..8e59310c1a94 100644
--- a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/postprocess.py
+++ b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/postprocess.py
@@ -1,10 +1,12 @@
 # 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("of_grad.csv")
+DAgrad = pd.read_csv("DOE/DOT/of_grad.csv")
 FDvals = pd.read_csv("doe.his")
 
 # additional values
@@ -13,16 +15,12 @@
 DAstring = 'AVG_TEMPERATURE gradient '
 
 # create FD gradient
-FDgrad = (FDvals[FDstring].iloc[1:] - FDvals[FDstring].iloc[0]) / FDstep
-# The above operation creates a pd.series that starts at 1 because the first FDval is the baseline.
-# This is reset to match the DAgrad
-FDgrad = FDgrad.reset_index(inplace=False)
+FDgrad = (FDvals[FDstring].iloc[:18] - FDvals[FDstring].iloc[18]) / FDstep
 
 # absolute difference
-absoluteDiff = DAgrad[DAstring] - FDgrad[FDstring]
+absoluteDiff = DAgrad[DAstring] - FDgrad
 print("DAgrad - FDgrad\n", absoluteDiff)
 
 # relative difference in percent
-relDiffPercent = abs(DAgrad[DAstring] - FDgrad[FDstring])/DAgrad[DAstring] * 100
+relDiffPercent = abs(DAgrad[DAstring] - FDgrad)/DAgrad[DAstring] * 100
 print("(DAgrad - FDgrad) / DAgrad * 100\n", relDiffPercent)
-
diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py
index d69d01f9b9dc..95c5e6ae2045 100644
--- a/TestCases/parallel_regression_AD.py
+++ b/TestCases/parallel_regression_AD.py
@@ -340,7 +340,7 @@ def main():
     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]
+    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

From f20947b20cf00c410d06ccdc25fa3558c1ee5f6a Mon Sep 17 00:00:00 2001
From: TobiKattmann <31306376+TobiKattmann@users.noreply.github.com>
Date: Fri, 14 May 2021 16:18:19 +0200
Subject: [PATCH 4/4] Add comment on accuracy

---
 TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/README.md b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/README.md
index 7440d8d4896e..7023e9819430 100644
--- a/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/README.md
+++ b/TestCases/coupled_cht/disc_adj_unsteadyCHT_cylinder/README.md
@@ -27,4 +27,4 @@ Note that the script can be used with the files as they are. Necessary adaptions
 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).