diff --git a/dafoam/mphys/mphys_dafoam.py b/dafoam/mphys/mphys_dafoam.py index 31f53da0..ab046831 100644 --- a/dafoam/mphys/mphys_dafoam.py +++ b/dafoam/mphys/mphys_dafoam.py @@ -209,16 +209,24 @@ def setup(self): if self.prop_coupling is not None: if self.prop_coupling == "Prop": + activeProp = None + nActiveProps = 0 + wingPropNameList = list(self.DASolver.getOption("wingProp").keys()) - if len(wingPropNameList) != 1: - raise RuntimeError("if prop_coupling = Prop, we should set only one propeller for wingProp key.") + for propName in wingPropNameList: + if self.DASolver.getOption("wingProp")[propName]["active"]: + activeProp = propName + nActiveProps += 1 - propName = wingPropNameList[0] + if nActiveProps != 1: + raise RuntimeError( + "if prop_coupling = Prop, we should set only one active propeller for wingProp key." + ) self.add_subsystem( - propName, - DAFoamPropForce(solver=self.DASolver, propName=propName), + activeProp, + DAFoamPropForce(solver=self.DASolver, propName=activeProp), promotes_inputs=["*"], promotes_outputs=["*"], ) @@ -947,6 +955,12 @@ def solve_linear(self, d_outputs, d_residuals, mode): DASolver.ksp = PETSc.KSP().create(self.comm) DASolver.solverAD.createMLRKSPMatrixFree(DASolver.dRdWTPC, DASolver.ksp) + # if useNonZeroInitGuess is False, we will manually reset self.psi to zero + # this is important because we need the correct psi to update the KSP tolerance + # in the next line + if not self.DASolver.getOption("adjEqnOption")["useNonZeroInitGuess"]: + self.psi.set(0) + if self.DASolver.getOption("adjEqnOption")["dynAdjustTol"]: # if we want to dynamically adjust the tolerance, call this function. This is mostly used # in the block Gauss-Seidel method in two discipline coupling @@ -1707,7 +1721,7 @@ class DAFoamPropForce(ExplicitComponent): """ DAFoam component that computes the propeller force and radius profile based on the CFD surface states """ - + def initialize(self): self.options.declare("solver", recordable=False) self.options.declare("propName", recordable=False) @@ -1726,9 +1740,9 @@ def setup(self): self.add_input("prop_center", distributed=False, shape=3, tags=["mphys_coupling"]) # outputs - self.add_output("axial_force_profile", distributed=False, shape=self.nForceSections, tags=["mphys_coupling"]) - self.add_output("tangential_force_profile", distributed=False, shape=self.nForceSections, tags=["mphys_coupling"]) - self.add_output("radial_location", distributed=False, shape=self.nForceSections, tags=["mphys_coupling"]) + self.add_output("axial_force", distributed=False, shape=self.nForceSections, tags=["mphys_coupling"]) + self.add_output("tangential_force", distributed=False, shape=self.nForceSections, tags=["mphys_coupling"]) + self.add_output("radial_location", distributed=False, shape=self.nForceSections, tags=["mphys_coupling"]) def compute(self, inputs, outputs): @@ -1742,10 +1756,6 @@ def compute(self, inputs, outputs): radialLocationVec = PETSc.Vec().createSeq(self.nForceSections, bsize=1, comm=PETSc.COMM_SELF) radialLocationVec.zeroEntries() - outputs["axial_force_profile"] = DASolver.vec2ArraySeq(axialForceProfileVec) - outputs["tangential_force_profile"] = DASolver.vec2ArraySeq(tangentialForceProfileVec) - outputs["radial_location"] = DASolver.vec2ArraySeq(radialLocationVec) - prop_center = inputs["prop_center"] prop_center_vec = DASolver.array2VecSeq(prop_center) @@ -1757,6 +1767,10 @@ def compute(self, inputs, outputs): radialLocationVec, ) + outputs["axial_force"] = DASolver.vec2ArraySeq(axialForceProfileVec) + outputs["tangential_force"] = DASolver.vec2ArraySeq(tangentialForceProfileVec) + outputs["radial_location"] = DASolver.vec2ArraySeq(radialLocationVec) + def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): DASolver = self.DASolver @@ -1768,48 +1782,178 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode): category=om.OpenMDAOWarning, ) return - + dafoam_states = inputs["%s_states" % self.discipline] dafoam_xv = inputs["%s_vol_coords" % self.discipline] + dafoam_center = inputs["prop_center"] stateVec = DASolver.array2Vec(dafoam_states) xvVec = DASolver.array2Vec(dafoam_xv) + centerVec = DASolver.array2Vec(dafoam_center) - if "axial_force_profile" in d_outputs: - aFBar = d_outputs["force_profile"] + if "axial_force" in d_outputs: + aFBar = d_outputs["axial_force"] aFBarVec = DASolver.array2VecSeq(aFBar) if "%s_states" % self.discipline in d_inputs: prodVec = self.DASolver.wVec.duplicate() prodVec.zeroEntries() - DASolver.solverAD.calcdForceProfiledXvWAD(self.propName.encode(), "aForce".encode(), "state".encode, xvVec, stateVec, aFBarVec, prodVec) - aBar = DASolver.vec2ArraySeq(prodVec) - aBar = self.comm.allreduce(aBar, op=MPI.SUM) - d_inputs["%s_states" % self.discipline] += aBar + DASolver.solverAD.calcdForceProfiledXvWAD( + self.propName.encode(), + "aForce".encode(), + "state".encode(), + xvVec, + stateVec, + centerVec, + aFBarVec, + prodVec, + ) + sBar = DASolver.vec2ArraySeq(prodVec) + # sBar = self.comm.allreduce(sBar, op=MPI.SUM) + d_inputs["%s_states" % self.discipline] += sBar if "%s_vol_coords" % self.discipline in d_inputs: prodVec = self.DASolver.xvVec.duplicate() prodVec.zeroEntries() - aBar = DASolver.vec2Array(prodVec) - d_inputs["%s_vol_coords" % self.discipline] += aBar + DASolver.solverAD.calcdForceProfiledXvWAD( + self.propName.encode(), + "aForce".encode(), + "mesh".encode(), + xvVec, + stateVec, + centerVec, + aFBarVec, + prodVec, + ) + vBar = DASolver.vec2ArraySeq(prodVec) + # vBar = self.comm.allreduce(vBar, op=MPI.SUM) + d_inputs["%s_vol_coords" % self.discipline] += vBar + + if "prop_center" in d_inputs: + prodVec = PETSc.Vec().createSeq(3, bsize=1, comm=PETSc.COMM_SELF) + prodVec.zeroEntries() + DASolver.solverAD.calcdForceProfiledXvWAD( + self.propName.encode(), + "aForce".encode(), + "center".encode(), + xvVec, + stateVec, + centerVec, + aFBarVec, + prodVec, + ) + cBar = DASolver.vec2ArraySeq(prodVec) + # cBar = self.comm.allreduce(cBar, op=MPI.SUM) + d_inputs["prop_center"] += cBar - if "tangential_force_profile" in d_outputs: - tFBar = d_outputs["force_profile"] + if "tangential_force" in d_outputs: + tFBar = d_outputs["tangential_force"] tFBarVec = DASolver.array2VecSeq(tFBar) if "%s_states" % self.discipline in d_inputs: prodVec = self.DASolver.wVec.duplicate() prodVec.zeroEntries() - DASolver.solverAD.calcdForceProfiledXvWAD(self.propName.encode(), "tForce".encode(), "state".encode, xvVec, stateVec, tFBarVec, prodVec) - tBar = DASolver.vec2ArraySeq(prodVec) - tBar = self.comm.allreduce(aBar, op=MPI.SUM) - d_inputs["%s_states" % self.discipline] += tBar + DASolver.solverAD.calcdForceProfiledXvWAD( + self.propName.encode(), + "tForce".encode(), + "state".encode(), + xvVec, + stateVec, + centerVec, + tFBarVec, + prodVec, + ) + sBar = DASolver.vec2ArraySeq(prodVec) + # sBar = self.comm.allreduce(sBar, op=MPI.SUM) + d_inputs["%s_states" % self.discipline] += sBar + + if "%s_vol_coords" % self.discipline in d_inputs: + prodVec = self.DASolver.xvVec.duplicate() + prodVec.zeroEntries() + DASolver.solverAD.calcdForceProfiledXvWAD( + self.propName.encode(), + "tForce".encode(), + "mesh".encode(), + xvVec, + stateVec, + centerVec, + tFBarVec, + prodVec, + ) + vBar = DASolver.vec2ArraySeq(prodVec) + # vBar = self.comm.allreduce(vBar, op=MPI.SUM) + d_inputs["%s_vol_coords" % self.discipline] += vBar + + if "prop_center" in d_inputs: + prodVec = PETSc.Vec().createSeq(3, bsize=1, comm=PETSc.COMM_SELF) + prodVec.zeroEntries() + DASolver.solverAD.calcdForceProfiledXvWAD( + self.propName.encode(), + "tForce".encode(), + "center".encode(), + xvVec, + stateVec, + centerVec, + tFBarVec, + prodVec, + ) + cBar = DASolver.vec2ArraySeq(prodVec) + # cBar = self.comm.allreduce(cBar, op=MPI.SUM) + d_inputs["prop_center"] += cBar + + if "radial_location" in d_inputs: + rDBar = d_outputs["radial_location"] + rDBarVec = DASolver.array2VecSeq(rDBar) + if "%s_states" % self.discipline in d_inputs: + prodVec = self.DASolver.wVec.duplicate() + prodVec.zeroEntries() + DASolver.solverAD.calcdForceProfiledXvWAD( + self.propName.encode(), + "rDist".encode(), + "state".encode(), + xvVec, + stateVec, + centerVec, + rDBarVec, + prodVec, + ) + sBar = DASolver.vec2ArraySeq(prodVec) + # sBar = self.comm.allreduce(sBar, op=MPI.SUM) + d_inputs["%s_states" % self.discipline] += sBar if "%s_vol_coords" % self.discipline in d_inputs: prodVec = self.DASolver.xvVec.duplicate() prodVec.zeroEntries() - tBar = DASolver.vec2Array(prodVec) - d_inputs["%s_vol_coords" % self.discipline] += tBar + DASolver.solverAD.calcdForceProfiledXvWAD( + self.propName.encode(), + "rDist".encode(), + "mesh".encode(), + xvVec, + stateVec, + centerVec, + rDBarVec, + prodVec, + ) + vBar = DASolver.vec2ArraySeq(prodVec) + # vBar = self.comm.allreduce(vBar, op=MPI.SUM) + d_inputs["%s_vol_coords" % self.discipline] += vBar + + if "prop_center" in d_inputs: + prodVec = PETSc.Vec().createSeq(3, bsize=1, comm=PETSc.COMM_SELF) + prodVec.zeroEntries() + DASolver.solverAD.calcdForceProfiledXvWAD( + self.propName.encode(), + "rDist".encode(), + "center".encode()(), + xvVec, + stateVec, + centerVec, + rDBarVec, + prodVec, + ) + cBar = DASolver.vec2ArraySeq(prodVec) + # cBar = self.comm.allreduce(cBar, op=MPI.SUM) + d_inputs["prop_center"] += cBar class DAFoamFvSource(ExplicitComponent): diff --git a/src/adjoint/DASolver/DASolver.C b/src/adjoint/DASolver/DASolver.C index 803b5362..ca0cc8d6 100644 --- a/src/adjoint/DASolver/DASolver.C +++ b/src/adjoint/DASolver/DASolver.C @@ -1800,11 +1800,12 @@ void DASolver::calcForceProfileInternal( void DASolver::calcdForceProfiledXvWAD( const word propName, - const word inputMode, const word outputMode, + const word inputMode, const Vec xvVec, const Vec wVec, Vec center, + Vec psi, Vec dForcedXvW) { /* @@ -1817,7 +1818,6 @@ void DASolver::calcdForceProfiledXvWAD( Info << "Calculating [dForceProfile/dInputs]^T*Psi using reverse-mode AD. PropName: " << propName << endl; - // Buradan basliyorum VecZeroEntries(dForcedXvW); this->updateOFField(wVec); @@ -1832,6 +1832,7 @@ void DASolver::calcdForceProfiledXvWAD( vector centerVector = vector::zero; pointField meshPoints = meshPtr_->points(); + const PetscScalar* vecArrayPsi; PetscScalar* vecArrayCenter; VecGetArray(center, &vecArrayCenter); centerVector[0] = vecArrayCenter[0]; @@ -1844,18 +1845,7 @@ void DASolver::calcdForceProfiledXvWAD( this->globalADTape_.setActive(); // Step 2 - /* - forAll(meshPoints, i) - { - for (label j = 0; j < 3; j++) - { - this->globalADTape_.registerInput(meshPoints[i][j]); - } - } - meshPtr_->movePoints(meshPoints); - meshPtr_->moving(false); - */ - if (outputMode == "mesh") + if (inputMode == "mesh") { forAll(meshPoints, i) { @@ -1867,13 +1857,20 @@ void DASolver::calcdForceProfiledXvWAD( meshPtr_->movePoints(meshPoints); meshPtr_->moving(false); } - else if (outputMode == "state") + else if (inputMode == "state") { this->registerStateVariableInput4AD(); } + else if (inputMode == "center") + { + for (label i = 0; i < 3; i++) + { + this->globalADTape_.registerInput(centerVector[i]); + } + } else { - FatalErrorIn("calcdFvSourcedInputsTPsiAD") << "outputMode not valid" + FatalErrorIn("calcdFvSourcedInputsTPsiAD") << "inputMode not valid" << abort(FatalError); } daResidualPtr_->correctBoundaryConditions(); @@ -1885,23 +1882,30 @@ void DASolver::calcdForceProfiledXvWAD( this->calcForceProfileInternal(propName, centerVector, aForce, tForce, rDist); // Step 4 - if (inputMode == "aForce") + if (outputMode == "aForce") { for(label i = 0; iglobalADTape_.registerOutput(aForce[i]); } } - else if (inputMode == "tForce") + else if (outputMode == "tForce") { for(label i = 0; iglobalADTape_.registerOutput(tForce[i]); } } + else if (outputMode == "rDist") + { + for(label i = 0; iglobalADTape_.registerOutput(rDist[i]); + } + } else { - FatalErrorIn("calcdFvSourcedInputsTPsiAD") << "inputMode not valid" + FatalErrorIn("calcdFvSourcedInputsTPsiAD") << "outputMode not valid" << abort(FatalError); } @@ -1909,37 +1913,36 @@ void DASolver::calcdForceProfiledXvWAD( this->globalADTape_.setPassive(); // Step 6 - if (inputMode == "aForce") + VecGetArrayRead(psi, &vecArrayPsi); + if (outputMode == "aForce") { forAll(aForce, i) { - aForce[i].setGradient(1.0); + aForce[i].setGradient(vecArrayPsi[i]); + } + } + else if (outputMode == "tForce") + { + forAll(tForce, i) + { + tForce[i].setGradient(vecArrayPsi[i]); } } - else if (inputMode == "tForce") + else if (outputMode == "rDist") { forAll(tForce, i) { - tForce[i].setGradient(1.0); + rDist[i].setGradient(vecArrayPsi[i]); } + } + VecRestoreArrayRead(psi, &vecArrayPsi); // Step 7 this->globalADTape_.evaluate(); // Step 8 - /* - forAll(meshPoints, i) - { - for (label j = 0; j < 3; j++) - { - label rowI = daIndexPtr_->getGlobalXvIndex(i, j); - PetscScalar val = meshPoints[i][j].getGradient(); - VecSetValue(dForcedXvW, rowI, val, INSERT_VALUES); - } - } - */ - if (outputMode == "mesh") + if (inputMode == "mesh") { forAll(meshPoints, i) { @@ -1950,17 +1953,30 @@ void DASolver::calcdForceProfiledXvWAD( VecSetValue(dForcedXvW, rowI, val, INSERT_VALUES); } } + VecAssemblyBegin(dForcedXvW); + VecAssemblyEnd(dForcedXvW); } - else if (outputMode == "state") + else if (inputMode == "state") { this->assignStateGradient2Vec(dForcedXvW); + VecAssemblyBegin(dForcedXvW); + VecAssemblyEnd(dForcedXvW); + } + else if (inputMode == "center") + { + PetscScalar* vecArrayCenter; + VecGetArray(dForcedXvW, &vecArrayCenter); + forAll(centerVector, i) + { + vecArrayCenter[i] = centerVector[i].getGradient(); + } + VecRestoreArray(dForcedXvW, &vecArrayCenter); } - VecAssemblyBegin(dForcedXvW); - VecAssemblyEnd(dForcedXvW); // Step 9 this->globalADTape_.clearAdjoints(); this->globalADTape_.reset(); + #endif } diff --git a/src/adjoint/DASolver/DASolver.H b/src/adjoint/DASolver/DASolver.H index b249f020..60b23698 100644 --- a/src/adjoint/DASolver/DASolver.H +++ b/src/adjoint/DASolver/DASolver.H @@ -818,6 +818,7 @@ public: const Vec xvVec, const Vec wVec, Vec center, + Vec psi, Vec dForcedXvW); void calcdForcedStateTPsiAD( diff --git a/src/pyDASolvers/DASolvers.H b/src/pyDASolvers/DASolvers.H index 2d9e1d09..4efbe5c4 100755 --- a/src/pyDASolvers/DASolvers.H +++ b/src/pyDASolvers/DASolvers.H @@ -786,9 +786,10 @@ public: const Vec xvVec, const Vec wVec, Vec center, + Vec psi, Vec dForcedXvW) { - DASolverPtr_->calcdForceProfiledXvWAD(propName, inputMode, outputMode, xvVec, wVec, center, dForcedXvW); + DASolverPtr_->calcdForceProfiledXvWAD(propName, inputMode, outputMode, xvVec, wVec, center, psi, dForcedXvW); } void calcdForcedStateTPsiAD( diff --git a/src/pyDASolvers/pyDASolvers.pyx b/src/pyDASolvers/pyDASolvers.pyx index e457c75f..fcb4d77c 100755 --- a/src/pyDASolvers/pyDASolvers.pyx +++ b/src/pyDASolvers/pyDASolvers.pyx @@ -126,7 +126,7 @@ cdef extern from "DASolvers.H" namespace "Foam": void calcFvSource(char *, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec) void calcdFvSourcedInputsTPsiAD(char *, char *, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec) void calcForceProfile(char *, PetscVec, PetscVec, PetscVec, PetscVec) - void calcdForceProfiledXvWAD(char *, char *, char *, PetscVec, PetscVec, PetscVec, PetscVec) + void calcdForceProfiledXvWAD(char *, char *, char *, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec) void calcdForcedStateTPsiAD(char *, PetscVec, PetscVec, PetscVec, PetscVec) int runFPAdj(PetscVec, PetscVec, PetscVec, PetscVec) void initTensorFlowFuncs(pyComputeInterface, void *, pyJacVecProdInterface, void *) @@ -512,8 +512,8 @@ cdef class pyDASolvers: def calcForceProfile(self, propName, Vec center, Vec aForce, Vec tForce, Vec rDist): self._thisptr.calcForceProfile(propName, center.vec, aForce.vec, tForce.vec, rDist.vec) - def calcdForceProfiledXvWAD(self, propName, inputMode, outputMode, Vec xvVec, Vec wVec, Vec center, Vec dForcedXvW): - self._thisptr.calcdForceProfiledXvWAD(propName, inputMode, outputMode, xvVec.vec, wVec.vec, center.vec, dForcedXvW.vec) + def calcdForceProfiledXvWAD(self, propName, inputMode, outputMode, Vec xvVec, Vec wVec, Vec center, Vec psi, Vec dForcedXvW): + self._thisptr.calcdForceProfiledXvWAD(propName, inputMode, outputMode, xvVec.vec, wVec.vec, center.vec, psi.vec, dForcedXvW.vec) def calcdForcedStateTPsiAD(self, mode, Vec xv, Vec state, Vec psi, Vec prod): self._thisptr.calcdForcedStateTPsiAD(mode, xv.vec, state.vec, psi.vec, prod.vec) diff --git a/tests/Allrun b/tests/Allrun index 48887943..49d1c333 100755 --- a/tests/Allrun +++ b/tests/Allrun @@ -139,6 +139,7 @@ case $argm in runTests MphysAeroThermal runTests MphysAeroAcoustic runTests MphysAeroProp + runTestsSerial MphysAeroPropCoupled fi echo " " echo "************************************************************" @@ -243,6 +244,7 @@ case $argm in runTests MphysAeroThermal runTests MphysAeroAcoustic runTests MphysAeroProp + runTestsSerial MphysAeroPropCoupled fi echo " " echo "************************************************************" diff --git a/tests/refs/DAFoam_Test_MphysAeroPropCoupledRef.txt b/tests/refs/DAFoam_Test_MphysAeroPropCoupledRef.txt new file mode 100755 index 00000000..ace935e6 --- /dev/null +++ b/tests/refs/DAFoam_Test_MphysAeroPropCoupledRef.txt @@ -0,0 +1,84 @@ +Dictionary Key: force_balance +Dictionary Key: prop_shape +@value 41.76664173305865 0.0001 1e-06 +@value -558.5423009863406 0.0001 1e-06 +@value 468.9352557829879 0.0001 1e-06 +@value -159.9533579028434 0.0001 1e-06 +@value -85.83173476500616 0.0001 1e-06 +@value -627.2594281956608 0.0001 1e-06 +@value 139.8394731154231 0.0001 1e-06 +@value 1.641078539768139 0.0001 1e-06 +@value -508.437483312048 0.0001 1e-06 +@value -540.9345413274771 0.0001 1e-06 +@value -103.7505310591076 0.0001 1e-06 +@value -160.0809170878194 0.0001 1e-06 +@value 59.59297327588049 0.0001 1e-06 +@value 115.0314926447795 0.0001 1e-06 +@value 355.4713791986433 0.0001 1e-06 +@value 337.3777404218209 0.0001 1e-06 +@value -348.2919525815905 0.0001 1e-06 +@value -128.1976413071988 0.0001 1e-06 +@value 526.8054649421039 0.0001 1e-06 +@value 297.380661951973 0.0001 1e-06 +@value -109.8985893597628 0.0001 1e-06 +@value -63.39374501729827 0.0001 1e-06 +@value 343.0495751681397 0.0001 1e-06 +@value 265.3078158195785 0.0001 1e-06 +Dictionary Key: wing_twist +@value -0.5310089984568579 0.0001 1e-06 +Dictionary Key: lift +Dictionary Key: prop_shape +@value -0.8024086546636775 0.0001 1e-06 +@value 0.2512220982632938 0.0001 1e-06 +@value -0.5993622560496008 0.0001 1e-06 +@value 0.559818053679464 0.0001 1e-06 +@value -0.1608821470899166 0.0001 1e-06 +@value 0.7204955855958434 0.0001 1e-06 +@value -0.2781452271888569 0.0001 1e-06 +@value 0.3060357594902175 0.0001 1e-06 +@value -0.1604939168535242 0.0001 1e-06 +@value -0.06051870356633195 0.0001 1e-06 +@value -0.0318493451378556 0.0001 1e-06 +@value 0.5191412468734024 0.0001 1e-06 +@value -0.1273093746931327 0.0001 1e-06 +@value -0.5509556229383084 0.0001 1e-06 +@value 0.2865302728076517 0.0001 1e-06 +@value 0.5460179164873511 0.0001 1e-06 +@value 0.5997199965573426 0.0001 1e-06 +@value 0.3005645753458265 0.0001 1e-06 +@value -0.9429833818751224 0.0001 1e-06 +@value -0.743581567217245 0.0001 1e-06 +@value -0.2083466786540431 0.0001 1e-06 +@value -0.5687393911954824 0.0001 1e-06 +@value 0.2166622744547796 0.0001 1e-06 +@value 0.7030345493081527 0.0001 1e-06 +Dictionary Key: wing_twist +@value 64.2479085997736 0.0001 1e-06 +Dictionary Key: power +Dictionary Key: prop_shape +@value 21.6445351162704 0.0001 1e-06 +@value 9.725675429407039 0.0001 1e-06 +@value -16.45985162956367 0.0001 1e-06 +@value -66.37815016921245 0.0001 1e-06 +@value 11.59270649428473 0.0001 1e-06 +@value 29.12956440705982 0.0001 1e-06 +@value -29.55343893978312 0.0001 1e-06 +@value -4.154151715622334 0.0001 1e-06 +@value 22.22751437099393 0.0001 1e-06 +@value 25.97478895167304 0.0001 1e-06 +@value 19.63489794537844 0.0001 1e-06 +@value 19.12132585323177 0.0001 1e-06 +@value -50.0379006566045 0.0001 1e-06 +@value -32.61051693461268 0.0001 1e-06 +@value -32.06355269971075 0.0001 1e-06 +@value -38.10472255014534 0.0001 1e-06 +@value 0.5728692414433036 0.0001 1e-06 +@value 23.22198660570084 0.0001 1e-06 +@value 33.42775034342064 0.0001 1e-06 +@value 10.06941189040855 0.0001 1e-06 +@value 13.88100830884428 0.0001 1e-06 +@value 16.7880028383877 0.0001 1e-06 +@value 27.70287108517477 0.0001 1e-06 +@value 24.32984391195807 0.0001 1e-06 +Dictionary Key: wing_twist +@value -0 0.0001 1e-06 \ No newline at end of file diff --git a/tests/runTests_MphysAeroPropCoupled.py b/tests/runTests_MphysAeroPropCoupled.py new file mode 100755 index 00000000..5c081fb2 --- /dev/null +++ b/tests/runTests_MphysAeroPropCoupled.py @@ -0,0 +1,291 @@ +#!/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("./input/WingProp") + +if gcomm.rank == 0: + os.system("rm -rf 0 processor*") + +# aero setup +daOptionsWing = { + "designSurfaces": ["wing"], + "solverName": "DARhoSimpleFoam", + "primalMinResTol": 1.0e-10, + "primalBC": { + "U0": {"variable": "U", "patches": ["inout"], "value": [100.0, 0.0, 0.0]}, + "useWallFunction": True, + }, + "objFunc": { + "drag": { + "part1": { + "type": "force", + "source": "patchToFace", + "patches": ["wing"], + "directionMode": "fixedDirection", + "direction": [1.0, 0.0, 0.0], + "scale": 1.0, + "addToAdjoint": True, + } + }, + "lift": { + "part1": { + "type": "force", + "source": "patchToFace", + "patches": ["wing"], + "directionMode": "fixedDirection", + "direction": [0.0, 1.0, 0.0], + "scale": 1.0, + "addToAdjoint": True, + } + }, + }, + "adjEqnOption": { + "gmresRelTol": 1.0e-8, + "pcFillLevel": 1, + "jacMatReOrdering": "rcm", + }, + "normalizeStates": { + "U": 100, + "p": 101325, + "T": 300, + "nuTilda": 1e-3, + "phi": 1.0, + }, + "designVar": { + "twist": {"designVarType": "FFD"}, + "fvSource": {"designVarType": "Field", "fieldName": "fvSource", "fieldType": "vector"}, + }, + "wingProp": { + "prop1": { + "active": True, + "nForceSections": 5, + "axis": [0.0, 0.0, 1.0], + "actEps": 0.2, + "rotDir": "right", + "interpScheme": "gauss", + }, + "prop2": { + "active": True, + "nForceSections": 5, + "axis": [0.0, 0.0, 1.0], + "actEps": 0.2, + "rotDir": "left", + "interpScheme": "gauss", + }, + }, +} + +meshOptionsWing = { + "gridFile": os.getcwd(), + "fileType": "OpenFOAM", + "useRotations": False, + # 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]]], +} + +daOptionsProp = { + "solverName": "DARhoSimpleFoam", + "designSurfaces": ["blade"], + "primalMinResTol": 1e-12, + "primalBC": { + "MRF": -500.0, + }, + "objFunc": { + "thrust": { + "part1": { + "type": "force", + "source": "patchToFace", + "patches": ["blade"], + "directionMode": "fixedDirection", + "direction": [0.0, 0.0, 1.0], + "scale": 1.0, + "addToAdjoint": True, + } + }, + "power": { + "part1": { + "type": "power", + "source": "patchToFace", + "patches": ["blade"], + "axis": [0.0, 0.0, 1.0], + "center": [0.0, 0.0, 0.0], + "scale": -2.0 / 1000.0, # kW + "addToAdjoint": True, + } + }, + }, + "normalizeStates": {"U": 10, "p": 50.0, "nuTilda": 1e-3, "phi": 1.0}, + "adjEqnOption": {"gmresRelTol": 1.0e-10, "pcFillLevel": 1, "jacMatReOrdering": "rcm"}, + # Design variable setup + "designVar": { + "shapey": {"designVarType": "FFD"}, + }, + "decomposeParDict": {"preservePatches": ["per1", "per2"]}, + "wingProp": { + "prop": { + "active": True, + "nForceSections": 5, + "axis": [0.0, 0.0, 1.0], + "actEps": 0.2, + "rotDir": "right", + "interpScheme": "gauss", + }, + }, +} + +meshOptionsProp = { + "gridFile": os.getcwd(), + "fileType": "OpenFOAM", + # point and normal for the symmetry plane + "symmetryPlanes": [], +} + + +class Top(Multipoint): + def setup(self): + dafoam_builder_wing = DAFoamBuilder( + daOptionsWing, meshOptionsWing, scenario="aerodynamic", prop_coupling="Wing", run_directory="Wing" + ) + dafoam_builder_wing.initialize(self.comm) + + dafoam_builder_prop = DAFoamBuilder( + daOptionsProp, meshOptionsProp, scenario="aerodynamic", prop_coupling="Prop", run_directory="Prop" + ) + dafoam_builder_prop.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_wing", dafoam_builder_wing.get_mesh_coordinate_subsystem()) + self.add_subsystem("mesh_prop", dafoam_builder_prop.get_mesh_coordinate_subsystem()) + + # add the geometry component, we dont need a builder because we do it here. + self.add_subsystem("geometry_wing", OM_DVGEOCOMP(file="Wing/FFD/wingFFD.xyz", type="ffd")) + self.add_subsystem("geometry_prop", OM_DVGEOCOMP(file="Prop/FFD/localFFD.xyz", type="ffd")) + + self.mphys_add_scenario("cruise_prop", ScenarioAerodynamic(aero_builder=dafoam_builder_prop)) + self.mphys_add_scenario("cruise_wing", ScenarioAerodynamic(aero_builder=dafoam_builder_wing)) + + self.connect("mesh_wing.x_aero0", "geometry_wing.x_aero_in") + self.connect("geometry_wing.x_aero0", "cruise_wing.x_aero") + self.connect("mesh_prop.x_aero0", "geometry_prop.x_aero_in") + self.connect("geometry_prop.x_aero0", "cruise_prop.x_aero") + + # add an exec comp to average two drags + self.add_subsystem("force_balance", om.ExecComp("value=2*thrust-drag")) + + def configure(self): + + self.cruise_wing.aero_post.mphys_add_funcs() + self.cruise_prop.aero_post.mphys_add_funcs() + + # create geometric DV setup + points_wing = self.mesh_wing.mphys_get_surface_mesh() + points_prop = self.mesh_prop.mphys_get_surface_mesh() + + # add pointset + self.geometry_wing.nom_add_discipline_coords("aero", points_wing) + self.geometry_prop.nom_add_discipline_coords("aero", points_prop) + + # geometry setup + + # WING + def fvSource(val, DASolver): + for idxI, v in enumerate(val): + cellI = idxI // 3 + compI = idxI % 3 + DASolver.setFieldValue4LocalCellI(b"fvSource", v, cellI, compI) + # DASolver.updateBoundaryConditions(b"fvSource", b"vector") + + self.cruise_wing.coupling.solver.add_dv_func("fvSource", fvSource) + # no need to give fvSource to aero_post because we don't need its derivs + # self.cruise.aero_post.add_dv_func("fvSource", fvSource) + + nRefAxPtsWing = self.geometry_wing.nom_addRefAxis(name="wingAxis", xFraction=0.25, alignIndex="k") + + # Select all points + def twist_wing(val, geo): + for i in range(nRefAxPtsWing): + geo.rot_z["wingAxis"].coef[i] = -val[0] + + self.geometry_wing.nom_addGlobalDV(dvName="twist_wing", value=np.array([2]), func=twist_wing) + + # PROP + nShapesProp = self.geometry_prop.nom_addLocalDV(dvName="shape_prop") + + # add dvs to ivc and connect + self.dvs.add_output("twist_wing", val=np.array([2])) + self.dvs.add_output("shape_prop", val=np.array([0] * nShapesProp)) + self.dvs.add_output("prop1_center", val=np.array([-0.2, 0.2, 0.05])) + self.dvs.add_output("prop2_center", val=np.array([-0.2, -0.2, 0.05])) + self.dvs.add_output("prop_rot_center", val=np.array([0, 0, 0])) + self.dvs.add_output("prop_integral_force", val=np.array([5, 2])) + + for i in [1, 2]: + self.connect("cruise_prop.axial_force", "cruise_wing.prop%d_axial_force" % i) + self.connect("cruise_prop.tangential_force", "cruise_wing.prop%d_tangential_force" % i) + self.connect("cruise_prop.radial_location", "cruise_wing.prop%d_radial_location" % i) + self.connect("prop%d_center" % i, "cruise_wing.prop%d_prop_center" % i) + self.connect("prop_integral_force", "cruise_wing.prop%d_integral_force" % i) + + self.connect("twist_wing", "geometry_wing.twist_wing") + self.connect("shape_prop", "geometry_prop.shape_prop") + self.connect("prop_rot_center", "cruise_prop.prop_center") + + # define the design variables + self.add_design_var("twist_wing", lower=-10.0, upper=10.0, scaler=1.0) + self.add_design_var("shape_prop", lower=-0.1, upper=0.1, scaler=1.0) + + # add constraints and the objective + self.add_objective("cruise_prop.aero_post.power", scaler=1.0) + self.add_constraint("cruise_wing.aero_post.lift", equals=130, scaler=1.0) + + self.connect("cruise_wing.aero_post.drag", "force_balance.drag") + self.connect("cruise_prop.aero_post.thrust", "force_balance.thrust") + self.add_constraint("force_balance.value", equals=0, scaler=1.0) + + +prob = om.Problem(reports=None) +prob.model = Top() + +prob.setup(mode="rev") +om.n2(prob, show_browser=False, outfile="mphys_wing_prop.html") + +prob.run_model() +totals = prob.compute_totals() + +print(totals) + +if gcomm.rank == 0: + derivDict = {} + derivDict["force_balance"] = {} + derivDict["force_balance"]["prop_shape"] = totals[("force_balance.value", "dvs.shape_prop")][0] + derivDict["force_balance"]["wing_twist"] = totals[("force_balance.value", "dvs.twist_wing")][0] + derivDict["lift"] = {} + derivDict["lift"]["prop_shape"] = totals[("cruise_wing.aero_post.functionals.lift", "dvs.shape_prop")][0] + derivDict["lift"]["wing_twist"] = totals[("cruise_wing.aero_post.functionals.lift", "dvs.twist_wing")][0] + derivDict["power"] = {} + derivDict["power"]["prop_shape"] = totals[("cruise_prop.aero_post.functionals.power", "dvs.shape_prop")][0] + derivDict["power"]["wing_twist"] = totals[("cruise_prop.aero_post.functionals.power", "dvs.twist_wing")][0] + reg_write_dict(derivDict, 1e-4, 1e-6)