Add adjoint for thermal radiation

Author: bburke38

funtofem/driver/funtofem_nlbgs_driver.py

@@ -289,6 +289,12 @@ def _solve_steady_adjoint(self, scenario):
                         print("Flow solver returned fail flag")
                     return fail
 
+                # Add contribution from thermal radiation
+                if scenario.thermal_rad:
+                    self.solvers.thermal_rad.iterate_adjoint(
+                        scenario, self.model.bodies, step
+                    )
+
                 # Get the structural adjoint rhs
                 for body in self.model.bodies:
                     body.transfer_disps_adjoint(scenario)

funtofem/interface/radiation_interface.py

@@ -20,15 +20,20 @@
 limitations under the License.
 """
 
-from __future__ import print_function
+from __future__ import print_function, annotations
 
 __all__ = ["RadiationInterface"]
 
 import numpy as np
 import os
+from typing import TYPE_CHECKING
 from funtofem import TransferScheme
 from ._solver_interface import SolverInterface
 
+if TYPE_CHECKING:
+    from ..model.body import Body
+    from ..model.scenario import Scenario
+
 
 class RadiationInterface(SolverInterface):
     """
@@ -38,7 +43,7 @@ class RadiationInterface(SolverInterface):
 
     """
 
-    def __init__(self, comm, model, conv_hist=False):
+    def __init__(self, comm, model, conv_hist=False, complex_mode=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.
@@ -54,6 +59,7 @@ def __init__(self, comm, model, conv_hist=False):
         """
         self.comm = comm
         self.conv_hist = conv_hist
+        self.complex_mode = complex_mode
 
         # setup forward and adjoint tolerances
         super().__init__()
@@ -68,6 +74,74 @@ def __init__(self, comm, model, conv_hist=False):
         # Get the initial aerodynamic surface meshes
         self.initialize(model.scenarios[0], model.bodies)
 
+    def set_variables(self, scenario, bodies):
+        """
+        Set the design variables into the solver.
+        The scenario and bodies objects have dictionaries of :class:`~variable.Variable` objects.
+        The interface class should pick out which type of variables it needs and pass them into the solver
+
+        Parameters
+        ----------
+        scenario: :class:`~scenario.Scenario`
+            The current scenario
+        bodies: list of :class:`~body.Body` objects
+            The bodies in the model
+        """
+        pass
+
+    def set_functions(self, scenario, bodies):
+        """
+        Set the function definitions into the solver.
+        The scenario has a list of function objects.
+
+        Parameters
+        ----------
+        scenario: :class:`~scenario.Scenario`
+            The current scenario
+        bodies: list of :class:`~body.Body` objects
+            The bodies in the model
+        """
+        pass
+
+    def get_functions(self, scenario, bodies):
+        """
+        Put the function values from the solver in the value attribute of the scneario's functions.
+        The scenario has the list of function objects where the functions owned by this solver will be set.
+        You can evaluate the functions based on the name or based on the functions set during
+        :func:`~solver_interface.SolverInterface.set_functions`.
+        The solver is only responsible for returning the values of functions it owns.
+
+        Parameters
+        ----------
+        scenario: :class:`~scenario.Scenario`
+            The current scenario
+        bodies: list of :class:`~body.Body` objects
+            The bodies in the model
+        """
+        pass
+
+    def get_function_gradients(self, scenario, bodies):
+        """
+        Get the derivatives of all the functions with respect to design variables associated with this solver.
+
+        Each solver sets the function gradients for its own variables into the function objects using either
+        ``function.set_gradient(var, value)`` or ``function.add_gradient(var, vaule)``. Note that before
+        this function is called, all gradient components are zeroed.
+
+        The derivatives are stored in a dictionary in each function class. As a result, the gradients are
+        stored in an unordered format. The gradients returned by ``model.get_function_gradients()`` are
+        flattened into a list of lists whose order is determined by the variable list stored in the model
+        class.
+
+        Parameters
+        ----------
+        scenario: :class:`~scenario.Scenario`
+            The current scenario
+        bodies: list of :class:`~body.Body` objects
+            The bodies in the model
+        """
+        pass
+
     def initialize(self, scenario, bodies):
         """
         Initialize the thermal radiation interface and solver.
@@ -109,26 +183,7 @@ def initialize(self, scenario, bodies):
 
         return 0
 
-    def get_functions(self, scenario, bodies):
-        """
-        Populate the scenario with the aerodynamic function values.
-
-        Parameters
-        ----------
-        scenario: :class:`~scenario.Scenario`
-            The scenario
-        bodies: :class:`~body.Body`
-            list of FUNtoFEM bodies
-        """
-        pass
-        # for function in scenario.functions:
-        #    if function.analysis_type=='aerodynamic':
-        #        # the [6] index returns the value
-        #        if self.comm.Get_rank() == 0:
-        #            function.value = interface.design_pull_composite_func(function.id)[6]
-        #        function.value = self.comm.bcast(function.value,root=0)
-
-    def iterate(self, scenario, bodies, step):
+    def iterate(self, scenario: Scenario, bodies: list[Body], step):
         """
         Forward iteration of thermal radiation.
 
@@ -156,19 +211,93 @@ def iterate(self, scenario, bodies, step):
             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
 
         return 0
 
+    def post(self, scenario, bodies):
+        """
+        Perform any tasks the solver needs to do after the forward steps are complete, e.g., evaluate functions,
+        post-process, deallocate unneeded memory.
+
+        Parameters
+        ----------
+        scenario: :class:`~scenario.Scenario`
+            The current scenario
+        bodies: list of :class:`~body.Body` objects
+            The bodies in the model
+        """
+        pass
+
+    def initialize_adjoint(self, scenario, bodies):
+        """
+        Perform any tasks the solver needs to do before taking adjoint steps
+
+        Parameters
+        ----------
+        scenario: :class:`~scenario.Scenario`
+            The current scenario
+        bodies: list of :class:`~body.Body` objects
+            The bodies in the model
+        """
+        return 0
+
+    def iterate_adjoint(self, scenario: Scenario, bodies: list[Body], step):
+        """
+        Adjoint iteration of thermal radiation.
+
+        Add contribution to aerodynamic temperature adjoint.
+        """
+
+        nfuncs = scenario.count_adjoint_functions()
+        for ibody, body in enumerate(bodies, 1):
+            # Get the adjoint-Jacobian product for the heat flux
+            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
+                psi_R = aero_flux_ajp
+
+                rad_heat_deriv = self.calc_heat_flux_deriv(aero_temps, scenario)
+
+                dtype = TransferScheme.dtype
+                lam = np.zeros((aero_nnodes, nfuncs), dtype=dtype)
+
+                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]
+
+        return
+
+    def post_adjoint(self, scenario, bodies):
+        """
+        Any actions that need to be performed after completing the adjoint solve, e.g., evaluating gradients, deallocating memory, etc.
+        """
+        pass
+
     def calc_heat_flux(self, temps, scenario=None):
         """
         Implementation of thermal radiation from a surface to a low-temperature environment.
         Calculate the heat flux per area as a function of temperature.
 
-        Paramters
+        Parameters
         ---------
         temps:
             Array of temperatures.
@@ -191,3 +320,27 @@ def calc_heat_flux(self, temps, scenario=None):
             rad_heat[indx] = -sigma_sb * emis * F_v * (temp_i**4 - T_v**4)
 
         return rad_heat
+
+    def calc_heat_flux_deriv(self, temps, scenario=None):
+        """
+        Calculate the derivative of radiative heat flux residual dR/dtA
+
+        This forms a diagonal matrix, since the heat flux depends only on the temperature
+        at its corresponding node.
+        """
+        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
+
+        return rad_heat_deriv