Add unsteady cht adjoint testcase

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

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).

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 )
+