black reformat + throw error if FUN3D 13.6 uses coupling frequency

Author: sean-engelstad

funtofem/driver/funtofem_nlbgs_driver.py

@@ -201,8 +201,8 @@ def _solve_steady_forward(self, scenario, steps=None):
                     if forward_resid > forward_tol:
                         all_converged = False
                         break
-                
-                if all_converged: 
+
+                if all_converged:
                     if self.comm.rank == 0:
                         print(
                             f"F2F Steady Forward analysis of scenario {scenario.name} exited early"
@@ -283,19 +283,26 @@ def _solve_steady_adjoint(self, scenario):
 
             # check for early stopping criterion, exit if meets criterion
             exit_early = False
-            print(f"scenario min adjoint steps = {scenario.min_adjoint_steps}", flush=True)
+            print(
+                f"scenario min adjoint steps = {scenario.min_adjoint_steps}", flush=True
+            )
             if scenario.early_stopping and step > scenario.min_adjoint_steps:
-                all_converged = True # assume all converged until proven otherwise (then when one isn't exit for loop)
-                for isolver,solver in enumerate(self.solvers.solver_list):
+                all_converged = True  # assume all converged until proven otherwise (then when one isn't exit for loop)
+                for isolver, solver in enumerate(self.solvers.solver_list):
                     adjoint_resid = abs(solver.get_adjoint_residual(step=step))
                     adjoint_tol = solver.adjoint_tolerance
                     if isolver == 0:
-                        adjoint_resids = solver.get_adjoint_residual(step=step, all=True)
+                        adjoint_resids = solver.get_adjoint_residual(
+                            step=step, all=True
+                        )
                         print(f"adjoint residuals = {adjoint_resids}", flush=True)
-                    print(f"adjoint step {step} solver {isolver} resid = {adjoint_resid}, tol = {adjoint_tol}", flush=True)
+                    print(
+                        f"adjoint step {step} solver {isolver} resid = {adjoint_resid}, tol = {adjoint_tol}",
+                        flush=True,
+                    )
                     if adjoint_resid > adjoint_tol:
                         all_converged = False
-                     
+
                 if all_converged:
                     if self.comm.rank == 0:
                         print(
@@ -307,7 +314,7 @@ def _solve_steady_adjoint(self, scenario):
                         )
                     exit_early = True
                     break
-            
+
             if exit_early:
                 break
 

funtofem/interface/fun3d_14_interface.py

@@ -315,7 +315,9 @@ def set_functions(self, scenario, bodies):
                     start = 1
                     stop = 1
                     if ct == 1 and scenario.early_stopping:
-                        raise AssertionError("Need to register an aerodynamic function first otherwise the Adjoint early stopping criterion fails")
+                        raise AssertionError(
+                            "Need to register an aerodynamic function first otherwise the Adjoint early stopping criterion fails"
+                        )
                 else:
                     start = 1 if function.start is None else function.start
                     if unsteady:
@@ -674,7 +676,10 @@ def post(self, scenario, bodies):
 
         scalar_resid = abs(np.linalg.norm(resid).real)
         if self.comm.rank == 0:
-            print(f"scalar forward resid in post scenario {scenario.name} = {scalar_resid}", flush=True)
+            print(
+                f"scalar forward resid in post scenario {scenario.name} = {scalar_resid}",
+                flush=True,
+            )
 
         # throw a runtime error if adjoint didn't converge sufficiently
         if scalar_resid > self.forward_tolerance:
@@ -850,7 +855,7 @@ def iterate_adjoint(self, scenario, bodies, step):
             rstep = step
 
         # update the last successful adjoint step in case FUN3D exits early
-        #self._last_adjoint_step = step
+        # self._last_adjoint_step = step
 
         nfuncs = scenario.count_adjoint_functions()
         for ibody, body in enumerate(bodies, 1):
@@ -1001,8 +1006,8 @@ def iterate_adjoint(self, scenario, bodies, step):
         # in FUN3D)
         if self.comm.rank == 0:
             print(f"iterate fun3d adjoint step {rstep}", flush=True)
-        for i_coupled in range(1, scenario.coupling_frequency+1):
-            adj_step = scenario.coupling_frequency * (rstep-1) + i_coupled
+        for i_coupled in range(1, scenario.coupling_frequency + 1):
+            adj_step = scenario.coupling_frequency * (rstep - 1) + i_coupled
             self.fun3d_adjoint.iterate(adj_step)
             self._last_adjoint_step = adj_step
 
@@ -1094,7 +1099,9 @@ def post_adjoint(self, scenario, bodies):
         """
 
         # report warning if flow residual too large
-        resid = self.get_adjoint_residual(step=self._last_adjoint_step, outer=False, all=True)
+        resid = self.get_adjoint_residual(
+            step=self._last_adjoint_step, outer=False, all=True
+        )
         if self.comm.rank == 0:
             print(f"Adjoint residuals = {resid}")
         self._adjoint_done = True
@@ -1106,7 +1113,10 @@ def post_adjoint(self, scenario, bodies):
 
         scalar_resid = abs(np.linalg.norm(resid).real)
         if self.comm.rank == 0:
-            print(f"scalar resid in adjoint post scenario {scenario.name} = {scalar_resid}",flush = True)
+            print(
+                f"scalar resid in adjoint post scenario {scenario.name} = {scalar_resid}",
+                flush=True,
+            )
 
         # throw a runtime error if adjoint didn't converge sufficiently
         if abs(np.linalg.norm(resid).real) > self.adjoint_tolerance:
@@ -1150,7 +1160,9 @@ def get_adjoint_residual(self, step=0, outer=True, all=False):
         all: bool
             whether to return a list of all residuals or a scalar
         """
-        if outer: # just overwrite to the saved last adjoint step since we don't have the scenario data
+        if (
+            outer
+        ):  # just overwrite to the saved last adjoint step since we don't have the scenario data
             step = self._last_adjoint_step
 
         if not self._adjoint_done:

funtofem/interface/fun3d_interface.py

@@ -760,6 +760,14 @@ def initialize_adjoint(self, scenario, bodies):
             if self.model.flow is not None:
                 self.fun3d_adjoint.set_mesh_morph(self.model.flow.mesh_morph)
 
+            # set the funtofem coupling frequency
+            # would need to be implemented in FUN3D adjoint
+            # self.fun3d_adjoint.set_coupling_frequency(scenario.coupling_frequency)
+            if scenario.coupling_frequency != 1:
+                raise AssertionError(
+                    "Funtofem has not implemented funtofem coupling in FUN3D 13.6"
+                )
+
             # Deform the aero mesh before finishing FUN3D initialization
             for ibody, body in enumerate(bodies, 1):
                 aero_disps = body.get_aero_disps(
@@ -945,9 +953,9 @@ def iterate_adjoint(self, scenario, bodies, step):
 
         # Update the aerodynamic and grid adjoint variables (Note: step starts at 1
         # in FUN3D)
-        for i_coupled in range(1,scenario.coupling_frequency+1):
-           adj_step = scenario.coupling_frequency * (rstep-1) + i_coupled
-           self.fun3d_adjoint.iterate(adj_step)
+        for i_coupled in range(1, scenario.coupling_frequency + 1):
+            adj_step = scenario.coupling_frequency * (rstep - 1) + i_coupled
+            self.fun3d_adjoint.iterate(adj_step)
 
         for ibody, body in enumerate(bodies, 1):
             # Extract aero_disps_ajp = dG/du_A^T psi_G from FUN3D

funtofem/model/scenario.py

@@ -153,7 +153,9 @@ def __init__(
 
         # early stopping criterion
         self.min_forward_steps = (
-            min_forward_steps if min_forward_steps is not None else uncoupled_steps+self.UNCOUPLED_STEP_BUFFER
+            min_forward_steps
+            if min_forward_steps is not None
+            else uncoupled_steps + self.UNCOUPLED_STEP_BUFFER
         )
         self.min_adjoint_steps = (
             min_adjoint_steps if min_adjoint_steps is not None else 0
@@ -180,7 +182,13 @@ def __init__(
             self.add_variable("aerodynamic", zrate)
 
     @classmethod
-    def steady(cls, name: str, steps: int, coupling_frequency:int=1, uncoupled_steps: int = 0):
+    def steady(
+        cls,
+        name: str,
+        steps: int,
+        coupling_frequency: int = 1,
+        uncoupled_steps: int = 0,
+    ):
         return cls(
             name=name,
             steady=True,
@@ -206,20 +214,20 @@ def unsteady(
         )
 
     @property
-    def adjoint_steps(self) -> int: 
+    def adjoint_steps(self) -> int:
         """
-        in the steady case it's best to choose the 
+        in the steady case it's best to choose the
         adjoint steps based on the funtofem coupling frequency
         """
         if self._adjoint_steps is not None and self.steady:
             return self._adjoint_steps
         elif not self.steady:
-            return None # defaults to number of steps in unsteady case
-        else: # choose it based on funtofem coupling frequency in steady case
-            return int(np.ceil(self.steps/self.coupling_frequency))
+            return None  # defaults to number of steps in unsteady case
+        else:  # choose it based on funtofem coupling frequency in steady case
+            return int(np.ceil(self.steps / self.coupling_frequency))
 
     @adjoint_steps.setter
-    def adjoint_steps(self, new_steps:int):
+    def adjoint_steps(self, new_steps: int):
         assert self.steady
         self._adjoint_steps = new_steps
 

tests/fun3d_tests/fully_coupled_disc_deriv/turb_AE/test_adjoint_speed.py

@@ -1,6 +1,7 @@
 """
 Unittest for FUN3D 14.0.2 finite-difference test
 """
+
 import numpy as np, unittest, importlib, os
 from mpi4py import MPI
 import time
@@ -39,15 +40,17 @@ class TestFun3dTacs(unittest.TestCase):
     FILEPATH = os.path.join(results_folder, FILENAME)
 
     def test_alpha_turbulent_aeroelastic(self):
-        start_time = time.time()        
+        start_time = time.time()
 
         # build the funtofem model with one body and scenario
         model = FUNtoFEMmodel("plate")
         plate = Body.aeroelastic("plate", boundary=2)
         plate.register_to(model)
 
         # build the scenario
-        test_scenario = Scenario.steady("turbulent_beta", steps=500, coupling_frequency=20, uncoupled_steps=10)
+        test_scenario = Scenario.steady(
+            "turbulent_beta", steps=500, coupling_frequency=20, uncoupled_steps=10
+        )
         test_scenario.set_stop_criterion(early_stopping=True, min_forward_steps=50)
         test_scenario.set_temperature(T_ref=300.0, T_inf=300.0)
         Function.lift().register_to(test_scenario)
@@ -60,7 +63,13 @@ def test_alpha_turbulent_aeroelastic(self):
 
         # build the solvers and coupled driver
         solvers = SolverManager(comm)
-        solvers.flow = Fun3d14Interface(comm, model, fun3d_dir="meshes", forward_tolerance=1e-11, adjoint_tolerance=1e-9)
+        solvers.flow = Fun3d14Interface(
+            comm,
+            model,
+            fun3d_dir="meshes",
+            forward_tolerance=1e-11,
+            adjoint_tolerance=1e-9,
+        )
 
         solvers.structural = TacsSteadyInterface.create_from_bdf(
             model, comm, nprocs=1, bdf_file=bdf_filename, prefix=output_dir
@@ -78,25 +87,26 @@ def test_alpha_turbulent_aeroelastic(self):
 
         # run the complex step test on the model and driver
         # max_rel_error = TestResult.complex_step(
-        #max_rel_error = TestResult.finite_difference(
+        # max_rel_error = TestResult.finite_difference(
         #    "fun3d+tacs-turbulent-aeroelastic-flow",
         #    model,
         #    driver,
         #    TestFun3dTacs.FILEPATH,
-        #)
-        #self.assertTrue(max_rel_error < 1e-7)
- 
+        # )
+        # self.assertTrue(max_rel_error < 1e-7)
+
         driver.solve_forward()
         time2 = time.time()
         driver.solve_adjoint()
 
         _end_time = time.time()
         dt_full = _end_time - start_time
         dt_adjoint = _end_time - time2
-        
+
         print(f"full time = {dt_full} sec", flush=True)
         print(f"adjoint time = {dt_adjoint} sec", flush=True)
 
+
 if __name__ == "__main__":
     # open and close the file to reset it
     if comm.rank == 0:

tests/fun3d_tests/fully_coupled_disc_deriv/turb_AE/test_early_stop_FD.py

@@ -1,6 +1,7 @@
 """
 Unittest for FUN3D 14.0.2 finite-difference test
 """
+
 import numpy as np, unittest, importlib, os
 from mpi4py import MPI
 import time
@@ -39,21 +40,26 @@
 forward_tol = 1e-11
 adjoint_tol = 1e-9
 
+
 class TestFun3dTacs(unittest.TestCase):
     FILENAME = "fun3d-tacs-driver.txt"
     FILEPATH = os.path.join(results_folder, FILENAME)
 
     def test_alpha_turbulent_aeroelastic(self):
-        start_time = time.time()        
+        start_time = time.time()
 
         # build the funtofem model with one body and scenario
         model = FUNtoFEMmodel("plate")
         plate = Body.aeroelastic("plate", boundary=2)
         plate.register_to(model)
 
         # build the scenario
-        test_scenario = Scenario.steady("turbulent_beta", steps=500, coupling_frequency=20, uncoupled_steps=10)
-        test_scenario.set_stop_criterion(early_stopping=early_stopping, min_forward_steps=50)
+        test_scenario = Scenario.steady(
+            "turbulent_beta", steps=500, coupling_frequency=20, uncoupled_steps=10
+        )
+        test_scenario.set_stop_criterion(
+            early_stopping=early_stopping, min_forward_steps=50
+        )
         test_scenario.set_temperature(T_ref=300.0, T_inf=300.0)
         Function.lift().register_to(test_scenario)
         Function.ksfailure(ks_weight=10.0).register_to(test_scenario)
@@ -65,7 +71,13 @@ def test_alpha_turbulent_aeroelastic(self):
 
         # build the solvers and coupled driver
         solvers = SolverManager(comm)
-        solvers.flow = Fun3d14Interface(comm, model, fun3d_dir="meshes", forward_tolerance=forward_tol, adjoint_tolerance=adjoint_tol)
+        solvers.flow = Fun3d14Interface(
+            comm,
+            model,
+            fun3d_dir="meshes",
+            forward_tolerance=forward_tol,
+            adjoint_tolerance=adjoint_tol,
+        )
 
         solvers.structural = TacsSteadyInterface.create_from_bdf(
             model, comm, nprocs=1, bdf_file=bdf_filename, prefix=output_dir
@@ -102,18 +114,28 @@ def test_thick_turbulent_aeroelastic(self):
         plate.register_to(model)
 
         # build the scenario
-        test_scenario = Scenario.steady("turbulent_beta", steps=500, coupling_frequency=20, uncoupled_steps=10)
-        test_scenario.set_stop_criterion(early_stopping=early_stopping, min_forward_steps=50)
+        test_scenario = Scenario.steady(
+            "turbulent_beta", steps=500, coupling_frequency=20, uncoupled_steps=10
+        )
+        test_scenario.set_stop_criterion(
+            early_stopping=early_stopping, min_forward_steps=50
+        )
         test_scenario.set_temperature(T_ref=300.0, T_inf=300.0)
-        Function.lift().register_to(test_scenario)        
+        Function.lift().register_to(test_scenario)
         Function.ksfailure(ks_weight=10.0).register_to(test_scenario)
         Function.drag().register_to(test_scenario)
         test_scenario.set_flow_ref_vals(qinf=1.05e5)
         test_scenario.register_to(model)
 
         # build the solvers and coupled driver
         solvers = SolverManager(comm)
-        solvers.flow = Fun3d14Interface(comm, model, fun3d_dir="meshes", forward_tolerance=forward_tol, adjoint_tolerance=adjoint_tol)
+        solvers.flow = Fun3d14Interface(
+            comm,
+            model,
+            fun3d_dir="meshes",
+            forward_tolerance=forward_tol,
+            adjoint_tolerance=adjoint_tol,
+        )
 
         solvers.structural = TacsSteadyInterface.create_from_bdf(
             model, comm, nprocs=1, bdf_file=bdf_filename, prefix=output_dir
@@ -139,6 +161,7 @@ def test_thick_turbulent_aeroelastic(self):
         )
         self.assertTrue(max_rel_error < 2e-3)
 
+
 if __name__ == "__main__":
     # open and close the file to reset it
     if comm.rank == 0: