From 582f7a964f2b76165e2eeb10e1a3711185b1d132 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Gerhard=20Br=C3=A4unlich?= Date: Fri, 17 Oct 2025 13:28:02 +0200 Subject: [PATCH 1/9] problems.power_electronics: Fix the interface to optimize Even though optimize is not implemented --- engibench/problems/power_electronics/v0.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/engibench/problems/power_electronics/v0.py b/engibench/problems/power_electronics/v0.py index b1dfc98..88fb9c9 100644 --- a/engibench/problems/power_electronics/v0.py +++ b/engibench/problems/power_electronics/v0.py @@ -192,7 +192,7 @@ def simulate(self, design: npt.NDArray, config: dict[str, Any] | None = None) -> DcGain, VoltageRipple = process_log_file(self.config.log_file_path) return np.array([DcGain, VoltageRipple]) - def optimize(self, _starting_point: npt.NDArray, _config: dict[str, Any] | None = None) -> NoReturn: + def optimize(self, starting_point: npt.NDArray, config: dict[str, Any] | None = None) -> NoReturn: """Optimize the design variable. Not applicable for this problem.""" raise NotImplementedError("Not yet implemented") From 7e608aa9b39f801e480a896f8dc82c9f1de020a9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Gerhard=20Br=C3=A4unlich?= Date: Fri, 17 Oct 2025 13:29:17 +0200 Subject: [PATCH 2/9] utils.container: Remove uninteresting frame from the stack trace in case of an error --- engibench/utils/container.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/engibench/utils/container.py b/engibench/utils/container.py index 3f603e4..a4d2fd8 100644 --- a/engibench/utils/container.py +++ b/engibench/utils/container.py @@ -53,7 +53,7 @@ def run( # noqa: PLR0913 Command: {" ".join(command)} stdout: {result.stdout.decode() if result.stdout else "No output"} stderr: {result.stderr.decode() if result.stderr else "No output"}""" - raise RuntimeError(msg) from e + raise RuntimeError(msg) from None class ContainerRuntime: From 399ecfb08de7c4628370ec82395a52da4506c084 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Gerhard=20Br=C3=A4unlich?= Date: Fri, 17 Oct 2025 13:34:17 +0200 Subject: [PATCH 3/9] airfoil: Replace sigint/sigterm cleanup handling by passing -m mpi4py to python --- .../airfoil/templates/airfoil_analysis.py | 394 +++++++-------- .../problems/airfoil/templates/airfoil_opt.py | 477 ++++++++---------- engibench/problems/airfoil/v0.py | 102 ++-- 3 files changed, 434 insertions(+), 539 deletions(-) diff --git a/engibench/problems/airfoil/templates/airfoil_analysis.py b/engibench/problems/airfoil/templates/airfoil_analysis.py index 5589c4d..418b8b8 100644 --- a/engibench/problems/airfoil/templates/airfoil_analysis.py +++ b/engibench/problems/airfoil/templates/airfoil_analysis.py @@ -1,212 +1,182 @@ -# mypy: ignore-errors -"""This file is based on the MACHAero tutorials. - -https://github.com/mdolab/MACH-Aero/blob/main/tutorial/ - -TEMPLATED VARS: -- $mesh_fname: Path to the mesh file. -- $output_dir: Path to the output directory. -- $task: The task to perform: "analysis" or "polar". -- $mach: The Mach number (float). -- $reynolds: The Reynolds number (float). -- $temperature: The temperature (float). -- $alpha: The angle of attack (float). -""" - -import numpy as np -import os -from adflow import ADFLOW -from baseclasses import AeroProblem -from mpi4py import MPI -import signal -import atexit -import sys - -_exit_code = 0 - -def set_exit_code(code): - global _exit_code - _exit_code = code - -def cleanup_handler(signum=None, frame=None): - # Clean termination code for MPI - global _exit_code - try: - from mpi4py import MPI - if MPI.Is_initialized(): - if signum is not None: - _exit_code = 1 - MPI.COMM_WORLD.Abort(_exit_code) - except: - pass - sys.exit(_exit_code) - -if __name__ == "__main__": - - # Register signal handlers - signal.signal(signal.SIGTERM, cleanup_handler) - signal.signal(signal.SIGINT, cleanup_handler) - atexit.register(cleanup_handler) - - try: - mesh_fname = $mesh_fname - output_dir = $output_dir - task = $task - - # mach number - mach = $mach - # Reynolds number - reynolds = $reynolds - # altitude - altitude = $altitude - # temperature - T = $temperature - # Whether to use altitude - use_altitude = $use_altitude - # Reynold's Length - reynoldsLength = 1.0 - - comm = MPI.COMM_WORLD - print(f"Processor {comm.rank} of {comm.size} is running") - if not os.path.exists(output_dir): - if comm.rank == 0: - os.mkdir(output_dir) - - # rst ADflow options - aeroOptions = { - # I/O Parameters - "gridFile": mesh_fname, - "outputDirectory": output_dir, - "monitorvariables": ["cl", "cd","resrho","resrhoe"], - "writeTecplotSurfaceSolution": True, - # Physics Parameters - "equationType": "RANS", - "smoother": "DADI", - "rkreset": True, - "nrkreset": 10, - # Solver Parameters - "MGCycle": "sg", - # ANK Solver Parameters - "useANKSolver": True, - "ankswitchtol": 1e-1, - "liftIndex": 2, - "nsubiterturb": 10, - # NK Solver Parameters - "useNKSolver": True, - "NKSwitchTol": 1e-4, - # Termination Criteria - "L2Convergence": 1e-9, - "L2ConvergenceCoarse": 1e-4, - "nCycles": 5000, - } - print("ADflow options:") - # rst Start ADflow - # Create solver - CFDSolver = ADFLOW(options=aeroOptions) - - # Add features - # CFDSolver.addLiftDistribution(150, "z") - span = 1.0 - pos = np.array([0.5]) * span - CFDSolver.addSlices("z", pos, sliceType="absolute") - - # rst Create AeroProblem - alpha = $alpha - - if use_altitude: - ap = AeroProblem( - name="fc", - alpha=alpha, - mach=mach, - altitude=altitude, - areaRef=1.0, - chordRef=1.0, - evalFuncs=["cl", "cd"], - ) - else: - ap = AeroProblem( - name="fc", - alpha=alpha, - mach=mach, - T=T, - reynolds=reynolds, - reynoldsLength=reynoldsLength, - areaRef=1.0, - chordRef=1.0, - evalFuncs=["cl", "cd"], - ) - - # rst Run ADflow - if task == "analysis": - print("Running analysis") - # Solve - CFDSolver(ap) - # rst Evaluate and print - funcs = {} - CFDSolver.evalFunctions(ap, funcs) - # Print the evaluated functions - if comm.rank == 0: - CL = funcs[f"{ap.name}_cl"] - CD = funcs[f"{ap.name}_cd"] - # Save the lift and drag coefficients to a file - outputs = np.array([mach, reynolds, alpha, CL, CD]) - np.save(os.path.join(output_dir, "outputs.npy"), outputs) - - # rst Create polar arrays - elif task == "polar": - print("Running polar") - # Create an array of alpha values. - # In this case we create 5 random alpha values between 0 and 10 - # from numpy.random import uniform - alphaList = np.linspace(0, 20, 50) - # Sort the alpha values - alphaList.sort() - - # Create storage for the evaluated lift and drag coefficients - CLList = [] - CDList = [] - reslist = [] - # rst Start loop - # Loop over the alpha values and evaluate the polar - for alpha in alphaList: - # rst update AP - # Update the name in the AeroProblem. This allows us to modify the - # output file names with the current alpha. - ap.name = f"fc_{alpha:4.2f}" - - # Update the alpha in aero problem and print it to the screen. - ap.alpha = alpha - if comm.rank == 0: - print(f"current alpha: {ap.alpha}") - - # rst Run ADflow polar - # Solve the flow - CFDSolver(ap) - - # Evaluate functions - funcs = {} - CFDSolver.evalFunctions(ap, funcs) - - # Store the function values in the output list - CLList.append(funcs[f"{ap.name}_cl"]) - CDList.append(funcs[f"{ap.name}_cd"]) - reslist.append(CFDSolver.getFreeStreamResidual(ap)) - if comm.rank == 0: - print(f"CL: {CLList[-1]}, CD: {CDList[-1]}") - - # rst Print polar - # Print the evaluated functions in a table - if comm.rank == 0: - outputs = [] - for alpha_v, cl, cd, res in zip(alphaList, CLList, CDList, reslist): - print(f"{alpha_v:6.1f} {cl:8.4f} {cd:8.4f}") - outputs.append([mach, reynolds, alpha_v, cl, cd, res]) - # Save the lift and drag coefficients to a file - np.save(os.path.join(output_dir, "M_Re_alpha_CL_CD_res.npy"), outputs) - - MPI.COMM_WORLD.Barrier() - set_exit_code(0) - - except Exception as e: - print(f"Error: {e}") - set_exit_code(1) +# mypy: ignore-errors +"""This file is based on the MACHAero tutorials. + +https://github.com/mdolab/MACH-Aero/blob/main/tutorial/ + +TEMPLATED VARS: +- $mesh_fname: Path to the mesh file. +- $output_dir: Path to the output directory. +- $task: The task to perform: "analysis" or "polar". +- $mach: The Mach number (float). +- $reynolds: The Reynolds number (float). +- $temperature: The temperature (float). +- $alpha: The angle of attack (float). +""" + +import numpy as np +import os +from adflow import ADFLOW +from baseclasses import AeroProblem +from mpi4py import MPI + + +def main() -> None: + mesh_fname = $mesh_fname + output_dir = $output_dir + task = $task + + # mach number + mach = $mach + # Reynolds number + reynolds = $reynolds + # altitude + altitude = $altitude + # temperature + T = $temperature + # Whether to use altitude + use_altitude = $use_altitude + # Reynold's Length + reynoldsLength = 1.0 + + comm = MPI.COMM_WORLD + print(f"Processor {comm.rank} of {comm.size} is running") + if not os.path.exists(output_dir): + if comm.rank == 0: + os.mkdir(output_dir) + + # rst ADflow options + aeroOptions = { + # I/O Parameters + "gridFile": mesh_fname, + "outputDirectory": output_dir, + "monitorvariables": ["cl", "cd","resrho","resrhoe"], + "writeTecplotSurfaceSolution": True, + # Physics Parameters + "equationType": "RANS", + "smoother": "DADI", + "rkreset": True, + "nrkreset": 10, + # Solver Parameters + "MGCycle": "sg", + # ANK Solver Parameters + "useANKSolver": True, + "ankswitchtol": 1e-1, + "liftIndex": 2, + "nsubiterturb": 10, + # NK Solver Parameters + "useNKSolver": True, + "NKSwitchTol": 1e-4, + # Termination Criteria + "L2Convergence": 1e-9, + "L2ConvergenceCoarse": 1e-4, + "nCycles": 5000, + } + print("ADflow options:") + # rst Start ADflow + # Create solver + CFDSolver = ADFLOW(options=aeroOptions) + + # Add features + # CFDSolver.addLiftDistribution(150, "z") + span = 1.0 + pos = np.array([0.5]) * span + CFDSolver.addSlices("z", pos, sliceType="absolute") + + # rst Create AeroProblem + alpha = $alpha + + if use_altitude: + ap = AeroProblem( + name="fc", + alpha=alpha, + mach=mach, + altitude=altitude, + areaRef=1.0, + chordRef=1.0, + evalFuncs=["cl", "cd"], + ) + else: + ap = AeroProblem( + name="fc", + alpha=alpha, + mach=mach, + T=T, + reynolds=reynolds, + reynoldsLength=reynoldsLength, + areaRef=1.0, + chordRef=1.0, + evalFuncs=["cl", "cd"], + ) + + # rst Run ADflow + if task == "analysis": + print("Running analysis") + # Solve + CFDSolver(ap) + # rst Evaluate and print + funcs = {} + CFDSolver.evalFunctions(ap, funcs) + # Print the evaluated functions + if comm.rank == 0: + CL = funcs[f"{ap.name}_cl"] + CD = funcs[f"{ap.name}_cd"] + # Save the lift and drag coefficients to a file + outputs = np.array([mach, reynolds, alpha, CL, CD]) + np.save(os.path.join(output_dir, "outputs.npy"), outputs) + + # rst Create polar arrays + elif task == "polar": + print("Running polar") + # Create an array of alpha values. + # In this case we create 5 random alpha values between 0 and 10 + # from numpy.random import uniform + alphaList = np.linspace(0, 20, 50) + # Sort the alpha values + alphaList.sort() + + # Create storage for the evaluated lift and drag coefficients + CLList = [] + CDList = [] + reslist = [] + # rst Start loop + # Loop over the alpha values and evaluate the polar + for alpha in alphaList: + # rst update AP + # Update the name in the AeroProblem. This allows us to modify the + # output file names with the current alpha. + ap.name = f"fc_{alpha:4.2f}" + + # Update the alpha in aero problem and print it to the screen. + ap.alpha = alpha + if comm.rank == 0: + print(f"current alpha: {ap.alpha}") + + # rst Run ADflow polar + # Solve the flow + CFDSolver(ap) + + # Evaluate functions + funcs = {} + CFDSolver.evalFunctions(ap, funcs) + + # Store the function values in the output list + CLList.append(funcs[f"{ap.name}_cl"]) + CDList.append(funcs[f"{ap.name}_cd"]) + reslist.append(CFDSolver.getFreeStreamResidual(ap)) + if comm.rank == 0: + print(f"CL: {CLList[-1]}, CD: {CDList[-1]}") + + # rst Print polar + # Print the evaluated functions in a table + if comm.rank == 0: + outputs = [] + for alpha_v, cl, cd, res in zip(alphaList, CLList, CDList, reslist): + print(f"{alpha_v:6.1f} {cl:8.4f} {cd:8.4f}") + outputs.append([mach, reynolds, alpha_v, cl, cd, res]) + # Save the lift and drag coefficients to a file + np.save(os.path.join(output_dir, "M_Re_alpha_CL_CD_res.npy"), outputs) + + MPI.COMM_WORLD.Barrier() + +if __name__ == "__main__": + main() diff --git a/engibench/problems/airfoil/templates/airfoil_opt.py b/engibench/problems/airfoil/templates/airfoil_opt.py index a48023e..b2e514d 100644 --- a/engibench/problems/airfoil/templates/airfoil_opt.py +++ b/engibench/problems/airfoil/templates/airfoil_opt.py @@ -31,29 +31,6 @@ from pyoptsparse import OPT from pyoptsparse import Optimization -import signal -import atexit -import sys - -_exit_code = 0 - -def set_exit_code(code): - global _exit_code - _exit_code = code - -def cleanup_handler(signum=None, frame=None): - # Clean termination code for MPI - global _exit_code - try: - from mpi4py import MPI - if MPI.Is_initialized(): - if signum is not None: - _exit_code = 1 - MPI.COMM_WORLD.Abort(_exit_code) - except: - pass - sys.exit(_exit_code) - # ====================================================================== # Functions: # ====================================================================== @@ -94,237 +71,229 @@ def objCon(funcs, printOK): return funcs -if __name__ == "__main__": - - # Register signal handlers - signal.signal(signal.SIGTERM, cleanup_handler) - signal.signal(signal.SIGINT, cleanup_handler) - atexit.register(cleanup_handler) - - try: - # ====================================================================== - # Specify parameters for optimization - # ====================================================================== - # cL constraint - mycl = $cl_target - # angle of attack - alpha = $alpha - # mach number - mach = $mach - # Reynolds number - reynolds = $reynolds - # cruising altitude - altitude = $altitude - # temperature - T = $temperature - # Whether to use altitude - use_altitude = $use_altitude - # Reynold's Length - reynoldsLength = 1.0 - # volume constraint ratio - area_ratio_min = $area_ratio_min - area_initial = $area_initial - area_input_design = $area_input_design - - # Optimization parameters - opt = $opt - opt_options = $opt_options - # ====================================================================== - # Create multipoint communication object - # ====================================================================== - MP = multiPointSparse(MPI.COMM_WORLD) - MP.addProcessorSet("cruise", nMembers=1, memberSizes=MPI.COMM_WORLD.size) - comm, setComm, setFlags, groupFlags, ptID = MP.createCommunicators() - if not os.path.exists($output_dir): - if comm.rank == 0: - os.mkdir($output_dir) - - # ====================================================================== - # ADflow Set-up - # ====================================================================== - aeroOptions = { - # Common Parameters - "gridFile": $mesh_fname, - "outputDirectory": $output_dir, - #"writeSurfaceSolution": False, - "writeVolumeSolution": False, - "writeTecplotSurfaceSolution": True, - "monitorvariables": ["cl", "cd", "yplus"], - # Physics Parameters - "equationType": "RANS", - "smoother": "DADI", - "nCycles": 5000, - "rkreset": True, - "nrkreset": 10, - # NK Options - "useNKSolver": True, - "nkswitchtol": 1e-8, - # ANK Options - "useanksolver": True, - "ankswitchtol": 1e-1, - # "ANKCoupledSwitchTol": 1e-6, - # "ANKSecondOrdSwitchTol": 1e-5, - "liftIndex": 2, - "infchangecorrection": True, - "nsubiterturb": 10, - # Convergence Parameters - "L2Convergence": 1e-8, - "L2ConvergenceCoarse": 1e-4, - # Adjoint Parameters - "adjointSolver": "GMRES", - "adjointL2Convergence": 1e-8, - "ADPC": True, - "adjointMaxIter": 1000, - "adjointSubspaceSize": 200, +def main(): + # ====================================================================== + # Specify parameters for optimization + # ====================================================================== + # cL constraint + mycl = $cl_target + # angle of attack + alpha = $alpha + # mach number + mach = $mach + # Reynolds number + reynolds = $reynolds + # cruising altitude + altitude = $altitude + # temperature + T = $temperature + # Whether to use altitude + use_altitude = $use_altitude + # Reynold's Length + reynoldsLength = 1.0 + # volume constraint ratio + area_ratio_min = $area_ratio_min + area_initial = $area_initial + area_input_design = $area_input_design + + # Optimization parameters + opt = $opt + opt_options = $opt_options + # ====================================================================== + # Create multipoint communication object + # ====================================================================== + MP = multiPointSparse(MPI.COMM_WORLD) + MP.addProcessorSet("cruise", nMembers=1, memberSizes=MPI.COMM_WORLD.size) + comm, setComm, setFlags, groupFlags, ptID = MP.createCommunicators() + if not os.path.exists($output_dir): + if comm.rank == 0: + os.mkdir($output_dir) + + # ====================================================================== + # ADflow Set-up + # ====================================================================== + aeroOptions = { + # Common Parameters + "gridFile": $mesh_fname, + "outputDirectory": $output_dir, + #"writeSurfaceSolution": False, + "writeVolumeSolution": False, + "writeTecplotSurfaceSolution": True, + "monitorvariables": ["cl", "cd", "yplus"], + # Physics Parameters + "equationType": "RANS", + "smoother": "DADI", + "nCycles": 5000, + "rkreset": True, + "nrkreset": 10, + # NK Options + "useNKSolver": True, + "nkswitchtol": 1e-8, + # ANK Options + "useanksolver": True, + "ankswitchtol": 1e-1, + # "ANKCoupledSwitchTol": 1e-6, + # "ANKSecondOrdSwitchTol": 1e-5, + "liftIndex": 2, + "infchangecorrection": True, + "nsubiterturb": 10, + # Convergence Parameters + "L2Convergence": 1e-8, + "L2ConvergenceCoarse": 1e-4, + # Adjoint Parameters + "adjointSolver": "GMRES", + "adjointL2Convergence": 1e-8, + "ADPC": True, + "adjointMaxIter": 1000, + "adjointSubspaceSize": 200, + } + + # Create solver + CFDSolver = ADFLOW(options=aeroOptions, comm=comm) + # ====================================================================== + # Set up flow conditions with AeroProblem + # ====================================================================== + + if use_altitude: + ap = AeroProblem( + name="fc", + alpha=alpha, + mach=mach, + altitude=altitude, + areaRef=1.0, + chordRef=1.0, + evalFuncs=["cl", "cd"], + ) + else: + ap = AeroProblem( + name="fc", + alpha=alpha, + mach=mach, + T=T, + reynolds=reynolds, + reynoldsLength=reynoldsLength, + areaRef=1.0, + chordRef=1.0, + evalFuncs=["cl", "cd"], + ) + + # Add angle of attack variable + ap.addDV("alpha", value=alpha, lower=0.0, upper=10.0, scale=1.0) + # ====================================================================== + # Geometric Design Variable Set-up + # ====================================================================== + # Create DVGeometry object + DVGeo = DVGeometry($ffd_fname) + DVGeo.addLocalDV("shape", lower=-0.025, upper=0.025, axis="y", scale=1.0) + + span = 1.0 + pos = np.array([0.5]) * span + CFDSolver.addSlices("z", pos, sliceType="absolute") + + # Add DVGeo object to CFD solver + CFDSolver.setDVGeo(DVGeo) + # ====================================================================== + # DVConstraint Setup + # ====================================================================== + + DVCon = DVConstraints() + DVCon.setDVGeo(DVGeo) + + # Only ADflow has the getTriangulatedSurface Function + DVCon.setSurface(CFDSolver.getTriangulatedMeshSurface()) + + # Le/Te constraints + lIndex = DVGeo.getLocalIndex(0) + indSetA = [] + indSetB = [] + for k in range(0, 1): + indSetA.append(lIndex[0, 0, k]) # all DV for upper and lower should be same but different sign + indSetB.append(lIndex[0, 1, k]) + for k in range(0, 1): + indSetA.append(lIndex[-1, 0, k]) + indSetB.append(lIndex[-1, 1, k]) + DVCon.addLeTeConstraints(0, indSetA=indSetA, indSetB=indSetB) + + # DV should be same along spanwise + lIndex = DVGeo.getLocalIndex(0) + indSetA = [] + indSetB = [] + for i in range(lIndex.shape[0]): + indSetA.append(lIndex[i, 0, 0]) + indSetB.append(lIndex[i, 0, 1]) + for i in range(lIndex.shape[0]): + indSetA.append(lIndex[i, 1, 0]) + indSetB.append(lIndex[i, 1, 1]) + DVCon.addLinearConstraintsShape(indSetA, indSetB, factorA=1.0, factorB=-1.0, lower=0, upper=0) + + le = 0.010001 + leList = [[le, 0, le], [le, 0, 1.0 - le]] + teList = [[1.0 - le, 0, le], [1.0 - le, 0, 1.0 - le]] + + DVCon.addVolumeConstraint(leList, teList, 2, 100, lower=area_ratio_min*area_initial/area_input_design, upper=1.2*area_initial/area_input_design, scaled=True) + DVCon.addThicknessConstraints2D(leList, teList, 2, 100, lower=0.15, upper=3.0) + # Final constraint to keep TE thickness at original or greater + DVCon.addThicknessConstraints1D(ptList=teList, nCon=2, axis=[0, 1, 0], lower=1.0, scaled=True) + + if comm.rank == 0: + fileName = os.path.join($output_dir, "constraints.dat") + DVCon.writeTecplot(fileName) + # ====================================================================== + # Mesh Warping Set-up + # ====================================================================== + meshOptions = {"gridFile": $mesh_fname} + + mesh = USMesh(options=meshOptions, comm=comm) + CFDSolver.setMesh(mesh) + + # ====================================================================== + # Optimization Problem Set-up + # ====================================================================== + # Create optimization problem + optProb = Optimization("opt", MP.obj, comm=MPI.COMM_WORLD) + + # Add objective + optProb.addObj("obj", scale=1e4) + + # Add variables from the AeroProblem + ap.addVariablesPyOpt(optProb) + + # Add DVGeo variables + DVGeo.addVariablesPyOpt(optProb) + + # Add constraints + DVCon.addConstraintsPyOpt(optProb) + optProb.addCon("cl_con_" + ap.name, lower=0.0, upper=0.0, scale=1.0) + + # The MP object needs the 'obj' and 'sens' function for each proc set, + # the optimization problem and what the objcon function is: + MP.setProcSetObjFunc("cruise", cruiseFuncs) + MP.setProcSetSensFunc("cruise", cruiseFuncsSens) + MP.setObjCon(objCon) + MP.setOptProb(optProb) + optProb.printSparsity() + # Set up optimizer + if opt == "SLSQP": + optOptions = {"IFILE": os.path.join($output_dir, "SLSQP.out")} + elif opt == "SNOPT": + optOptions = { + "Major feasibility tolerance": 1e-4, + "Major optimality tolerance": 1e-4, + "Hessian full memory": None, + "Function precision": 1e-8, + "Print file": os.path.join($output_dir, "SNOPT_print.out"), + "Summary file": os.path.join($output_dir, "SNOPT_summary.out"), } + optOptions.update(opt_options) + opt = OPT(opt, options=optOptions) - # Create solver - CFDSolver = ADFLOW(options=aeroOptions, comm=comm) - # ====================================================================== - # Set up flow conditions with AeroProblem - # ====================================================================== - - if use_altitude: - ap = AeroProblem( - name="fc", - alpha=alpha, - mach=mach, - altitude=altitude, - areaRef=1.0, - chordRef=1.0, - evalFuncs=["cl", "cd"], - ) - else: - ap = AeroProblem( - name="fc", - alpha=alpha, - mach=mach, - T=T, - reynolds=reynolds, - reynoldsLength=reynoldsLength, - areaRef=1.0, - chordRef=1.0, - evalFuncs=["cl", "cd"], - ) - - # Add angle of attack variable - ap.addDV("alpha", value=alpha, lower=0.0, upper=10.0, scale=1.0) - # ====================================================================== - # Geometric Design Variable Set-up - # ====================================================================== - # Create DVGeometry object - DVGeo = DVGeometry($ffd_fname) - DVGeo.addLocalDV("shape", lower=-0.025, upper=0.025, axis="y", scale=1.0) - - span = 1.0 - pos = np.array([0.5]) * span - CFDSolver.addSlices("z", pos, sliceType="absolute") - - # Add DVGeo object to CFD solver - CFDSolver.setDVGeo(DVGeo) - # ====================================================================== - # DVConstraint Setup - # ====================================================================== - - DVCon = DVConstraints() - DVCon.setDVGeo(DVGeo) - - # Only ADflow has the getTriangulatedSurface Function - DVCon.setSurface(CFDSolver.getTriangulatedMeshSurface()) - - # Le/Te constraints - lIndex = DVGeo.getLocalIndex(0) - indSetA = [] - indSetB = [] - for k in range(0, 1): - indSetA.append(lIndex[0, 0, k]) # all DV for upper and lower should be same but different sign - indSetB.append(lIndex[0, 1, k]) - for k in range(0, 1): - indSetA.append(lIndex[-1, 0, k]) - indSetB.append(lIndex[-1, 1, k]) - DVCon.addLeTeConstraints(0, indSetA=indSetA, indSetB=indSetB) - - # DV should be same along spanwise - lIndex = DVGeo.getLocalIndex(0) - indSetA = [] - indSetB = [] - for i in range(lIndex.shape[0]): - indSetA.append(lIndex[i, 0, 0]) - indSetB.append(lIndex[i, 0, 1]) - for i in range(lIndex.shape[0]): - indSetA.append(lIndex[i, 1, 0]) - indSetB.append(lIndex[i, 1, 1]) - DVCon.addLinearConstraintsShape(indSetA, indSetB, factorA=1.0, factorB=-1.0, lower=0, upper=0) - - le = 0.010001 - leList = [[le, 0, le], [le, 0, 1.0 - le]] - teList = [[1.0 - le, 0, le], [1.0 - le, 0, 1.0 - le]] - - DVCon.addVolumeConstraint(leList, teList, 2, 100, lower=area_ratio_min*area_initial/area_input_design, upper=1.2*area_initial/area_input_design, scaled=True) - DVCon.addThicknessConstraints2D(leList, teList, 2, 100, lower=0.15, upper=3.0) - # Final constraint to keep TE thickness at original or greater - DVCon.addThicknessConstraints1D(ptList=teList, nCon=2, axis=[0, 1, 0], lower=1.0, scaled=True) + # Run Optimization + sol = opt(optProb, MP.sens, sensMode='pgc', sensStep=1e-6, storeHistory=os.path.join($output_dir, "opt.hst")) + if MPI.COMM_WORLD.rank == 0: + print(sol) - if comm.rank == 0: - fileName = os.path.join($output_dir, "constraints.dat") - DVCon.writeTecplot(fileName) - # ====================================================================== - # Mesh Warping Set-up - # ====================================================================== - meshOptions = {"gridFile": $mesh_fname} - - mesh = USMesh(options=meshOptions, comm=comm) - CFDSolver.setMesh(mesh) - - # ====================================================================== - # Optimization Problem Set-up - # ====================================================================== - # Create optimization problem - optProb = Optimization("opt", MP.obj, comm=MPI.COMM_WORLD) - - # Add objective - optProb.addObj("obj", scale=1e4) - - # Add variables from the AeroProblem - ap.addVariablesPyOpt(optProb) - - # Add DVGeo variables - DVGeo.addVariablesPyOpt(optProb) - - # Add constraints - DVCon.addConstraintsPyOpt(optProb) - optProb.addCon("cl_con_" + ap.name, lower=0.0, upper=0.0, scale=1.0) - - # The MP object needs the 'obj' and 'sens' function for each proc set, - # the optimization problem and what the objcon function is: - MP.setProcSetObjFunc("cruise", cruiseFuncs) - MP.setProcSetSensFunc("cruise", cruiseFuncsSens) - MP.setObjCon(objCon) - MP.setOptProb(optProb) - optProb.printSparsity() - # Set up optimizer - if opt == "SLSQP": - optOptions = {"IFILE": os.path.join($output_dir, "SLSQP.out")} - elif opt == "SNOPT": - optOptions = { - "Major feasibility tolerance": 1e-4, - "Major optimality tolerance": 1e-4, - "Hessian full memory": None, - "Function precision": 1e-8, - "Print file": os.path.join($output_dir, "SNOPT_print.out"), - "Summary file": os.path.join($output_dir, "SNOPT_summary.out"), - } - optOptions.update(opt_options) - opt = OPT(opt, options=optOptions) - - # Run Optimization - sol = opt(optProb, MP.sens, sensMode='pgc', sensStep=1e-6, storeHistory=os.path.join($output_dir, "opt.hst")) - if MPI.COMM_WORLD.rank == 0: - print(sol) - - MPI.COMM_WORLD.Barrier() - set_exit_code(0) - - except Exception as e: - print(f"Error: {e}") - set_exit_code(1) + MPI.COMM_WORLD.Barrier() + + +if __name__ == "__main__": + main() diff --git a/engibench/problems/airfoil/v0.py b/engibench/problems/airfoil/v0.py index 6594442..4a110d9 100644 --- a/engibench/problems/airfoil/v0.py +++ b/engibench/problems/airfoil/v0.py @@ -427,35 +427,15 @@ def __design_to_simulator_input(self, design: DesignType, config: dict[str, Any] # Launches a docker container with the pre_process.py script # The script generates the mesh and FFD files - try: - bash_command = ( - "source /home/mdolabuser/.bashrc_mdolab && cd /home/mdolabuser/mount/engibench && python " - + self.__docker_study_dir - + "/pre_process.py" - + " > " - + self.__docker_study_dir - + "/output_preprocess.log" - ) - assert self.container_id is not None, "Container ID is not set" - container.run( - command=["/bin/bash", "-c", bash_command], - image=self.container_id, - name="machaero", - mounts=[(self.__local_base_directory, self.__docker_base_dir)], - sync_uid=True, - ) - - except Exception as e: - # Verify output files exist - mesh_file = self.__local_study_dir + "/" + filename + ".cgns" - ffd_file = self.__local_study_dir + "/" + filename + "_ffd.xyz" - msg = "" - - if not os.path.exists(mesh_file): - msg += f"Mesh file not generated: {mesh_file}." - if not os.path.exists(ffd_file): - msg += f"FFD file not generated: {ffd_file}." - raise RuntimeError(f"Pre-processing failed: {e!s}. {msg} Check logs in {self.__local_study_dir}") from e + bash_command = f"source /home/mdolabuser/.bashrc_mdolab && cd {self.__docker_base_dir} && python {self.__docker_study_dir}/pre_process.py" + assert self.container_id is not None, "Container ID is not set" + container.run( + command=["/bin/bash", "-c", bash_command], + image=self.container_id, + name="machaero", + mounts=[(self.__local_base_directory, self.__docker_base_dir)], + sync_uid=True, + ) return filename @@ -542,28 +522,15 @@ def simulate(self, design: DesignType, config: dict[str, Any] | None = None, mpi # Launches a docker container with the airfoil_analysis.py script # The script takes a mesh and ffd and performs an optimization - try: - bash_command = ( - "source /home/mdolabuser/.bashrc_mdolab && cd /home/mdolabuser/mount/engibench && mpirun -np " - + str(mpicores) - + " python " - + self.__docker_study_dir - + "/airfoil_analysis.py > " - + self.__docker_study_dir - + "/output.log" - ) - assert self.container_id is not None, "Container ID is not set" - container.run( - command=["/bin/bash", "-c", bash_command], - image=self.container_id, - name="machaero", - mounts=[(self.__local_base_directory, self.__docker_base_dir)], - sync_uid=True, - ) - except Exception as e: - raise RuntimeError( - f"Failed to run airfoil analysis: {e!s}. Please check logs in {self.__local_study_dir}." - ) from e + bash_command = f"source /home/mdolabuser/.bashrc_mdolab && cd {self.__docker_base_dir} && mpirun -np {mpicores} python -m mpi4py {self.__docker_study_dir}/airfoil_analysis.py" + assert self.container_id is not None, "Container ID is not set" + container.run( + command=["/bin/bash", "-c", bash_command], + image=self.container_id, + name="machaero", + mounts=[(self.__local_base_directory, self.__docker_base_dir)], + sync_uid=True, + ) outputs = np.load(self.__local_study_dir + "/output/outputs.npy") lift = float(outputs[3]) @@ -615,28 +582,17 @@ def optimize( base_config, ) - try: - # Launches a docker container with the optimize_airfoil.py script - # The script takes a mesh and ffd and performs an optimization - bash_command = ( - "source /home/mdolabuser/.bashrc_mdolab && cd /home/mdolabuser/mount/engibench && mpirun -np " - + str(mpicores) - + " python " - + self.__docker_study_dir - + "/airfoil_opt.py > " - + self.__docker_study_dir - + "/airfoil_opt.log" - ) - assert self.container_id is not None, "Container ID is not set" - container.run( - command=["/bin/bash", "-c", bash_command], - image=self.container_id, - name="machaero", - mounts=[(self.__local_base_directory, self.__docker_base_dir)], - sync_uid=True, - ) - except Exception as e: - raise RuntimeError(f"Optimization failed: {e!s}. Check logs in {self.__local_study_dir}") from e + # Launches a docker container with the optimize_airfoil.py script + # The script takes a mesh and ffd and performs an optimization + bash_command = f"source /home/mdolabuser/.bashrc_mdolab && cd {self.__docker_base_dir} && mpirun -np {mpicores} python -m mpi4py {self.__docker_study_dir}/airfoil_opt.py" + assert self.container_id is not None, "Container ID is not set" + container.run( + command=["/bin/bash", "-c", bash_command], + image=self.container_id, + name="machaero", + mounts=[(self.__local_base_directory, self.__docker_base_dir)], + sync_uid=True, + ) # post process -- extract the shape and objective values optisteps_history = [] From 1df0e978c64b7720b625127768893b3ff8096e83 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Gerhard=20Br=C3=A4unlich?= Date: Fri, 17 Oct 2025 16:05:36 +0200 Subject: [PATCH 4/9] airfoil: Pass arguments to container scripts instead of using templating --- .../problems/airfoil/templates/__init__.py | 0 .../airfoil/templates/airfoil_analysis.py | 110 +++---- .../problems/airfoil/templates/airfoil_opt.py | 269 +++++++++--------- .../airfoil/templates/cli_interface.py | 129 +++++++++ .../problems/airfoil/templates/pre_process.py | 61 ++-- engibench/problems/airfoil/v0.py | 146 +++++----- pyproject.toml | 6 +- 7 files changed, 408 insertions(+), 313 deletions(-) create mode 100644 engibench/problems/airfoil/templates/__init__.py create mode 100644 engibench/problems/airfoil/templates/cli_interface.py diff --git a/engibench/problems/airfoil/templates/__init__.py b/engibench/problems/airfoil/templates/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/engibench/problems/airfoil/templates/airfoil_analysis.py b/engibench/problems/airfoil/templates/airfoil_analysis.py index 418b8b8..684ffbf 100644 --- a/engibench/problems/airfoil/templates/airfoil_analysis.py +++ b/engibench/problems/airfoil/templates/airfoil_analysis.py @@ -1,55 +1,53 @@ -# mypy: ignore-errors """This file is based on the MACHAero tutorials. https://github.com/mdolab/MACH-Aero/blob/main/tutorial/ - -TEMPLATED VARS: -- $mesh_fname: Path to the mesh file. -- $output_dir: Path to the output directory. -- $task: The task to perform: "analysis" or "polar". -- $mach: The Mach number (float). -- $reynolds: The Reynolds number (float). -- $temperature: The temperature (float). -- $alpha: The angle of attack (float). """ -import numpy as np +import itertools +import json import os +import sys +from typing import Any + from adflow import ADFLOW from baseclasses import AeroProblem +from cli_interface import AnalysisParameters +from cli_interface import Task from mpi4py import MPI +import numpy as np -def main() -> None: - mesh_fname = $mesh_fname - output_dir = $output_dir - task = $task +def main() -> None: # noqa: C901, PLR0915 + """Entry point of the script.""" + args = AnalysisParameters.from_dict(json.loads(sys.argv[1])) + mesh_fname = args.mesh_fname + output_dir = args.output_dir + task = args.task # mach number - mach = $mach + mach = args.mach # Reynolds number - reynolds = $reynolds + reynolds = args.reynolds # altitude - altitude = $altitude + altitude = args.altitude # temperature - T = $temperature + temperature = args.temperature # Whether to use altitude - use_altitude = $use_altitude + use_altitude = args.use_altitude # Reynold's Length - reynoldsLength = 1.0 + reynolds_length = 1.0 comm = MPI.COMM_WORLD print(f"Processor {comm.rank} of {comm.size} is running") - if not os.path.exists(output_dir): - if comm.rank == 0: - os.mkdir(output_dir) + if not os.path.exists(output_dir) and comm.rank == 0: + os.mkdir(output_dir) # rst ADflow options - aeroOptions = { + aero_options = { # I/O Parameters "gridFile": mesh_fname, "outputDirectory": output_dir, - "monitorvariables": ["cl", "cd","resrho","resrhoe"], + "monitorvariables": ["cl", "cd", "resrho", "resrhoe"], "writeTecplotSurfaceSolution": True, # Physics Parameters "equationType": "RANS", @@ -74,16 +72,15 @@ def main() -> None: print("ADflow options:") # rst Start ADflow # Create solver - CFDSolver = ADFLOW(options=aeroOptions) + cfd_solver = ADFLOW(options=aero_options) # Add features - # CFDSolver.addLiftDistribution(150, "z") span = 1.0 pos = np.array([0.5]) * span - CFDSolver.addSlices("z", pos, sliceType="absolute") + cfd_solver.addSlices("z", pos, sliceType="absolute") # rst Create AeroProblem - alpha = $alpha + alpha = args.alpha if use_altitude: ap = AeroProblem( @@ -100,47 +97,46 @@ def main() -> None: name="fc", alpha=alpha, mach=mach, - T=T, + T=temperature, reynolds=reynolds, - reynoldsLength=reynoldsLength, + reynoldsLength=reynolds_length, areaRef=1.0, chordRef=1.0, evalFuncs=["cl", "cd"], ) # rst Run ADflow - if task == "analysis": + if task == Task.ANALYSIS: print("Running analysis") # Solve - CFDSolver(ap) + cfd_solver(ap) # rst Evaluate and print - funcs = {} - CFDSolver.evalFunctions(ap, funcs) + funcs: dict[str, Any] = {} + cfd_solver.evalFunctions(ap, funcs) # Print the evaluated functions if comm.rank == 0: - CL = funcs[f"{ap.name}_cl"] - CD = funcs[f"{ap.name}_cd"] + cl = funcs[f"{ap.name}_cl"] + cd = funcs[f"{ap.name}_cd"] # Save the lift and drag coefficients to a file - outputs = np.array([mach, reynolds, alpha, CL, CD]) + outputs = np.array([mach, reynolds, alpha, cl, cd]) np.save(os.path.join(output_dir, "outputs.npy"), outputs) # rst Create polar arrays - elif task == "polar": + elif task == Task.POLAR: print("Running polar") # Create an array of alpha values. # In this case we create 5 random alpha values between 0 and 10 - # from numpy.random import uniform - alphaList = np.linspace(0, 20, 50) + alphas = np.linspace(0, 20, 50) # Sort the alpha values - alphaList.sort() + alphas.sort() # Create storage for the evaluated lift and drag coefficients - CLList = [] - CDList = [] + cl_list = [] + cd_list = [] reslist = [] # rst Start loop # Loop over the alpha values and evaluate the polar - for alpha in alphaList: + for alpha in alphas: # rst update AP # Update the name in the AeroProblem. This allows us to modify the # output file names with the current alpha. @@ -153,30 +149,34 @@ def main() -> None: # rst Run ADflow polar # Solve the flow - CFDSolver(ap) + cfd_solver(ap) # Evaluate functions funcs = {} - CFDSolver.evalFunctions(ap, funcs) + cfd_solver.evalFunctions(ap, funcs) # Store the function values in the output list - CLList.append(funcs[f"{ap.name}_cl"]) - CDList.append(funcs[f"{ap.name}_cd"]) - reslist.append(CFDSolver.getFreeStreamResidual(ap)) + cl_list.append(funcs[f"{ap.name}_cl"]) + cd_list.append(funcs[f"{ap.name}_cd"]) + reslist.append(cfd_solver.getFreeStreamResidual(ap)) if comm.rank == 0: - print(f"CL: {CLList[-1]}, CD: {CDList[-1]}") + print(f"CL: {cl_list[-1]}, CD: {cd_list[-1]}") # rst Print polar # Print the evaluated functions in a table if comm.rank == 0: - outputs = [] - for alpha_v, cl, cd, res in zip(alphaList, CLList, CDList, reslist): + outputs = np.array( + list( + zip(itertools.repeat(mach), itertools.repeat(reynolds), alphas, cl_list, cd_list, reslist, strict=False) + ) + ) + for *_, alpha_v, cl, cd, _res in outputs: print(f"{alpha_v:6.1f} {cl:8.4f} {cd:8.4f}") - outputs.append([mach, reynolds, alpha_v, cl, cd, res]) # Save the lift and drag coefficients to a file np.save(os.path.join(output_dir, "M_Re_alpha_CL_CD_res.npy"), outputs) MPI.COMM_WORLD.Barrier() + if __name__ == "__main__": main() diff --git a/engibench/problems/airfoil/templates/airfoil_opt.py b/engibench/problems/airfoil/templates/airfoil_opt.py index b2e514d..97a5cfe 100644 --- a/engibench/problems/airfoil/templates/airfoil_opt.py +++ b/engibench/problems/airfoil/templates/airfoil_opt.py @@ -1,27 +1,19 @@ -# mypy: ignore-errors """This file is largely based on the MACHAero tutorials. https://github.com/mdolab/MACH-Aero/blob/main/tutorial/ - -TEMPLATED VARS: -- $cl_target: The lift coefficient constraint (float). -- $alpha: The angle of attack (float). -- $mach: The Mach number (float). -- $altitude: The cruising altitude (int). -- $ffd_fname: Path to the FFD file. -- $mesh_fname: Path to the mesh file. -- $output_dir: Path to the output directory. -- $opt: The optimization algorithm: SLSQP or SNOPT. -- $opt_options: The optimization options (dict). """ # ====================================================================== # Import modules # ====================================================================== +import json import os +import sys from adflow import ADFLOW from baseclasses import AeroProblem +from cli_interface import Algorithm +from cli_interface import OptimizeParameters from idwarp import USMesh from mpi4py import MPI from multipoint import multiPointSparse @@ -31,92 +23,53 @@ from pyoptsparse import OPT from pyoptsparse import Optimization -# ====================================================================== -# Functions: -# ====================================================================== -def cruiseFuncs(x): - if MPI.COMM_WORLD.rank == 0: - print(x) - # Set design vars - DVGeo.setDesignVars(x) - ap.setDesignVars(x) - # Run CFD - CFDSolver(ap) - # Evaluate functions - funcs = {} - DVCon.evalFunctions(funcs) - CFDSolver.evalFunctions(ap, funcs) - CFDSolver.checkSolutionFailure(ap, funcs) - if MPI.COMM_WORLD.rank == 0: - print(funcs) - return funcs - -def cruiseFuncsSens(x, funcs): - funcsSens = {} - DVCon.evalFunctionsSens(funcsSens) - CFDSolver.evalFunctionsSens(ap, funcsSens) - CFDSolver.checkAdjointFailure(ap, funcsSens) - if MPI.COMM_WORLD.rank == 0: - print(funcsSens) - return funcsSens - - -def objCon(funcs, printOK): - # Assemble the objective and any additional constraints: - funcs["obj"] = funcs[ap["cd"]] - funcs["cl_con_" + ap.name] = funcs[ap["cl"]] - mycl - if printOK: - print("funcs in obj:", funcs) - return funcs - - -def main(): +def main() -> None: # noqa: C901, PLR0915 + """Entry point of the script.""" + args = OptimizeParameters.from_dict(json.loads(sys.argv[1])) # ====================================================================== # Specify parameters for optimization # ====================================================================== # cL constraint - mycl = $cl_target + mycl = args.cl_target # angle of attack - alpha = $alpha + alpha = args.alpha # mach number - mach = $mach + mach = args.mach # Reynolds number - reynolds = $reynolds + reynolds = args.reynolds # cruising altitude - altitude = $altitude + altitude = args.altitude # temperature - T = $temperature + temperature = args.temperature # Whether to use altitude - use_altitude = $use_altitude + use_altitude = args.use_altitude # Reynold's Length - reynoldsLength = 1.0 + reynolds_length = 1.0 # volume constraint ratio - area_ratio_min = $area_ratio_min - area_initial = $area_initial - area_input_design = $area_input_design + area_ratio_min = args.area_ratio_min + area_initial = args.area_initial + area_input_design = args.area_input_design # Optimization parameters - opt = $opt - opt_options = $opt_options + opt = args.opt + opt_options = args.opt_options # ====================================================================== # Create multipoint communication object # ====================================================================== - MP = multiPointSparse(MPI.COMM_WORLD) - MP.addProcessorSet("cruise", nMembers=1, memberSizes=MPI.COMM_WORLD.size) - comm, setComm, setFlags, groupFlags, ptID = MP.createCommunicators() - if not os.path.exists($output_dir): - if comm.rank == 0: - os.mkdir($output_dir) + mp = multiPointSparse(MPI.COMM_WORLD) + mp.addProcessorSet("cruise", nMembers=1, memberSizes=MPI.COMM_WORLD.size) + comm, *_ = mp.createCommunicators() + if not os.path.exists(args.output_dir) and comm.rank == 0: + os.mkdir(args.output_dir) # ====================================================================== # ADflow Set-up # ====================================================================== - aeroOptions = { + aero_options = { # Common Parameters - "gridFile": $mesh_fname, - "outputDirectory": $output_dir, - #"writeSurfaceSolution": False, + "gridFile": args.mesh_fname, + "outputDirectory": args.output_dir, "writeVolumeSolution": False, "writeTecplotSurfaceSolution": True, "monitorvariables": ["cl", "cd", "yplus"], @@ -132,8 +85,6 @@ def main(): # ANK Options "useanksolver": True, "ankswitchtol": 1e-1, - # "ANKCoupledSwitchTol": 1e-6, - # "ANKSecondOrdSwitchTol": 1e-5, "liftIndex": 2, "infchangecorrection": True, "nsubiterturb": 10, @@ -149,7 +100,7 @@ def main(): } # Create solver - CFDSolver = ADFLOW(options=aeroOptions, comm=comm) + cfd_solver = ADFLOW(options=aero_options, comm=comm) # ====================================================================== # Set up flow conditions with AeroProblem # ====================================================================== @@ -169,9 +120,9 @@ def main(): name="fc", alpha=alpha, mach=mach, - T=T, + T=temperature, reynolds=reynolds, - reynoldsLength=reynoldsLength, + reynoldsLength=reynolds_length, areaRef=1.0, chordRef=1.0, evalFuncs=["cl", "cd"], @@ -183,112 +134,154 @@ def main(): # Geometric Design Variable Set-up # ====================================================================== # Create DVGeometry object - DVGeo = DVGeometry($ffd_fname) - DVGeo.addLocalDV("shape", lower=-0.025, upper=0.025, axis="y", scale=1.0) + dv_geo = DVGeometry(args.ffd_fname) + dv_geo.addLocalDV("shape", lower=-0.025, upper=0.025, axis="y", scale=1.0) span = 1.0 pos = np.array([0.5]) * span - CFDSolver.addSlices("z", pos, sliceType="absolute") + cfd_solver.addSlices("z", pos, sliceType="absolute") # Add DVGeo object to CFD solver - CFDSolver.setDVGeo(DVGeo) + cfd_solver.setDVGeo(dv_geo) # ====================================================================== # DVConstraint Setup # ====================================================================== - DVCon = DVConstraints() - DVCon.setDVGeo(DVGeo) + dv_con = DVConstraints() + dv_con.setDVGeo(dv_geo) # Only ADflow has the getTriangulatedSurface Function - DVCon.setSurface(CFDSolver.getTriangulatedMeshSurface()) + dv_con.setSurface(cfd_solver.getTriangulatedMeshSurface()) # Le/Te constraints - lIndex = DVGeo.getLocalIndex(0) - indSetA = [] - indSetB = [] - for k in range(0, 1): - indSetA.append(lIndex[0, 0, k]) # all DV for upper and lower should be same but different sign - indSetB.append(lIndex[0, 1, k]) - for k in range(0, 1): - indSetA.append(lIndex[-1, 0, k]) - indSetB.append(lIndex[-1, 1, k]) - DVCon.addLeTeConstraints(0, indSetA=indSetA, indSetB=indSetB) + l_index = dv_geo.getLocalIndex(0) + ind_set_a = [] + ind_set_b = [] + for k in range(1): + ind_set_a.append(l_index[0, 0, k]) # all DV for upper and lower should be same but different sign + ind_set_b.append(l_index[0, 1, k]) + for k in range(1): + ind_set_a.append(l_index[-1, 0, k]) + ind_set_b.append(l_index[-1, 1, k]) + dv_con.addLeTeConstraints(0, indSetA=ind_set_a, indSetB=ind_set_b) # DV should be same along spanwise - lIndex = DVGeo.getLocalIndex(0) - indSetA = [] - indSetB = [] - for i in range(lIndex.shape[0]): - indSetA.append(lIndex[i, 0, 0]) - indSetB.append(lIndex[i, 0, 1]) - for i in range(lIndex.shape[0]): - indSetA.append(lIndex[i, 1, 0]) - indSetB.append(lIndex[i, 1, 1]) - DVCon.addLinearConstraintsShape(indSetA, indSetB, factorA=1.0, factorB=-1.0, lower=0, upper=0) + l_index = dv_geo.getLocalIndex(0) + ind_set_a = [] + ind_set_b = [] + for i in range(l_index.shape[0]): + ind_set_a.append(l_index[i, 0, 0]) + ind_set_b.append(l_index[i, 0, 1]) + for i in range(l_index.shape[0]): + ind_set_a.append(l_index[i, 1, 0]) + ind_set_b.append(l_index[i, 1, 1]) + dv_con.addLinearConstraintsShape(ind_set_a, ind_set_b, factorA=1.0, factorB=-1.0, lower=0, upper=0) le = 0.010001 - leList = [[le, 0, le], [le, 0, 1.0 - le]] - teList = [[1.0 - le, 0, le], [1.0 - le, 0, 1.0 - le]] - - DVCon.addVolumeConstraint(leList, teList, 2, 100, lower=area_ratio_min*area_initial/area_input_design, upper=1.2*area_initial/area_input_design, scaled=True) - DVCon.addThicknessConstraints2D(leList, teList, 2, 100, lower=0.15, upper=3.0) + le_list = [[le, 0, le], [le, 0, 1.0 - le]] + te_list = [[1.0 - le, 0, le], [1.0 - le, 0, 1.0 - le]] + + dv_con.addVolumeConstraint( + le_list, + te_list, + 2, + 100, + lower=area_ratio_min * area_initial / area_input_design, + upper=1.2 * area_initial / area_input_design, + scaled=True, + ) + dv_con.addThicknessConstraints2D(le_list, te_list, 2, 100, lower=0.15, upper=3.0) # Final constraint to keep TE thickness at original or greater - DVCon.addThicknessConstraints1D(ptList=teList, nCon=2, axis=[0, 1, 0], lower=1.0, scaled=True) + dv_con.addThicknessConstraints1D(ptList=te_list, nCon=2, axis=[0, 1, 0], lower=1.0, scaled=True) if comm.rank == 0: - fileName = os.path.join($output_dir, "constraints.dat") - DVCon.writeTecplot(fileName) + file_name = os.path.join(args.output_dir, "constraints.dat") + dv_con.writeTecplot(file_name) # ====================================================================== # Mesh Warping Set-up # ====================================================================== - meshOptions = {"gridFile": $mesh_fname} + mesh_options = {"gridFile": args.mesh_fname} - mesh = USMesh(options=meshOptions, comm=comm) - CFDSolver.setMesh(mesh) + mesh = USMesh(options=mesh_options, comm=comm) + cfd_solver.setMesh(mesh) # ====================================================================== # Optimization Problem Set-up # ====================================================================== # Create optimization problem - optProb = Optimization("opt", MP.obj, comm=MPI.COMM_WORLD) + opt_prob = Optimization("opt", mp.obj, comm=MPI.COMM_WORLD) # Add objective - optProb.addObj("obj", scale=1e4) + opt_prob.addObj("obj", scale=1e4) # Add variables from the AeroProblem - ap.addVariablesPyOpt(optProb) + ap.addVariablesPyOpt(opt_prob) # Add DVGeo variables - DVGeo.addVariablesPyOpt(optProb) + dv_geo.addVariablesPyOpt(opt_prob) # Add constraints - DVCon.addConstraintsPyOpt(optProb) - optProb.addCon("cl_con_" + ap.name, lower=0.0, upper=0.0, scale=1.0) + dv_con.addConstraintsPyOpt(opt_prob) + opt_prob.addCon("cl_con_" + ap.name, lower=0.0, upper=0.0, scale=1.0) # The MP object needs the 'obj' and 'sens' function for each proc set, # the optimization problem and what the objcon function is: - MP.setProcSetObjFunc("cruise", cruiseFuncs) - MP.setProcSetSensFunc("cruise", cruiseFuncsSens) - MP.setObjCon(objCon) - MP.setOptProb(optProb) - optProb.printSparsity() + def cruise_funcs(x): + if MPI.COMM_WORLD.rank == 0: + print(x) + # Set design vars + dv_geo.setDesignVars(x) + ap.setDesignVars(x) + # Run CFD + cfd_solver(ap) + # Evaluate functions + funcs = {} + dv_con.evalFunctions(funcs) + cfd_solver.evalFunctions(ap, funcs) + cfd_solver.checkSolutionFailure(ap, funcs) + if MPI.COMM_WORLD.rank == 0: + print(funcs) + return funcs + + def cruise_funcs_sens(_x, _funcs): + funcs_sens = {} + dv_con.evalFunctionsSens(funcs_sens) + cfd_solver.evalFunctionsSens(ap, funcs_sens) + cfd_solver.checkAdjointFailure(ap, funcs_sens) + if MPI.COMM_WORLD.rank == 0: + print(funcs_sens) + return funcs_sens + + def obj_con(funcs, print_ok): + # Assemble the objective and any additional constraints: + funcs["obj"] = funcs[ap["cd"]] + funcs["cl_con_" + ap.name] = funcs[ap["cl"]] - mycl + if print_ok: + print("funcs in obj:", funcs) + return funcs + + mp.setProcSetObjFunc("cruise", cruise_funcs) + mp.setProcSetSensFunc("cruise", cruise_funcs_sens) + mp.setObjCon(obj_con) + mp.setOptProb(opt_prob) + opt_prob.printSparsity() # Set up optimizer - if opt == "SLSQP": - optOptions = {"IFILE": os.path.join($output_dir, "SLSQP.out")} - elif opt == "SNOPT": - optOptions = { + if opt == Algorithm.SLSQP: + opt_options = {"IFILE": os.path.join(args.output_dir, "SLSQP.out")} + elif opt == Algorithm.SNOPT: + opt_options = { "Major feasibility tolerance": 1e-4, "Major optimality tolerance": 1e-4, "Hessian full memory": None, "Function precision": 1e-8, - "Print file": os.path.join($output_dir, "SNOPT_print.out"), - "Summary file": os.path.join($output_dir, "SNOPT_summary.out"), + "Print file": os.path.join(args.output_dir, "SNOPT_print.out"), + "Summary file": os.path.join(args.output_dir, "SNOPT_summary.out"), } - optOptions.update(opt_options) - opt = OPT(opt, options=optOptions) + opt_options.update(opt_options) + opt = OPT(opt.name, options=opt_options) # Run Optimization - sol = opt(optProb, MP.sens, sensMode='pgc', sensStep=1e-6, storeHistory=os.path.join($output_dir, "opt.hst")) + sol = opt(opt_prob, mp.sens, sensMode="pgc", sensStep=1e-6, storeHistory=os.path.join(args.output_dir, "opt.hst")) if MPI.COMM_WORLD.rank == 0: print(sol) diff --git a/engibench/problems/airfoil/templates/cli_interface.py b/engibench/problems/airfoil/templates/cli_interface.py new file mode 100644 index 0000000..7372098 --- /dev/null +++ b/engibench/problems/airfoil/templates/cli_interface.py @@ -0,0 +1,129 @@ +"""Dataclasses sent from the main script to scripts inside the airfoil container.""" + +from dataclasses import asdict +from dataclasses import dataclass +from enum import auto +from enum import Enum +from typing import Any + + +class Task(Enum): + """The task to perform by analysis.""" + + ANALYSIS = auto() + POLAR = auto() + + +@dataclass +class AnalysisParameters: + """Parameters for airfoil_analyse.py.""" + + mesh_fname: str + """Path to the mesh file""" + output_dir: str + """Path to the output directory""" + task: Task + """The task to perform: "analysis" or "polar""" + mach: float + """mach number""" + reynolds: float + """Reynolds number""" + altitude: float + """altitude""" + temperature: float + """temperature""" + use_altitude: bool + """Whether to use altitude""" + alpha: float + """Angle of attack""" + + def to_dict(self) -> dict[str, Any]: + """Serialize to python primitives.""" + return {**asdict(self), "task": self.task.name} + + @classmethod + def from_dict(cls, data: dict[str, Any]) -> "AnalysisParameters": + """Deserialize from python primitives.""" + return cls(task=Task[data.pop("task")], **data) + + +class Algorithm(Enum): + """Algorithm to be used by optimize.""" + + SLSQP = auto() + SNOPT = auto() + + +@dataclass +class OptimizeParameters: + """Parameters for airfoil_opt.py.""" + + cl_target: float + """The lift coefficient constraint""" + alpha: float + """The angle of attack""" + mach: float + """The Mach number""" + reynolds: float + """Reynolds number""" + temperature: float + """Temerature""" + altitude: int + """The cruising altitude""" + area_ratio_min: float + area_initial: float + area_input_design: float + ffd_fname: str + """Path to the FFD file""" + mesh_fname: str + """Path to the mesh file""" + output_dir: str + """Path to the output directory""" + opt: Algorithm + """The optimization algorithm: SLSQP or SNOPT""" + opt_options: dict[str, Any] + """The optimization options""" + use_altitude: bool + """Whether to use altitude""" + + def to_dict(self) -> dict[str, Any]: + """Serialize to python primitives.""" + return {**asdict(self), "opt": self.opt.name} + + @classmethod + def from_dict(cls, data: dict[str, Any]) -> "OptimizeParameters": + """Deserialize from python primitives.""" + return cls(opt=Algorithm[data.pop("opt")], **data) + + +@dataclass +class PreprocessParameters: + """Parameters for pre_process.py.""" + + design_fname: str + """Path to the design file""" + N_sample: int + """Number of points to sample on the airfoil surface. Defines part of the mesh resolution""" + n_tept_s: int + """Number of points on the trailing edge""" + x_cut: float + """Blunt edge dimensionless cut location""" + tmp_xyz_fname: str + """Path to the temporary xyz file""" + mesh_fname: str + """Path to the generated mesh file""" + ffd_fname: str + """Path to the generated FFD file""" + ffd_ymarginu: float + """Upper (y-axis) margin for the fitted FFD cage""" + ffd_ymarginl: float + """Lower (y-axis) margin for the fitted FFD cage""" + ffd_pts: int + """Number of FFD points""" + N_grid: int + """Number of grid levels to march from the airfoil surface. Defines part of the mesh resolution""" + s0: float + """Off-the-wall spacing for the purpose of modeling the boundary layer""" + input_blunted: bool + march_dist: float + """Distance to march the grid from the airfoil surface""" diff --git a/engibench/problems/airfoil/templates/pre_process.py b/engibench/problems/airfoil/templates/pre_process.py index 3bdd541..ae0a978 100644 --- a/engibench/problems/airfoil/templates/pre_process.py +++ b/engibench/problems/airfoil/templates/pre_process.py @@ -3,63 +3,53 @@ https://github.com/mdolab/MACH-Aero/blob/main/tutorial/ -TEMPLATED VARS: -- $design_fname: Path to the design file. -- $N_sample: Number of points to sample on the airfoil surface. Defines part of the mesh resolution. -- $nTEPts: Number of points on the trailing edge. -- $xCut: Blunt edge dimensionless cut location. -- $tmp_xyz_fname: Path to the temporary xyz file. -- $mesh_fname: Path to the generated mesh file. -- $ffd_fname: Path to the generated FFD file. -- $ffd_ymarginu: Upper (y-axis) margin for the fitted FFD cage. -- $ffd_ymarginl: Lower (y-axis) margin for the fitted FFD cage. -- $ffd_pts: Number of FFD points. -- $N_grid: Number of grid levels to march from the airfoil surface. Defines part of the mesh resolution. -- $s0: Off-the-wall spacing for the purpose of modeling the boundary layer. # TODO: Add the automatic grid spacing calculation. -- $marchDist: Distance to march the grid from the airfoil surface. - +TODO: Add the automatic grid spacing calculation. """ -import numpy as np -from pyhyp import pyHyp +import json +import sys + +from cli_interface import PreprocessParameters import prefoil +from pyhyp import pyHyp if __name__ == "__main__": + args = PreprocessParameters(**json.loads(sys.argv[1])) - coords = prefoil.utils.readCoordFile($design_fname) # type: ignore + coords = prefoil.utils.readCoordFile(args.design_fname) airfoil = prefoil.Airfoil(coords) print("Running pre-process.py") - input_blunted = $input_blunted + input_blunted = args.input_blunted if not input_blunted: airfoil.normalizeAirfoil() - airfoil.makeBluntTE(xCut=$xCut) + airfoil.makeBluntTE(xCut=args.x_cut) - N_sample = $N_sample - nTEPts = $nTEPts + N_sample = args.N_sample + n_tept_s = args.n_tept_s coords = airfoil.getSampledPts( N_sample, spacingFunc=prefoil.sampling.conical, func_args={"coeff": 1.2}, - nTEPts=nTEPts, + nTEPts=n_tept_s, ) # Write a fitted FFD with 10 chordwise points - ffd_ymarginu = $ffd_ymarginu - ffd_ymarginl = $ffd_ymarginl - ffd_fname = $ffd_fname - ffd_pts = $ffd_pts + ffd_ymarginu = args.ffd_ymarginu + ffd_ymarginl = args.ffd_ymarginl + ffd_fname = args.ffd_fname + ffd_pts = args.ffd_pts airfoil.generateFFD(ffd_pts, ffd_fname, ymarginu=ffd_ymarginu, ymarginl=ffd_ymarginl) # write out plot3d - airfoil.writeCoords($tmp_xyz_fname, file_format="plot3d") + airfoil.writeCoords(args.tmp_xyz_fname, file_format="plot3d") # GenOptions options = { # --------------------------- # Input Parameters # --------------------------- - "inputFile": $tmp_xyz_fname + ".xyz", + "inputFile": args.tmp_xyz_fname + ".xyz", "unattachedEdgesAreSymmetry": False, "outerFaceBC": "farfield", "autoConnect": True, @@ -68,19 +58,18 @@ # --------------------------- # Grid Parameters # --------------------------- - "N": $N_grid, + "N": args.N_grid, "nConstantStart": 8, - "s0": $s0, - "marchDist": $marchDist, + "s0": args.s0, + "marchDist": args.march_dist, # Smoothing parameters "volSmoothIter": 150, "volCoef": 0.25, - "volBlend": 0.001 - # "volSmoothSchedule": [[0, 0], [0.2, 2], [0.5, 200], [1.0, 1000]], + "volBlend": 0.001, } hyp = pyHyp(options=options) hyp.run() - hyp.writeCGNS($mesh_fname) + hyp.writeCGNS(args.mesh_fname) - print(f"Generated files FFD and mesh in ${ffd_fname}, ${mesh_fname}") + print(f"Generated files FFD and mesh in {ffd_fname}, {args.mesh_fname}") diff --git a/engibench/problems/airfoil/v0.py b/engibench/problems/airfoil/v0.py index 4a110d9..72459a6 100644 --- a/engibench/problems/airfoil/v0.py +++ b/engibench/problems/airfoil/v0.py @@ -19,6 +19,7 @@ from dataclasses import dataclass from dataclasses import field import importlib +import json import os import shutil import sys @@ -38,13 +39,13 @@ from engibench.core import OptiStep from engibench.core import Problem from engibench.problems.airfoil.pyopt_history import History +from engibench.problems.airfoil.templates import cli_interface from engibench.problems.airfoil.utils import calc_area from engibench.problems.airfoil.utils import calc_off_wall_distance from engibench.problems.airfoil.utils import reorder_coords from engibench.problems.airfoil.utils import scale_coords from engibench.utils import container from engibench.utils.files import clone_dir -from engibench.utils.files import replace_template_values # Allow loading pyoptsparse histories even if pyoptsparse is not installed: if importlib.util.find_spec("pyoptsparse") is None: @@ -289,7 +290,9 @@ class Conditions: float, bounded(lower=0.0).category(IMPL), bounded(lower=1e5, upper=1e9).warning().category(IMPL) ] = 1e6 area_initial: float | None = None + """actual initial airfoil area""" area_ratio_min: Annotated[float, bounded(lower=0.0, upper=1.2).category(THEORY)] = 0.7 + """Minimum ratio the initial area is allowed to decrease to i.e minimum_area = area_initial*area_target""" cl_target: float = 0.5 conditions = Conditions() @@ -304,8 +307,8 @@ class Config(Conditions): use_altitude: bool = False output_dir: str | None = None mesh_fname: str | None = None - task: str = "'analysis'" - opt: str = "'SLSQP'" + task: str = "analysis" + opt: str = "SLSQP" opt_options: dict = field(default_factory=dict) ffd_fname: str | None = None area_input_design: float | None = None @@ -364,7 +367,9 @@ def reset(self, seed: int | None = None, *, cleanup: bool = False) -> None: self.__local_study_dir = self.__local_target_dir + "/" + self.current_study self.__docker_study_dir = self.__docker_target_dir + "/" + self.current_study - def __design_to_simulator_input(self, design: DesignType, config: dict[str, Any], filename: str = "design") -> str: + def __design_to_simulator_input( + self, design: DesignType, mach: float, reynolds: float, temperature: float, filename: str = "design" + ) -> str: """Converts a design to a simulator input. The simulator inputs are two files: a mesh file (.cgns) and a FFD file (.xyz). This function generates these files from the design. @@ -372,7 +377,9 @@ def __design_to_simulator_input(self, design: DesignType, config: dict[str, Any] Args: design (dict): The design to convert. - config (dict): A dictionary with configuration (e.g., boundary conditions) for the simulation. + mach: mach number + reynolds: reynolds number + temperature: temperature filename (str): The filename to save the design to. """ # Creates the study directory @@ -382,52 +389,38 @@ def __design_to_simulator_input(self, design: DesignType, config: dict[str, Any] # Calculate the off-the-wall distance estimate_s0 = True - if estimate_s0: - s0 = calc_off_wall_distance( - mach=config["mach"], reynolds=config["reynolds"], freestreamTemp=config["temperature"] - ) - else: - s0 = 1e-5 - base_config = { - "design_fname": f"'{self.__docker_study_dir}/{filename}.dat'", - "tmp_xyz_fname": f"'{tmp}'", - "mesh_fname": "'" + self.__docker_study_dir + "/" + filename + ".cgns'", - "ffd_fname": "'" + self.__docker_study_dir + "/" + filename + "_ffd'", - "marchDist": 100.0, # Distance to march the grid from the airfoil surface - "N_sample": 180, - "nTEPts": 4, - "xCut": 0.99, - "ffd_ymarginu": 0.05, - "ffd_ymarginl": 0.05, - "ffd_pts": 10, - "N_grid": 100, - "estimate_s0": estimate_s0, - "make_input_design_blunt": True, - "input_blunted": False, - "s0": s0, - **dataclasses.asdict(self.Conditions()), - } + s0 = calc_off_wall_distance(mach=mach, reynolds=reynolds, freestreamTemp=temperature) if estimate_s0 else 1e-5 # Scale the design to fit in the design space + x_cut = 0.99 scaled_design, input_blunted = scale_coords( design["coords"], - blunted=bool(base_config["input_blunted"]), - xcut=base_config["xCut"], + blunted=False, + xcut=x_cut, + ) + args = cli_interface.PreprocessParameters( + design_fname=f"{self.__docker_study_dir}/{filename}.dat", + tmp_xyz_fname=tmp, + mesh_fname=self.__docker_study_dir + "/" + filename + ".cgns", + ffd_fname=self.__docker_study_dir + "/" + filename + "_ffd", + N_sample=180, + n_tept_s=4, + x_cut=x_cut, + ffd_ymarginu=0.05, + ffd_ymarginl=0.05, + ffd_pts=10, + N_grid=100, + s0=s0, + input_blunted=input_blunted, + march_dist=100.0, ) - base_config["input_blunted"] = input_blunted # Save the design to a temporary file. Format to 1e-6 rounding np.savetxt(self.__local_study_dir + "/" + filename + ".dat", scaled_design.transpose()) - # Prepares the preprocess.py script with the design - replace_template_values( - self.__local_study_dir + "/pre_process.py", - base_config, - ) - # Launches a docker container with the pre_process.py script # The script generates the mesh and FFD files - bash_command = f"source /home/mdolabuser/.bashrc_mdolab && cd {self.__docker_base_dir} && python {self.__docker_study_dir}/pre_process.py" + bash_command = f"source /home/mdolabuser/.bashrc_mdolab && cd {self.__docker_base_dir} && python {self.__docker_study_dir}/pre_process.py '{json.dumps(dataclasses.asdict(args))}'" assert self.container_id is not None, "Container ID is not set" container.run( command=["/bin/bash", "-c", bash_command], @@ -503,26 +496,24 @@ def simulate(self, design: DesignType, config: dict[str, Any] | None = None, mpi # pre-process the design and run the simulation # Prepares the airfoil_analysis.py script with the simulation configuration - base_config = { - "alpha": design["angle_of_attack"], - "altitude": 10000, - "temperature": 300, - "use_altitude": False, - "output_dir": "'" + self.__docker_study_dir + "/output/'", - "mesh_fname": "'" + self.__docker_study_dir + "/design.cgns'", - "task": "'analysis'", - **dataclasses.asdict(self.Conditions()), - **(config or {}), - } - self.__design_to_simulator_input(design, base_config) - replace_template_values( - self.__local_study_dir + "/airfoil_analysis.py", - base_config, + conditions = self.Conditions() + config = config or {} + args = cli_interface.AnalysisParameters( + alpha=design["angle_of_attack"], + altitude=config.get("altitude", 10000), + temperature=config.get("temperature", 300), + reynolds=config.get("reynolds", conditions.reynolds), + mach=config.get("mach", conditions.mach), + use_altitude=config.get("use_altitude", False), + output_dir=config.get("output_dir", self.__docker_study_dir + "/output/"), + mesh_fname=config.get("mesh_fname", self.__docker_study_dir + "/design.cgns"), + task=cli_interface.Task[config["task"]] if "task" in config else cli_interface.Task.ANALYSIS, ) + self.__design_to_simulator_input(design, mach=args.mach, reynolds=args.reynolds, temperature=args.temperature) # Launches a docker container with the airfoil_analysis.py script # The script takes a mesh and ffd and performs an optimization - bash_command = f"source /home/mdolabuser/.bashrc_mdolab && cd {self.__docker_base_dir} && mpirun -np {mpicores} python -m mpi4py {self.__docker_study_dir}/airfoil_analysis.py" + bash_command = f"source /home/mdolabuser/.bashrc_mdolab && cd {self.__docker_base_dir} && mpirun -np {mpicores} python -m mpi4py {self.__docker_study_dir}/airfoil_analysis.py '{json.dumps(args.to_dict())}'" assert self.container_id is not None, "Container ID is not set" container.run( command=["/bin/bash", "-c", bash_command], @@ -557,34 +548,31 @@ def optimize( filename = "candidate_design" # Prepares the optimize_airfoil.py script with the optimization configuration - base_config = { - "cl_target": 0.5, - "alpha": starting_point["angle_of_attack"], - "mach": 0.75, - "reynolds": 1e6, - "altitude": 10000, - "temperature": 300, # should specify either mach + altitude or mach + reynolds + reynoldsLength (default to 1) + temperature - "use_altitude": False, - "area_initial": None, # actual initial airfoil area - "area_ratio_min": 0.7, # Minimum ratio the initial area is allowed to decrease to i.e minimum_area = area_initial*area_target - "opt": "'SLSQP'", - "opt_options": {}, - "output_dir": "'" + self.__docker_study_dir + "/output/'", - "ffd_fname": "'" + self.__docker_study_dir + "/" + filename + "_ffd.xyz'", - "mesh_fname": "'" + self.__docker_study_dir + "/" + filename + ".cgns'", - "area_input_design": calc_area(starting_point["coords"]), + fields = {f.name for f in dataclasses.fields(cli_interface.OptimizeParameters)} + config = {key: val for key, val in (config or {}).items() if key in fields} + if "opt" in config: + config["opt"] = cli_interface.Algorithm[config["opt"]] + args = cli_interface.OptimizeParameters( **dataclasses.asdict(self.Conditions()), - **(config or {}), - } - self.__design_to_simulator_input(starting_point, base_config, filename) - replace_template_values( - self.__local_study_dir + "/airfoil_opt.py", - base_config, + alpha=starting_point["angle_of_attack"], + altitude=10000, + temperature=300, # should specify either mach + altitude or mach + reynolds + reynoldsLength (default to 1) + temperature + use_altitude=False, + opt=cli_interface.Algorithm.SLSQP, + opt_options={}, + output_dir=self.__docker_study_dir + "/output", + ffd_fname=self.__docker_study_dir + "/" + filename + "_ffd.xyz", + mesh_fname=self.__docker_study_dir + "/" + filename + ".cgns", + area_input_design=calc_area(starting_point["coords"]), + **config, + ) + self.__design_to_simulator_input( + starting_point, reynolds=args.reynolds, mach=args.reynolds, temperature=args.temperature, filename=filename ) # Launches a docker container with the optimize_airfoil.py script # The script takes a mesh and ffd and performs an optimization - bash_command = f"source /home/mdolabuser/.bashrc_mdolab && cd {self.__docker_base_dir} && mpirun -np {mpicores} python -m mpi4py {self.__docker_study_dir}/airfoil_opt.py" + bash_command = f"source /home/mdolabuser/.bashrc_mdolab && cd {self.__docker_base_dir} && mpirun -np {mpicores} python -m mpi4py {self.__docker_study_dir}/airfoil_opt.py '{json.dumps(args.to_dict())}'" assert self.container_id is not None, "Container ID is not set" container.run( command=["/bin/bash", "-c", bash_command], diff --git a/pyproject.toml b/pyproject.toml index 0cf7c71..ebaad33 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -106,7 +106,6 @@ exclude = [ "dist", "node_modules", "site-packages", - "templates", "venv", "docs", ] @@ -190,7 +189,7 @@ order-by-type = false "TRY", # tryceratops ] -"**/templates/**/*.py" = ["ALL"] +"**/heatconduction*/templates/**/*.py" = ["ALL"] ###################################### FORMAT ######################################## @@ -244,9 +243,6 @@ warn_unused_ignores = true files = ["engibench/", "tests/"] exclude = [ "^build/", - "^templates/.*", - "templates/.*", - ".*/templates/.*", "^engibench_studies/problems/airfoil/study_[^/]*/", "^engibench/problems/airfoil/pyopt_history.py", "^docs/", From 35e9ba4b992fda539d9a6ac99f564da18b5298df Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Gerhard=20Br=C3=A4unlich?= Date: Mon, 3 Nov 2025 10:17:53 +0100 Subject: [PATCH 5/9] code quality: Also lint scripts in template folders --- .../airfoil/templates/airfoil_analysis.py | 2 +- .../problems/airfoil/templates/airfoil_opt.py | 2 +- .../airfoil/templates/cli_interface.py | 2 +- .../templates/initialize_design_2d.py | 15 +++-- .../templates/optimize_heat_conduction_2d.py | 66 ++++++++++++------- .../templates/simulate_heat_conduction_2d.py | 48 ++++++++++---- .../templates/initialize_design_3d.py | 17 +++-- .../templates/optimize_heat_conduction_3d.py | 63 ++++++++++++------ .../templates/simulate_heat_conduction_3d.py | 47 ++++++++++--- engibench/utils/files.py | 32 --------- pyproject.toml | 10 ++- tests/templates/__init__.py | 0 tests/templates/template_1.py | 4 +- tests/templates/template_2.yml | 5 +- tests/test_utils.py | 21 ++++-- 15 files changed, 215 insertions(+), 119 deletions(-) create mode 100644 tests/templates/__init__.py diff --git a/engibench/problems/airfoil/templates/airfoil_analysis.py b/engibench/problems/airfoil/templates/airfoil_analysis.py index 684ffbf..e73b1eb 100644 --- a/engibench/problems/airfoil/templates/airfoil_analysis.py +++ b/engibench/problems/airfoil/templates/airfoil_analysis.py @@ -11,7 +11,7 @@ from adflow import ADFLOW from baseclasses import AeroProblem -from cli_interface import AnalysisParameters +from cli_interface import AnalysisParameters # type:ignore[import-not-found] from cli_interface import Task from mpi4py import MPI import numpy as np diff --git a/engibench/problems/airfoil/templates/airfoil_opt.py b/engibench/problems/airfoil/templates/airfoil_opt.py index 97a5cfe..6a0b9c7 100644 --- a/engibench/problems/airfoil/templates/airfoil_opt.py +++ b/engibench/problems/airfoil/templates/airfoil_opt.py @@ -12,7 +12,7 @@ from adflow import ADFLOW from baseclasses import AeroProblem -from cli_interface import Algorithm +from cli_interface import Algorithm # type:ignore[import-not-found] from cli_interface import OptimizeParameters from idwarp import USMesh from mpi4py import MPI diff --git a/engibench/problems/airfoil/templates/cli_interface.py b/engibench/problems/airfoil/templates/cli_interface.py index 7372098..b4768ec 100644 --- a/engibench/problems/airfoil/templates/cli_interface.py +++ b/engibench/problems/airfoil/templates/cli_interface.py @@ -67,7 +67,7 @@ class OptimizeParameters: reynolds: float """Reynolds number""" temperature: float - """Temerature""" + """Temperature""" altitude: int """The cruising altitude""" area_ratio_min: float diff --git a/engibench/problems/heatconduction2d/templates/initialize_design_2d.py b/engibench/problems/heatconduction2d/templates/initialize_design_2d.py index 74bfc1d..6ee4418 100755 --- a/engibench/problems/heatconduction2d/templates/initialize_design_2d.py +++ b/engibench/problems/heatconduction2d/templates/initialize_design_2d.py @@ -1,13 +1,18 @@ #!/usr/bin/env python3 """This script sets up and initializes the design problem for a finite element analysis using dolfin adjoint based on SIMP method. + It defines the resolution, reads the design variables, and writes out the initial design to a file. """ -import os +from fenics import Constant +from fenics import FunctionSpace +from fenics import interpolate +from fenics import UnitSquareMesh +from fenics import XDMFFile import numpy as np -from fenics import * -from engibench.utils.cli import np_array_from_stdin, cast_argv + +from engibench.utils.cli import cast_argv # Extract parameters # NN: Resolution of the grid (arbitrary, affects performance) @@ -44,8 +49,8 @@ # Populate the results array with the mesh coordinates and the corresponding volume value ind = 0 - for xs in x_values: - for ys in y_values: + for _xs in x_values: + for _ys in y_values: results[ind, 0] = V # Store the volume value ind += 1 results = results.reshape(NN + 1, NN + 1) diff --git a/engibench/problems/heatconduction2d/templates/optimize_heat_conduction_2d.py b/engibench/problems/heatconduction2d/templates/optimize_heat_conduction_2d.py index ed8a412..dce194a 100755 --- a/engibench/problems/heatconduction2d/templates/optimize_heat_conduction_2d.py +++ b/engibench/problems/heatconduction2d/templates/optimize_heat_conduction_2d.py @@ -1,28 +1,49 @@ #!/usr/bin/env python3 """Topology optimization for heat conduction using the SIMP method with dolfin-adjoint. + The script reads initial design data, solves the heat conduction problem, and optimizes material distribution to minimize thermal complaicen under a volume constraint. """ -import sys +import glob +import importlib import os -from io import BytesIO import re + +from fenics import dof_to_vertex_map +from fenics import dx +from fenics import File +from fenics import FunctionSpace +from fenics import grad +from fenics import inner +from fenics import MPI +from fenics import SubDomain +from fenics import TestFunction +from fenics import XDMFFile +from fenics_adjoint import assemble +from fenics_adjoint import Constant +from fenics_adjoint import Control +from fenics_adjoint import DirichletBC +from fenics_adjoint import Function +from fenics_adjoint import InequalityConstraint +from fenics_adjoint import interpolate +from fenics_adjoint import IPOPTSolver +from fenics_adjoint import MinimizationProblem +from fenics_adjoint import ReducedFunctional +from fenics_adjoint import solve +from fenics_adjoint import UnitSquareMesh import numpy as np -from fenics import * -from fenics_adjoint import * -from typing import Any -from engibench.utils.cli import np_array_from_stdin, cast_argv +from pyadjoint.reduced_functional_numpy import set_local + +from engibench.utils.cli import cast_argv +from engibench.utils.cli import np_array_from_stdin # Ensure IPOPT is available -try: - from pyadjoint import ipopt -except ImportError: - print("""This example depends on IPOPT and Python ipopt bindings. \ +if importlib.util.find_spec("pyadjoint.ipopt") is None: + raise ImportError("""This example depends on IPOPT and Python ipopt bindings. \ When compiling IPOPT, make sure to link against HSL, as it \ is a necessity for practical problems.""") - raise # Extract parameters @@ -92,7 +113,8 @@ def k(a): class BoundaryConditions(SubDomain): """Defines Dirichlet boundary conditions on specific edges.""" - def inside(self, x, on_boundary): + def inside(self, x, _on_boundary): + """True if in the interior of the domain.""" return x[0] == 0.0 or x[1] == 1.0 or x[0] == 1.0 or (x[1] == 0.0 and (x[0] < lb_2 or x[0] > ub_2)) @@ -109,6 +131,7 @@ def inside(self, x, on_boundary): # ------------------------------- def forward(a): """Solve the heat conduction PDE given a material distribution 'a'.""" + # ruff: noqa: N806 T = Function(P, name="Temperature") v = TestFunction(P) @@ -146,6 +169,7 @@ def forward(a): class VolumeConstraint(InequalityConstraint): """Constraint to maintain volume fraction.""" + # ruff: noqa: N803 def __init__(self, V): self.V = float(V) self.smass = assemble(TestFunction(A) * Constant(1) * dx) @@ -153,17 +177,16 @@ def __init__(self, V): def function(self, m): """Compute volume constraint value.""" - from pyadjoint.reduced_functional_numpy import set_local - set_local(self.tmpvec, m) integral = self.smass.inner(self.tmpvec.vector()) return [self.V - integral] if MPI.rank(MPI.comm_world) == 0 else [] - def jacobian(self, m): + def jacobian(self, _m): """Compute Jacobian of volume constraint.""" return [-self.smass] def output_workspace(self): + """Return an object like the output of c(m) for calculations.""" return [0.0] def length(self): @@ -189,15 +212,13 @@ def length(self): objective_values = [] # Open and read the log file -with open(log_filename, "r") as f: +with open(log_filename) as f: for line in f: # Match lines that start with an iteration number followed by an objective value match = re.match(r"^\s*\d+\s+([-+]?\d*\.\d+e[-+]?\d+)", line) if match: objective_values.append(float(match.group(1))) # Extract and convert to float -# Convert to NumPy array -objective_values = np.array(objective_values) # Save optimized design mesh_output = UnitSquareMesh(NN, NN) V_output = FunctionSpace(mesh_output, "CG", 1) @@ -216,7 +237,8 @@ def length(self): np.save(output_npy, RES_OPTults) xdmf_filename = XDMFFile(MPI.comm_world, os.path.join(output_dir, f"final_solution_v={vol_f}_w={width}_.xdmf")) xdmf_filename.write(a_opt) -print("v=" + "{}".format(vol_f)) -print("w=" + "{}".format(width)) -np.savez(output_path, design=RES_OPTults, OptiStep=objective_values) -os.system("rm /home/fenics/shared/templates/RES_OPT/TEMP*") +print("v={vol_f}") +print("w={width}") +np.savez(output_path, design=RES_OPTults, OptiStep=np.array(objective_values)) +for f in glob.glob("/home/fenics/shared/templates/RES_OPT/TEMP*"): + os.remove(f) diff --git a/engibench/problems/heatconduction2d/templates/simulate_heat_conduction_2d.py b/engibench/problems/heatconduction2d/templates/simulate_heat_conduction_2d.py index e7e8e84..83c31d5 100755 --- a/engibench/problems/heatconduction2d/templates/simulate_heat_conduction_2d.py +++ b/engibench/problems/heatconduction2d/templates/simulate_heat_conduction_2d.py @@ -1,17 +1,40 @@ #!/usr/bin/env python3 """This script evaluates the design using finite element analysis with dolfin-adjoint based on the SIMP method. + It sets up the computational domain, reads the design variables, solves the forward heat conduction problem, and saves performance (thermal conductivity) metric. """ +import glob import os -from io import BytesIO -import sys + +from fenics import dof_to_vertex_map +from fenics import dx +from fenics import FunctionSpace +from fenics import grad +from fenics import inner +from fenics import MPI +from fenics import SubDomain +from fenics import TestFunction +from fenics import XDMFFile +from fenics_adjoint import assemble +from fenics_adjoint import Constant +from fenics_adjoint import Control +from fenics_adjoint import DirichletBC +from fenics_adjoint import Function +from fenics_adjoint import InequalityConstraint +from fenics_adjoint import interpolate +from fenics_adjoint import IPOPTSolver +from fenics_adjoint import MinimizationProblem +from fenics_adjoint import ReducedFunctional +from fenics_adjoint import solve +from fenics_adjoint import UnitSquareMesh import numpy as np -from fenics import * -from fenics_adjoint import * -from engibench.utils.cli import np_array_from_stdin, cast_argv +from pyadjoint.reduced_functional_numpy import set_local + +from engibench.utils.cli import cast_argv +from engibench.utils.cli import np_array_from_stdin # ------------------------------- # Initialization and Parameter Setup @@ -87,7 +110,8 @@ def k(a): class BoundaryConditions(SubDomain): """Defines Dirichlet boundary conditions on specific edges.""" - def inside(self, x, on_boundary): + def inside(self, x, _on_boundary): + """True if in the interior of the domain.""" return x[0] == 0.0 or x[1] == 1.0 or x[0] == 1.0 or (x[1] == 0.0 and (x[0] < lb_2 or x[0] > ub_2)) @@ -105,6 +129,7 @@ def inside(self, x, on_boundary): def forward(a): """Solve the heat conduction PDE given a material distribution 'a'.""" + # ruff: noqa: N806 T = Function(P, name="Temperature") v = TestFunction(P) @@ -142,6 +167,7 @@ def forward(a): class VolumeConstraint(InequalityConstraint): """Constraint to maintain volume fraction.""" + # ruff: noqa: N803 def __init__(self, V): self.V = float(V) self.smass = assemble(TestFunction(A) * Constant(1) * dx) @@ -149,17 +175,16 @@ def __init__(self, V): def function(self, m): """Compute volume constraint value.""" - from pyadjoint.reduced_functional_numpy import set_local - set_local(self.tmpvec, m) integral = self.smass.inner(self.tmpvec.vector()) return [self.V - integral] if MPI.rank(MPI.comm_world) == 0 else [] - def jacobian(self, m): + def jacobian(self, _m): """Compute Jacobian of volume constraint.""" return [-self.smass] def output_workspace(self): + """Return an object like the output of c(m) for calculations.""" return [0.0] def length(self): @@ -209,9 +234,10 @@ def length(self): # Save performance metric with open(output_path, "w") as f: - f.write("%.14f" % J_CONTROL.tape_value()) + f.write(f"{J_CONTROL.tape_value():.14f}") # Clean up temporary files -os.system("rm /home/fenics/shared/templates/RES_SIM/TEMP*") +for f in glob.glob("/home/fenics/shared/templates/RES_SIM/TEMP*"): + os.remove(f) print(f"Optimization complete: v={vol_f}, w={width}") diff --git a/engibench/problems/heatconduction3d/templates/initialize_design_3d.py b/engibench/problems/heatconduction3d/templates/initialize_design_3d.py index e3df877..407b159 100755 --- a/engibench/problems/heatconduction3d/templates/initialize_design_3d.py +++ b/engibench/problems/heatconduction3d/templates/initialize_design_3d.py @@ -1,13 +1,18 @@ #!/usr/bin/env python3 """This script sets up and initializes the design problem for a finite element analysis using dolfin adjoint based on SIMP method. + It defines the resolution, reads the design variables, and writes out the initial design to a file. """ -import os +from fenics import Constant +from fenics import FunctionSpace +from fenics import interpolate +from fenics import UnitCubeMesh +from fenics import XDMFFile import numpy as np -from fenics import * -from engibench.utils.cli import np_array_from_stdin, cast_argv + +from engibench.utils.cli import cast_argv # Extract parameters # NN: Resolution of the grid (arbitrary, affects performance) @@ -45,9 +50,9 @@ # Populate the results array with the mesh coordinates and the corresponding volume value ind = 0 - for xs in x_values: - for ys in y_values: - for zs in z_values: + for _xs in x_values: + for _ys in y_values: + for _zs in z_values: results[ind, 0] = V # Store the volume value ind += 1 results = results.reshape(NN + 1, NN + 1, NN + 1) diff --git a/engibench/problems/heatconduction3d/templates/optimize_heat_conduction_3d.py b/engibench/problems/heatconduction3d/templates/optimize_heat_conduction_3d.py index 120f0a9..780c49e 100755 --- a/engibench/problems/heatconduction3d/templates/optimize_heat_conduction_3d.py +++ b/engibench/problems/heatconduction3d/templates/optimize_heat_conduction_3d.py @@ -1,26 +1,50 @@ #!/usr/bin/env python3 """Topology optimization for heat conduction using the SIMP method with dolfin-adjoint. + The script reads initial design data, solves the heat conduction problem, and optimizes material distribution to minimize thermal complaicen under a volume constraint. """ +import glob +import importlib import os import re + +from fenics import dof_to_vertex_map +from fenics import dx +from fenics import File +from fenics import FunctionSpace +from fenics import grad +from fenics import inner +from fenics import MPI +from fenics import parameters +from fenics import SubDomain +from fenics import TestFunction +from fenics import XDMFFile +from fenics_adjoint import assemble +from fenics_adjoint import Constant +from fenics_adjoint import Control +from fenics_adjoint import DirichletBC +from fenics_adjoint import Function +from fenics_adjoint import InequalityConstraint +from fenics_adjoint import interpolate +from fenics_adjoint import IPOPTSolver +from fenics_adjoint import MinimizationProblem +from fenics_adjoint import ReducedFunctional +from fenics_adjoint import solve +from fenics_adjoint import UnitCubeMesh import numpy as np -from fenics import * -from fenics_adjoint import * -from engibench.utils.cli import np_array_from_stdin, cast_argv +from pyadjoint.reduced_functional_numpy import set_local +from engibench.utils.cli import cast_argv +from engibench.utils.cli import np_array_from_stdin # Ensure IPOPT is available -try: - from pyadjoint import ipopt -except ImportError: - print("""This example depends on IPOPT and Python ipopt bindings. \ +if importlib.util.find_spec("pyadjoint.ipopt") is None: + raise ImportError("""This example depends on IPOPT and Python ipopt bindings. \ When compiling IPOPT, make sure to link against HSL, as it \ is a necessity for practical problems.""") - raise # Extract parameters @@ -91,7 +115,8 @@ def k(a): class BoundaryConditions(SubDomain): """Defines Dirichlet boundary conditions on specific edges.""" - def inside(self, x, on_boundary): + def inside(self, x, _on_boundary): + """True if in the interior of the domain.""" return ( (x[2] > 0 and x[0] == 0) or (x[2] > 0 and x[0] == 1) @@ -119,6 +144,7 @@ def inside(self, x, on_boundary): def forward(a): """Solve the heat conduction PDE given a material distribution 'a'.""" + # ruff: noqa: N806 T = Function(P, name="Temperature") v = TestFunction(P) @@ -168,6 +194,7 @@ def forward(a): class VolumeConstraint(InequalityConstraint): """Constraint to maintain volume fraction.""" + # ruff: noqa: N803 def __init__(self, V): self.V = float(V) self.smass = assemble(TestFunction(A) * Constant(1) * dx) @@ -175,17 +202,16 @@ def __init__(self, V): def function(self, m): """Compute volume constraint value.""" - from pyadjoint.reduced_functional_numpy import set_local - set_local(self.tmpvec, m) integral = self.smass.inner(self.tmpvec.vector()) return [self.V - integral] if MPI.rank(MPI.comm_world) == 0 else [] - def jacobian(self, m): + def jacobian(self, _m): """Compute Jacobian of volume constraint.""" return [-self.smass] def output_workspace(self): + """Return an object like the output of c(m) for calculations.""" return [0.0] def length(self): @@ -211,15 +237,13 @@ def length(self): objective_values = [] # Open and read the log file -with open(log_filename, "r") as f: +with open(log_filename) as f: for line in f: # Match lines that start with an iteration number followed by an objective value match = re.match(r"^\s*\d+\s+([-+]?\d*\.\d+e[-+]?\d+)", line) if match: objective_values.append(float(match.group(1))) # Extract and convert to float -# Convert to NumPy array -objective_values = np.array(objective_values) # Save optimized design mesh_output = UnitCubeMesh(NN, NN, NN) V_output = FunctionSpace(mesh_output, "CG", 1) @@ -239,7 +263,8 @@ def length(self): np.save(output_npy, RES_OPTults) xdmf_filename = XDMFFile(MPI.comm_world, os.path.join(output_dir, f"final_solution_v={vol_f}_w={width}_.xdmf")) xdmf_filename.write(a_opt) -print("v=" + "{}".format(vol_f)) -print("w=" + "{}".format(width)) -np.savez(output_path, design=RES_OPTults, OptiStep=objective_values) -os.system("rm /home/fenics/shared/templates/RES_OPT/TEMP*") +print(f"v={vol_f}") +print(f"w={width}") +np.savez(output_path, design=RES_OPTults, OptiStep=np.array(objective_values)) +for f in glob.glob("/home/fenics/shared/templates/RES_OPT/TEMP*"): + os.remove(f) diff --git a/engibench/problems/heatconduction3d/templates/simulate_heat_conduction_3d.py b/engibench/problems/heatconduction3d/templates/simulate_heat_conduction_3d.py index 68e5546..25f143a 100755 --- a/engibench/problems/heatconduction3d/templates/simulate_heat_conduction_3d.py +++ b/engibench/problems/heatconduction3d/templates/simulate_heat_conduction_3d.py @@ -1,15 +1,41 @@ #!/usr/bin/env python3 """This script evaluates the design using finite element analysis with dolfin-adjoint based on the SIMP method. + It sets up the computational domain, reads the design variables, solves the forward heat conduction problem, and saves performance (thermal conductivity) metric. """ +import glob import os + +from fenics import dof_to_vertex_map +from fenics import dx +from fenics import FunctionSpace +from fenics import grad +from fenics import inner +from fenics import MPI +from fenics import parameters +from fenics import SubDomain +from fenics import TestFunction +from fenics import XDMFFile +from fenics_adjoint import assemble +from fenics_adjoint import Constant +from fenics_adjoint import Control +from fenics_adjoint import DirichletBC +from fenics_adjoint import Function +from fenics_adjoint import InequalityConstraint +from fenics_adjoint import interpolate +from fenics_adjoint import IPOPTSolver +from fenics_adjoint import MinimizationProblem +from fenics_adjoint import ReducedFunctional +from fenics_adjoint import solve +from fenics_adjoint import UnitCubeMesh import numpy as np -from fenics import * -from fenics_adjoint import * -from engibench.utils.cli import np_array_from_stdin, cast_argv +from pyadjoint.reduced_functional_numpy import set_local + +from engibench.utils.cli import cast_argv +from engibench.utils.cli import np_array_from_stdin # Extract parameters # NN: Grid size @@ -82,7 +108,8 @@ def k(a): class BoundaryConditions(SubDomain): """Defines Dirichlet boundary conditions on specific edges.""" - def inside(self, x, on_boundary): + def inside(self, x, _on_boundary): + """True if in the interior of the domain.""" return ( (x[2] > 0 and x[0] == 0) or (x[2] > 0 and x[0] == 1) @@ -111,6 +138,7 @@ def inside(self, x, on_boundary): def forward(a): """Solve the heat conduction PDE given a material distribution 'a'.""" + # ruff: noqa: N806 T = Function(P, name="Temperature") v = TestFunction(P) @@ -160,6 +188,7 @@ def forward(a): class VolumeConstraint(InequalityConstraint): """Constraint to maintain volume fraction.""" + # ruff: noqa: N803 def __init__(self, V): self.V = float(V) self.smass = assemble(TestFunction(A) * Constant(1) * dx) @@ -167,17 +196,16 @@ def __init__(self, V): def function(self, m): """Compute volume constraint value.""" - from pyadjoint.reduced_functional_numpy import set_local - set_local(self.tmpvec, m) integral = self.smass.inner(self.tmpvec.vector()) return [self.V - integral] if MPI.rank(MPI.comm_world) == 0 else [] - def jacobian(self, m): + def jacobian(self, _m): """Compute Jacobian of volume constraint.""" return [-self.smass] def output_workspace(self): + """Return an object like the output of c(m) for calculations.""" return [0.0] def length(self): @@ -229,9 +257,10 @@ def length(self): # Save performance metric with open(output_path, "w") as f: - f.write("%.14f" % J_CONTROL.tape_value()) + f.write(f"{J_CONTROL.tape_value():.14f}") # Clean up temporary files -os.system("rm /home/fenics/shared/templates/RES_SIM/TEMP*") +for f in glob.glob("/home/fenics/shared/templates/RES_SIM/TEMP*"): + os.remove(f) print(f"Optimization complete: v={vol_f}, w={width}") diff --git a/engibench/utils/files.py b/engibench/utils/files.py index bf642a4..0eda933 100644 --- a/engibench/utils/files.py +++ b/engibench/utils/files.py @@ -2,8 +2,6 @@ import os import shutil -from string import Template -from typing import Any def _create_study_dir(study_dir: str) -> None: @@ -51,33 +49,3 @@ def clone_dir(source_dir: str, target_dir: str) -> None: except Exception as e: msg = f"Failed to clone directory from {source_dir} to {target_dir}: {e!s}" raise RuntimeError(msg) from e - - -def replace_template_values(template_fname: str, values: dict[str, Any]) -> None: - """Replace values in a template file. - - Args: - template_fname (str): Path to the template file. - values (dict[str, Any]): Dictionary with the values to replace. - """ - if not os.path.exists(template_fname): - msg = f"Template file does not exist: {template_fname}" - raise FileNotFoundError(msg) - - try: - with open(template_fname) as f: - template = Template(f.read()) - try: - content = template.substitute(values) - except KeyError as e: - msg = f"Missing required template value: {e}" - raise ValueError(msg) from e - except ValueError as e: - msg = f"Invalid template value: {e}" - raise ValueError(msg) from e - - with open(template_fname, "w") as f: - f.write(content) - except OSError as e: - msg = f"Failed to process template file {template_fname}: {e!s}" - raise RuntimeError(msg) from e diff --git a/pyproject.toml b/pyproject.toml index ebaad33..85893e0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -189,7 +189,6 @@ order-by-type = false "TRY", # tryceratops ] -"**/heatconduction*/templates/**/*.py" = ["ALL"] ###################################### FORMAT ######################################## @@ -222,7 +221,7 @@ docstring-code-line-length = "dynamic" ###################################### PYRIGHT ######################################## [tool.pyright] include = ["engibench/**"] -exclude = ["**/node_modules", "**/__pycache__", "**/templates/**", "**/study*", "**/pyopt_history.py"] +exclude = ["**/node_modules", "**/__pycache__", "**/study*", "**/pyopt_history.py"] strict = [] typeCheckingMode = "basic" @@ -269,5 +268,12 @@ module = [ "ceviche.*", "sqlitedict", "sqlitedict.*", + "pyadjoint.*", + "adflow", + "baseclasses", + "idwarp", + "mpi4py", + "multipoint", + "pygeo", ] ignore_missing_imports = true diff --git a/tests/templates/__init__.py b/tests/templates/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/templates/template_1.py b/tests/templates/template_1.py index 7de8453..cb4172a 100644 --- a/tests/templates/template_1.py +++ b/tests/templates/template_1.py @@ -1,3 +1,3 @@ -hello = $hello -world = $world +hello = "hello" +world = "world" print(hello, world) diff --git a/tests/templates/template_2.yml b/tests/templates/template_2.yml index 388f634..da7e4ea 100644 --- a/tests/templates/template_2.yml +++ b/tests/templates/template_2.yml @@ -1,2 +1,3 @@ -hi: $hi -world: $world +--- +hi: hello +world: world diff --git a/tests/test_utils.py b/tests/test_utils.py index 92b1a53..5b435dc 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -1,7 +1,6 @@ import os from engibench.utils.files import clone_dir -from engibench.utils.files import replace_template_values def test_clone_template() -> None: @@ -14,20 +13,30 @@ def test_clone_template() -> None: # Cloning clone_dir(template_dir, study_dir) - # Replacement - replace_template_values(study_dir + "/template_1.py", {"hello": "hello", "world": "world"}) # Tests the replacement with open(study_dir + "/template_1.py") as f: content = f.read() - assert content == "hello = hello\nworld = world\nprint(hello, world)\n" + assert ( + content + == """hello = "hello" +world = "world" +print(hello, world) +""" + ) # Another file, in another format - replace_template_values(study_dir + "/template_2.yml", {"hi": "hello", "world": "world"}) with open(study_dir + "/template_2.yml") as f: content = f.read() - assert content == "hi: hello\nworld: world\n" + assert ( + content + == """--- +hi: hello +world: world +""" + ) # Cleanup os.remove(study_dir + "/template_1.py") os.remove(study_dir + "/template_2.yml") + os.remove(study_dir + "/__init__.py") os.rmdir(study_dir) From 02c7f1810520c2b3133e4442ec161deda84b0e25 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Gerhard=20Br=C3=A4unlich?= Date: Tue, 4 Nov 2025 09:24:14 +0100 Subject: [PATCH 6/9] airfoil: Fix duplicate argument bug in optimize --- engibench/problems/airfoil/v0.py | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/engibench/problems/airfoil/v0.py b/engibench/problems/airfoil/v0.py index 72459a6..e207119 100644 --- a/engibench/problems/airfoil/v0.py +++ b/engibench/problems/airfoil/v0.py @@ -553,18 +553,20 @@ def optimize( if "opt" in config: config["opt"] = cli_interface.Algorithm[config["opt"]] args = cli_interface.OptimizeParameters( - **dataclasses.asdict(self.Conditions()), - alpha=starting_point["angle_of_attack"], - altitude=10000, - temperature=300, # should specify either mach + altitude or mach + reynolds + reynoldsLength (default to 1) + temperature - use_altitude=False, - opt=cli_interface.Algorithm.SLSQP, - opt_options={}, - output_dir=self.__docker_study_dir + "/output", - ffd_fname=self.__docker_study_dir + "/" + filename + "_ffd.xyz", - mesh_fname=self.__docker_study_dir + "/" + filename + ".cgns", - area_input_design=calc_area(starting_point["coords"]), - **config, + **{ + **dataclasses.asdict(self.Conditions()), + "alpha": starting_point["angle_of_attack"], + "altitude": 10000, + "temperature": 300, # should specify either mach + altitude or mach + reynolds + reynoldsLength (default to 1) + temperature + "use_altitude": False, + "opt": cli_interface.Algorithm.SLSQP, + "opt_options": {}, + "output_dir": self.__docker_study_dir + "/output", + "ffd_fname": self.__docker_study_dir + "/" + filename + "_ffd.xyz", + "mesh_fname": self.__docker_study_dir + "/" + filename + ".cgns", + "area_input_design": calc_area(starting_point["coords"]), + **config, + }, ) self.__design_to_simulator_input( starting_point, reynolds=args.reynolds, mach=args.reynolds, temperature=args.temperature, filename=filename From 6cee3779262ee4d01a8ab43f288d586ce53a2642 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Gerhard=20Br=C3=A4unlich?= Date: Tue, 4 Nov 2025 09:25:46 +0100 Subject: [PATCH 7/9] airfoil: Raise an early error when area_initial is not specified --- engibench/problems/airfoil/v0.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/engibench/problems/airfoil/v0.py b/engibench/problems/airfoil/v0.py index e207119..7a96cf0 100644 --- a/engibench/problems/airfoil/v0.py +++ b/engibench/problems/airfoil/v0.py @@ -289,7 +289,7 @@ class Conditions: reynolds: Annotated[ float, bounded(lower=0.0).category(IMPL), bounded(lower=1e5, upper=1e9).warning().category(IMPL) ] = 1e6 - area_initial: float | None = None + area_initial: float = float("NAN") """actual initial airfoil area""" area_ratio_min: Annotated[float, bounded(lower=0.0, upper=1.2).category(THEORY)] = 0.7 """Minimum ratio the initial area is allowed to decrease to i.e minimum_area = area_initial*area_target""" @@ -315,12 +315,12 @@ class Config(Conditions): @constraint(categories=THEORY) @staticmethod - def area_ratio_bound(area_ratio_min: float, area_initial: float | None, area_input_design: float | None) -> None: + def area_ratio_bound(area_ratio_min: float, area_initial: float, area_input_design: float | None) -> None: """Constraint for area_ratio_min <= area_ratio <= 1.2.""" area_ratio_max = 1.2 if area_input_design is None: return - assert area_initial is not None + assert not np.isnan(area_initial) area_ratio = area_input_design / area_initial assert area_ratio_min <= area_ratio <= area_ratio_max, ( f"Config.area_ratio: {area_ratio} ∉ [area_ratio_min={area_ratio_min}, 1.2]" @@ -550,6 +550,8 @@ def optimize( # Prepares the optimize_airfoil.py script with the optimization configuration fields = {f.name for f in dataclasses.fields(cli_interface.OptimizeParameters)} config = {key: val for key, val in (config or {}).items() if key in fields} + if "area_initial" not in config: + raise ValueError("optimize(): config is missing the required parameter 'area_initial'") if "opt" in config: config["opt"] = cli_interface.Algorithm[config["opt"]] args = cli_interface.OptimizeParameters( From d836688baf009059a478ba21269e5616f89cb896 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Gerhard=20Br=C3=A4unlich?= Date: Fri, 14 Nov 2025 07:56:32 +0100 Subject: [PATCH 8/9] airfoil: include pid in the name of the study directory This prevents file collisions when running multiple simultations with the same seed in parallel. --- engibench/problems/airfoil/v0.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/engibench/problems/airfoil/v0.py b/engibench/problems/airfoil/v0.py index 7a96cf0..c96be38 100644 --- a/engibench/problems/airfoil/v0.py +++ b/engibench/problems/airfoil/v0.py @@ -363,7 +363,7 @@ def reset(self, seed: int | None = None, *, cleanup: bool = False) -> None: shutil.rmtree(self.__local_study_dir) super().reset(seed) - self.current_study = f"study_{self.seed}" + self.current_study = f"study_{self.seed}-pid{os.getpid()}" self.__local_study_dir = self.__local_target_dir + "/" + self.current_study self.__docker_study_dir = self.__docker_target_dir + "/" + self.current_study From fdc26af7d63f985681a4085ab9b3449cb0149d14 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Gerhard=20Br=C3=A4unlich?= Date: Fri, 14 Nov 2025 08:46:06 +0100 Subject: [PATCH 9/9] utils.container: Make CONTAINER_RUNTIME case insensitive --- engibench/utils/container.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/engibench/utils/container.py b/engibench/utils/container.py index a4d2fd8..f59e4cc 100644 --- a/engibench/utils/container.py +++ b/engibench/utils/container.py @@ -124,9 +124,9 @@ def runtime() -> type[ContainerRuntime] | None: Class object of the first available container runtime or the container runtime selected by the `CONTAINER_RUNTIME` environment variable if set. """ - runtimes_by_name = {rt.name: rt for rt in RUNTIMES} + runtimes_by_name = {rt.name.lower(): rt for rt in RUNTIMES} rt_name = os.environ.get("CONTAINER_RUNTIME") - rt = runtimes_by_name.get(rt_name) if rt_name is not None else None + rt = runtimes_by_name.get(rt_name.lower()) if rt_name is not None else None if rt is not None: return rt for rt in RUNTIMES: