From 19311048f18faf464edf3be5082edfa93f90e747 Mon Sep 17 00:00:00 2001 From: Ping He Date: Sun, 2 Feb 2025 12:36:25 -0600 Subject: [PATCH] Added aeroOpt test (#755) * Added unit test for primal. * Added the aeroOpt test. * Deleted the rhoSimple test. * Minor changes to the primal test. --- ....py => backup_runRegTests_DASimpleFoam.py} | 0 tests/refs/DAFoam_Test_AeroOptRef.txt | 4 + tests/refs/DAFoam_Test_DARhoSimpleFoamRef.txt | 14 -- tests/refs/DAFoam_Test_DASimpleFoamRef.txt | 39 --- tests/runRegTests_AeroOpt.py | 233 ++++++++++++++++++ tests/runRegTests_DARhoSimpleFoam.py | 168 ------------- tests/runUnitTests_primal.py | 67 +++++ 7 files changed, 304 insertions(+), 221 deletions(-) rename tests/{runRegTests_DASimpleFoam.py => backup_runRegTests_DASimpleFoam.py} (100%) create mode 100755 tests/refs/DAFoam_Test_AeroOptRef.txt delete mode 100755 tests/refs/DAFoam_Test_DARhoSimpleFoamRef.txt delete mode 100755 tests/refs/DAFoam_Test_DASimpleFoamRef.txt create mode 100755 tests/runRegTests_AeroOpt.py delete mode 100755 tests/runRegTests_DARhoSimpleFoam.py create mode 100755 tests/runUnitTests_primal.py diff --git a/tests/runRegTests_DASimpleFoam.py b/tests/backup_runRegTests_DASimpleFoam.py similarity index 100% rename from tests/runRegTests_DASimpleFoam.py rename to tests/backup_runRegTests_DASimpleFoam.py diff --git a/tests/refs/DAFoam_Test_AeroOptRef.txt b/tests/refs/DAFoam_Test_AeroOptRef.txt new file mode 100755 index 00000000..02aaf073 --- /dev/null +++ b/tests/refs/DAFoam_Test_AeroOptRef.txt @@ -0,0 +1,4 @@ +Dictionary Key: CD +@value 0.0167462707303083 1e-05 1e-10 +Dictionary Key: CL +@value 0.2992242625146929 1e-05 1e-10 \ No newline at end of file diff --git a/tests/refs/DAFoam_Test_DARhoSimpleFoamRef.txt b/tests/refs/DAFoam_Test_DARhoSimpleFoamRef.txt deleted file mode 100755 index 1f393a1e..00000000 --- a/tests/refs/DAFoam_Test_DARhoSimpleFoamRef.txt +++ /dev/null @@ -1,14 +0,0 @@ -Dictionary Key: CD -@value 0.0074800436855696 1e-10 1e-12 -Dictionary Key: CL -@value 0.2191086912905469 1e-10 1e-12 -Dictionary Key: CD -Dictionary Key: twist-Adjoint -@value 0.0017099075339183 1e-08 1e-12 -Dictionary Key: twist-FD -@value 0.0017099093516104 1e-08 1e-12 -Dictionary Key: CL -Dictionary Key: twist-Adjoint -@value 0.1497633960444238 1e-08 1e-12 -Dictionary Key: twist-FD -@value 0.1497632302527308 1e-08 1e-12 \ No newline at end of file diff --git a/tests/refs/DAFoam_Test_DASimpleFoamRef.txt b/tests/refs/DAFoam_Test_DASimpleFoamRef.txt deleted file mode 100755 index e7609ede..00000000 --- a/tests/refs/DAFoam_Test_DASimpleFoamRef.txt +++ /dev/null @@ -1,39 +0,0 @@ -Dictionary Key: CD -@value 0.0170618714819692 1e-10 1e-12 -Dictionary Key: CL -@value 0.2984067151382958 1e-10 1e-12 -Dictionary Key: nonOrtho -@value 23.8934356542776314 1e-10 1e-12 -Dictionary Key: skewness -@value 1.4671221512193253 1e-10 1e-12 -Dictionary Key: CL -Dictionary Key: patchV-Adjoint -@value 0.0602938548210897 1e-08 1e-12 -@value 0.0971078786360775 1e-08 1e-12 -Dictionary Key: patchV-FD -@value 0.0602938310007630 1e-08 1e-12 -@value 0.0971078612989515 1e-08 1e-12 -Dictionary Key: shape-Adjoint -@value 1.2435369168112151 1e-08 1e-12 -Dictionary Key: shape-FD -@value 1.2476798608176409 1e-08 1e-12 -Dictionary Key: LoD -Dictionary Key: patchV-Adjoint -@value 0.3116714610927169 1e-08 1e-12 -@value 4.4183903526843844 1e-08 1e-12 -Dictionary Key: patchV-FD -@value 0.3116693801530346 1e-08 1e-12 -@value 4.4183875452326902 1e-08 1e-12 -Dictionary Key: shape-Adjoint -@value 62.9167283705385358 1e-08 1e-12 -Dictionary Key: shape-FD -@value 62.9468890552234370 1e-08 1e-12 -Dictionary Key: nonOrtho -Dictionary Key: shape-Adjoint -@value -0.0178410856249778 1e-08 1e-12 -Dictionary Key: shape-FD -@value -0.0243671058851760 1e-08 1e-12 -Dictionary Key: skewness -@value -0.0001091622840756 1e-08 1e-12 -Dictionary Key: shape-FD -@value -0.0003604027958772 1e-08 1e-12 \ No newline at end of file diff --git a/tests/runRegTests_AeroOpt.py b/tests/runRegTests_AeroOpt.py new file mode 100755 index 00000000..afd85c20 --- /dev/null +++ b/tests/runRegTests_AeroOpt.py @@ -0,0 +1,233 @@ +#!/usr/bin/env python +""" +Run Python tests for optimization integration +""" + +from mpi4py import MPI +import os +import numpy as np +from testFuncs import * + +import openmdao.api as om +from mphys.multipoint import Multipoint +from dafoam.mphys import DAFoamBuilder, OptFuncs +from mphys.scenario_aerodynamic import ScenarioAerodynamic +from pygeo.mphys import OM_DVGEOCOMP +from pygeo import geo_utils + +gcomm = MPI.COMM_WORLD + +os.chdir("./reg_test_files-main/NACA0012") +if gcomm.rank == 0: + os.system("rm -rf 0 system processor* *.bin") + os.system("cp -r 0.incompressible 0") + os.system("cp -r system.incompressible system") + os.system("cp -r constant/turbulenceProperties.sa constant/turbulenceProperties") + replace_text_in_file("system/fvSchemes", "meshWave;", "meshWaveFrozen;") + +# aero setup +U0 = 10.0 +p0 = 0.0 +A0 = 0.1 +twist0 = 3.0 +LRef = 1.0 +nuTilda0 = 4.5e-5 + +daOptions = { + "designSurfaces": ["wing"], + "solverName": "DASimpleFoam", + "primalMinResTol": 1.0e-10, + "primalMinResTolDiff": 1e4, + "writeDeformedFFDs": True, + "writeDeformedConstraints": True, + "primalBC": { + "U0": {"variable": "U", "patches": ["inout"], "value": [U0, 0.0, 0.0]}, + "p0": {"variable": "p", "patches": ["inout"], "value": [p0]}, + "nuTilda0": {"variable": "nuTilda", "patches": ["inout"], "value": [nuTilda0]}, + "useWallFunction": True, + "transport:nu": 1.5e-5, + }, + "function": { + "CD": { + "type": "force", + "source": "patchToFace", + "patches": ["wing"], + "directionMode": "parallelToFlow", + "patchVelocityInputName": "patchV", + "scale": 1.0 / (0.5 * U0 * U0 * A0), + }, + "CL": { + "type": "force", + "source": "patchToFace", + "patches": ["wing"], + "directionMode": "normalToFlow", + "patchVelocityInputName": "patchV", + "scale": 1.0 / (0.5 * U0 * U0 * A0), + }, + "skewness": { + "type": "meshQualityKS", + "source": "allCells", + "coeffKS": 10.0, + "metric": "faceSkewness", + "scale": 1.0, + }, + "nonOrtho": { + "type": "meshQualityKS", + "source": "allCells", + "coeffKS": 1.0, + "metric": "nonOrthoAngle", + "scale": 1.0, + }, + }, + "adjEqnOption": {"gmresRelTol": 1.0e-10, "pcFillLevel": 1, "jacMatReOrdering": "rcm"}, + "normalizeStates": {"U": U0, "p": U0 * U0 / 2.0, "phi": 1.0, "nuTilda": 1e-3}, + "inputInfo": { + "aero_vol_coords": {"type": "volCoord", "components": ["solver", "function"]}, + "patchV": { + "type": "patchVelocity", + "patches": ["inout"], + "flowAxis": "x", + "normalAxis": "y", + "components": ["solver", "function"], + }, + }, +} + +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]], [[0.0, 0.0, 0.1], [0.0, 0.0, 1.0]]], +} + + +class Top(Multipoint): + def setup(self): + dafoam_builder = DAFoamBuilder(daOptions, meshOptions, scenario="aerodynamic") + dafoam_builder.initialize(self.comm) + + ################################################################################ + # MPHY setup + ################################################################################ + + # ivc to keep the top level DVs + self.add_subsystem("dvs", om.IndepVarComp(), promotes=["*"]) + + # create the mesh and cruise scenario because we only have one analysis point + self.add_subsystem("mesh", dafoam_builder.get_mesh_coordinate_subsystem()) + + # add the geometry component, we dont need a builder because we do it here. + self.add_subsystem("geometry", OM_DVGEOCOMP(file="FFD/wingFFD.xyz", type="ffd")) + + self.mphys_add_scenario("cruise", ScenarioAerodynamic(aero_builder=dafoam_builder)) + + self.connect("mesh.x_aero0", "geometry.x_aero_in") + self.connect("geometry.x_aero0", "cruise.x_aero") + + self.add_subsystem("LoD", om.ExecComp("val=CL/CD")) + + def configure(self): + + # create geometric DV setup + points = self.mesh.mphys_get_surface_mesh() + + # add pointset + self.geometry.nom_add_discipline_coords("aero", points) + + # create constraint DV setup + tri_points = self.mesh.mphys_get_triangulated_surface() + self.geometry.nom_setConstraintSurface(tri_points) + + # add the dv_geo object to the builder solver. This will be used to write deformed FFDs + self.cruise.coupling.solver.add_dvgeo(self.geometry.DVGeo) + + # add the dv_con object to the builder solver. This will be used to write deformed constraints + self.cruise.coupling.solver.add_dvcon(self.geometry.DVCon) + + # geometry setup + + # Select all points + pts = self.geometry.DVGeo.getLocalIndex(0) + indexList = pts[:, :, :].flatten() + PS = geo_utils.PointSelect("list", indexList) + nShapes = self.geometry.nom_addLocalDV(dvName="shape", pointSelect=PS) + + # setup the symmetry constraint to link the y displacement between k=0 and k=1 + nFFDs_x = pts.shape[0] + nFFDs_y = pts.shape[1] + indSetA = [] + indSetB = [] + for i in range(nFFDs_x): + for j in range(nFFDs_y): + indSetA.append(pts[i, j, 0]) + indSetB.append(pts[i, j, 1]) + self.geometry.nom_addLinearConstraintsShape("linearcon", indSetA, indSetB, factorA=1.0, factorB=-1.0) + + # setup the volume and thickness constraints + leList = [[1e-4, 0.0, 1e-4], [1e-4, 0.0, 0.1 - 1e-4]] + teList = [[0.998 - 1e-4, 0.0, 1e-4], [0.998 - 1e-4, 0.0, 0.1 - 1e-4]] + self.geometry.nom_addThicknessConstraints2D("thickcon", leList, teList, nSpan=2, nChord=10) + self.geometry.nom_addVolumeConstraint("volcon", leList, teList, nSpan=2, nChord=10) + # add the LE/TE constraints + self.geometry.nom_add_LETEConstraint("lecon", volID=0, faceID="iLow", topID="k") + self.geometry.nom_add_LETEConstraint("tecon", volID=0, faceID="iHigh", topID="k") + + # add the design variables to the dvs component's output + self.dvs.add_output("shape", val=np.zeros(nShapes)) + self.dvs.add_output("patchV", val=np.array([10.0, 3.0])) + # manually connect the dvs output to the geometry and cruise + self.connect("shape", "geometry.shape") + self.connect("patchV", "cruise.patchV") + + # define the design variables to the top level + self.add_design_var("shape", lower=-0.1, upper=0.1, scaler=1.0) + self.add_design_var("patchV", lower=-50.0, upper=50.0, scaler=1.0, indices=[1]) + + # add constraints and the objective + self.connect("cruise.aero_post.CD", "LoD.CD") + self.connect("cruise.aero_post.CL", "LoD.CL") + self.add_objective("LoD.val", scaler=-1.0) + self.add_constraint("cruise.aero_post.CL", equals=0.3) + self.add_constraint("cruise.aero_post.skewness", upper=4.0) + self.add_constraint("cruise.aero_post.nonOrtho", upper=70.0) + self.add_constraint("geometry.thickcon", lower=0.5, upper=3.0, scaler=1.0) + self.add_constraint("geometry.volcon", lower=1.0, scaler=1.0) + self.add_constraint("geometry.tecon", equals=0.0, scaler=1.0, linear=True) + self.add_constraint("geometry.lecon", equals=0.0, scaler=1.0, linear=True) + self.add_constraint("geometry.linearcon", equals=0.0, scaler=1.0, linear=True) + + +prob = om.Problem() +prob.model = Top() + +prob.driver = om.pyOptSparseDriver() +prob.driver.options["optimizer"] = "IPOPT" +prob.driver.opt_settings = { + "tol": 1.0e-5, + "constr_viol_tol": 1.0e-5, + "max_iter": 2, + "print_level": 5, + "mu_strategy": "adaptive", + "limited_memory_max_history": 10, + "nlp_scaling_method": "none", + "alpha_for_y": "full", + "recalc_y": "yes", +} +prob.driver.options["debug_print"] = ["nl_cons", "objs", "desvars"] + +prob.setup(mode="rev") +om.n2(prob, show_browser=False, outfile="mphys_aero.html") + +optFuncs = OptFuncs(daOptions, prob) + +optFuncs.findFeasibleDesign( + ["cruise.aero_post.CL"], ["patchV"], designVarsComp=[1], targets=[0.3] +) + +prob.run_driver() + +if gcomm.rank == 0: + funcDict = {} + funcDict["CD"] = prob.get_val("cruise.aero_post.CD") + funcDict["CL"] = prob.get_val("cruise.aero_post.CL") + reg_write_dict(funcDict, 1e-5, 1e-10) diff --git a/tests/runRegTests_DARhoSimpleFoam.py b/tests/runRegTests_DARhoSimpleFoam.py deleted file mode 100755 index 559f2680..00000000 --- a/tests/runRegTests_DARhoSimpleFoam.py +++ /dev/null @@ -1,168 +0,0 @@ -#!/usr/bin/env python -""" -Run Python tests for optimization integration -""" - -from mpi4py import MPI -import os -import numpy as np -from testFuncs import * - -import openmdao.api as om -from mphys.multipoint import Multipoint -from dafoam.mphys import DAFoamBuilder, OptFuncs -from mphys.scenario_aerodynamic import ScenarioAerodynamic -from pygeo.mphys import OM_DVGEOCOMP - -gcomm = MPI.COMM_WORLD - -os.chdir("./reg_test_files-main/NACA0012") -if gcomm.rank == 0: - os.system("rm -rf 0 system processor* *.bin") - os.system("cp -r 0.compressible 0") - os.system("cp -r system.subsonic system") - os.system("cp -r constant/turbulenceProperties.sst constant/turbulenceProperties") - replace_text_in_file("system/fvSchemes", "meshWave;", "meshWaveFrozen;") - -# aero setup -U0 = 100.0 -p0 = 101325.0 -T0 = 300.0 -A0 = 0.1 -twist0 = 2.0 -LRef = 1.0 -nuTilda0 = 4.5e-5 - -daOptions = { - "designSurfaces": ["wing"], - "solverName": "DARhoSimpleFoam", - "primalMinResTol": 1.0e-11, - "primalMinResTolDiff": 1e4, - "primalBC": { - "U0": {"variable": "U", "patches": ["inout"], "value": [U0, 0.0, 0.0]}, - "T0": {"variable": "T", "patches": ["inout"], "value": [T0]}, - "p0": {"variable": "p", "patches": ["inout"], "value": [p0]}, - "nuTilda0": {"variable": "nuTilda", "patches": ["inout"], "value": [nuTilda0]}, - "useWallFunction": True, - }, - "function": { - "CD": { - "type": "force", - "source": "patchToFace", - "patches": ["wing"], - "directionMode": "fixedDirection", - "direction": [1.0, 0.0, 0.0], - "scale": 1.0 / (0.5 * U0 * U0 * A0), - }, - "CL": { - "type": "force", - "source": "patchToFace", - "patches": ["wing"], - "directionMode": "fixedDirection", - "direction": [0.0, 1.0, 0.0], - "scale": 1.0 / (0.5 * U0 * U0 * A0), - }, - }, - "adjEqnOption": {"gmresRelTol": 1.0e-11, "pcFillLevel": 1, "jacMatReOrdering": "rcm"}, - "normalizeStates": {"U": U0, "p": p0, "phi": 1.0, "T": T0, "nuTilda": 1e-3}, - "inputInfo": { - "aero_vol_coords": {"type": "volCoord", "components": ["solver", "function"]}, - }, -} - -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]], [[0.0, 0.0, 0.1], [0.0, 0.0, 1.0]]], -} - - -class Top(Multipoint): - def setup(self): - dafoam_builder = DAFoamBuilder(daOptions, meshOptions, scenario="aerodynamic") - dafoam_builder.initialize(self.comm) - - ################################################################################ - # MPHY setup - ################################################################################ - - # ivc to keep the top level DVs - self.add_subsystem("dvs", om.IndepVarComp(), promotes=["*"]) - - # create the mesh and cruise scenario because we only have one analysis point - self.add_subsystem("mesh", dafoam_builder.get_mesh_coordinate_subsystem()) - - # add the geometry component, we dont need a builder because we do it here. - self.add_subsystem("geometry", OM_DVGEOCOMP(file="FFD/wingFFD.xyz", type="ffd")) - - self.mphys_add_scenario("cruise", ScenarioAerodynamic(aero_builder=dafoam_builder)) - - self.connect("mesh.x_aero0", "geometry.x_aero_in") - self.connect("geometry.x_aero0", "cruise.x_aero") - - self.add_subsystem("LoD", om.ExecComp("val=CL/CD")) - - def configure(self): - - # create geometric DV setup - points = self.mesh.mphys_get_surface_mesh() - - # add pointset - self.geometry.nom_add_discipline_coords("aero", points) - - # geometry setup - self.geometry.nom_addRefAxis(name="wingAxis", xFraction=0.25, alignIndex="k") - - # Set up global design variables. We dont change the root twist - def twist(val, geo): - for i in range(2): - geo.rot_z["wingAxis"].coef[i] = -val[0] - - self.geometry.nom_addGlobalDV(dvName="twist", value=np.ones(1) * twist0, func=twist) - - # add the design variables to the dvs component's output - self.dvs.add_output("twist", val=np.ones(1) * twist0) - # manually connect the dvs output to the geometry and cruise - self.connect("twist", "geometry.twist") - - # define the design variables to the top level - self.add_design_var("twist", lower=-10.0, upper=10.0, scaler=1.0) - - # add constraints and the objective - self.add_objective("cruise.aero_post.CD", scaler=1.0) - self.add_constraint("cruise.aero_post.CL", equals=0.3) - - -prob = om.Problem() -prob.model = Top() - -prob.setup(mode="rev") -om.n2(prob, show_browser=False, outfile="mphys_aero.html") - -# optFuncs = OptFuncs(daOptions, prob) - -# verify the total derivatives against the finite-difference -prob.run_model() -results = prob.check_totals( - of=["cruise.aero_post.CD", "cruise.aero_post.CL"], - wrt=["twist"], - compact_print=True, - step=1e-3, - form="central", - step_calc="abs", -) - -if gcomm.rank == 0: - funcDict = {} - funcDict["CD"] = prob.get_val("cruise.aero_post.CD") - funcDict["CL"] = prob.get_val("cruise.aero_post.CL") - derivDict = {} - derivDict["CD"] = {} - derivDict["CD"]["twist-Adjoint"] = results[("cruise.aero_post.CD", "twist")]["J_fwd"][0] - derivDict["CD"]["twist-FD"] = results[("cruise.aero_post.CD", "twist")]["J_fd"][0] - derivDict["CL"] = {} - derivDict["CL"]["twist-Adjoint"] = results[("cruise.aero_post.CL", "twist")]["J_fwd"][0] - derivDict["CL"]["twist-FD"] = results[("cruise.aero_post.CL", "twist")]["J_fd"][0] - reg_write_dict(funcDict, 1e-10, 1e-12) - reg_write_dict(derivDict, 1e-8, 1e-12) diff --git a/tests/runUnitTests_primal.py b/tests/runUnitTests_primal.py new file mode 100755 index 00000000..25a57ce2 --- /dev/null +++ b/tests/runUnitTests_primal.py @@ -0,0 +1,67 @@ +#!/usr/bin/env python +""" +Run Python tests for optimization integration +""" + +from mpi4py import MPI +from dafoam import PYDAFOAM +import os + +gcomm = MPI.COMM_WORLD + +os.chdir("./reg_test_files-main/ConvergentChannel") +if gcomm.rank == 0: + os.system("rm -rf processor* *.bin") + os.system("cp -r 0.compressible/* 0/") + os.system("cp -r system.subsonic/* system/") + os.system("cp -r constant/turbulenceProperties.ke constant/turbulenceProperties") + +# aero setup +U0 = 10.0 + +daOptions = {"solverName": "DARhoSimpleFoam", "primalMinResTol": 1e-12, "printDAOptions": False} + +DASolver = PYDAFOAM(options=daOptions, comm=gcomm) +DASolver() +states = DASolver.getStates() +norm = np.linalg.norm(states) +norm = gcomm.allreduce(norm, op=MPI.SUM) +if ((3787018.3272404578 - norm) / norm) > 1e-10: + print("ke test failed!") + exit(1) +else: + print("ke test passed!") + +if gcomm.rank == 0: + os.system("rm -rf processor* *.bin") + os.system("cp -r 0.compressible/* 0/") + os.system("cp -r system.subsonic/* system/") + os.system("cp -r constant/turbulenceProperties.kw constant/turbulenceProperties") + +DASolver = PYDAFOAM(options=daOptions, comm=gcomm) +DASolver() +states = DASolver.getStates() +norm = np.linalg.norm(states) +norm = gcomm.allreduce(norm, op=MPI.SUM) +if ((3787032.628925756 - norm) / norm) > 1e-10: + print("kw test failed!") + exit(1) +else: + print("kw test passed!") + +if gcomm.rank == 0: + os.system("rm -rf processor* *.bin") + os.system("cp -r 0.compressible/* 0/") + os.system("cp -r system.subsonic/* system/") + os.system("cp -r constant/turbulenceProperties.sstlm constant/turbulenceProperties") + +DASolver = PYDAFOAM(options=daOptions, comm=gcomm) +DASolver() +states = DASolver.getStates() +norm = np.linalg.norm(states) +norm = gcomm.allreduce(norm, op=MPI.SUM) +if ((3787217.019831411 - norm) / norm) > 1e-10: + print("SSTLM test failed!") + exit(1) +else: + print("SSTLM test passed!")