Working thermal radiation adjoint (with test solvers)

Author: bburke38

funtofem/interface/radiation_interface.py

@@ -43,7 +43,7 @@ class RadiationInterface(SolverInterface):
 
     """
 
-    def __init__(self, comm, model, conv_hist=False, complex_mode=False):
+    def __init__(self, comm, model, conv_hist=False, complex_mode=False, debug=False):
         """
         The instantiation of the thermal radiation interface class will populate the model
         with the aerodynamic surface mesh, body.aero_X and body.aero_nnodes.
@@ -60,6 +60,9 @@ def __init__(self, comm, model, conv_hist=False, complex_mode=False):
         self.comm = comm
         self.conv_hist = conv_hist
         self.complex_mode = complex_mode
+        self.debug = debug
+
+        self.sigma_sb = 5.670374419e-8
 
         # setup forward and adjoint tolerances
         super().__init__()
@@ -206,16 +209,15 @@ def iterate(self, scenario: Scenario, bodies: list[Body], step):
                 f.write("{0:03d} ".format(step))
 
         for ibody, body in enumerate(bodies, 1):
-            aero_X = body.get_aero_nodes()
-            aero_id = body.get_aero_node_ids()
             aero_nnodes = body.get_num_aero_nodes()
 
             aero_temps = body.get_aero_temps(scenario, time_index=step)
-            print(f"aero_temps: {aero_temps}")
-            heat_rad = self.calc_heat_flux(aero_temps, scenario)
-
             heat_flux = body.get_aero_heat_flux(scenario, time_index=step)
-            heat_flux += heat_rad
+
+            if aero_temps is not None and aero_nnodes > 0:
+                heat_rad = self.calc_heat_flux(aero_temps, scenario)
+
+                heat_flux += heat_rad
 
         return 0
 
@@ -255,17 +257,16 @@ def iterate_adjoint(self, scenario: Scenario, bodies: list[Body], step):
 
         nfuncs = scenario.count_adjoint_functions()
         for ibody, body in enumerate(bodies, 1):
-            # Get the adjoint-Jacobian product for the heat flux
+            # Get the adjoint-Jacobian product for the aero temperature
             aero_flux_ajp = body.get_aero_heat_flux_ajp(scenario)
             aero_nnodes = body.get_num_aero_nodes()
-            aero_flux = body.get_aero_heat_flux(scenario, time_index=step)
             aero_temps = body.get_aero_temps(scenario, time_index=step)
 
             aero_temps_ajp = body.get_aero_temps_ajp(scenario)
 
-            if aero_flux_ajp is not None and aero_nnodes > 0:
-                # Solve the aero heat flux computation adjoint
-                # dR/dhA^{T} * psi_R = - dQ/dhA^{T} * psi_Q = - aero_flux_ajp
+            if aero_temps_ajp is not None and aero_nnodes > 0:
+                # Add contribution to aero_temps_ajp from radiation
+                # dR/dhR^{T} * psi_R = dA_dhA^{T} * psi_A = - dQ/dhA^{T} * psi_Q = - aero_flux_ajp
                 psi_R = aero_flux_ajp
 
                 rad_heat_deriv = self.calc_heat_flux_deriv(aero_temps, scenario)
@@ -276,13 +277,11 @@ def iterate_adjoint(self, scenario: Scenario, bodies: list[Body], step):
                 for func in range(nfuncs):
                     lam[:, func] = psi_R[:, func] * rad_heat_deriv[:]
 
-                if not self.complex_mode:
-                    lam = lam.astype(np.double)
-
-                for func in range(nfuncs):
-                    print(f"aero_flux_ajp: {aero_flux_ajp[:, func]}")
-                    print(f"lam: {lam[:, func]}")
-                    aero_flux_ajp[:, func] += lam[:, func]
+                for ifunc in range(nfuncs):
+                    if self.debug and self.comm.rank == 0:
+                        print(f"aero_temps_ajp: {aero_temps_ajp[:, func]}")
+                        print(f"lam: {lam[:, func]}")
+                    aero_temps_ajp[:, ifunc] += lam[:, ifunc]
 
         return
 
@@ -313,11 +312,10 @@ def calc_heat_flux(self, temps, scenario=None):
             F_v = scenario.F_v
             T_v = scenario.T_v
 
-        sigma_sb = 5.670374419e-8
         rad_heat = np.zeros_like(temps, dtype=TransferScheme.dtype)
 
         for indx, temp_i in enumerate(temps):
-            rad_heat[indx] = -sigma_sb * emis * F_v * (temp_i**4 - T_v**4)
+            rad_heat[indx] = -self.sigma_sb * emis * F_v * (temp_i**4 - T_v**4)
 
         return rad_heat
 
@@ -331,16 +329,13 @@ def calc_heat_flux_deriv(self, temps, scenario=None):
         if scenario is None:
             emis = 0.8
             F_v = 1.0
-            T_v = 0.0
         else:
             emis = scenario.emis
             F_v = scenario.F_v
-            T_v = scenario.T_v
 
-        sigma_sb = 5.670374419e-8
         rad_heat_deriv = np.zeros_like(temps, dtype=TransferScheme.dtype)
 
         for indx, temp_i in enumerate(temps):
-            rad_heat_deriv[indx] = -4 * sigma_sb * emis * F_v * temp_i**3
+            rad_heat_deriv[indx] = -4 * self.sigma_sb * emis * F_v * temp_i**3
 
         return rad_heat_deriv

funtofem/interface/test_interfaces/_test_struct_solver.py

@@ -92,7 +92,7 @@ def __init__(self, npts, struct_dvs):
                     * (np.random.rand(self.npts, 3 * self.npts) - 0.5*0)
                 )
                 self.c2 = (
-                    0.01
+                    0.5
                     * thermal_scale
                     * (np.random.rand(self.npts, len(self.struct_dvs)) - 0.5*0)
                 )

tests/unit_tests/framework/test_radiationFlux_adjoint.py

@@ -23,9 +23,9 @@ def _setup_model_and_driver(self):
 
         # Create a structural variable
         for i in range(5):
-            thickness = np.random.rand()
+            thickness = np.random.rand() * 200
             svar = Variable(
-                "thickness %d" % (i), value=thickness, lower=0.01, upper=0.1
+                "thickness %d" % (i), value=thickness, lower=0.01, upper=10.0
             )
             plate.add_variable("structural", svar)
 
@@ -44,7 +44,7 @@ def _setup_model_and_driver(self):
         temp = Function("temperature", analysis_type="structural")
         steady.add_function(temp)
         steady.set_temperature(T_ref=200, T_inf=300)
-        steady.set_therm_rad_vals()
+        steady.set_therm_rad_vals(emis=0.8)
 
         # Add the steady-state scenario
         model.add_scenario(steady)
@@ -64,46 +64,9 @@ def _setup_model_and_driver(self):
             solvers, transfer_settings=transfer_settings, model=model
         )
 
-        return model, driver
+        model.print_summary()
 
-    # def test_model_derivatives(self):
-    #     model, driver = self._setup_model_and_driver()
-
-    #     # Check whether to use the complex-step method or now
-    #     complex_step = False
-    #     epsilon = 1e-6
-    #     rtol = 1e-6
-    #     if TransferScheme.dtype == complex:
-    #         complex_step = True
-    #         epsilon = 1e-30
-    #         rtol = 1e-9
-
-    #     # Manual test of the disciplinary solvers
-    #     scenario = model.scenarios[0]
-    #     bodies = model.bodies
-    #     solvers = driver.solvers
-
-    #     fail = solvers.flow.test_adjoint(
-    #         "flow",
-    #         scenario,
-    #         bodies,
-    #         epsilon=epsilon,
-    #         complex_step=complex_step,
-    #         rtol=rtol,
-    #     )
-    #     assert fail == False
-
-    #     solvers.structural.test_adjoint(
-    #         "structural",
-    #         scenario,
-    #         bodies,
-    #         epsilon=epsilon,
-    #         complex_step=complex_step,
-    #         rtol=rtol,
-    #     )
-    #     assert fail == False
-
-    #     return
+        return model, driver
 
     def test_coupled_derivatives(self):
         model, driver = self._setup_model_and_driver()
@@ -119,7 +82,7 @@ def test_coupled_derivatives(self):
 
         # Solve the forward analysis
         driver.solve_forward()
-        driver.solve_adjoint()
+        # driver.solve_adjoint()
 
         # Get the functions
         functions = model.get_functions()
@@ -134,6 +97,8 @@ def test_coupled_derivatives(self):
         driver.solve_adjoint()
         grads = model.get_function_gradients()
 
+        print(f"Gradients: {grads[0][:]}")
+
         # Set the new variable values
         if complex_step:
             variables[0].value = variables[0].value + 1j * epsilon
@@ -156,7 +121,7 @@ def test_coupled_derivatives(self):
             pass_ = False
             if driver.comm.rank == 0:
                 pass_ = abs(rel_error) < rtol
-                print("Approximate gradient  = ", deriv.real)
+                print("Approx. gradient (CS) = ", deriv.real)
                 print("Adjoint gradient      = ", grads[0][0].real)
                 print("Relative error        = ", rel_error.real)
                 print("Pass flag             = ", pass_)
@@ -170,7 +135,7 @@ def test_coupled_derivatives(self):
             pass_ = False
             if driver.comm.rank == 0:
                 pass_ = abs(rel_error) < rtol
-                print("Approximate gradient  = ", deriv)
+                print("Approx. gradient (FD) = ", deriv)
                 print("Adjoint gradient      = ", grads[0][0])
                 print("Relative error        = ", rel_error)
                 print("Pass flag             = ", pass_)