diff --git a/funtofem/driver/funtofem_nlbgs_driver.py b/funtofem/driver/funtofem_nlbgs_driver.py index e08106b6..efd5aea7 100644 --- a/funtofem/driver/funtofem_nlbgs_driver.py +++ b/funtofem/driver/funtofem_nlbgs_driver.py @@ -19,6 +19,7 @@ See the License for the specific language governing permissions and limitations under the License. """ +from __future__ import annotations __all__ = ["FUNtoFEMnlbgs"] @@ -158,6 +159,10 @@ def _solve_steady_forward(self, scenario, steps=None): print("Flow solver returned fail flag") return fail + # Add contribution from thermal radiation + if scenario.thermal_rad: + self.solvers.thermal_rad.iterate(scenario, self.model.bodies, step) + # Transfer the loads and heat flux for body in self.model.bodies: body.transfer_loads(scenario) @@ -284,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) diff --git a/funtofem/interface/__init__.py b/funtofem/interface/__init__.py index 3539243f..e65514a9 100644 --- a/funtofem/interface/__init__.py +++ b/funtofem/interface/__init__.py @@ -39,6 +39,9 @@ # caps2fun utils from .caps2fun import * +# Radiation interface +from .radiation_interface import * + # FUN3D interface, with python package "fun3d" fun3d_loader = importlib.util.find_spec("fun3d") if fun3d_loader is not None: diff --git a/funtofem/interface/radiation_interface.py b/funtofem/interface/radiation_interface.py new file mode 100644 index 00000000..9807a601 --- /dev/null +++ b/funtofem/interface/radiation_interface.py @@ -0,0 +1,341 @@ +#!/usr/bin/env python +""" +This file is part of the package FUNtoFEM for coupled aeroelastic simulation +and design optimization. + +Copyright (C) 2015 Georgia Tech Research Corporation. +Additional copyright (C) 2015 Kevin Jacobson, Jan Kiviaho and Graeme Kennedy. +All rights reserved. + +FUNtoFEM is licensed under the Apache License, Version 2.0 (the "License"); +you may not use this software except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +""" + +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): + """ + FUNtoFEM interface class for radiative heating. Works for steady analysis only. + Unlike other interfaces in FUNtoFEM, this interface also implements the physical model; + no external radiative heating discipline solver is used here. + + """ + + 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. + + Parameters + ---------- + comm: MPI.comm + MPI communicator + model : :class:`FUNtoFEMmodel` + FUNtoFEM model + conv_hist : bool + output convergence history to file + """ + 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__() + + # Store previous iteration's aerodynamic forces and displacements for + # each body in order to compute RMS error at each iteration + if self.conv_hist: + self.temps_prev = {} + self.heat_prev = {} + self.conv_hist_file = "conv_hist.dat" + + # 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. + + Parameters + ---------- + scenario: :class:`~scenario.Scenario` + The scenario that needs to be initialized + bodies: :class:`~body.Body` + list of FUNtoFEM bodies to either get new surface meshes from or to + set the original mesh in + + Returns + ------- + fail: int + Returns zero for successful completion of initialization + + """ + + # Touch a file to record the RMS error output for convergence study + if self.conv_hist: + with open(self.conv_hist_file, "w") as f: + pass + + # Extract the relevant components' node locations to body object + 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() + + # Initialize the state values used for convergence study as well + if self.conv_hist: + self.temps_prev[body.id] = np.zeros( + 3 * aero_nnodes, dtype=TransferScheme.dtype + ) + self.heat_prev[body.id] = np.zeros( + 3 * aero_nnodes, dtype=TransferScheme.dtype + ) + + return 0 + + def iterate(self, scenario: Scenario, bodies: list[Body], step): + """ + Forward iteration of thermal radiation. + + Add heat flux contribution from thermal radiation to aero heat flux. + + Parameters + ---------- + scenario: :class:`~scenario.Scenario` + The scenario + bodies: :class:`~body.Body` + list of FUNtoFEM bodies. + step: int + the time step number + + """ + + # Write step number to file output + if self.conv_hist: + with open(self.conv_hist_file, "a") as f: + f.write("{0:03d} ".format(step)) + + for ibody, body in enumerate(bodies, 1): + aero_nnodes = body.get_num_aero_nodes() + + aero_temps = body.get_aero_temps(scenario, time_index=step) + heat_flux = body.get_aero_heat_flux(scenario, time_index=step) + + 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 + + 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 aero temperature + aero_flux_ajp = body.get_aero_heat_flux_ajp(scenario) + aero_nnodes = body.get_num_aero_nodes() + aero_temps = body.get_aero_temps(scenario, time_index=step) + + aero_temps_ajp = body.get_aero_temps_ajp(scenario) + + 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) + + dtype = TransferScheme.dtype + lam = np.zeros((aero_nnodes, nfuncs), dtype=dtype) + + for func in range(nfuncs): + lam[:, func] = psi_R[:, func] * rad_heat_deriv[:] + + 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 + + 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. + + Parameters + --------- + temps: + Array of temperatures. + scenario: + Used to set the model parameters. + """ + 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 + + rad_heat = np.zeros_like(temps, dtype=TransferScheme.dtype) + + for indx, temp_i in enumerate(temps): + rad_heat[indx] = -self.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 + else: + emis = scenario.emis + F_v = scenario.F_v + + rad_heat_deriv = np.zeros_like(temps, dtype=TransferScheme.dtype) + + for indx, temp_i in enumerate(temps): + rad_heat_deriv[indx] = -4 * self.sigma_sb * emis * F_v * temp_i**3 + + return rad_heat_deriv diff --git a/funtofem/interface/solver_manager.py b/funtofem/interface/solver_manager.py index 9f579a4a..e327eef3 100644 --- a/funtofem/interface/solver_manager.py +++ b/funtofem/interface/solver_manager.py @@ -13,6 +13,8 @@ if fun3d_loader is not None: from .fun3d_14_interface import Fun3d14Interface +from .radiation_interface import RadiationInterface + class CommManager: def __init__( @@ -55,7 +57,13 @@ def __str__(self): class SolverManager: - def __init__(self, comm, use_flow: bool = True, use_struct: bool = True): + def __init__( + self, + comm, + use_flow: bool = True, + use_struct: bool = True, + use_thermal_rad: bool = False, + ): """ Create a solver manager object which holds flow, struct solvers and in the future might be expanded to hold dynamics, etc. @@ -68,13 +76,17 @@ def __init__(self, comm, use_flow: bool = True, use_struct: bool = True): whether to require flow solvers like Fun3dInterface use_struct: bool whether to require structural solvers like TacsInterface + use_thermal_rad: bool + whether to require thermal radiation solvers like RadiationInterface """ self.comm = comm self._use_flow = use_flow self._use_struct = use_struct + self._use_thermal_rad = use_thermal_rad self._flow = None self._structural = None + self._thermal_rad = None @property def use_flow(self) -> bool: @@ -84,6 +96,10 @@ def use_flow(self) -> bool: def use_struct(self) -> bool: return self._use_struct + @property + def use_thermal_rad(self) -> bool: + return self._use_thermal_rad + @property def uses_fun3d(self) -> bool: if fun3d_loader is None or self.flow is None: @@ -91,6 +107,13 @@ def uses_fun3d(self) -> bool: else: return isinstance(self.flow, Fun3d14Interface) + @property + def uses_radiation(self) -> bool: + if self.thermal_rad is None: + return False + else: + return isinstance(self.thermal_rad, RadiationInterface) + @property def solver_list(self): """ @@ -101,6 +124,8 @@ def solver_list(self): mlist.append(self.flow) if self.use_struct: mlist.append(self.structural) + if self.use_thermal_rad: + mlist.append(self.thermal_rad) return mlist @property @@ -142,6 +167,14 @@ def structural(self): def structural(self, new_structural_solver): self._structural = new_structural_solver + @property + def thermal_rad(self): + return self._thermal_rad + + @thermal_rad.setter + def thermal_rad(self, thermal_rad_solver): + self._thermal_rad = thermal_rad_solver + @property def comm_manager(self) -> CommManager: """ diff --git a/funtofem/interface/test_interfaces/_test_struct_solver.py b/funtofem/interface/test_interfaces/_test_struct_solver.py index c4fd041f..600028a9 100644 --- a/funtofem/interface/test_interfaces/_test_struct_solver.py +++ b/funtofem/interface/test_interfaces/_test_struct_solver.py @@ -7,7 +7,15 @@ class TestStructuralSolver(SolverInterface): - def __init__(self, comm, model, elastic_k=1.0, thermal_k=1.0, default_mesh=True): + def __init__( + self, + comm, + model, + elastic_k=1.0, + thermal_k=1.0, + default_mesh=True, + hot_temps=False, + ): """ A test solver that provides the functionality that FUNtoFEM expects from a structural solver. @@ -32,6 +40,8 @@ def __init__(self, comm, model, elastic_k=1.0, thermal_k=1.0, default_mesh=True) self.npts = 25 np.random.seed(54321) + self.hot_temps = hot_temps + # setup forward and adjoint tolerances super().__init__() @@ -58,9 +68,10 @@ def __init__(self, comm, model, elastic_k=1.0, thermal_k=1.0, default_mesh=True) # scenario data for the multi scenario case class ScenarioData: - def __init__(self, npts, struct_dvs): + def __init__(self, npts, struct_dvs, hot_temps=False): self.npts = npts self.struct_dvs = struct_dvs + self.hot_temps = hot_temps # choose time step self.dt = 0.01 @@ -82,20 +93,37 @@ def __init__(self, npts, struct_dvs): * (np.random.rand(3 * self.npts, len(self.struct_dvs)) - 0.5) ) - # Struct temps = Jac2 * struct_flux + b2 * struct_X + c2 * struct_dvs + omega2 * time - self.Jac2 = ( - 0.05 * thermal_scale * (np.random.rand(self.npts, self.npts) - 0.5) - ) - self.b2 = ( - 0.1 - * thermal_scale - * (np.random.rand(self.npts, 3 * self.npts) - 0.5) - ) - self.c2 = ( - 0.01 - * thermal_scale - * (np.random.rand(self.npts, len(self.struct_dvs)) - 0.5) - ) + if not self.hot_temps: + # Struct temps = Jac2 * struct_flux + b2 * struct_X + c2 * struct_dvs + omega2 * time + self.Jac2 = ( + 0.05 + * thermal_scale + * (np.random.rand(self.npts, self.npts) - 0.5) + ) + self.b2 = ( + 0.1 + * thermal_scale + * (np.random.rand(self.npts, 3 * self.npts) - 0.5) + ) + self.c2 = ( + 0.01 + * thermal_scale + * (np.random.rand(self.npts, len(self.struct_dvs)) - 0.5) + ) + elif self.hot_temps: + self.Jac2 = ( + 0.05 + * thermal_scale + * (np.random.rand(self.npts, self.npts) + 0.5) + ) + self.b2 = ( + 0.1 * thermal_scale * (np.random.rand(self.npts, 3 * self.npts)) + ) + self.c2 = ( + 0.5 + * thermal_scale + * (np.random.rand(self.npts, len(self.struct_dvs)) - 0.5 * 0) + ) # Data for output functional values self.func_coefs1 = np.random.rand(3 * self.npts) @@ -109,7 +137,9 @@ def __init__(self, npts, struct_dvs): # create scenario data for each scenario self.scenario_data = {} for scenario in model.scenarios: - self.scenario_data[scenario.id] = ScenarioData(self.npts, self.struct_dvs) + self.scenario_data[scenario.id] = ScenarioData( + self.npts, self.struct_dvs, self.hot_temps + ) if default_mesh: # Set random initial node locations diff --git a/funtofem/model/scenario.py b/funtofem/model/scenario.py index a19de897..87ad46a2 100644 --- a/funtofem/model/scenario.py +++ b/funtofem/model/scenario.py @@ -68,6 +68,10 @@ def __init__( gamma=1.4, R_specific=287.058, Pr=0.72, + emis=0.8, + F_v=1.0, + T_v=0.0, + thermal_rad=False, ): """ Parameters @@ -133,6 +137,14 @@ def __init__( Specific gas constant of the working fluid (assumed air). Units of J/kg-K Pr: double Prandtl number. + emis: float + Surface emissivity (for thermal radiation model). + F_v: float + View factor (for thermal radiation model). + T_v: float + Temperature of background environment (vacuum) (for thermal radiation model). + thermal_rad: bool + Whether to use thermal radiation model in this scenario. See Also -------- :mod:`base` : Scenario inherits from Base @@ -187,6 +199,12 @@ def __init__( cp = self.R_specific * self.gamma / (self.gamma - 1) self.cp = cp + # Thermal radiation parameters + self.emis = emis + self.F_v = F_v + self.T_v = T_v + self.thermal_rad = thermal_rad + if fun3d: mach = Variable("Mach", id=1, upper=5.0, active=False) aoa = Variable("AOA", id=2, lower=-15.0, upper=15.0, active=False) @@ -456,6 +474,33 @@ def set_flow_ref_vals(self, qinf: float = 1.0, flow_dt: float = 1.0): raise ValueError("For steady cases, flow_dt must be set to 1.") return self + def set_therm_rad_vals( + self, + emis: float = 0.8, + F_v: float = 1.0, + T_v: float = 0.0, + thermal_rad: bool = True, + ): + """ + Set parameters for thermal radiation model. + + Parameters + ---------- + emis: float + Surface emissivity. + F_v: float + View factor. + T_v: float + Temperature of background environment (vacuum). + thermal_rad: bool + Whether to use thermal radiation model in this scenario. + """ + self.emis = emis + self.F_v = F_v + self.T_v = T_v + self.thermal_rad = thermal_rad + return self + def set_id(self, id): """ **[model call]** diff --git a/tests/unit_tests/framework/test_radiationFlux_adjoint.py b/tests/unit_tests/framework/test_radiationFlux_adjoint.py new file mode 100644 index 00000000..667d9446 --- /dev/null +++ b/tests/unit_tests/framework/test_radiationFlux_adjoint.py @@ -0,0 +1,150 @@ +import numpy as np +from mpi4py import MPI +from funtofem import TransferScheme + +from funtofem.model import FUNtoFEMmodel, Variable, Scenario, Body, Function +from funtofem.interface import ( + TestAerodynamicSolver, + TestStructuralSolver, + RadiationInterface, + SolverManager, +) +from funtofem.driver import FUNtoFEMnlbgs, TransferSettings + +import unittest + + +class CoupledRadiationTest(unittest.TestCase): + + def _setup_model_and_driver(self): + # Build the model + model = FUNtoFEMmodel("model") + plate = Body("plate", "aerothermal", group=0, boundary=1) + + # Create a structural variable + for i in range(5): + thickness = np.random.rand() * 200 + svar = Variable( + "thickness %d" % (i), value=thickness, lower=0.01, upper=10.0 + ) + plate.add_variable("structural", svar) + + model.add_body(plate) + + # Create a scenario to run + steady = Scenario("steady", group=0, steps=100) + + # Add the aerodynamic variables to the scenario + for i in range(4): + value = np.random.rand() + avar = Variable("aero var %d" % (i), value=value, lower=-10.0, upper=10.0) + steady.add_variable("aerodynamic", avar) + + # Add a function to the scenario + temp = Function("temperature", analysis_type="structural") + steady.add_function(temp) + steady.set_temperature(T_ref=200, T_inf=300) + steady.set_therm_rad_vals(emis=0.8) + + # Add the steady-state scenario + model.add_scenario(steady) + + # Instantiate a test solver for the flow and structures + comm = MPI.COMM_WORLD + solvers = SolverManager(comm) + solvers.flow = TestAerodynamicSolver(comm, model) + solvers.structural = TestStructuralSolver(comm, model, hot_temps=True) + solvers.thermal_rad = RadiationInterface(comm, model) + + # L&D transfer options + transfer_settings = TransferSettings(npts=5) + + # instantiate the driver + driver = FUNtoFEMnlbgs( + solvers, transfer_settings=transfer_settings, model=model + ) + + model.print_summary() + + return model, driver + + def test_coupled_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-5 + if TransferScheme.dtype == complex: + complex_step = True + epsilon = 1e-30 + rtol = 1e-9 + + # Solve the forward analysis + driver.solve_forward() + # driver.solve_adjoint() + + # Get the functions + functions = model.get_functions() + variables = model.get_variables() + + # Store the function values + fvals_init = [] + for func in functions: + fvals_init.append(func.value) + + # Solve the adjoint and get the function gradients + 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 + model.set_variables(variables) + else: + variables[0].value = variables[0].value + epsilon + model.set_variables(variables) + + driver.solve_forward() + + # Store the function values + fvals = [] + for func in functions: + fvals.append(func.value) + + if complex_step: + deriv = fvals[0].imag / epsilon + + rel_error = (deriv - grads[0][0]) / deriv + pass_ = False + if driver.comm.rank == 0: + pass_ = abs(rel_error) < rtol + print("Approx. gradient (CS) = ", deriv.real) + print("Adjoint gradient = ", grads[0][0].real) + print("Relative error = ", rel_error.real) + print("Pass flag = ", pass_) + + pass_ = driver.comm.bcast(pass_, root=0) + assert pass_ + else: + deriv = (fvals[0] - fvals_init[0]) / epsilon + + rel_error = (deriv - grads[0][0]) / deriv + pass_ = False + if driver.comm.rank == 0: + pass_ = abs(rel_error) < rtol + print("Approx. gradient (FD) = ", deriv) + print("Adjoint gradient = ", grads[0][0]) + print("Relative error = ", rel_error) + print("Pass flag = ", pass_) + + pass_ = driver.comm.bcast(pass_, root=0) + assert pass_ + + return + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/unit_tests/framework/test_radiationFlux_deriv.py b/tests/unit_tests/framework/test_radiationFlux_deriv.py new file mode 100644 index 00000000..e5687e04 --- /dev/null +++ b/tests/unit_tests/framework/test_radiationFlux_deriv.py @@ -0,0 +1,62 @@ +import numpy as np +from mpi4py import MPI +from funtofem import TransferScheme +from funtofem.model import Scenario, FUNtoFEMmodel +from funtofem.interface.radiation_interface import RadiationInterface +import unittest + +np.random.seed(343) + +nprocs = 1 +comm = MPI.COMM_WORLD + + +class RadiationFluxTest(unittest.TestCase): + def test_radiative_flux_deriv(self): + """ + Perform a finite difference check to test the implementation of the thermal radiation + heat flux derivative. + + hR = q_rad(tA) + dhR/dtA = dq_rad/dtA + + fd = (q_rad(tA+h*p)-q_rad(tA))/h + """ + + myType = TransferScheme.dtype + model = FUNtoFEMmodel("radiation") + steady = Scenario("steady", group=0, steps=1) + steady.register_to(model) + + aero_temps = np.array(np.random.rand(100) * 400, dtype=myType) + rtol = 1e-6 + + p = np.ones(aero_temps.shape, dtype=TransferScheme.dtype) + h = 1e-5 + + temps_pert = aero_temps + h * p + + thermal_rad_solver = RadiationInterface(comm=comm, model=model) + + q_rad0 = thermal_rad_solver.calc_heat_flux(aero_temps) + + q_rad1 = thermal_rad_solver.calc_heat_flux(temps_pert) + + fd = (q_rad1 - q_rad0) / h + + dq_dtA = thermal_rad_solver.calc_heat_flux_deriv(aero_temps) + + fd_scalar = np.dot(fd, p) + dq_dtA_scalar = np.dot(dq_dtA, p) + + rel_err = (fd_scalar - dq_dtA_scalar) / dq_dtA_scalar + + print( + f"Finite difference check for thermal radiation heat flux derivative: {rel_err}", + flush=True, + ) + assert abs(rel_err) < rtol + + +if __name__ == "__main__": + unittest.main()