diff --git a/adflow/__init__.py b/adflow/__init__.py index 065b095d5..1320315f0 100644 --- a/adflow/__init__.py +++ b/adflow/__init__.py @@ -6,15 +6,3 @@ from .pyADflow_C import ADFLOW_C from .oversetCheck import OversetCheck from .checkZipper import checkZipper - - -try: - import openmdao -except ImportError as err: - openmdao = None - - if MPI.COMM_WORLD.rank == 0: - print("Warning: OpenMDAO dependency is not installed. OM_ADFLOW wrapper will not be active") - -if openmdao is not None: - from .om_adflow import OM_ADFLOW diff --git a/adflow/mphys/__init__.py b/adflow/mphys/__init__.py new file mode 100644 index 000000000..7cb7b067b --- /dev/null +++ b/adflow/mphys/__init__.py @@ -0,0 +1,3 @@ +from .mphys_adflow import ADflowBuilder + +__all__ = ["ADflowBuilder"] diff --git a/adflow/mphys/mphys_adflow.py b/adflow/mphys/mphys_adflow.py new file mode 100644 index 000000000..c9f6679d7 --- /dev/null +++ b/adflow/mphys/mphys_adflow.py @@ -0,0 +1,1270 @@ +import numpy as np +from pprint import pprint as pp + +from adflow import ADFLOW +from idwarp import USMesh + +from openmdao.api import Group, ImplicitComponent, ExplicitComponent, AnalysisError +from mphys.builder import Builder + +from .om_utils import get_dvs_and_cons + + +class ADflowMesh(ExplicitComponent): + """ + Component to get the partitioned initial surface mesh coordinates + + """ + + def initialize(self): + self.options.declare("aero_solver", recordable=False) + + def setup(self): + + self.aero_solver = self.options["aero_solver"] + + self.x_a0 = self.aero_solver.getSurfaceCoordinates(includeZipper=False).flatten(order="C") + + coord_size = self.x_a0.size + self.add_output( + "x_aero0", + distributed=True, + shape=coord_size, + desc="initial aerodynamic surface node coordinates", + tags=["mphys_coordinates"], + ) + + def mphys_add_coordinate_input(self): + self.add_input( + "x_aero0_points", distributed=True, shape_by_conn=True, desc="aerodynamic surface with geom changes" + ) + + # return the promoted name and coordinates + return "x_aero0_points", self.x_a0 + + def mphys_get_surface_mesh(self): + return self.x_a0 + + def mphys_get_triangulated_surface(self, groupName=None): + # this is a list of lists of 3 points + # p0, v1, v2 + + return self._getTriangulatedMeshSurface() + + def _getTriangulatedMeshSurface(self, groupName=None, **kwargs): + """ + This function returns a trianguled verision of the surface + mesh on all processors. The intent is to use this for doing + constraints in DVConstraints. + + Returns + ------- + surf : list + List of points and vectors describing the surface. This may + be passed directly to DVConstraint setSurface() function. + """ + + if groupName is None: + groupName = self.aero_solver.allWallsGroup + + # Obtain the points and connectivity for the specified + # groupName + pts = self.aero_solver.comm.allgather(self.aero_solver.getSurfaceCoordinates(groupName, **kwargs)) + conn, faceSizes = self.aero_solver.getSurfaceConnectivity(groupName) + conn = np.array(conn).flatten() + conn = self.aero_solver.comm.allgather(conn) + faceSizes = self.aero_solver.comm.allgather(faceSizes) + # Triangle info...point and two vectors + p0 = [] + v1 = [] + v2 = [] + + # loop over the faces + for iProc in range(len(faceSizes)): + + connCounter = 0 + for iFace in range(len(faceSizes[iProc])): + # Get the number of nodes on this face + faceSize = faceSizes[iProc][iFace] + faceNodes = conn[iProc][connCounter : connCounter + faceSize] + # Start by getting the centerpoint + ptSum = [0, 0, 0] + for i in range(faceSize): + # idx = ptCounter+i + idx = faceNodes[i] + ptSum += pts[iProc][idx] + + avgPt = ptSum / faceSize + + # Now go around the face and add a triangle for each adjacent pair + # of points. This assumes an ordered connectivity from the + # meshwarping + for i in range(faceSize): + idx = faceNodes[i] + p0.append(avgPt) + v1.append(pts[iProc][idx] - avgPt) + if i < (faceSize - 1): + idxp1 = faceNodes[i + 1] + v2.append(pts[iProc][idxp1] - avgPt) + else: + # wrap back to the first point for the last element + idx0 = faceNodes[0] + v2.append(pts[iProc][idx0] - avgPt) + + # Now increment the connectivity + connCounter += faceSize + return [p0, v1, v2] + + def compute(self, inputs, outputs): + if "x_aero0_points" in inputs: + outputs["x_aero0"] = inputs["x_aero0_points"] + else: + outputs["x_aero0"] = self.x_a0 + + def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): + if mode == "fwd": + if "x_aero0_points" in d_inputs: + d_outputs["x_aero0"] += d_inputs["x_aero0_points"] + elif mode == "rev": + if "x_aero0_points" in d_inputs: + d_inputs["x_aero0_points"] += d_outputs["x_aero0"] + + +class ADflowWarper(ExplicitComponent): + """ + OpenMDAO component that wraps the warping. + + """ + + def initialize(self): + self.options.declare("aero_solver", recordable=False) + # self.options.declare('use_OM_KSP', default=False, types=bool, + # desc="uses OpenMDAO's PestcKSP linear solver with ADflow's preconditioner to solve the adjoint.") + + def setup(self): + # self.set_check_partial_options(wrt='*',directional=True) + + self.solver = self.options["aero_solver"] + solver = self.solver + + # self.ap_vars,_ = get_dvs_and_cons(ap=ap) + + # state inputs and outputs + local_volume_coord_size = solver.mesh.getSolverGrid().size + + self.add_input("x_aero", distributed=True, shape_by_conn=True, tags=["mphys_coupling"]) + self.add_output("adflow_vol_coords", distributed=True, shape=local_volume_coord_size, tags=["mphys_coupling"]) + + # self.declare_partials(of='adflow_vol_coords', wrt='x_aero') + + def compute(self, inputs, outputs): + + solver = self.solver + + x_a = inputs["x_aero"].reshape((-1, 3)) + solver.setSurfaceCoordinates(x_a) + solver.updateGeometryInfo() + outputs["adflow_vol_coords"] = solver.mesh.getSolverGrid() + + def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): + + if mode == "fwd": + if "adflow_vol_coords" in d_outputs: + if "x_aero" in d_inputs: + dxS = d_inputs["x_aero"] + dxV = self.solver.mesh.warpDerivFwd(dxS) + d_outputs["adflow_vol_coords"] += dxV + + elif mode == "rev": + if "adflow_vol_coords" in d_outputs: + if "x_aero" in d_inputs: + dxV = d_outputs["adflow_vol_coords"] + self.solver.mesh.warpDeriv(dxV) + dxS = self.solver.mesh.getdXs() + dxS = self.solver.mapVector( + dxS, self.solver.meshFamilyGroup, self.solver.designFamilyGroup, includeZipper=False + ) + d_inputs["x_aero"] += dxS.flatten() + + +class ADflowSolver(ImplicitComponent): + """ + OpenMDAO component that wraps the ADflow flow solver + + """ + + def initialize(self): + self.options.declare("aero_solver", recordable=False) + # self.options.declare('use_OM_KSP', default=False, types=bool, + # desc="uses OpenMDAO's PestcKSP linear solver with ADflow's preconditioner to solve the adjoint.") + self.options.declare("restart_failed_analysis", default=False) + self.options.declare("err_on_convergence_fail", default=False) + + # testing flag used for unit-testing to prevent the call to actually solve + # NOT INTENDED FOR USERS!!! FOR TESTING ONLY + self._do_solve = True + self.analysis_error_on_failure = True + + def setup(self): + # self.set_check_partial_options(wrt='*',directional=True) + + self.restart_failed_analysis = self.options["restart_failed_analysis"] + self.err_on_convergence_fail = self.options["err_on_convergence_fail"] + self.solver = self.options["aero_solver"] + solver = self.solver + + # this is the solution counter for failed solution outputs. + # the converged solutions are written by the adflow functionals group + self.solution_counter = 0 + + # flag to keep track if the current solution started from a clean restart, + # or it was restarted from the previous converged state. + self.cleanRestart = True + + # state inputs and outputs + local_state_size = solver.getStateSize() + + self.add_input("adflow_vol_coords", distributed=True, shape_by_conn=True, tags=["mphys_coupling"]) + self.add_output("adflow_states", distributed=True, shape=local_state_size, tags=["mphys_coupling"]) + + # self.declare_partials(of='adflow_states', wrt='*') + + def _set_ap(self, inputs, print_dict=True): + tmp = {} + for (args, _kwargs) in self.ap_vars: + name = args[0] + tmp[name] = inputs[name] + + # enable if you want to print all aero dv inputs + # if self.comm.rank == 0: + # print('aero dv inputs:') + # pp(tmp) + + self.ap.setDesignVars(tmp) + if self.comm.rank == 0 and print_dict: + pp(tmp) + + def set_ap(self, ap): + # this is the external function to set the ap to this component + self.ap = ap + + self.ap_vars, _ = get_dvs_and_cons(ap=ap) + + # parameter inputs + if self.comm.rank == 0: + print("adding ap var inputs") + for (args, kwargs) in self.ap_vars: + name = args[0] + size = args[1] + self.add_input( + name, distributed=False, shape=size, val=kwargs["value"], units=kwargs["units"], tags=["mphys_input"] + ) + if self.comm.rank == 0: + print("%s (%s)" % (name, kwargs["units"])) + + def _set_states(self, outputs): + self.solver.setStates(outputs["adflow_states"]) + + def apply_nonlinear(self, inputs, outputs, residuals): + + solver = self.solver + + self._set_states(outputs) + self._set_ap(inputs, print_dict=False) + + ap = self.ap + + # Set the warped mesh + # solver.mesh.setSolverGrid(inputs['adflow_vol_coords']) + # ^ This call does not exist. Assume the mesh hasn't changed since the last call to the warping comp for now + # TODO we must fix this. we need to put in the modified volume coordinates + + # flow residuals + residuals["adflow_states"] = solver.getResidual(ap) + + def solve_nonlinear(self, inputs, outputs): + solver = self.solver + ap = self.ap + if self._do_solve: + + # Set the warped mesh + # solver.mesh.setSolverGrid(inputs['adflow_vol_coords']) + # ^ This call does not exist. Assume the mesh hasn't changed since the last call to the warping comp for now + + self._set_ap(inputs) + ap.solveFailed = False # might need to clear this out? + ap.fatalFail = False + + # do not write solution files inside the solver loop + solver(ap, writeSolution=False) + + if ap.fatalFail: + if self.comm.rank == 0: + print("###############################################################") + print("# Solve Fatal Fail. Analysis Error") + print("###############################################################") + + raise AnalysisError("ADFLOW Solver Fatal Fail") + + if ap.solveFailed: + if self.restart_failed_analysis: # the mesh was fine, but it didn't converge + # if the previous iteration was already a clean restart, dont try again + if self.cleanRestart: + if self.comm.rank == 0: + print("###############################################################") + print("# This was a clean restart. Will not try another one.") + print("###############################################################") + + # write the solution so that we can diagnose + solver.writeSolution(baseName="analysis_fail", number=self.solution_counter) + self.solution_counter += 1 + + if self.analysis_error_on_failure: + solver.resetFlow(ap) + self.cleanRestart = True + raise AnalysisError("ADFLOW Solver Fatal Fail") + + # the previous iteration restarted from another solution, so we can try again + # with a re-set flowfield for the initial guess. + else: + if self.comm.rank == 0: + print("###############################################################") + print("# Solve Failed, attempting a clean restart!") + print("###############################################################") + + # write the solution so that we can diagnose + solver.writeSolution(baseName="analysis_fail", number=self.solution_counter) + self.solution_counter += 1 + + ap.solveFailed = False + ap.fatalFail = False + solver.resetFlow(ap) + solver(ap, writeSolution=False) + + if ap.solveFailed or ap.fatalFail: # we tried, but there was no saving it + if self.comm.rank == 0: + print("###############################################################") + print("# Clean Restart failed. There is no saving this one!") + print("###############################################################") + + # write the solution so that we can diagnose + solver.writeSolution(baseName="analysis_fail", number=self.solution_counter) + self.solution_counter += 1 + + if self.analysis_error_on_failure: + # re-set the flow for the next iteration: + solver.resetFlow(ap) + # set the reset flow flag + self.cleanRestart = True + raise AnalysisError("ADFLOW Solver Fatal Fail") + + # see comment for the same flag below + else: + self.cleanRestart = False + + elif self.err_on_convergence_fail: + # the solve failed but we dont want to re-try. We also want to raise an analysis error + if self.comm.rank == 0: + print("###############################################################") + print("# Solve Failed, not attempting a clean restart") + print("###############################################################") + raise AnalysisError("ADFLOW Solver Fatal Fail") + + else: + # the solve failed, but we dont want to restart or raise an error, we will just let this one pass + pass + + # solve did not fail, therefore we will re-use this converged flowfield for the next iteration. + # change the flag so that if the next iteration fails with current initial guess, it can retry + # with a clean restart + else: + self.cleanRestart = False + + outputs["adflow_states"] = solver.getStates() + + def linearize(self, inputs, outputs, residuals): + solver = self.solver + ap = self.ap + + self._set_ap(inputs, print_dict=False) + self._set_states(outputs) + + # check if we changed APs, then we have to do a bunch of updates + if ap != solver.curAP: + # AP is changed, so we have to update the AP and + # run a residual to make sure all intermediate vairables are up to date + # we assume the AP has the last converged state information, + # which is automatically set in the getResidual call + solver.getResidual(ap) + + def apply_linear(self, inputs, outputs, d_inputs, d_outputs, d_residuals, mode): + + solver = self.solver + ap = self.ap + + self._set_ap(inputs, print_dict=False) + self._set_states(outputs) + + # check if we changed APs, then we have to do a bunch of updates + if ap != solver.curAP: + # AP is changed, so we have to update the AP and + # run a residual to make sure all intermediate vairables are up to date + # we assume the AP has the last converged state information, + # which is automatically set in the getResidual call + solver.getResidual(ap) + + if mode == "fwd": + if "adflow_states" in d_residuals: + xDvDot = {} + for var_name in d_inputs: + xDvDot[var_name] = d_inputs[var_name] + if "adflow_vol_coords" in d_inputs: + xVDot = d_inputs["adflow_vol_coords"] + else: + xVDot = None + if "adflow_states" in d_outputs: + wDot = d_outputs["adflow_states"] + else: + wDot = None + + dwdot = solver.computeJacobianVectorProductFwd( + xDvDot=xDvDot, xVDot=xVDot, wDot=wDot, residualDeriv=True + ) + d_residuals["adflow_states"] += dwdot + + elif mode == "rev": + if "adflow_states" in d_residuals: + resBar = d_residuals["adflow_states"] + + wBar, xVBar, xDVBar = solver.computeJacobianVectorProductBwd( + resBar=resBar, wDeriv=True, xVDeriv=True, xDvDeriv=False, xDvDerivAero=True + ) + + if "adflow_states" in d_outputs: + d_outputs["adflow_states"] += wBar + + if "adflow_vol_coords" in d_inputs: + d_inputs["adflow_vol_coords"] += xVBar + + for dv_name, dv_bar in xDVBar.items(): + if dv_name in d_inputs: + d_inputs[dv_name] += dv_bar.flatten() + + def solve_linear(self, d_outputs, d_residuals, mode): + solver = self.solver + ap = self.ap + + # check if we changed APs, then we have to do a bunch of updates + if ap != solver.curAP: + # AP is changed, so we have to update the AP and + # run a residual to make sure all intermediate vairables are up to date + # we assume the AP has the last converged state information, + # which is automatically set in the getResidual call + solver.getResidual(ap) + + # the adjoint might not be set up regardless if we changed APs + # this is because the first call with any AP will not have this set up, so we have to check + # if we changed APs, then we also freed adjoint memory, + # and then again we would need to setup adjoint again + # finally, we generally want to avoid extra calls here + # because this routine can be call multiple times back to back in a LBGS solver. + if not solver.adjointSetup: + solver._setupAdjoint() + + if self.comm.rank == 0: + print("Solving linear in mphys_adflow", flush=True) + if mode == "fwd": + d_outputs["adflow_states"] = solver.solveDirectForRHS(d_residuals["adflow_states"]) + elif mode == "rev": + # d_residuals['adflow_states'] = solver.solveAdjointForRHS(d_outputs['adflow_states']) + solver.adflow.adjointapi.solveadjoint(d_outputs["adflow_states"], d_residuals["adflow_states"], True) + + return True, 0, 0 + + +class ADflowForces(ExplicitComponent): + """ + OpenMDAO component that wraps force integration + + """ + + def initialize(self): + self.options.declare("aero_solver", recordable=False) + + def setup(self): + # self.set_check_partial_options(wrt='*',directional=True) + + self.solver = self.options["aero_solver"] + solver = self.solver + + self.add_input("adflow_vol_coords", distributed=True, shape_by_conn=True, tags=["mphys_coupling"]) + self.add_input("adflow_states", distributed=True, shape_by_conn=True, tags=["mphys_coupling"]) + + local_surface_coord_size = solver.mesh.getSurfaceCoordinates().size + self.add_output("f_aero", distributed=True, shape=local_surface_coord_size, tags=["mphys_coupling"]) + + # self.declare_partials(of='f_aero', wrt='*') + + def _set_ap(self, inputs, print_dict=True): + tmp = {} + for (args, _kwargs) in self.ap_vars: + name = args[0] + tmp[name] = inputs[name] + + self.ap.setDesignVars(tmp) + if self.comm.rank == 0 and print_dict: + pp(tmp) + + def set_ap(self, ap): + # this is the external function to set the ap to this component + self.ap = ap + + self.ap_vars, _ = get_dvs_and_cons(ap=ap) + + # parameter inputs + # if self.comm.rank == 0: + # print('adding ap var inputs:') + for (args, kwargs) in self.ap_vars: + name = args[0] + size = args[1] + self.add_input(name, distributed=False, shape=size, units=kwargs["units"], tags=["mphys_input"]) + # if self.comm.rank == 0: + # print('%s (%s)'%(name, kwargs['units'])) + + def _set_states(self, inputs): + self.solver.setStates(inputs["adflow_states"]) + + def compute(self, inputs, outputs): + + solver = self.solver + + self._set_ap(inputs) + + # Set the warped mesh + # solver.mesh.setSolverGrid(inputs['adflow_vol_coords']) + # ^ This call does not exist. Assume the mesh hasn't changed since the last call to the warping comp for now + # TODO we must fix this. we need to put in the modified volume coordinates + self._set_states(inputs) + + outputs["f_aero"] = solver.getForces().flatten(order="C") + + def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): + + solver = self.solver + ap = self.ap + + self._set_ap(inputs, print_dict=False) + + # check if we changed APs, then we have to do a bunch of updates + if ap != solver.curAP: + # AP is changed, so we have to update the AP and + # run a residual to make sure all intermediate vairables are up to date + # we assume the AP has the last converged state information, + # which is automatically set in the getResidual call + solver.getResidual(ap) + + if mode == "fwd": + if "f_aero" in d_outputs: + xDvDot = {} + for var_name in d_inputs: + xDvDot[var_name] = d_inputs[var_name] + if "adflow_states" in d_inputs: + wDot = d_inputs["adflow_states"] + else: + wDot = None + if "adflow_vol_coords" in d_inputs: + xVDot = d_inputs["adflow_vol_coords"] + else: + xVDot = None + if not (xVDot is None and wDot is None): + dfdot = solver.computeJacobianVectorProductFwd(xDvDot=xDvDot, xVDot=xVDot, wDot=wDot, fDeriv=True) + d_outputs["f_aero"] += dfdot.flatten() + + elif mode == "rev": + if "f_aero" in d_outputs: + fBar = d_outputs["f_aero"] + + wBar, xVBar, xDVBar = solver.computeJacobianVectorProductBwd( + fBar=fBar, wDeriv=True, xVDeriv=True, xDvDeriv=False, xDvDerivAero=True + ) + + if "adflow_vol_coords" in d_inputs: + d_inputs["adflow_vol_coords"] += xVBar + if "adflow_states" in d_inputs: + d_inputs["adflow_states"] += wBar + + for dv_name, dv_bar in xDVBar.items(): + if dv_name in d_inputs: + d_inputs[dv_name] += dv_bar.flatten() + + +class AdflowHeatTransfer(ExplicitComponent): + """ + OpenMDAO component that wraps heat transfer integration + + """ + + def initialize(self): + self.options.declare("aero_solver") + + def setup(self): + # self.set_check_partial_options(wrt='*',directional=True) + + self.solver = self.options["aero_solver"] + solver = self.solver + + local_nodes, _ = solver._getSurfaceSize(solver.allIsothermalWallsGroup) + + self.add_input("adflow_vol_coords", distributed=True, shape_by_conn=True, tags=["mphys_coupling"]) + self.add_input("adflow_states", distributed=True, shape_by_conn=True, tags=["mphys_coupling"]) + + self.add_output( + "q_convect", + distributed=True, + val=np.ones(local_nodes) * -499, + shape=local_nodes, + units="W/m**2", + tags=["mphys_coupling"], + ) + + # self.declare_partials(of='f_aero', wrt='*') + + def _set_ap(self, inputs): + tmp = {} + for (args, _kwargs) in self.ap_vars: + name = args[0] + tmp[name] = inputs[name] + + self.ap.setDesignVars(tmp) + + def set_ap(self, ap): + # this is the external function to set the ap to this component + self.ap = ap + + self.ap_vars, _ = get_dvs_and_cons(ap=ap) + + # parameter inputs + if self.comm.rank == 0: + print("adding ap var inputs") + for (args, kwargs) in self.ap_vars: + name = args[0] + size = args[1] + self.add_input(name, distributed=False, shape=size, units=kwargs["units"], tags=["mphys_input"]) + if self.comm.rank == 0: + print(name) + + def _set_states(self, inputs): + self.solver.setStates(inputs["adflow_states"]) + + def compute(self, inputs, outputs): + + solver = self.solver + + ## already done by solver + self._set_ap(inputs) + + # Set the warped mesh + # solver.mesh.setSolverGrid(inputs['adflow_vol_coords']) + # ^ This call does not exist. Assume the mesh hasn't changed since the last call to the warping comp for now + # TODO we must fix this. we need to put in the modified volume coordinates + + # + # self._set_states(inputs) + + outputs["q_convect"] = solver.getHeatFluxes().flatten(order="C") + # print() + + def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): + + solver = self.solver + ap = self.options["ap"] + + self._set_ap(inputs) + + # check if we changed APs, then we have to do a bunch of updates + if ap != solver.curAP: + # AP is changed, so we have to update the AP and + # run a residual to make sure all intermediate vairables are up to date + # we assume the AP has the last converged state information, + # which is automatically set in the getResidual call + solver.getResidual(ap) + + if mode == "fwd": + if "q_convect" in d_outputs: + xDvDot = {} + for var_name in d_inputs: + xDvDot[var_name] = d_inputs[var_name] + if "adflow_states" in d_inputs: + wDot = d_inputs["adflow_states"] + else: + wDot = None + if "adflow_vol_coords" in d_inputs: + xVDot = d_inputs["adflow_vol_coords"] + else: + xVDot = None + if not (xVDot is None and wDot is None): + dhfdot = solver.computeJacobianVectorProductFwd(xDvDot=xDvDot, xVDot=xVDot, wDot=wDot, hfDeriv=True) + dhfdot_map = np.zeros((dhfdot.size, 3)) + dhfdot_map[:, 0] = dhfdot.flatten() + dhfdot_map = self.solver.mapVector( + dhfdot_map, self.solver.allWallsGroup, self.solver.allIsothermalWallsGroup + ) + dhfdot = dhfdot_map[:, 0] + d_outputs["q_convect"] += dhfdot + + elif mode == "rev": + if "q_convect" in d_outputs: + hfBar = d_outputs["q_convect"] + + hfBar_map = np.zeros((hfBar.size, 3)) + hfBar_map[:, 0] = hfBar.flatten() + hfBar_map = self.solver.mapVector( + hfBar_map, self.solver.allIsothermalWallsGroup, self.solver.allWallsGroup + ) + hfBar = hfBar_map[:, 0] + + wBar, xVBar, xDVBar = solver.computeJacobianVectorProductBwd( + hfBar=hfBar, wDeriv=True, xVDeriv=True, xDvDeriv=True + ) + + if "adflow_vol_coords" in d_inputs: + d_inputs["adflow_vol_coords"] += xVBar + if "adflow_states" in d_inputs: + d_inputs["adflow_states"] += wBar + + for dv_name, dv_bar in xDVBar.items(): + if dv_name in d_inputs: + d_inputs[dv_name] += dv_bar.flatten() + + +FUNCS_UNITS = { + "mdot": "kg/s", + "mavgptot": "Pa", + "mavgps": "Pa", + "aavgptot": "Pa", + "aavgps": "Pa", + "mavgttot": "degK", + "mavgvx": "m/s", + "mavgvy": "m/s", + "mavgvz": "m/s", + "drag": "N", + "lift": "N", + "dragpressure": "N", + "dragviscous": "N", + "dragmomentum": "N", + "fx": "N", + "fy": "N", + "fz": "N", + "forcexpressure": "N", + "forceypressure": "N", + "forcezpressure": "N", + "forcexviscous": "N", + "forceyviscous": "N", + "forcezviscous": "N", + "forcexmomentum": "N", + "forceymomentum": "N", + "forcezmomentum": "N", + "flowpower": "W", + "area": "m**2", +} + + +class ADflowFunctions(ExplicitComponent): + def initialize(self): + self.options.declare("aero_solver", recordable=False) + # flag to automatically add the AP functions as output + self.options.declare("ap_funcs", default=True) + self.options.declare("write_solution", default=True) + + # testing flag used for unit-testing to prevent the call to actually solve + # NOT INTENDED FOR USERS!!! FOR TESTING ONLY + self._do_solve = True + + self.extra_funcs = None + + def setup(self): + + self.solver = self.options["aero_solver"] + self.ap_funcs = self.options["ap_funcs"] + self.write_solution = self.options["write_solution"] + # self.set_check_partial_options(wrt='*',directional=True) + self.solution_counter = 0 + + self.add_input("adflow_vol_coords", distributed=True, shape_by_conn=True, tags=["mphys_coupling"]) + self.add_input("adflow_states", distributed=True, shape_by_conn=True, tags=["mphys_coupling"]) + + # self.declare_partials(of=f_name, wrt='*') + + def _set_ap(self, inputs, print_dict=True): + tmp = {} + for (args, _kwargs) in self.ap_vars: + name = args[0] + tmp[name] = inputs[name] + + self.ap.setDesignVars(tmp) + # self.options['solver'].setAeroProblem(self.options['ap']) + if self.comm.rank == 0 and print_dict: + pp(tmp) + + def mphys_set_ap(self, ap): + # this is the external function to set the ap to this component + self.ap = ap + self.ap_vars, _ = get_dvs_and_cons(ap=ap) + + # parameter inputs + # if self.comm.rank == 0: + # print('adding ap var inputs:') + for (args, kwargs) in self.ap_vars: + name = args[0] + size = args[1] + self.add_input(name, distributed=False, shape=size, units=kwargs["units"], tags=["mphys_input"]) + # if self.comm.rank == 0: + # print('%s with units %s'%(name, kwargs['units'])) + + if self.ap_funcs: + if self.comm.rank == 0: + print("adding adflow funcs as output:") + + for f_name in sorted(list(self.ap.evalFuncs)): + # get the function type. this is the first word before the first underscore + f_type = f_name.split("_")[0] + + # check if we have a unit defined for this + if f_type in FUNCS_UNITS: + units = FUNCS_UNITS[f_type] + else: + units = None + + # print the function name and units + if self.comm.rank == 0: + print("%s (%s)" % (f_name, units)) + + self.add_output(f_name, distributed=False, shape=1, units=units, tags=["mphys_result"]) + + # def mphys_add_prop_funcs(self, prop_funcs): + # save this list + # self.extra_funcs = prop_funcs + + # if self.comm.rank == 0: + # print("adding adflow funcs as propulsion output:", prop_funcs) + + # call the add_funcs function + # self.mphys_add_funcs(prop_funcs) + + def mphys_add_funcs(self, funcs): + + self.extra_funcs = funcs + + # loop over the functions here and create the output + for f_name in funcs: + # get the function type. this is the first word before the first underscore + f_type = f_name.split("_")[0] + + # check if we have a unit defined for this + if f_type in FUNCS_UNITS: + units = FUNCS_UNITS[f_type] + else: + units = None + + # print the function name and units + # if self.comm.rank == 0: + # print("%s (%s)" % (f_name, units)) + + self.add_output(f_name, distributed=False, shape=1, units=units, tags=["mphys_result"]) + + def _set_states(self, inputs): + self.solver.setStates(inputs["adflow_states"]) + + def _get_func_name(self, name): + return "%s_%s" % (self.ap.name, name.lower()) + + def nom_write_solution(self, **kwargs): + # this writes the solution files and is callable from outside openmdao call routines + solver = self.solver + ap = self.ap + + # re-set the AP so that we are sure state is updated + solver.setAeroProblem(ap) + + # write the solution files. Internally, this checks the + # types of solution files specified in the options and + # only outsputs these + solver.writeSolution(number=self.solution_counter, **kwargs) + self.solution_counter += 1 + + def compute(self, inputs, outputs): + solver = self.solver + # print('funcs compute') + # actually setting things here triggers some kind of reset, so we only do it if you're actually solving + if self._do_solve: + self._set_ap(inputs) + # Set the warped mesh + # solver.mesh.setSolverGrid(inputs['adflow_vol_coords']) + # ^ This call does not exist. Assume the mesh hasn't changed since the last call to the warping comp for now + self._set_states(inputs) + + if self.write_solution: + # write the solution files. Internally, this checks the + # types of solution files specified in the options and + # only outsputs these + solver.writeSolution(number=self.solution_counter) + self.solution_counter += 1 + + funcs = {} + + if self.ap_funcs: + # without the sorted, each proc might get a different order... + eval_funcs = sorted(list(self.ap.evalFuncs)) + solver.evalFunctions(self.ap, funcs, evalFuncs=eval_funcs) + + for name in self.ap.evalFuncs: + f_name = self._get_func_name(name) + if f_name in funcs: + outputs[name.lower()] = funcs[f_name] + + if self.extra_funcs is not None: + # also do the prop + solver.evalFunctions(self.ap, funcs, evalFuncs=self.extra_funcs) + for name in self.extra_funcs: + f_name = self._get_func_name(name) + if f_name in funcs: + outputs[name.lower()] = funcs[f_name] + + def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): + solver = self.solver + ap = self.ap + self._set_ap(inputs, print_dict=False) + + # check if we changed APs, then we have to do a bunch of updates + if ap != solver.curAP: + # AP is changed, so we have to update the AP and + # run a residual to make sure all intermediate vairables are up to date + # we assume the AP has the last converged state information, + # which is automatically set in the getResidual call + solver.getResidual(ap) + + if mode == "fwd": + xDvDot = {} + for key in ap.DVs: + if key in d_inputs: + mach_name = key.split("_")[0] + xDvDot[mach_name] = d_inputs[key] + + if "adflow_states" in d_inputs: + wDot = d_inputs["adflow_states"] + else: + wDot = None + + if "adflow_vol_coords" in d_inputs: + xVDot = d_inputs["adflow_vol_coords"] + else: + xVDot = None + + funcsdot = solver.computeJacobianVectorProductFwd(xDvDot=xDvDot, xVDot=xVDot, wDot=wDot, funcDeriv=True) + + for name in funcsdot: + func_name = name.lower() + if name in d_outputs: + d_outputs[name] += funcsdot[func_name] + + elif mode == "rev": + funcsBar = {} + + if self.ap_funcs: + for name in self.ap.evalFuncs: + func_name = name.lower() + + # we have to check for 0 here, so we don't include any unnecessary variables in funcsBar + # becasue it causes ADflow to do extra work internally even if you give it extra variables, even if the seed is 0 + if func_name in d_outputs and d_outputs[func_name] != 0.0: + funcsBar[func_name] = d_outputs[func_name][0] + # this stuff is fixed now. no need to divide + # funcsBar[func_name] = d_outputs[func_name][0] / self.comm.size + # print(self.comm.rank, func_name, funcsBar[func_name]) + + # also do the same for prop functions + if self.extra_funcs is not None: + for name in self.extra_funcs: + func_name = name.lower() + if func_name in d_outputs and d_outputs[func_name] != 0.0: + funcsBar[func_name] = d_outputs[func_name][0] + + # print(funcsBar, flush=True) + + wBar = None + xVBar = None + xDVBar = None + + wBar, xVBar, xDVBar = solver.computeJacobianVectorProductBwd( + funcsBar=funcsBar, wDeriv=True, xVDeriv=True, xDvDeriv=False, xDvDerivAero=True + ) + if "adflow_states" in d_inputs: + d_inputs["adflow_states"] += wBar + if "adflow_vol_coords" in d_inputs: + d_inputs["adflow_vol_coords"] += xVBar + + for dv_name, dv_bar in xDVBar.items(): + if dv_name in d_inputs: + d_inputs[dv_name] += dv_bar.flatten() + + +class ADflowGroup(Group): + def initialize(self): + self.options.declare("solver", recordable=False) + self.options.declare("struct_coupling", default=False) + self.options.declare("prop_coupling", default=False) + self.options.declare("heat_transfer", default=False) + self.options.declare("use_warper", default=True) + self.options.declare("restart_failed_analysis", default=False) + self.options.declare("err_on_convergence_fail", default=False) + self.options.declare("balance_group", default=None, recordable=False) + + def setup(self): + + self.aero_solver = self.options["solver"] + self.struct_coupling = self.options["struct_coupling"] + self.prop_coupling = self.options["prop_coupling"] + self.restart_failed_analysis = self.options["restart_failed_analysis"] + self.err_on_convergence_fail = self.options["err_on_convergence_fail"] + + self.use_warper = self.options["use_warper"] + + balance_group = self.options["balance_group"] + self.heat_transfer = self.options["heat_transfer"] + + if self.use_warper: + # if we dont have geo_disp, we also need to promote the x_a as x_a0 from the deformer component + self.add_subsystem( + "deformer", + ADflowWarper( + aero_solver=self.aero_solver, + ), + promotes_inputs=["x_aero"], + promotes_outputs=["adflow_vol_coords"], + ) + + self.add_subsystem( + "solver", + ADflowSolver( + aero_solver=self.aero_solver, + restart_failed_analysis=self.restart_failed_analysis, + err_on_convergence_fail=self.err_on_convergence_fail, + ), + promotes_inputs=["adflow_vol_coords"], + promotes_outputs=["adflow_states"], + ) + + if self.struct_coupling: + self.add_subsystem( + "force", + ADflowForces(aero_solver=self.aero_solver), + promotes_inputs=["adflow_vol_coords", "adflow_states"], + promotes_outputs=["f_aero"], + ) + if self.prop_coupling: + self.add_subsystem( + "prop", + ADflowFunctions( + aero_solver=self.aero_solver, + ap_funcs=False, + write_solution=False, + ), + promotes_inputs=["adflow_vol_coords", "adflow_states"], + ) + + if self.heat_transfer: + self.add_subsystem( + "heat_xfer", AdflowHeatTransfer(aero_solver=self.aero_solver), promotes_outputs=["q_convect"] + ) + + if balance_group is not None: + self.add_subsystem("balance", balance_group) + + def mphys_set_ap(self, ap): + # set the ap, add inputs and outputs, promote? + self.solver.set_ap(ap) + # self.funcs.set_ap(ap) + if self.struct_coupling: + self.force.set_ap(ap) + if self.prop_coupling: + self.prop.mphys_set_ap(ap) + + if self.heat_transfer: + self.heat_xfer.set_ap(ap) + + # promote the DVs for this ap + ap_vars, _ = get_dvs_and_cons(ap=ap) + + for (args, _kwargs) in ap_vars: + name = args[0] + self.promotes("solver", inputs=[name]) + # self.promotes('funcs', inputs=[name]) + if self.struct_coupling: + self.promotes("force", inputs=[name]) + if self.prop_coupling: + self.promotes("prop", inputs=[name]) + + def mphys_add_prop_funcs(self, prop_funcs): + # this is the main routine to enable outputs from the propulsion element + + # call the method of the prop element + self.prop.mphys_add_funcs(prop_funcs) + + # promote these variables to the aero group level + self.promotes("prop", outputs=prop_funcs) + + +class ADflowMeshGroup(Group): + def initialize(self): + self.options.declare("aero_solver", recordable=False) + + def setup(self): + aero_solver = self.options["aero_solver"] + + self.add_subsystem("surface_mesh", ADflowMesh(aero_solver=aero_solver), promotes=["*"]) + self.add_subsystem( + "volume_mesh", + ADflowWarper(aero_solver=aero_solver), + promotes_inputs=[("x_aero", "x_aero0")], + promotes_outputs=["adflow_vol_coords"], + ) + + def mphys_add_coordinate_input(self): + # just pass through the call + return self.surface_mesh.mphys_add_coordinate_input() + + def mphys_get_triangulated_surface(self): + # just pass through the call + return self.surface_mesh.mphys_get_triangulated_surface() + + +class ADflowBuilder(Builder): + def __init__( + self, + options, # adflow options + mesh_options=None, # idwarp options + scenario="aerodynamic", # scenario type to configure the groups + restart_failed_analysis=False, # retry after failed analysis + err_on_convergence_fail=False, # raise an analysis error if the solver stalls + balance_group=None, + ): + + # options dictionary for ADflow + self.options = options + + # MACH tools require separate option dictionaries for solver and mesh + # if user did not provide a separate mesh_options dictionary, we just use + # the grid file option from the aero options. + if mesh_options is None: + if "gridFile" in options: + self.mesh_options = { + "gridFile": options["gridFile"], + } + elif "gridfile" in options: + self.mesh_options = { + "gridFile": options["gridfile"], + } + else: + self.mesh_options = mesh_options + + # defaults: + + # flag to determine if the mesh warping component is added + # in the nonlinear solver loop (e.g. for aerostructural) + # or as a preprocessing step like the surface mesh coordinates + # (e.g. for aeropropulsive). This will avoid doing extra work + # for mesh deformation when the volume mesh does not change + # during nonlinear iterations + self.warp_in_solver = False + # flag for aerostructural coupling variables + self.struct_coupling = False + # flag to enable propulsion coupling variables + self.prop_coupling = False + # flag to enable heat transfer coupling variables + # TODO AY-JA: Can you rename heat_transfer to thermal_coupling to be consistent with other flags? + self.heat_transfer = False + + # depending on the scenario we are building for, we adjust a few internal parameters: + if scenario.lower() == "aerodynamic": + # default + pass + + elif scenario.lower() == "aerostructural": + # volume mesh warping needs to be inside the coupling loop for aerostructural + self.warp_in_solver = True + self.struct_coupling = True + + elif scenario.lower() == "aeropropulsive": + self.prop_coupling = True + + elif scenario.lower() == "aerothermal": + self.heat_transfer = True + + # flag to determine if we want to restart a failed solution from free stream + self.restart_failed_analysis = restart_failed_analysis + + # flag for raising an error on convergence stall + self.err_on_convergence_fail = err_on_convergence_fail + + # balance group for propulsion + self.balance_group = balance_group + + # TODO AY-JA: Check if we still need this + # if self.heat_transfer: + # self.promotes('heat_xfer', inputs=[name]) + + # api level method for all builders + def initialize(self, comm): + self.solver = ADFLOW(options=self.options, comm=comm) + mesh = USMesh(options=self.mesh_options, comm=comm) + self.solver.setMesh(mesh) + + def get_solver(self): + # this method is only used by the RLT transfer scheme + return self.solver + + # api level method for all builders + def get_coupling_group_subsystem(self, scenario_name=None): + adflow_group = ADflowGroup( + solver=self.solver, + use_warper=self.warp_in_solver, + struct_coupling=self.struct_coupling, + prop_coupling=self.prop_coupling, + restart_failed_analysis=self.restart_failed_analysis, + err_on_convergence_fail=self.err_on_convergence_fail, + balance_group=self.balance_group, + ) + return adflow_group + + def get_mesh_coordinate_subsystem(self, scenario_name=None): + + # TODO modify this so that we can move the volume mesh warping to the top level + # we need this to do mesh warping only once for all serial points. + # for now, volume warping is duplicated on all scenarios, which is not efficient + + # use_warper = not self.warp_in_solver + # # if we do warper in the mesh element, we will do a group thing + # if use_warper: + # return ADflowMeshGroup(aero_solver=self.solver) + # else: + + # just return the component that outputs the surface mesh. + return ADflowMesh(aero_solver=self.solver) + + def get_pre_coupling_subsystem(self, scenario_name=None): + if self.warp_in_solver: + # if we warp in the solver, then we wont have any pre-coupling systems + return None + else: + # we warp as a pre-processing step + return ADflowWarper(aero_solver=self.solver) + + def get_post_coupling_subsystem(self, scenario_name=None): + return ADflowFunctions(aero_solver=self.solver) + + # TODO the get_nnodes is deprecated. will remove + def get_nnodes(self, groupName=None): + return int(self.solver.getSurfaceCoordinates(groupName=groupName).size / 3) + + def get_number_of_nodes(self, groupName=None): + return int(self.solver.getSurfaceCoordinates(groupName=groupName).size / 3) diff --git a/adflow/om_utils.py b/adflow/mphys/om_utils.py similarity index 92% rename from adflow/om_utils.py rename to adflow/mphys/om_utils.py index 29d4b89b1..6d82b40b8 100644 --- a/adflow/om_utils.py +++ b/adflow/mphys/om_utils.py @@ -21,13 +21,13 @@ def addObj(self, name, *args, **kwargs): def get_dvs_and_cons(ap=None, geo=None, con=None): - vars = [] + dvs = [] cons = [] if ap is not None: ap_vars = DummyOptProb() ap.addVariablesPyOpt(ap_vars) - vars.extend(ap_vars.variables) + dvs.extend(ap_vars.variables) # for long_name, dv in ap.DVs.items(): # dv_data = {'scalar': True, # 'value': dv.value, @@ -41,11 +41,11 @@ def get_dvs_and_cons(ap=None, geo=None, con=None): if geo is not None: geo_vars = DummyOptProb() geo.addVariablesPyOpt(geo_vars) - vars.extend(geo_vars.variables) + dvs.extend(geo_vars.variables) if con is not None: dv_cons = DummyOptProb() con.addConstraintsPyOpt(dv_cons) cons.extend(dv_cons.constraints) - return vars, cons + return dvs, cons diff --git a/adflow/om_adflow.py b/adflow/om_adflow.py deleted file mode 100644 index f9c7d829d..000000000 --- a/adflow/om_adflow.py +++ /dev/null @@ -1,90 +0,0 @@ -import types - -from baseclasses import AeroProblem - -from openmdao.api import Group, PETScKrylov, LinearUserDefined, IndepVarComp - -from .om_states_comp import OM_STATES_COMP -from .om_func_comp import OM_FUNC_COMP -from .om_geocon_comp import OM_GEOCON_COMP -from .om_utils import get_dvs_and_cons - - -class OM_ADFLOW(Group): - """Group integrating the states, funcs, and geometry constraints into a single system""" - - def initialize(self): - self.options.declare("ap", types=AeroProblem) - self.options.declare("setup_cb", types=types.FunctionType, allow_none=True, default=None) - - self.options.declare( - "use_OM_KSP", - default=False, - types=bool, - desc="uses OpenMDAO's PestcKSP linear solver with ADflow's preconditioner to solve the adjoint.", - ) - - self.options.declare( - "owns_indeps", default=False, desc="when True will create IndepVarComp with outputs for all des vars" - ) - self.options.declare("debug", default=False) - - def setup(self): - ap = self.options["ap"] - - self.solver, self.mesh, self.dvgeo, self.dvcon = self.options["setup_cb"](self.comm) - - # necessary to get some memory properly allocated - # self.solver.getResidual(ap) - - des_vars, constraints = get_dvs_and_cons(ap, self.dvgeo, self.dvcon) - geo_vars, _ = get_dvs_and_cons(geo=self.dvgeo) - - self.des_var_names = des_vars - self.geo_var_names = geo_vars - self.constraint_names = constraints - - des_var_names = [dv[0][0] for dv in des_vars] - geo_var_names = [dv[0][0] for dv in geo_vars] - - if self.options["owns_indeps"]: - indeps = self.add_subsystem("indeps", IndepVarComp(), promotes=["*"]) - for (args, kwargs) in des_vars: - name = args[0] - value = kwargs["value"] - if "units" in kwargs: - indeps.add_output(name, value, units=kwargs["units"]) - else: - indeps.add_output(name, value) - - states = OM_STATES_COMP(ap=ap, dvgeo=self.dvgeo, solver=self.solver, use_OM_KSP=self.options["use_OM_KSP"]) - - self.add_subsystem("states", states, promotes_inputs=des_var_names) - - # this lets the OpenMDAO KSPsolver converge the adjoint using the ADflow PC - if self.options["use_OM_KSP"]: - states.linear_solver = PETScKrylov(iprint=2, atol=1e-8, rtol=1e-8, maxiter=300, ksp_type="gmres") - states.linear_solver.precon = LinearUserDefined() - - functionals = OM_FUNC_COMP(ap=ap, dvgeo=self.dvgeo, solver=self.solver) - self.add_subsystem("functionals", functionals, promotes_inputs=des_var_names, max_procs=self.comm.size) - self.connect("states.states", "functionals.states") - - if self.dvcon is not None: - - geocon = OM_GEOCON_COMP(dvgeo=self.dvgeo, dvcon=self.dvcon) - self.add_subsystem("geocon", geocon, promotes_inputs=geo_var_names) - - for cons in self.dvcon.constraints.values(): - for name, con in cons.items(): - if self.comm.rank == 0: - print("adding nonlinear dvcon constarint: ", name) - self.add_constraint( - "geocon.%s" % name, lower=con.lower, upper=con.upper, ref=1 / con.scale, vectorize_derivs=True - ) - for name, con in self.dvcon.linearCon.items(): - if self.comm.rank == 0: - print("adding linear dvcon constraint: ", name) - self.add_constraint( - "geocon.%s" % name, lower=con.lower, upper=con.upper, linear=True, vectorize_derivs=True - ) diff --git a/adflow/om_func_comp.py b/adflow/om_func_comp.py deleted file mode 100644 index 92c5bce2b..000000000 --- a/adflow/om_func_comp.py +++ /dev/null @@ -1,244 +0,0 @@ -import numpy as np - -from openmdao.api import ExplicitComponent -from .om_utils import get_dvs_and_cons - -from baseclasses import AeroProblem -from pygeo import DVGeometry - - -DVGEO_CLASSES = (DVGeometry,) -try: - from pygeo import DVGeometryVSP - - DVGEO_CLASSES = (DVGeometry, DVGeometryVSP) -except ImportError: - pass - -FUNCS_UNITS = { - "mdot": "kg/s", - "mavgptot": "Pa", - "mavgps": "Pa", - "aavgptot": "Pa", - "aavgps": "Pa", - "mavgttot": "degK", - "mavgvx": "m/s", - "mavgvy": "m/s", - "mavgvz": "m/s", - "drag": "N", - "lift": "N", - "dragpressure": "N", - "dragviscous": "N", - "dragmomentum": "N", - "fx": "N", - "fy": "N", - "fz": "N", - "forcexpressure": "N", - "forceypressure": "N", - "forcezpressure": "N", - "forcexviscous": "N", - "forceyviscous": "N", - "forcezviscous": "N", - "forcexmomentum": "N", - "forceymomentum": "N", - "forcezmomentum": "N", - "flowpower": "W", - "area": "m**2", -} - - -class OM_FUNC_COMP(ExplicitComponent): - def initialize(self): - self.options.declare( - "ap", - types=AeroProblem, - ) - self.options.declare("dvgeo", types=DVGEO_CLASSES, allow_none=True, default=None) - self.options.declare("solver") - - # testing flag used for unit-testing to prevent the call to actually solve - # NOT INTENDED FOR USERS!!! FOR TESTING ONLY - self._do_solve = True - - # self.distributed=True - - def setup(self): - solver = self.options["solver"] - ap = self.options["ap"] - geo = self.options["dvgeo"] - - self.ap_vars, _ = get_dvs_and_cons(ap=ap) - self.geo_vars, _ = get_dvs_and_cons(geo=geo) - - for (args, kwargs) in self.geo_vars: - name = args[0] - size = args[1] - self.add_input(name, shape=size) - - for (args, kwargs) in self.ap_vars: - name = args[0] - size = args[1] - - self.add_input(name, shape=size, units=kwargs["units"]) - - local_state_size = self.options["solver"].getStateSize() - local_state_sizes = self.comm.allgather(local_state_size) - iproc = self.comm.rank - - ind1 = np.sum(local_state_sizes[:iproc]) - ind2 = np.sum(local_state_sizes[: iproc + 1]) - - self.add_input("states", src_indices=np.arange(ind1, ind2, dtype=int)) - - for f_name, f_meta in solver.adflowCostFunctions.items(): - f_type = f_meta[1] - units = None - if f_type in FUNCS_UNITS: - units = FUNCS_UNITS[f_type] - - if self.comm.rank == 0: - print("adding adflow func as output: {}".format(f_name)) - self.add_output(f_name, shape=1, units=units) - - # self.declare_partials(of=f_name, wrt='*') - - def _set_ap(self, inputs): - tmp = {} - for (args, kwargs) in self.ap_vars: - name = args[0] - tmp[name] = inputs[name][0] - - self.options["ap"].setDesignVars(tmp) - # self.options['solver'].setAeroProblem(self.options['ap']) - - def _set_geo(self, inputs, update_jacobian=True): - dvgeo = self.options["dvgeo"] - if dvgeo is None: - return - - tmp = {} - for (args, kwargs) in self.geo_vars: - name = args[0] - tmp[name] = inputs[name] - - # if self.comm.rank == 0: - # import pprint - # pprint.pprint(tmp) - try: - self.options["dvgeo"].setDesignVars(tmp, update_jacobian) - except TypeError: # this is needed because dvGeo and dvGeoVSP have different APIs - self.options["dvgeo"].setDesignVars(tmp) - - def _set_states(self, inputs): - self.options["solver"].setStates(inputs["states"]) - - def _get_func_name(self, name): - return "%s_%s" % (self.options["ap"].name, name.lower()) - - def compute(self, inputs, outputs): - solver = self.options["solver"] - ap = self.options["ap"] - # print('funcs compute') - # actually setting things here triggers some kind of reset, so we only do it if you're actually solving - if self._do_solve: - self._set_ap(inputs) - self._set_geo(inputs, update_jacobian=False) - self._set_states(inputs) - - funcs = {} - - eval_funcs = [f_name for f_name, f_meta in solver.adflowCostFunctions.items()] - solver.evalFunctions(ap, funcs, eval_funcs) - # solver.evalFunctions(ap, funcs) - - # for name in ap.evalFuncs: - for name in solver.adflowCostFunctions.keys(): - f_name = self._get_func_name(name) - if f_name in funcs: - outputs[name.lower()] = funcs[f_name] - - def _compute_partials(self, inputs, J): - - # self._set_ap(inputs) - # self._set_geo(inputs) - # self._set_states(inputs) - - # print('om_funcs linearize') - - pass - - def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): - solver = self.options["solver"] - ap = self.options["ap"] - - # self.options['solver'].setAeroProblem(ap) - # print('func matvec') - - # if self._do_solve: - # self._set_ap(inputs) - # self._set_geo(inputs) - # self._set_states(inputs) - - if mode == "fwd": - xDvDot = {} - for key in ap.DVs: - if key in d_inputs: - mach_name = key.split("_")[0] - xDvDot[mach_name] = d_inputs[key] - for (args, kwargs) in self.geo_vars: - name = args[0] - xDvDot[name] = d_inputs[name] - - if "states" in d_inputs: - wDot = d_inputs["states"] - else: - wDot = None - - funcsdot = solver.computeJacobianVectorProductFwd(xDvDot=xDvDot, wDot=wDot, funcDeriv=True) - - for name, meta in ap.evalFuncs: - func_name = name.lower() - if name in d_outputs: - d_outputs[name] += funcsdot[func_name] - - elif mode == "rev": - funcsBar = {} - - for name, meta in solver.adflowCostFunctions.items(): - func_name = name.lower() - - # we have to check for 0 here, so we don't include any unnecessary variables in funcsBar - # becasue it causes ADflow to do extra work internally even if you give it extra variables, even if the seed is 0 - if func_name in d_outputs and d_outputs[func_name] != 0.0: - funcsBar[func_name] = d_outputs[func_name][0] - - # because of the 0 checking, the funcsBar is now only correct on the root proc, - # so we need to broadcast it to everyone. its not actually imporant that the seeds are the same everywhere, - # but the keys in the dictionary need to be the same. - funcsBar = self.comm.bcast(funcsBar, root=0) - - # print(funcsBar, flush=True) - - d_input_vars = list(d_inputs.keys()) - n_input_vars = len(d_input_vars) - - wBar = None - xDVBar = None - - if "states" in d_inputs and n_input_vars == 1: - wBar = solver.computeJacobianVectorProductBwd(funcsBar=funcsBar, wDeriv=True) - elif ("states" not in d_input_vars) and n_input_vars: - xDVBar = solver.computeJacobianVectorProductBwd(funcsBar=funcsBar, xDvDeriv=True) - - elif ("states" in d_input_vars) and n_input_vars: - wBar, xDVBar = solver.computeJacobianVectorProductBwd(funcsBar=funcsBar, wDeriv=True, xDvDeriv=True) - else: # nothing to do, so why is this being called? - return - - if wBar is not None: - d_inputs["states"] += wBar - - if xDVBar is not None: - for dv_name, dv_bar in xDVBar.items(): - if dv_name in d_inputs: - d_inputs[dv_name] += dv_bar.flatten() diff --git a/adflow/om_geocon_comp.py b/adflow/om_geocon_comp.py deleted file mode 100644 index 7e00c5c48..000000000 --- a/adflow/om_geocon_comp.py +++ /dev/null @@ -1,98 +0,0 @@ -from six import iteritems, itervalues - -from openmdao.api import ExplicitComponent -from .om_utils import get_dvs_and_cons - -from pygeo import DVGeometry, DVConstraints - -DVGEO_CLASSES = (DVGeometry,) -try: - from pygeo import DVGeometryVSP - - DVGEO_CLASSES = (DVGeometry, DVGeometryVSP) -except ImportError: - pass - - -class OM_GEOCON_COMP(ExplicitComponent): - """OpenMDAO Component that wraps the geometry constraints calculation""" - - def initialize(self): - self.options.declare("dvgeo", types=DVGEO_CLASSES) - self.options.declare("dvcon", types=DVConstraints) - - # testing flag used for unit-testing to prevent the call to actually solve - # NOT INTENDED FOR USERS!!! FOR TESTING ONLY - self._do_solve = True - - def setup(self): - self.dvs, _ = get_dvs_and_cons(self.options["dvgeo"]) - for (args, kwargs) in self.dvs: - name = args[0] - size = args[1] - self.add_input(name, shape=size) - - funcsSens = {} - self.options["dvcon"].evalFunctionsSens(funcsSens, includeLinear=True) - - for cons in itervalues(self.options["dvcon"].constraints): - for name, con in iteritems(cons): - if self.comm.rank == 0: - print("nonlinear con", name) - self.add_output(name, shape=con.nCon) - jac = funcsSens[name] - for wrt_var, subjac in iteritems(jac): - self.declare_partials(of=name, wrt=wrt_var) - - for name, con in iteritems(self.options["dvcon"].linearCon): - if self.comm.rank == 0: - print("linear_con", name) - self.add_output(name, shape=con.ncon) - jac = funcsSens[name] - for wrt_var, subjac in iteritems(jac): - self.declare_partials(of=name, wrt=wrt_var, val=subjac) - - def _set_geo(self, inputs, updateJacobian=True): - tmp = {} - for (args, kwargs) in self.dvs: - name = args[0] - tmp[name] = inputs[name] - - try: - self.options["dvgeo"].setDesignVars(tmp, updateJacobian) - except TypeError: # this is needed because dvGeo and dvGeoVSP have different APIs - self.options["dvgeo"].setDesignVars(tmp) - - def compute(self, inputs, outputs): - - if self._do_solve: - # self._set_geo(inputs, updateJacobian=False) - pass - - funcs = {} - self.options["dvcon"].evalFunctions(funcs, includeLinear=True) - - for cons in itervalues(self.options["dvcon"].constraints): - for name, con in iteritems(cons): - # print('nonlinear foobar', name, funcs[name]) - outputs[name] = funcs[name] - - for name, con in iteritems(self.options["dvcon"].linearCon): - # print('linear foobar', name, funcs[name]) - outputs[name] = funcs[name] - - def compute_partials(self, inputs, J): - # print('om_geocon linearize') - if self._do_solve: - # self._set_geo(inputs) - pass - - funcsSens = {} - - self.options["dvcon"].evalFunctionsSens(funcsSens) - - for of_var, jac in funcsSens.items(): - # print('of: ', of_var) - for wrt_var, subjac in jac.items(): - # print(' wrt: ', wrt_var) - J[of_var, wrt_var] = subjac diff --git a/adflow/om_states_comp.py b/adflow/om_states_comp.py deleted file mode 100644 index 7cdddbfd5..000000000 --- a/adflow/om_states_comp.py +++ /dev/null @@ -1,241 +0,0 @@ -from pygeo import DVGeometry -from baseclasses import AeroProblem - -from openmdao.api import ImplicitComponent -from openmdao.core.analysis_error import AnalysisError -from .om_utils import get_dvs_and_cons - -DVGEO_CLASSES = (DVGeometry,) -try: - from pygeo import DVGeometryVSP - - DVGEO_CLASSES = (DVGeometry, DVGeometryVSP) -except ImportError: - pass - - -class OM_STATES_COMP(ImplicitComponent): - """OpenMDAO component that wraps the flow solve""" - - def initialize(self): - self.options.declare("ap", types=AeroProblem) - self.options.declare("dvgeo", types=DVGEO_CLASSES, allow_none=True, default=None) - self.options.declare("solver") - self.options.declare( - "use_OM_KSP", - default=False, - types=bool, - desc="uses OpenMDAO's PestcKSP linear solver with ADflow's preconditioner to solve the adjoint.", - ) - - # self.options.declare('max_procs', default=64, types=int) - - self.options["distributed"] = True - - # testing flag used for unit-testing to prevent the call to actually solve - # NOT INTENDED FOR USERS!!! FOR TESTING ONLY - self._do_solve = True - - def setup(self): - solver = self.options["solver"] - ap = self.options["ap"] - geo = self.options["dvgeo"] - - # flag to keep track of clean restarts. We need this because sometimes - # the current solution fails even if it already was a clean restart because - # the previous one also failed. In these cases, we dont want to retry - # with a clean restart because there is no point. - self.cleanRestart = True - - self.ap_vars, _ = get_dvs_and_cons(ap=ap) - self.geo_vars, _ = get_dvs_and_cons(geo=geo) - - if self.comm.rank == 0: - print("adding ap var inputs") - for (args, kwargs) in self.ap_vars: - name = args[0] - size = args[1] - self.add_input(name, shape=size, units=kwargs["units"]) - if self.comm.rank == 0: - print(name) - - if self.comm.rank == 0: - print("adding geo var inputs") - - for (args, kwargs) in self.geo_vars: - name = args[0] - size = args[1] - self.add_input(name, shape=size) - if self.comm.rank == 0: - print(name) - - local_state_size = solver.getStateSize() - self.add_output("states", shape=local_state_size) - - # apparently needed to initialize the state arrays - # but also somehow borks up the restarts - # solver.getResidual(ap) - - # self.declare_partials(of='states', wrt='*') - - def _set_ap(self, inputs): - tmp = {} - for (args, kwargs) in self.ap_vars: - name = args[0] - tmp[name] = inputs[name] - - self.options["ap"].setDesignVars(tmp) - - def _set_geo(self, inputs, update_jacobian=True): - dvgeo = self.options["dvgeo"] - if dvgeo is None: - return - - tmp = {} - for (args, kwargs) in self.geo_vars: - name = args[0] - tmp[name] = inputs[name] - try: - self.options["dvgeo"].setDesignVars(tmp, update_jacobian) - except TypeError: # this is needed because dvGeo and dvGeoVSP have different APIs - self.options["dvgeo"].setDesignVars(tmp) - - def _set_states(self, outputs): - self.options["solver"].setStates(outputs["states"]) - - def apply_nonlinear(self, inputs, outputs, residuals): - - self._set_states(outputs) - self._set_ap(inputs) - self._set_geo(inputs, update_jacobian=False) - - ap = self.options["ap"] - residuals["states"] = self.options["solver"].getResidual(ap) - - def solve_nonlinear(self, inputs, outputs): - - solver = self.options["solver"] - ap = self.options["ap"] - - if self._do_solve: - - self._set_ap(inputs) - self._set_geo(inputs, update_jacobian=False) - ap.solveFailed = False # might need to clear this out? - ap.fatalFail = False - - solver(ap) - - if ap.fatalFail: - if self.comm.rank == 0: - print("###############################################################") - print("# Solve Fatal Fail. Analysis Error") - print("###############################################################") - - raise AnalysisError("ADFLOW Solver Fatal Fail") - - if ap.solveFailed: # the mesh was fine, but it didn't converge - # if the previous iteration was already a clean restart, dont try again - if self.cleanRestart: - print("###############################################################") - print("# This was a clean restart. Will not try another one.") - print("###############################################################") - solver.resetFlow(ap) - self.cleanRestart = True - raise AnalysisError("ADFLOW Solver Fatal Fail") - - # the previous iteration restarted from another solution, so we can try again - # with a re-set flowfield for the initial guess. - else: - if self.comm.rank == 0: - print("###############################################################") - print("# Solve Failed, attempting a clean restart!") - print("###############################################################") - - ap.solveFailed = False - ap.fatalFail = False - solver.resetFlow(ap) - solver(ap) - - if ap.solveFailed or ap.fatalFail: # we tried, but there was no saving it - print("###############################################################") - print("# Clean Restart failed. There is no saving this one!") - print("###############################################################") - - # re-set the flow for the next iteration: - solver.resetFlow(ap) - # set the reset flow flag - self.cleanRestart = True - raise AnalysisError("ADFLOW Solver Fatal Fail") - - # see comment for the same flag below - else: - self.cleanRestart = False - - # solve did not fail, therefore we will re-use this converged flowfield for the next iteration. - # change the flag so that if the next iteration fails with current initial guess, it can retry - # with a clean restart - else: - self.cleanRestart = False - - outputs["states"] = solver.getStates() - - def linearize(self, inputs, outputs, residuals): - - self.options["solver"]._setupAdjoint() - - self._set_ap(inputs) - self._set_geo(inputs, update_jacobian=False) - self._set_states(outputs) - - # print('om_states linearize') - - def apply_linear(self, inputs, outputs, d_inputs, d_outputs, d_residuals, mode): - - solver = self.options["solver"] - - # self._set_ap(inputs) - # self._set_geo(inputs) - # self._set_states(outputs) - - if mode == "fwd": - if "states" in d_residuals: - xDvDot = {} - for var_name in d_inputs: - xDvDot[var_name] = d_inputs[var_name] - if "states" in d_outputs: - wDot = d_outputs["states"] - else: - wDot = None - - dwdot = solver.computeJacobianVectorProductFwd(xDvDot=xDvDot, wDot=wDot, residualDeriv=True) - d_residuals["states"] += dwdot - - elif mode == "rev": - if "states" in d_residuals: - resBar = d_residuals["states"] - - wBar, xDVBar = solver.computeJacobianVectorProductBwd(resBar=resBar, wDeriv=True, xDvDeriv=True) - - if "states" in d_outputs: - d_outputs["states"] += wBar - - for dv_name, dv_bar in xDVBar.items(): - if dv_name in d_inputs: - d_inputs[dv_name] += dv_bar.flatten() - - def solve_linear(self, d_outputs, d_residuals, mode): - solver = self.options["solver"] - if self.options["use_OM_KSP"]: - if mode == "fwd": - d_outputs["states"] = solver.globalNKPreCon(d_residuals["states"], d_outputs["states"]) - elif mode == "rev": - d_residuals["states"] = solver.globalAdjointPreCon(d_outputs["states"], d_residuals["states"]) - else: - if mode == "fwd": - d_outputs["states"] = solver.solveDirectForRHS(d_residuals["states"]) - elif mode == "rev": - # d_residuals['states'] = solver.solveAdjointForRHS(d_outputs['states']) - solver.adflow.adjointapi.solveadjoint(d_outputs["states"], d_residuals["states"], True) - - return True, 0, 0 diff --git a/setup.py b/setup.py index 915bcc309..7382d2beb 100644 --- a/setup.py +++ b/setup.py @@ -1,4 +1,4 @@ -from setuptools import setup +from setuptools import setup, find_packages import re import os @@ -22,9 +22,7 @@ author_email="", url="https://github.com/mdolab/adflow", license="LGPL version 2.1", - packages=[ - "adflow", - ], + packages=find_packages(include=["adflow*"]), package_data={"adflow": ["*.so"]}, install_requires=[ "numpy>=1.16", @@ -32,6 +30,9 @@ "mpi4py>=3.0", "petsc4py>=3.11", ], - extras_require={"testing": ["parameterized", "testflo"]}, + extras_require={ + "testing": ["parameterized", "testflo"], + "mphys": ["openmdao", "mphys", "idwarp"], + }, classifiers=["Operating System :: Linux", "Programming Language :: Python, Fortran"], )