From 6df852eab9398983485913d5156addfbc1ccfc4f Mon Sep 17 00:00:00 2001 From: Sean Engelstad Date: Mon, 6 Nov 2023 11:37:21 -0500 Subject: [PATCH] fix aero loads (were almost zero before) (#250) --- funtofem/driver/oneway_struct_driver.py | 8 +- funtofem/model/body.py | 31 ++- funtofem/model/funtofem_model.py | 206 +++++++++--------- ..._loads_file.py => test_aero_loads_file.py} | 112 +++++----- .../framework/test_struct_loads_file.py | 76 +++++++ 5 files changed, 261 insertions(+), 172 deletions(-) rename tests/unit_tests/framework/{test_loads_file.py => test_aero_loads_file.py} (65%) create mode 100644 tests/unit_tests/framework/test_struct_loads_file.py diff --git a/funtofem/driver/oneway_struct_driver.py b/funtofem/driver/oneway_struct_driver.py index a7d62870..e6039423 100644 --- a/funtofem/driver/oneway_struct_driver.py +++ b/funtofem/driver/oneway_struct_driver.py @@ -221,6 +221,7 @@ def prime_loads_from_file( nprocs, transfer_settings, external_shape=False, + init_transfer=False, ): """ Used to prime aero loads for optimization over tacs analysis with shape change and tacs aim @@ -242,6 +243,8 @@ def prime_loads_from_file( 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() @@ -271,12 +274,15 @@ def prime_loads_from_file( body.initialize_variables(scenario) body._distribute_aero_loads(loads_data) - return cls( + tacs_driver = cls( solvers, model, nprocs=nprocs, external_shape=external_shape, ) + if init_transfer: + tacs_driver._transfer_fixed_aero_loads() + return tacs_driver @property def manager(self, hot_start: bool = False): diff --git a/funtofem/model/body.py b/funtofem/model/body.py index c82e44ef..d1f2fa37 100644 --- a/funtofem/model/body.py +++ b/funtofem/model/body.py @@ -1598,18 +1598,31 @@ def _distribute_aero_loads(self, data): """ distribute the aero loads and heat flux from a loads file """ + print(f"F2F - starting to distribute loads") + for scenario_id in data: scenario_data = data[scenario_id] + + # create a dict for this entry + scenario_entry_dict = {} for entry in scenario_data: - for ind, aero_id in enumerate(self.aero_id): - if entry["aeroID"] == aero_id and entry["bodyName"] == self.name: - if self.transfer is not None: - self.aero_loads[scenario_id][3 * ind : 3 * ind + 3] = entry[ - "load" - ] - if self.thermal_transfer is not None: - self.aero_heat_flux[scenario_id][ind] = entry["hflux"] - break + if entry["bodyName"] == self.name: + scenario_entry_dict[entry["aeroID"]] = { + "load": entry["load"], + "hflux": entry["hflux"], + } + + 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"] + + print(f"F2F - done distribute loads") def _collect_aero_mesh(self, comm, root=0): """ diff --git a/funtofem/model/funtofem_model.py b/funtofem/model/funtofem_model.py index afee76ca..06bf1f1a 100644 --- a/funtofem/model/funtofem_model.py +++ b/funtofem/model/funtofem_model.py @@ -400,9 +400,9 @@ def evaluate_composite_functions(self, compute_grad=True): composite_func.evaluate_gradient() return - def write_aero_loads(self, comm, filename, root=0): + def read_aero_loads(self, comm, filename, root=0): """ - Write the aerodynamic loads file for the OnewayStructDriver. + Read the aerodynamic loads file for the OnewayStructDriver. This file contains the following information: @@ -429,50 +429,81 @@ def write_aero_loads(self, comm, filename, root=0): root: int The rank of the processor that will write the file """ + loads_data = None + mesh_data = None + if comm.rank == root: - data = "" - # Specify the number of scenarios in file - data += f"{len(self.bodies)} {len(self.scenarios)} \n" - data += "aeromesh" + "\n" + scenario_data = None + loads_data = {} + mesh_data = {} + + with open(filename, "r") as fp: + for line in fp.readlines(): + entries = line.strip().split(" ") + # print("==> entries: ", entries) + if len(entries) == 2: + assert int(entries[1]) == len(self.scenarios) + assert int(entries[0]) == len(self.bodies) + + elif len(entries) == 3 and entries[0] == "scenario": + matching_scenario = False + for scenario in self.scenarios: + if str(scenario.name).strip() == str(entries[2]).strip(): + matching_scenario = True + break + assert matching_scenario + if scenario_data is not None: + loads_data[scenario.id] = scenario_data + scenario_data = [] + elif len(entries) == 4 and entries[0] == "body_mesh": + body_name = entries[2] + mesh_data[body_name] = {"aeroID": [], "aeroX": []} + elif len(entries) == 4 and entries[0] != "body": + mesh_data[body_name]["aeroID"] += [entries[0]] + mesh_data[body_name]["aeroX"] += entries[1:4] + elif len(entries) == 5: + entry = { + "bodyName": body_name, + "aeroID": entries[0], + "load": entries[1:4], + "hflux": entries[4], + } + scenario_data.append(entry) + + loads_data[scenario.id] = scenario_data + + loads_data = comm.bcast(loads_data, root=root) + mesh_data = comm.bcast(mesh_data, root=root) + + # initialize the mesh data for body in self.bodies: + global_aero_x = np.array(mesh_data[body.name]["aeroX"]) + global_aero_ids = np.array(mesh_data[body.name]["aeroID"]) + + body_ind = np.array([_ for _ in range(len(global_aero_ids))]) if comm.rank == root: - data += f"body_mesh {body.id} {body.name} {body.aero_nnodes} \n" + split_body_ind = np.array_split(body_ind, comm.Get_size()) + else: + split_body_ind = None - id, aeroX = body._collect_aero_mesh(comm, root=root) + local_body_ind = comm.scatter(split_body_ind, root=root) - if comm.rank == root: - for i in range(len(id)): - data += "{} {} {} {} \n".format( - int(id[i]), - aeroX[3 * i + 0].real, - aeroX[3 * i + 1].real, - aeroX[3 * i + 2].real, - ) - if comm.rank == root: - data += f"aeroloads \n" + local_aero_ids = global_aero_ids[local_body_ind] - for scenario in self.scenarios: - if comm.rank == root: - data += f"scenario {scenario.id} {scenario.name} \n" + aero_x_ind = ( + [3 * i for i in local_body_ind] + + [3 * i + 1 for i in local_body_ind] + + [3 * i + 2 for i in local_body_ind] + ) + aero_x_ind = sorted(aero_x_ind) - for body in self.bodies: - id, hflux, load = body._collect_aero_loads(comm, scenario, root=root) + local_aero_x = list(global_aero_x[aero_x_ind]) - if comm.rank == root: - data += f"body {body.id} {body.name} {body.aero_nnodes} \n" - for i in range(len(id)): - data += "{} {} {} {} {} \n".format( - int(id[i]), - load[3 * i + 0].real, - load[3 * i + 1].real, - load[3 * i + 2].real, - float(hflux[i].real), - ) + body.initialize_aero_nodes(local_aero_x, local_aero_ids) - with open(filename, "w") as fp: - fp.write(data) - return + # return the loads data + return loads_data def write_struct_loads(self, comm, filename, root=0): """ @@ -528,9 +559,9 @@ def write_struct_loads(self, comm, filename, root=0): fp.write(data) return - def read_aero_loads(self, comm, filename, root=0): + def write_aero_loads(self, comm, filename, root=0): """ - Read the aerodynamic loads file for the OnewayStructDriver. + Write the aerodynamic loads file for the OnewayStructDriver. This file contains the following information: @@ -557,81 +588,50 @@ def read_aero_loads(self, comm, filename, root=0): root: int The rank of the processor that will write the file """ - loads_data = None - mesh_data = None - if comm.rank == root: - scenario_data = None - loads_data = {} - mesh_data = {} - - with open(filename, "r") as fp: - for line in fp.readlines(): - entries = line.strip().split(" ") - # print("==> entries: ", entries) - if len(entries) == 2: - assert int(entries[1]) == len(self.scenarios) - assert int(entries[0]) == len(self.bodies) - - elif len(entries) == 3 and entries[0] == "scenario": - matching_scenario = False - for scenario in self.scenarios: - if str(scenario.name).strip() == str(entries[2]).strip(): - matching_scenario = True - break - assert matching_scenario - if scenario_data is not None: - loads_data[scenario.id] = scenario_data - scenario_data = [] - elif len(entries) == 4 and entries[0] == "body_mesh": - body_name = entries[2] - mesh_data[body_name] = {"aeroID": [], "aeroX": []} - elif len(entries) == 4 and entries[0] != "body": - mesh_data[body_name]["aeroID"] += [entries[0]] - mesh_data[body_name]["aeroX"] += entries[1:4] - - elif len(entries) == 5: - entry = { - "bodyName": body_name, - "aeroID": entries[0], - "load": entries[1:4], - "hflux": entries[4], - } - scenario_data.append(entry) - - loads_data[scenario.id] = scenario_data - - loads_data = comm.bcast(loads_data, root=root) - mesh_data = comm.bcast(mesh_data, root=root) + data = "" + # Specify the number of scenarios in file + data += f"{len(self.bodies)} {len(self.scenarios)} \n" + data += "aeromesh" + "\n" - # initialize the mesh data for body in self.bodies: - global_aero_x = np.array(mesh_data[body.name]["aeroX"]) - global_aero_ids = np.array(mesh_data[body.name]["aeroID"]) - - body_ind = np.array([_ for _ in range(len(global_aero_ids))]) if comm.rank == root: - split_body_ind = np.array_split(body_ind, comm.Get_size()) - else: - split_body_ind = None + data += f"body_mesh {body.id} {body.name} {body.aero_nnodes} \n" - local_body_ind = comm.scatter(split_body_ind, root=root) + id, aeroX = body._collect_aero_mesh(comm, root=root) - local_aero_ids = global_aero_ids[local_body_ind] + if comm.rank == root: + for i in range(len(id)): + data += "{} {} {} {} \n".format( + int(id[i]), + aeroX[3 * i + 0].real, + aeroX[3 * i + 1].real, + aeroX[3 * i + 2].real, + ) + if comm.rank == root: + data += f"aeroloads \n" - aero_x_ind = ( - [3 * i for i in local_body_ind] - + [3 * i + 1 for i in local_body_ind] - + [3 * i + 2 for i in local_body_ind] - ) - aero_x_ind = sorted(aero_x_ind) + for scenario in self.scenarios: + if comm.rank == root: + data += f"scenario {scenario.id} {scenario.name} \n" - local_aero_x = list(global_aero_x[aero_x_ind]) + for body in self.bodies: + id, hflux, load = body._collect_aero_loads(comm, scenario, root=root) - body.initialize_aero_nodes(local_aero_x, local_aero_ids) + if comm.rank == root: + data += f"body {body.id} {body.name} {body.aero_nnodes} \n" + for i in range(len(id)): + data += "{} {} {} {} {} \n".format( + int(id[i]), + load[3 * i + 0].real, + load[3 * i + 1].real, + load[3 * i + 2].real, + float(hflux[i].real), + ) - # return the loads data - return loads_data + with open(filename, "w") as fp: + fp.write(data) + return def write_sensitivity_file( self, comm, filename, discipline="aerodynamic", root=0, write_dvs: bool = True diff --git a/tests/unit_tests/framework/test_loads_file.py b/tests/unit_tests/framework/test_aero_loads_file.py similarity index 65% rename from tests/unit_tests/framework/test_loads_file.py rename to tests/unit_tests/framework/test_aero_loads_file.py index ab3af93e..bd66ec4a 100644 --- a/tests/unit_tests/framework/test_loads_file.py +++ b/tests/unit_tests/framework/test_aero_loads_file.py @@ -29,9 +29,9 @@ struct_loads_file = os.path.join(output_dir, "struct_loads.txt") -class TestLoadsFile(unittest.TestCase): - # N_PROCS = 2 - FILENAME = "test_oneway_loads_file.txt" +class TestAeroLoadsFile(unittest.TestCase): + N_PROCS = 2 + FILENAME = "test_aero_loads_file.txt" FILEPATH = os.path.join(results_folder, FILENAME) def test_loads_file_aeroelastic(self): @@ -66,25 +66,35 @@ def test_loads_file_aeroelastic(self): FUNtoFEMnlbgs( solvers, transfer_settings=transfer_settings, model=f2f_model ).solve_forward() + + # save initial loads in an array + orig_aero_loads = plate.aero_loads[scenario.id] * 1.0 f2f_model.write_aero_loads(comm, aero_loads_file, root=0) + # zero the aero loads + plate.aero_loads[scenario.id] *= 0.0 + # ----------------------------------------------- # Read the loads file and test the oneway driver # ----------------------------------------------- solvers.flow = None - oneway_driver = OnewayStructDriver.prime_loads_from_file( + OnewayStructDriver.prime_loads_from_file( aero_loads_file, solvers, f2f_model, 1, transfer_settings ) - max_rel_error = TestResult.derivative_test( - "testaero=>tacs_loads-aeroelastic", - f2f_model, - oneway_driver, - TestLoadsFile.FILEPATH, - complex_mode=complex_mode, - ) - rtol = 1e-7 if complex_mode else 1e-3 - self.assertTrue(max_rel_error < rtol) + # verify the aero loads are the same + new_aero_loads = plate.aero_loads[scenario.id] + diff_aero_loads = new_aero_loads - orig_aero_loads + orig_norm = np.max(np.abs(orig_aero_loads)) + abs_err_norm = np.max(np.abs(diff_aero_loads)) + rel_err_norm = abs_err_norm / orig_norm + + print("aeroelastic aero loads test:") + print(f"\taero load norm = {orig_norm:.5f}") + print(f"\tabs error = {abs_err_norm:.5f}") + print(f"\trel error = {rel_err_norm:.5f}") + rtol = 1e-7 + self.assertTrue(rel_err_norm < rtol) return def test_loads_file_aerothermoelastic(self): @@ -119,64 +129,48 @@ def test_loads_file_aerothermoelastic(self): FUNtoFEMnlbgs( solvers, transfer_settings=transfer_settings, model=f2f_model ).solve_forward() + + # save initial loads in an array + orig_aero_loads = plate.aero_loads[scenario.id] * 1.0 + orig_aero_hflux = plate.aero_heat_flux[scenario.id] * 1.0 f2f_model.write_aero_loads(comm, aero_loads_file, root=0) + # zero the aero loads + plate.aero_loads[scenario.id] *= 0.0 + # ----------------------------------------------- # Read the loads file and test the oneway driver # ----------------------------------------------- solvers.flow = None - oneway_driver = OnewayStructDriver.prime_loads_from_file( + OnewayStructDriver.prime_loads_from_file( aero_loads_file, solvers, f2f_model, 1, transfer_settings ) - max_rel_error = TestResult.derivative_test( - "testaero=>tacs_loads-aerothermoelastic", - f2f_model, - oneway_driver, - TestLoadsFile.FILEPATH, - complex_mode=complex_mode, - ) - rtol = 1e-7 if complex_mode else 1e-3 - self.assertTrue(max_rel_error < rtol) - return - - def test_struct_loads_file_aerothermoelastic(self): - # --------------------------- - # Write the loads file - # --------------------------- - # build the model and driver - f2f_model = FUNtoFEMmodel("wedge") - plate = Body.aerothermoelastic("plate", boundary=1) - Variable.structural("thickness").set_bounds( - lower=0.01, value=0.1, upper=1.0 - ).register_to(plate) - plate.register_to(f2f_model) - - # build the scenario - scenario = Scenario.steady("test", steps=150) - Function.ksfailure().register_to(scenario) - scenario.register_to(f2f_model) - - # make the solvers for a CFD analysis to store and write the loads file - solvers = SolverManager(comm) - solvers.flow = TestAerodynamicSolver(comm, f2f_model) - solvers.structural = TacsInterface.create_from_bdf( - f2f_model, - comm, - 1, - bdf_filename, - callback=thermoelasticity_callback, - output_dir=output_dir, - ) - transfer_settings = TransferSettings(npts=5) - FUNtoFEMnlbgs( - solvers, transfer_settings=transfer_settings, model=f2f_model - ).solve_forward() - f2f_model.write_struct_loads(comm, struct_loads_file, root=0) + # verify the aero loads are the same + new_aero_loads = plate.aero_loads[scenario.id] + diff_aero_loads = new_aero_loads - orig_aero_loads + abs_err_norm = np.max(np.abs(diff_aero_loads)) + rel_err_norm = abs_err_norm / np.max(np.abs(orig_aero_loads)) + print("aero loads + aero hflux test:") + print("aero loads comparison") + print(f"\tabs error = {abs_err_norm:.5f}") + print(f"\trel error = {rel_err_norm:.5f}") + + # verify the aero heat flux is the same + new_aero_hflux = plate.aero_heat_flux[scenario.id] + diff_aero_hflux = new_aero_hflux - orig_aero_hflux + abs_err_norm2 = np.max(np.abs(diff_aero_hflux)) + rel_err_norm2 = abs_err_norm / np.max(np.abs(orig_aero_hflux)) + print("aero hflux comparison") + print(f"\tabs error = {abs_err_norm2:.5f}") + print(f"\trel error = {rel_err_norm2:.5f}") + rtol = 1e-7 + self.assertTrue(rel_err_norm < rtol) + self.assertTrue(rel_err_norm2 < rtol) return if __name__ == "__main__": if comm.rank == 0: - open(TestLoadsFile.FILEPATH, "w").close() + open(TestAeroLoadsFile.FILEPATH, "w").close() unittest.main() diff --git a/tests/unit_tests/framework/test_struct_loads_file.py b/tests/unit_tests/framework/test_struct_loads_file.py new file mode 100644 index 00000000..25be216e --- /dev/null +++ b/tests/unit_tests/framework/test_struct_loads_file.py @@ -0,0 +1,76 @@ +import unittest, os, numpy as np, sys +from _bdf_test_utils import elasticity_callback, thermoelasticity_callback +from mpi4py import MPI +from tacs import TACS + +np.set_printoptions(threshold=sys.maxsize) + +from funtofem import TransferScheme +from funtofem.model import FUNtoFEMmodel, Variable, Scenario, Body, Function +from funtofem.interface import ( + TestAerodynamicSolver, + TacsInterface, + SolverManager, + TestResult, + make_test_directories, +) +from funtofem.driver import FUNtoFEMnlbgs, TransferSettings, OnewayStructDriver + +np.random.seed(1234567) + +base_dir = os.path.dirname(os.path.abspath(__file__)) +bdf_filename = os.path.join(base_dir, "input_files", "test_bdf_file.bdf") +comm = MPI.COMM_WORLD + +complex_mode = TransferScheme.dtype == complex and TACS.dtype == complex + +results_folder, output_dir = make_test_directories(comm, base_dir) +aero_loads_file = os.path.join(output_dir, "aero_loads.txt") +struct_loads_file = os.path.join(output_dir, "struct_loads.txt") + + +class TestStructLoadsFile(unittest.TestCase): + N_PROCS = 2 + FILENAME = "test_struct_loads_file.txt" + FILEPATH = os.path.join(results_folder, FILENAME) + + def test_struct_loads_file_aerothermoelastic(self): + # --------------------------- + # Write the loads file + # --------------------------- + # build the model and driver + f2f_model = FUNtoFEMmodel("wedge") + plate = Body.aerothermoelastic("plate", boundary=1) + Variable.structural("thickness").set_bounds( + lower=0.01, value=0.1, upper=1.0 + ).register_to(plate) + plate.register_to(f2f_model) + + # build the scenario + scenario = Scenario.steady("test", steps=150) + Function.ksfailure().register_to(scenario) + scenario.register_to(f2f_model) + + # make the solvers for a CFD analysis to store and write the loads file + solvers = SolverManager(comm) + solvers.flow = TestAerodynamicSolver(comm, f2f_model) + solvers.structural = TacsInterface.create_from_bdf( + f2f_model, + comm, + 1, + bdf_filename, + callback=thermoelasticity_callback, + output_dir=output_dir, + ) + transfer_settings = TransferSettings(npts=5) + FUNtoFEMnlbgs( + solvers, transfer_settings=transfer_settings, model=f2f_model + ).solve_forward() + f2f_model.write_struct_loads(comm, struct_loads_file, root=0) + return + + +if __name__ == "__main__": + if comm.rank == 0: + open(TestStructLoadsFile.FILEPATH, "w").close() + unittest.main()