FUNtoFEM read and write unsteady aero loads files (#289)

* working unsteady aero loads files

* black reformat

Author: sean-engelstad

funtofem/driver/oneway_struct_driver.py

@@ -257,6 +257,86 @@ def prime_loads(
             timing_file=timing_file,
         )
 
+    @classmethod
+    def prime_loads_from_unsteady_files(
+        cls,
+        files: list,
+        solvers,
+        model,
+        nprocs,
+        transfer_settings,
+        external_shape=False,
+        init_transfer=False,
+        timing_file=None,
+    ):
+        """
+        Used to prime aero loads for optimization over tacs analysis with shape change and tacs aim
+        Built from an aero loads file of a previous CFD analysis in FuntofemNlbgs driver (TODO : make uncoupled fluid solver features)
+        The loads file is written from the FUNtoFEM model after a forward analysis of a flow solver
+
+        Parameters
+        ----------
+        filename : str or path to aero loads file
+            the filepath of the aerodynamic loads file from a previous CFD analysis (written from FUNtoFEM model)
+        solvers: :class:`~interface.SolverManager`
+            Solver Interface for CFD such as Fun3dInterface, Su2Interface, or test aero solver
+        model : :class:`~model.FUNtoFEMmodel
+        nprocs: int
+            Number of processes that TACS is running on.
+        transfer_settings: :class:`~transfer_settings.TransferSettings`
+            Settings for transfer of state variables across aero and struct meshes
+        struct_aim: `caps2tacs.TacsAim`
+            Interface object from TACS to ESP/CAPS, wraps the tacsAIM object.
+        external_shape: bool
+            whether the tacs aim shape analysis is performed outside this class
+        timing_file: str or path
+            path to funtofem timing file statistics
+        """
+        comm = solvers.comm
+        world_rank = comm.Get_rank()
+        if world_rank < nprocs:
+            color = 1
+        else:
+            color = MPI.UNDEFINED
+        tacs_comm = comm.Split(color, world_rank)
+
+        # initialize transfer settings
+        comm_manager = solvers.comm_manager
+
+        for itime, file in enumerate(files):
+            # read in the loads from the file for each time step
+            loads_data = model._read_aero_loads(comm, file)
+
+            # initialize the transfer scheme then distribute aero loads
+            for body in model.bodies:
+                if (
+                    itime == 0
+                ):  # only initialize transfer and scenario data on the first time step after the aero mesh is loaded
+                    body.initialize_transfer(
+                        comm=comm,
+                        struct_comm=tacs_comm,
+                        struct_root=comm_manager.struct_root,
+                        aero_comm=comm_manager.aero_comm,
+                        aero_root=comm_manager.aero_root,
+                        transfer_settings=transfer_settings,
+                    )
+                    for scenario in model.scenarios:
+                        assert not scenario.steady
+                        body.initialize_variables(scenario)
+
+                body._distribute_aero_loads(loads_data, steady=False, itime=itime)
+
+        tacs_driver = cls(
+            solvers,
+            model,
+            nprocs=nprocs,
+            external_shape=external_shape,
+            timing_file=timing_file,
+        )
+        if init_transfer:
+            tacs_driver._transfer_fixed_aero_loads()
+        return tacs_driver
+
     @classmethod
     def prime_loads_from_file(
         cls,
@@ -304,7 +384,7 @@ def prime_loads_from_file(
         comm_manager = solvers.comm_manager
 
         # read in the loads from the file
-        loads_data = model.read_aero_loads(comm, filename)
+        loads_data = model._read_aero_loads(comm, filename)
 
         # initialize the transfer scheme then distribute aero loads
         for body in model.bodies:
@@ -318,7 +398,8 @@ def prime_loads_from_file(
             )
             for scenario in model.scenarios:
                 body.initialize_variables(scenario)
-            body._distribute_aero_loads(loads_data)
+                assert scenario.steady
+            body._distribute_aero_loads(loads_data, steady=True)
 
         tacs_driver = cls(
             solvers,
@@ -418,16 +499,36 @@ def _transfer_fixed_aero_loads(self):
                 )
 
                 # initialize new elastic struct vectors
-                if body.transfer is not None:
-                    body.struct_loads[scenario.id] = np.zeros(3 * ns, dtype=dtype)
-                    body.struct_disps[scenario.id] = np.zeros(3 * ns, dtype=dtype)
-
-                # initialize new struct heat flux
-                if body.thermal_transfer is not None:
-                    body.struct_heat_flux[scenario.id] = np.zeros(ns, dtype=dtype)
-                    body.struct_temps[scenario.id] = (
-                        np.ones(ns, dtype=dtype) * scenario.T_ref
-                    )
+                if scenario.steady:
+                    if body.transfer is not None:
+                        body.struct_loads[scenario.id] = np.zeros(3 * ns, dtype=dtype)
+                        body.struct_disps[scenario.id] = np.zeros(3 * ns, dtype=dtype)
+
+                    # initialize new struct heat flux
+                    if body.thermal_transfer is not None:
+                        body.struct_heat_flux[scenario.id] = np.zeros(ns, dtype=dtype)
+                        body.struct_temps[scenario.id] = (
+                            np.ones(ns, dtype=dtype) * scenario.T_ref
+                        )
+
+                else:  # unsteady
+                    if body.transfer is not None:
+                        body.struct_loads[scenario.id] = [
+                            np.zeros(3 * ns, dtype=dtype) for _ in range(scenario.steps)
+                        ]
+                        body.struct_disps[scenario.id] = [
+                            np.zeros(3 * ns, dtype=dtype) for _ in range(scenario.steps)
+                        ]
+
+                    # initialize new struct heat flux
+                    if body.thermal_transfer is not None:
+                        body.struct_heat_flux[scenario.id] = [
+                            np.zeros(ns, dtype=dtype) for _ in range(scenario.steps)
+                        ]
+                        body.struct_temps[scenario.id] = [
+                            (np.ones(ns, dtype=dtype) * scenario.T_ref)
+                            for _ in range(scenario.steps)
+                        ]
 
                 # transfer disps to prevent seg fault if coming from OnewayAeroDriver
                 body.transfer_disps(scenario)

funtofem/model/body.py

@@ -1594,11 +1594,11 @@ def aitken_adjoint_relax(self, comm, scenario, tol=1e-16):
 
         return
 
-    def _distribute_aero_loads(self, data):
+    def _distribute_aero_loads(self, data, steady: bool = True, itime: int = 0):
         """
         distribute the aero loads and heat flux from a loads file
         """
-        print(f"F2F - starting to distribute loads")
+        print(f"F2F - starting to distribute loads Time Step {itime}")
 
         for scenario_id in data:
             scenario_data = data[scenario_id]
@@ -1613,16 +1613,31 @@ def _distribute_aero_loads(self, data):
                     }
 
             for ind, aero_id in enumerate(self.aero_id):
-                if self.transfer is not None:
-                    self.aero_loads[scenario_id][3 * ind : 3 * ind + 3] = (
-                        scenario_entry_dict[aero_id]["load"]
-                    )
-                if self.thermal_transfer is not None:
-                    self.aero_heat_flux[scenario_id][ind] = scenario_entry_dict[
-                        aero_id
-                    ]["hflux"]
+                if steady:
+                    if self.transfer is not None:
+                        self.aero_loads[scenario_id][3 * ind : 3 * ind + 3] = (
+                            scenario_entry_dict[aero_id]["load"]
+                        )
+                    if self.thermal_transfer is not None:
+                        self.aero_heat_flux[scenario_id][ind] = scenario_entry_dict[
+                            aero_id
+                        ]["hflux"]
+
+                else:  # unsteady
+                    if self.transfer is not None:
+                        # make sure this time index exists in the case of scenarios with different # of time steps
+                        if itime < len(self.aero_loads[scenario_id]):
+                            self.aero_loads[scenario_id][itime][
+                                3 * ind : 3 * ind + 3
+                            ] = scenario_entry_dict[aero_id]["load"]
+                    if self.thermal_transfer is not None:
+                        # make sure this time index exists in the case of scenarios with different # of time steps
+                        if itime < len(self.aero_heat_flux[scenario_id]):
+                            self.aero_heat_flux[scenario_id][itime][ind] = (
+                                scenario_entry_dict[aero_id]["hflux"]
+                            )
 
-        print(f"F2F - done distribute loads")
+        print(f"\tF2F - done distribute loads Time step {itime}")
 
     def _collect_aero_mesh(self, comm, root=0):
         """
@@ -1658,18 +1673,30 @@ def _collect_aero_mesh(self, comm, root=0):
 
         return aero_ids, aero_X
 
-    def _collect_aero_loads(self, comm, scenario, root=0):
+    def _collect_aero_loads(self, comm, scenario, itime: int = 0, root=0):
         """
         gather the aerodynamic load and heat flux from each MPI processor onto the root
         Then return the global aero ids, heat fluxes, and loads, which are later written to a file
         """
+        if itime >= scenario.steps:
+            return [], [], []
         all_aero_ids = comm.gather(self.aero_id, root=root)
         if self.transfer is not None:
-            all_aero_loads = comm.gather(self.aero_loads[scenario.id], root=root)
+            _loads = (
+                self.aero_loads[scenario.id]
+                if scenario.steady
+                else self.aero_loads[scenario.id][itime]
+            )
+            all_aero_loads = comm.gather(_loads, root=root)
         else:
             all_aero_loads = []
         if self.thermal_transfer is not None:
-            all_aero_hflux = comm.gather(self.aero_heat_flux[scenario.id], root=root)
+            _hflux = (
+                self.aero_heat_flux[scenario.id]
+                if scenario.steady
+                else self.aero_heat_flux[scenario.id][itime]
+            )
+            all_aero_hflux = comm.gather(_hflux, root=root)
         else:
             all_aero_hflux = []
 

funtofem/model/funtofem_model.py

@@ -414,7 +414,7 @@ def evaluate_composite_functions(self, compute_grad=True):
                 composite_func.evaluate_gradient()
         return
 
-    def read_aero_loads(self, comm, filename, root=0):
+    def _read_aero_loads(self, comm, filename, root=0):
         """
         Read the aerodynamic loads file for the OnewayStructDriver.
 
@@ -573,7 +573,20 @@ def write_struct_loads(self, comm, filename, root=0):
                         fp.write(data)
         return
 
-    def write_aero_loads(self, comm, filename, root=0):
+    @classmethod
+    def _get_loads_filename(cls, prefix, itime: int, suffix=".txt", padding=3):
+        """routine to get default padded 0 aero loads filenames"""
+        return prefix + "_" + f"%0{padding}d" % itime + suffix
+
+    def get_loads_files(self, prefix, suffix=".txt"):
+        """get a list of the loads files for this unsteady scenario"""
+        max_steps = max([scenario.steps for scenario in self.scenarios])
+        return [
+            FUNtoFEMmodel._get_loads_filename(prefix, itime=itime, suffix=suffix)
+            for itime in range(max_steps)
+        ]
+
+    def write_aero_loads(self, comm, filename, itime=0, root=0):
         """
         Write the aerodynamic loads file for the OnewayStructDriver.
 
@@ -599,6 +612,8 @@ def write_aero_loads(self, comm, filename, root=0):
             Global communicator across all FUNtoFEM processors
         filename: str
             The name of the file to be generated
+        itime: int
+            The time step of the loads to write  out (irrelevant if steady)
         root: int
             The rank of the processor that will write the file
         """
@@ -630,7 +645,9 @@ def write_aero_loads(self, comm, filename, root=0):
                 data += f"scenario {scenario.id} {scenario.name} \n"
 
             for body in self.bodies:
-                id, hflux, load = body._collect_aero_loads(comm, scenario, root=root)
+                id, hflux, load = body._collect_aero_loads(
+                    comm, scenario, itime=itime, root=root
+                )
 
                 if comm.rank == root:
                     data += f"body {body.id} {body.name} {body.aero_nnodes} \n"
@@ -647,6 +664,24 @@ def write_aero_loads(self, comm, filename, root=0):
                         fp.write(data)
         return
 
+    def write_unsteady_aero_loads(self, comm, prefix, suffix=".txt", root=0):
+        """
+        Write the aerodynamic loads files for unsteady scenarios for the OnewayStructDriver.
+
+        Parameters
+        ----------
+        comm: MPI communicator
+            Global communicator across all FUNtoFEM processors
+        prefix: str
+            file prefix
+        root: int
+            The rank of the processor that will write the file
+        """
+        loads_files = self.get_loads_files(prefix, suffix)
+        for itime, load_file in enumerate(loads_files):
+            self.write_aero_loads(comm, load_file, itime=itime, root=root)
+        return loads_files
+
     def write_sensitivity_file(
         self, comm, filename, discipline="aerodynamic", root=0, write_dvs: bool = True
     ):