From 88ae03e093e49226444fc5259228daabb671448b Mon Sep 17 00:00:00 2001 From: Chris Psenica Date: Fri, 19 Sep 2025 20:31:53 -0500 Subject: [PATCH] Updating the runscript for the square CS UBend case to V4 --- UBend_Channel/runScript.py | 303 +++++++++++++++++++++++++++++++ UBend_Channel/runScript_v2.py | 328 ---------------------------------- 2 files changed, 303 insertions(+), 328 deletions(-) create mode 100755 UBend_Channel/runScript.py delete mode 100755 UBend_Channel/runScript_v2.py diff --git a/UBend_Channel/runScript.py b/UBend_Channel/runScript.py new file mode 100755 index 00000000..20f644b0 --- /dev/null +++ b/UBend_Channel/runScript.py @@ -0,0 +1,303 @@ +#!/usr/bin/env python +""" +DAFoam run script for the U bend channel case v4 +""" + +# ============================================================================= +# Imports +# ============================================================================= +import argparse +from mpi4py import MPI +import os +import numpy as np +import openmdao.api as om +from mphys.multipoint import Multipoint +from dafoam.mphys import DAFoamBuilder +from mphys.scenario_aerodynamic import ScenarioAerodynamic +from pygeo.mphys import OM_DVGEOCOMP +from pygeo import geo_utils + +# ============================================================================= +# Input Parameters +# ============================================================================= +parser = argparse.ArgumentParser() +# which optimizer to use. Options are: IPOPT (default), SLSQP, and SNOPT +parser.add_argument("-optimizer", help="optimizer to use", type=str, default="SLSQP") +# which task to run. Options are: run_driver (default), run_model, compute_totals, check_totals +parser.add_argument("-task", help="type of run to do", type=str, default="run_driver") +args = parser.parse_args() +gcomm = MPI.COMM_WORLD + +# ============================================================================= +# daOptions Setup +# ============================================================================= +CPL_weight = 0.50 # Weight of pressure loss (PL) in Obj. Function +HFX_weight = CPL_weight - 1.0 # Weight of heat flux (HFX) in obj. Function +HFX0 = 305 # HFX value for baseline design +CPL0 = 85.23 - 35.62 # PL value for baseline design +U0 = 8.4 # Fluid flow (m/s) in x-direction + +daOptionsAero = { + "solverName": "DASimpleFoam", + "designSurfaces": ["ubend"], + "useAD": {"mode": "reverse"}, + "primalMinResTol": 1e-8, + "primalMinResTolDiff": 1e7, + "writeMinorIterations": True, + "wallDistanceMethod": "daCustom", + "primalBC" : { + "useWallFunction": True, + }, + + "function": { + + "TP1": { + "type": "totalPressure", + "source": "patchToFace", + "patches": ["inlet"], + "scale": 1.0, + "addToAdjoint": True, + }, + + "TP2": { + "type": "totalPressure", + "source": "patchToFace", + "patches": ["outlet"], + "scale": 1.0, + "addToAdjoint": True, + }, + + "HFX": { + "type": "wallHeatFlux", + "source": "patchToFace", + "patches": ["ubend"], + "scale": 1.0, + "addToAdjoint": True, + }, + + }, + + "adjStateOrdering": "cell", + + "adjEqnOption": {"gmresRelTol" : 1e-5, + "gmresTolDiff" : 1e4, + "pcFillLevel" : 2, + "jacMatReOrdering" : "natural", + "gmresMaxIters" : 3000, + "gmresRestart" : 3000}, + + "normalizeStates": {"U" : U0, + "p" : (U0 * U0) / 2., + "nuTilda" : 1e-3, + "phi" : 1.0, + "T" : 300}, + + "inputInfo": { + "aero_vol_coords": {"type": "volCoord", "components": ["solver" , "function"]}, + }, + + "outputInfo": { + "q_convect": { + "type" : "thermalCouplingOutput", + "patches" : ["ubend"], + "components" : ["thermalCoupling"], + }, + }, + +} + +# ============================================================================= +# Mesh Setup +# ============================================================================= +meshOptions = { + "gridFile" : os.getcwd(), + "fileType" : "OpenFOAM", + "symmetryPlanes" : [], +} + +# ============================================================================= +# Top class To Setup The Optimization Problem +# ============================================================================= +class Top(Multipoint): + def setup(self): + + # initialize builders + dafoam_builder = DAFoamBuilder(daOptionsAero , meshOptions , scenario = "aerodynamic") + dafoam_builder.initialize(self.comm) + + # add design variable component and promote to top level + self.add_subsystem("dvs" , om.IndepVarComp() , promotes = ["*"]) + + # add mesh component + self.add_subsystem("mesh_aero" , dafoam_builder.get_mesh_coordinate_subsystem()) + + # add geometry component + self.add_subsystem("geometry_aero" , OM_DVGEOCOMP(file = "FFD/testFFD.xyz" , type = "ffd")) + + # add a scenario (flow condition) for optimization. For no themal (solid) use ScenarioAerodynamic, for thermal (solid) use ScenarioAerothermal + self.mphys_add_scenario("scenario" , ScenarioAerodynamic(aero_builder = dafoam_builder)) + + # need to manually connect the x_aero0 between the mesh and geometry components + self.connect("mesh_aero.x_aero0" , "geometry_aero.x_aero_in") + self.connect("geometry_aero.x_aero0" , "scenario.x_aero") + + # add obj val for PL + self.add_subsystem("OBJ" , om.ExecComp("val = scalePL * (TP1 - TP2) + (scaleHFX * HFX)" , scalePL = {'val' : CPL_weight / CPL0 , 'constant' : True} , scaleHFX = {'val' : HFX_weight / HFX0 , 'constant' : True})) + + def configure(self): + + # initialize the optimization + super().configure() + + # get surface coordinates from mesh component + points_aero = self.mesh_aero.mphys_get_surface_mesh() + + # add pointset to the geometry component + self.geometry_aero.nom_add_discipline_coords("aero" , points_aero) + + # get FFD points + pts = self.geometry_aero.nom_getDVGeo().getLocalIndex(0) + + #---------- setup DVs ---------- + # shapex + indexList = [] + indexList.extend(pts[7:16 , 1 , :].flatten()) + PS = geo_utils.PointSelect("list" , indexList) + shapexUpper = self.geometry_aero.nom_addLocalDV(dvName = "shapexUpper" , pointSelect = PS , axis = "x") + + # shapey + indexList = [] + indexList.extend(pts[7:16 , 1 , :].flatten()) + PS = geo_utils.PointSelect("list" , indexList) + shapeyUpper = self.geometry_aero.nom_addLocalDV(dvName = "shapeyUpper" , pointSelect = PS , axis = "y") + + # shapez + indexList = [] + indexList.extend(pts[7:16 , 1 , :].flatten()) + PS = geo_utils.PointSelect("list", indexList) + shapezUpper = self.geometry_aero.nom_addLocalDV(dvName = "shapezUpper" , pointSelect = PS , axis = "z") + + # shapex + indexList = [] + indexList.extend(pts[7:16 , 0 , :].flatten()) + PS = geo_utils.PointSelect("list" , indexList) + shapexLower = self.geometry_aero.nom_addLocalDV(dvName = "shapexLower" , pointSelect = PS , axis = "x") + + # shapey + indexList = [] + indexList.extend(pts[7:16 , 0 , :].flatten()) + PS = geo_utils.PointSelect("list" , indexList) + shapeyLower = self.geometry_aero.nom_addLocalDV(dvName = "shapeyLower" , pointSelect = PS , axis = "y") + + # shapez + indexList = [] + indexList.extend(pts[7:16 , 0 , :].flatten()) + PS = geo_utils.PointSelect("list", indexList) + shapezLower = self.geometry_aero.nom_addLocalDV(dvName = "shapezLower" , pointSelect = PS , axis = "z") + + #---------- Add Outputs For The DVs ---------- + self.dvs.add_output("shapexUpper" , val = np.array([0]*shapexUpper)) + self.dvs.add_output("shapeyUpper" , val = np.array([0]*shapeyUpper)) + self.dvs.add_output("shapezUpper" , val = np.array([0]*shapezUpper)) + + self.dvs.add_output("shapexLower" , val = np.array([0]*shapexLower)) + self.dvs.add_output("shapeyLower" , val = np.array([0]*shapeyLower)) + self.dvs.add_output("shapezLower" , val = np.array([0]*shapezLower)) + + #---------- Connect The Design Variables To The Geometry ---------- + self.connect("shapexUpper" , "geometry_aero.shapexUpper") + self.connect("shapeyUpper" , "geometry_aero.shapeyUpper") + self.connect("shapezUpper" , "geometry_aero.shapezUpper") + + self.connect("shapexLower" , "geometry_aero.shapexLower") + self.connect("shapeyLower" , "geometry_aero.shapeyLower") + self.connect("shapezLower" , "geometry_aero.shapezLower") + + #---------- Define The Design Variables To The Top Level ---------- + self.add_design_var("shapexUpper" , lower = -0.04 , upper = 0.04 , scaler = 25.0) + self.add_design_var("shapeyUpper" , lower = -0.04 , upper = 0.04 , scaler = 25.0) + self.add_design_var("shapezUpper" , lower = -0.04 , upper = 0.04 , scaler = 25.0) + self.add_design_var("shapexLower" , lower = -0.04 , upper = 0.04 , scaler = 25.0) + self.add_design_var("shapeyLower" , lower = -0.04 , upper = 0.04 , scaler = 25.0) + self.add_design_var("shapezLower" , lower = -0.04 , upper = 0.04 , scaler = 25.0) + + # add objective and constraints + self.connect("scenario.aero_post.TP1" , "OBJ.TP1") + self.connect("scenario.aero_post.TP2" , "OBJ.TP2") + self.connect("scenario.aero_post.HFX" , "OBJ.HFX") + self.add_objective("OBJ.val" , scaler = 1.0) + +# ============================================================================= +# Problem Setup +# ============================================================================= +prob = om.Problem(reports = None) +prob.model = Top() +prob.setup(mode = "rev") +prob.driver = om.pyOptSparseDriver() +prob.driver.options["optimizer"] = args.optimizer + +if args.optimizer == "SNOPT": + prob.driver.opt_settings = { + "Major feasibility tolerance" : 1.0e-5, + "Major optimality tolerance" : 1.0e-5, + "Minor feasibility tolerance" : 1.0e-5, + "Verify level" : -1, + "Function precision" : 1.0e-5, + "Major iterations limit" : 100, + "Nonderivative linesearch" : None, + "Print file" : "opt_SNOPT_print.txt", + "Summary file" : "opt_SNOPT_summary.txt", + } +elif args.optimizer == "IPOPT": + prob.driver.opt_settings = { + "tol" : 1.0e-5, + "constr_viol_tol" : 1.0e-5, + "max_iter" : 100, + "print_level" : 5, + "output_file" : "opt_IPOPT.txt", + "mu_strategy" : "adaptive", + "limited_memory_max_history" : 10, + "nlp_scaling_method" : "none", + "alpha_for_y" : "full", + "recalc_y" : "yes", + } +elif args.optimizer == "SLSQP": + prob.driver.opt_settings = { + "ACC" : 1.0e-5, + "MAXIT" : 100, + "IFILE" : "opt_SLSQP.txt", + } +else: + print("optimizer arg not valid!") + exit(1) + +prob.driver.options["debug_print"] = ["nl_cons" , "objs" , "desvars"] +prob.driver.options["print_opt_prob"] = True +prob.driver.hist_file = "OptView.hst" + +if args.task == "run_driver": + prob.run_driver() + +elif args.task == "run_model": + prob.run_model() + +elif args.task == "compute_totals": + prob.run_model() + totals = prob.compute_totals() + + if MPI.COMM_WORLD.rank == 0: + print(totals) + +elif args.task == "check_totals": + prob.run_model() + + prob.check_totals( + compact_print = False, + step = 1e-4, + form = "central", + step_calc = "abs", + ) + +else: + print("task arg not found!") + exit(1) \ No newline at end of file diff --git a/UBend_Channel/runScript_v2.py b/UBend_Channel/runScript_v2.py deleted file mode 100755 index fb0edb5d..00000000 --- a/UBend_Channel/runScript_v2.py +++ /dev/null @@ -1,328 +0,0 @@ -#!/usr/bin/env python -""" -DAFoam run script for the U bend channel case -""" - -# ============================================================================= -# Imports -# ============================================================================= -import os -import argparse -from mpi4py import MPI -from dafoam import PYDAFOAM, optFuncs -from pygeo import * -from pyspline import * -from idwarp import * -from pyoptsparse import Optimization, OPT -import numpy as np - -# ============================================================================= -# Input Parameters -# ============================================================================= -parser = argparse.ArgumentParser() -# which optimizer to use. Options are: slsqp (default), snopt, or ipopt -parser.add_argument("--opt", help="optimizer to use", type=str, default="ipopt") -# which task to run. Options are: opt (default), run, testSensShape, or solveCL -parser.add_argument("--task", help="type of run to do", type=str, default="opt") -args = parser.parse_args() -gcomm = MPI.COMM_WORLD - -HFL0 = 241.3 -HFL_weight = -0.5 -CPL0 = 40.26 -CPL_weight = 0.5 - -# Set the parameters for optimization -daOptions = { - "solverName": "DASimpleTFoam", - "useAD": {"mode": "reverse"}, - "designSurfaces": ["ubend", "ubendup"], - "primalMinResTol": 1e-8, - "primalMinResTolDiff": 1e5, - "objFunc": { - "obj": { - "part1": { - "type": "totalPressure", - "source": "patchToFace", - "patches": ["inlet"], - "scale": 1.0 / CPL0 * CPL_weight, - "addToAdjoint": True, - }, - "part2": { - "type": "totalPressure", - "source": "patchToFace", - "patches": ["outlet"], - "scale": -1.0 / CPL0 * CPL_weight, - "addToAdjoint": True, - }, - "part3": { - "type": "wallHeatFlux", - "source": "patchToFace", - "patches": ["ubendup"], - "scale": 1.0 / HFL0 * HFL_weight, - "addToAdjoint": True, - } - }, - }, - "normalizeStates": {"U": 10.0, "p": 30.0, "nuTilda": 1e-3, "phi": 1.0, "T": 300.0}, - "adjPartDerivFDStep": {"State": 1e-6, "FFD": 1e-3}, - "adjEqnOption": {"gmresRelTol": 1.0e-6, "pcFillLevel": 1, "jacMatReOrdering": "rcm"}, - "adjPCLag": 5, - # Design variable setup - "designVar": { - "shapez": {"designVarType": "FFD"}, - "shapeyouter": {"designVarType": "FFD"}, - "shapeyinner": {"designVarType": "FFD"}, - "shapexinner": {"designVarType": "FFD"}, - "shapexouter1": {"designVarType": "FFD"}, - "shapexouter2": {"designVarType": "FFD"}, - }, -} - -# mesh warping parameters, users need to manually specify the symmetry plane -meshOptions = { - "gridFile": os.getcwd(), - "fileType": "OpenFOAM", - # point and normal for the symmetry plane - "symmetryPlanes": [[[0.0, 0.0, 0.0], [0.0, 0.0, -1.0]]], -} - -# options for optimizers -if args.opt == "snopt": - optOptions = { - "Major feasibility tolerance": 1.0e-7, - "Major optimality tolerance": 1.0e-7, - "Minor feasibility tolerance": 1.0e-7, - "Verify level": -1, - "Function precision": 1.0e-7, - "Major iterations limit": 50, - "Nonderivative linesearch": None, - "Print file": "opt_SNOPT_print.txt", - "Summary file": "opt_SNOPT_summary.txt", - } -elif args.opt == "ipopt": - optOptions = { - "tol": 1.0e-7, - "constr_viol_tol": 1.0e-7, - "max_iter": 50, - "print_level": 5, - "output_file": "opt_IPOPT.txt", - "mu_strategy": "adaptive", - "limited_memory_max_history": 10, - "nlp_scaling_method": "none", - "alpha_for_y": "full", - "recalc_y": "yes", - } -elif args.opt == "slsqp": - optOptions = { - "ACC": 1.0e-7, - "MAXIT": 50, - "IFILE": "opt_SLSQP.txt", - } -else: - print("opt arg not valid!") - exit(0) - - -# ============================================================================= -# Design variable setup -# ============================================================================= -DVGeo = DVGeometry("./FFD/UBendDuctFFDSym.xyz") -# Select points -pts = DVGeo.getLocalIndex(0) -# shapez -indexList = [] -indexList.extend(pts[7:16, :, -1].flatten()) -PS = geo_utils.PointSelect("list", indexList) -DVGeo.addLocalDV("shapez", lower=-0.04, upper=0.0077, axis="z", scale=1.0, pointSelect=PS, config="configz") -# shapeyouter -indexList = [] -indexList.extend(pts[7:16, -1, :].flatten()) -PS = geo_utils.PointSelect("list", indexList) -DVGeo.addLocalDV("shapeyouter", lower=-0.02, upper=0.02, axis="y", scale=1.0, pointSelect=PS, config="configyouter") -# shapeyinner -indexList = [] -indexList.extend(pts[7:16, 0, :].flatten()) -PS = geo_utils.PointSelect("list", indexList) -DVGeo.addLocalDV("shapeyinner", lower=-0.04, upper=0.04, axis="y", scale=1.0, pointSelect=PS, config="configyinner") -# shapexinner -indexList = [] -indexList.extend(pts[7:16, 0, :].flatten()) -PS = geo_utils.PointSelect("list", indexList) -DVGeo.addLocalDV("shapexinner", lower=-0.04, upper=0.04, axis="x", scale=1.0, pointSelect=PS, config="configxinner") - -# shapexouter1 -indexList = [] -indexList.extend(pts[9, -1, :].flatten()) -PS = geo_utils.PointSelect("list", indexList) -DVGeo.addLocalDV( - "shapexouter1", lower=-0.05, upper=0.05, axis="x", scale=1.0, pointSelect=PS, config="configxouter1" -) - -# shapexouter2 -indexList = [] -indexList.extend(pts[10, -1, :].flatten()) -PS = geo_utils.PointSelect("list", indexList) -DVGeo.addLocalDV("shapexouter2", lower=-0.05, upper=0.0, axis="x", scale=1.0, pointSelect=PS, config="configxouter2") - - -# ============================================================================= -# DAFoam initialization -# ============================================================================= -DASolver = PYDAFOAM(options=daOptions, comm=gcomm) -DASolver.setDVGeo(DVGeo) -mesh = USMesh(options=meshOptions, comm=gcomm) -DASolver.printFamilyList() -DASolver.setMesh(mesh) -evalFuncs = [] -DASolver.setEvalFuncs(evalFuncs) - -# ============================================================================= -# Constraint setup -# ============================================================================= -DVCon = DVConstraints() -DVCon.setDVGeo(DVGeo) -DVCon.setSurface(DASolver.getTriangulatedMeshSurface(groupName=DASolver.designSurfacesGroup)) - -if args.task == "opt": - # Create a linear constraint so that the curvature at the symmetry plane is zero - indSetA = [] - indSetB = [] - for i in range(7, 16, 1): - indSetA.append(pts[i, 0, 0]) - indSetB.append(pts[i, 0, 1]) - DVCon.addLinearConstraintsShape( - indSetA, indSetB, factorA=1.0, factorB=-1.0, lower=0.0, upper=0.0, config="configyinner" - ) - - indSetA = [] - indSetB = [] - for i in range(7, 16, 1): - indSetA.append(pts[i, -1, 0]) - indSetB.append(pts[i, -1, 1]) - DVCon.addLinearConstraintsShape( - indSetA, indSetB, factorA=1.0, factorB=-1.0, lower=0.0, upper=0.0, config="configyouter" - ) - - indSetA = [] - indSetB = [] - for i in range(7, 16, 1): - indSetA.append(pts[i, 0, 0]) - indSetB.append(pts[i, 0, 1]) - DVCon.addLinearConstraintsShape( - indSetA, indSetB, factorA=1.0, factorB=-1.0, lower=0.0, upper=0.0, config="configxinner" - ) - - indSetA = [] - indSetB = [] - for i in [9]: - indSetA.append(pts[i, -1, 0]) - indSetB.append(pts[i, -1, 1]) - DVCon.addLinearConstraintsShape( - indSetA, indSetB, factorA=1.0, factorB=-1.0, lower=0.0, upper=0.0, config="configxouter1" - ) - - indSetA = [] - indSetB = [] - for i in [10]: - indSetA.append(pts[i, -1, 0]) - indSetB.append(pts[i, -1, 1]) - DVCon.addLinearConstraintsShape( - indSetA, indSetB, factorA=1.0, factorB=-1.0, lower=0.0, upper=0.0, config="configxouter2" - ) - # linear constraint to make sure the inner bend does not intersect - minInnerBend = 0.005 - indSetA = [] - indSetB = [] - for k in range(3): - indSetA.append(pts[8, 0, k]) - indSetB.append(pts[14, 0, k]) - DVCon.addLinearConstraintsShape( - indSetA, - indSetB, - factorA=1.0, - factorB=-1.0, - lower=-0.03300 + minInnerBend, - upper=0.03300 - minInnerBend, - config="configyinner", - ) - indSetA = [] - indSetB = [] - for k in range(3): - indSetA.append(pts[9, 0, k]) - indSetB.append(pts[13, 0, k]) - DVCon.addLinearConstraintsShape( - indSetA, - indSetB, - factorA=1.0, - factorB=-1.0, - lower=-0.02853 + minInnerBend, - upper=0.02853 - minInnerBend, - config="configyinner", - ) - indSetA = [] - indSetB = [] - for k in range(3): - indSetA.append(pts[10, 0, k]) - indSetB.append(pts[12, 0, k]) - DVCon.addLinearConstraintsShape( - indSetA, - indSetB, - factorA=1.0, - factorB=-1.0, - lower=-0.01635 + minInnerBend, - upper=0.01635 - minInnerBend, - config="configyinner", - ) - -DVCon.writeTecplot("DVConstraints.dat") - -# ============================================================================= -# Initialize optFuncs for optimization -# ============================================================================= -optFuncs.DASolver = DASolver -optFuncs.DVGeo = DVGeo -optFuncs.DVCon = DVCon -optFuncs.evalFuncs = evalFuncs -optFuncs.gcomm = gcomm - -# ============================================================================= -# Task -# ============================================================================= -if args.task == "opt": - - optProb = Optimization("opt", objFun=optFuncs.calcObjFuncValues, comm=gcomm) - DVGeo.addVariablesPyOpt(optProb) - DVCon.addConstraintsPyOpt(optProb) - - # Add objective - optProb.addObj("obj", scale=1) - # Add physical constraints - #optProb.addCon("HFL", lower=HFL_target, upper=HFL_target, scale=1) - - if gcomm.rank == 0: - print(optProb) - - DASolver.runColoring() - - opt = OPT(args.opt, options=optOptions) - histFile = "./%s_hist.hst" % args.opt - sol = opt(optProb, sens=optFuncs.calcObjFuncSens, storeHistory=histFile) - if gcomm.rank == 0: - print(sol) - -elif args.task == "runPrimal": - - optFuncs.runPrimal() - -elif args.task == "runAdjoint": - - optFuncs.runAdjoint() - -elif args.task == "runForwardAD": - - optFuncs.runForwardAD("shapez", 0) - -else: - print("task arg not found!") - exit(0)