From b3ed666e6d7353ada91f539e56f1ffcf4577d9d6 Mon Sep 17 00:00:00 2001 From: Sabet Seraj <48863473+sseraj@users.noreply.github.com> Date: Wed, 27 Jul 2022 23:34:15 -0400 Subject: [PATCH] Overset update fixes (#225) * fixed cutCallBack signature for full update * allocated xCen for fast overset update in parallel * re-enable fast overset update in AD routines * print on root only for full update * update connectivity for fast mode only in blocketteRes * added overset update tests * renamed test_zipper * updated docs * add cutCallback to test * added fast update fraction check * modified fraction check to work with complexify * loosened fracTol --- adflow/pyADflow.py | 3 +- doc/options.yaml | 6 +- src/NKSolver/blockette.F90 | 4 +- src/adjoint/masterRoutines.F90 | 32 +++--- src/overset/finalOversetCommStructures.F90 | 4 +- src/overset/oversetAPI.F90 | 8 +- src/overset/oversetCommUtilites.F90 | 17 ++- tests/reg_tests/test_overset.py | 115 +++++++++++++++++++++ tests/reg_tests/test_zipper.py | 63 ----------- 9 files changed, 164 insertions(+), 88 deletions(-) create mode 100644 tests/reg_tests/test_overset.py delete mode 100644 tests/reg_tests/test_zipper.py diff --git a/adflow/pyADflow.py b/adflow/pyADflow.py index 6eda094f9..04dba6a67 100644 --- a/adflow/pyADflow.py +++ b/adflow/pyADflow.py @@ -3759,7 +3759,8 @@ def updateGeometryInfo(self, warpMesh=True): cutCallBack = self.getOption("cutCallBack") if cutCallBack is not None: xCen = self.adflow.utils.getcellcenters(1, n).T - cutCallBack(xCen, flag) + cellIDs = self.adflow.utils.getcellcgnsblockids(1, n) + cutCallBack(xCen, self.CGNSZoneNameIDs, cellIDs, flag) # Verify previous mesh failures self.adflow.killsignals.routinefailed = self.comm.allreduce( diff --git a/doc/options.yaml b/doc/options.yaml index 8a3481764..e025ccb85 100644 --- a/doc/options.yaml +++ b/doc/options.yaml @@ -491,11 +491,15 @@ oversetUpdateMode: This is the cheapest option to use when the entire overset mesh is warped together, such as with USMesh in IDWarp. fast: > Only the weights are updated and the donors are unchanged. - This is usually a fairly small correction and is fast to run. This is used when the entire overset mesh is warped together, such as with USMesh in IDWarp. + This is fast to run but only works for small design changes. + Large design changes can shift the position of interpolate cells relative to their donors enough to cause nonphysical interpolation weights. + In this case, an error is raised. full: > The overset connectivity is recomputed from stratch. This is usually only used when the component meshes are warped independently, such as with MultiUSMesh in IDWarp. + The hole cutting is likely to fail if there are large design changes. + Even if the hole cutting works, the aerodynamic derivatives will be inaccurate because the hole cutting routines are not differentiated. nRefine: desc: > diff --git a/src/NKSolver/blockette.F90 b/src/NKSolver/blockette.F90 index 9952a803f..8293db443 100644 --- a/src/NKSolver/blockette.F90 +++ b/src/NKSolver/blockette.F90 @@ -189,9 +189,7 @@ subroutine blocketteRes(useDissApprox, useViscApprox, useUpdateIntermed, useFlow do sps=1, nTimeIntervalsSpectral ! Update overset connectivity if necessary - if (oversetPresent .and. & - (oversetUpdateMode == updateFast .or. & - oversetUpdateMode == updateFull)) then + if (oversetPresent .and. oversetUpdateMode == updateFast) then call updateOversetConnectivity(1_intType, sps) end if end do diff --git a/src/adjoint/masterRoutines.F90 b/src/adjoint/masterRoutines.F90 index cfb166a21..1f01bfcd7 100644 --- a/src/adjoint/masterRoutines.F90 +++ b/src/adjoint/masterRoutines.F90 @@ -83,9 +83,7 @@ subroutine master(useSpatial, famLists, funcValues, forces, & do sps=1, nTimeIntervalsSpectral ! Update overset connectivity if necessary - if (oversetPresent .and. & - (oversetUpdateMode == updateFast .or. & - oversetUpdateMode == updateFull)) then + if (oversetPresent .and. oversetUpdateMode == updateFast) then call updateOversetConnectivity(1_intType, sps) end if end do @@ -241,7 +239,7 @@ subroutine master_d(wdot, xdot, forcesDot, dwDot, famLists, funcValues, funcValu use constants use diffsizes, only : ISIZE1OFDrfbcdata, ISIZE1OFDrfviscsubface - use communication, only : adflow_comm_world + use communication, only : adflow_comm_world, myID use iteration, only : currentLevel use BCExtra_d, only : applyAllBC_Block_d use inputAdjoint, only : viscPC @@ -353,11 +351,14 @@ subroutine master_d(wdot, xdot, forcesDot, dwDot, famLists, funcValues, funcValu do sps=1, nTimeIntervalsSpectral ! Update overset connectivity if necessary - if (oversetPresent .and. & - (oversetUpdateMode == updateFast .or. & - oversetUpdateMode == updateFull)) then - print *,'Full overset update derivatives not implemented' - !call updateOversetConnectivity_d(1_intType, sps) + if (oversetPresent) then + if (oversetUpdateMode == updateFast) then + call updateOversetConnectivity_d(1_intType, sps) + else if (oversetUpdateMode == updateFull) then + if (myID == 0) then + print *,'Full overset update derivatives not implemented' + end if + end if end if end do @@ -879,11 +880,14 @@ subroutine master_b(wbar, xbar, extraBar, forcesBar, dwBar, nState, famLists, & do sps=1, nTimeIntervalsSpectral ! Update overset connectivity if necessary - if (oversetPresent .and. & - (oversetUpdateMode == updateFast .or. & - oversetUpdateMode == updateFull)) then - print *,'Full overset update derivatives not implemented' - !call updateOversetConnectivity_b(1_intType, sps) + if (oversetPresent) then + if (oversetUpdateMode == updateFast) then + call updateOversetConnectivity_b(1_intType, sps) + else if (oversetUpdateMode == updateFull) then + if (myID == 0) then + print *,'Full overset update derivatives not implemented' + end if + end if end if end do ! Now the adjoint of the coordinate exhcange diff --git a/src/overset/finalOversetCommStructures.F90 b/src/overset/finalOversetCommStructures.F90 index ac621d922..9b6b03af1 100644 --- a/src/overset/finalOversetCommStructures.F90 +++ b/src/overset/finalOversetCommStructures.F90 @@ -303,7 +303,8 @@ subroutine finalOversetCommStructures(level, sps) commPatternOverset(level, sps)%sendList(iSendProc)%block(n), & commPatternOverset(level, sps)%sendList(iSendProc)%indices(n, 3), & commPatternOverset(level, sps)%sendList(iSendProc)%interp(n, 8), & - commPatternOverset(level, sps)%sendList(iSendProc)%interpd(n, 8)) + commPatternOverset(level, sps)%sendList(iSendProc)%interpd(n, 8), & + commPatternOverset(level, sps)%sendList(iSendProc)%xCen(n, 3)) ! Now set the data do i=1, n @@ -323,6 +324,7 @@ subroutine finalOversetCommStructures(level, sps) call fracToWeights(realRecvBuf(jj+1:jj+3), & commPatternOverset(level, sps)%sendList(iSendProc)%interp(i, :)) commPatternOverset(level, sps)%sendList(iSendProc)%interpd(i, :) = zero + commPatternOverset(level, sps)%sendList(iSendProc)%xCen(i, :) = zero jj = jj + 3 end do end if diff --git a/src/overset/oversetAPI.F90 b/src/overset/oversetAPI.F90 index fe667c1cf..952838238 100644 --- a/src/overset/oversetAPI.F90 +++ b/src/overset/oversetAPI.F90 @@ -2101,16 +2101,16 @@ end subroutine setExplicitHoleCut subroutine updateOverset(flag, n, closedFamList, nFam) ! This is the main gateway routine for updating the overset - ! connecitivty during a solution. It is wrapped and intended to be + ! connectivity during a solution. It is wrapped and intended to be ! called from Python. What this routine does depends on the value ! of oversetUpdateMode: ! updateFrozen: Nothing happens. The initial - ! connecitivty computed during initialzation is kept. + ! connectivity computed during initialization is kept. - ! updteFast: Update just the weight, but leave the donors + ! updateFast: Update just the weight, but leave the donors ! unchanged. This is only applicable when the entire mesh is - ! warped at the same time with pyWarpuStruct. + ! warped at the same time like with USMesh in IDWarp. ! updateFull: Complete from scratch update. Run the full ! oversetComm routine. diff --git a/src/overset/oversetCommUtilites.F90 b/src/overset/oversetCommUtilites.F90 index ae565166a..8c6dc4dc4 100644 --- a/src/overset/oversetCommUtilites.F90 +++ b/src/overset/oversetCommUtilites.F90 @@ -1740,7 +1740,7 @@ subroutine updateOversetConnectivity(level, sps) ! mesh. It does *not* completely redo the connectivity. Rather, a ! newton search on the existing donors are performed using the ! updated coordinates. This type of update is only applicable if the - ! entire volume mesh is warped as one a-la pyWarpUstruct. This + ! entire volume mesh is warped as one like with USMesh in IDWarp. This ! actually ends up being a fairly small correction most of the time, ! however, desipite looks to the contrary is actually quite fast to ! run. @@ -1767,6 +1767,9 @@ subroutine updateOversetConnectivity(level, sps) real(kind=realType) :: frac(3), frac0(3), xCen(3) integer(kind=intType), dimension(8), parameter :: indices=(/1,2,4,3,5,6,8,7/) + ! Set a tolerance for checking whether fractions are between 0 and 1 + real(kind=realType) :: fracTol=1e-4 + ! Pointers to the overset comms to make it easier to read commPattern => commPatternOverset(level, sps) internal => internalOverset(level, sps) @@ -1913,6 +1916,12 @@ subroutine updateOversetConnectivity(level, sps) call newtonUpdate(xCen, & flowDoms(d1, level, sps)%x(i1-1:i1+1, j1-1:j1+1, k1-1:k1+1, :), frac0, frac) + ! Check if the fractions are between 0 and 1 + if (MAXVAL(frac) > one + fracTol .or. MINVAL(frac) < zero - fracTol) then + print *, "Invalid overset connectivity update. Use 'frozen' or 'full' oversetUpdateMode instead." + error stop + end if + ! Set the new weights call fracToWeights(frac, internal%donorInterp(i, :)) enddo localInterp @@ -1951,6 +1960,12 @@ subroutine updateOversetConnectivity(level, sps) call newtonUpdate(xCen, & flowDoms(d2, level, sps)%x(i2-1:i2+1, j2-1:j2+1, k2-1:k2+1, :), frac0, frac) + ! Check if the fractions are between zero and one + if (MAXVAL(frac) > one + fracTol .or. MINVAL(frac) < zero - fracTol) then + print *, "Invalid overset connectivity update. Use 'frozen' or 'full' oversetUpdateMode instead." + error stop + end if + ! Set the new weights call fracToWeights(frac, commPattern%sendList(ii)%interp(j, :)) enddo diff --git a/tests/reg_tests/test_overset.py b/tests/reg_tests/test_overset.py new file mode 100644 index 000000000..7419e11d1 --- /dev/null +++ b/tests/reg_tests/test_overset.py @@ -0,0 +1,115 @@ +import unittest +import os +import copy +from parameterized import parameterized_class +from adflow import ADFLOW +from reg_default_options import adflowDefOpts, IDWarpDefOpts +from reg_aeroproblems import ap_simple_cart_cube +from idwarp import USMesh +from pygeo import DVGeometry + +baseDir = os.path.dirname(os.path.abspath(__file__)) + +# Tests for overset and zipper meshes using a mesh with two overlapping cubes + +# This does not affect the hole cutting for this case but tests that cutCallback is working +def cutCallback(xCen, CGNSZoneNameIDs, cellIDs, flags): + + # Blank cells in the negative y-axis + flags[xCen[:, 1] < 0] = 1 + + +gridFile = os.path.join(baseDir, "../../input_files/cube_overset.cgns") +commonTestOptions = { + "gridFile": gridFile, + "equationType": "RANS", + "writeTecplotSurfaceSolution": True, + "monitorVariables": ["cpu", "resrho", "resturb", "cd"], + "volumeVariables": ["resrho", "resturb", "cp", "mach", "blank"], + "mgcycle": "sg", + "L2Convergence": 1e-13, + "nCycles": 500, + "useANKSolver": True, + "useNKSolver": True, + "nearWallDist": 1.0, + "cutCallback": cutCallback, +} + +test_params = [ + { + "name": "frozen", + "options": { + "oversetUpdateMode": "frozen", + }, + }, + { + "name": "fast", + "options": { + "oversetUpdateMode": "fast", + }, + }, + { + "name": "full", + "options": { + "oversetUpdateMode": "full", + }, + }, +] + + +@parameterized_class(test_params) +class TestCubeOverset(unittest.TestCase): + N_PROCS = 2 + + def setUp(self): + + super().setUp() + + # Start with the default testing options dictionary + options = copy.copy(adflowDefOpts) + + # Set the output directory + options["outputDirectory"] = os.path.join(baseDir, options["outputDirectory"]) + + # These are the modified options common to these tests + options.update(commonTestOptions) + + # Add the parameterized options + options.update(self.options) + + # Define the AeroProblem + self.ap = copy.deepcopy(ap_simple_cart_cube) + + # Create the solver + self.CFDSolver = ADFLOW(options=options, debug=False) + + # Create mesh warping object + meshOptions = copy.copy(IDWarpDefOpts) + meshOptions.update({"gridFile": options["gridfile"]}) + self.CFDSolver.setMesh(USMesh(options=meshOptions)) + + # Create geometry object + ffdFile = os.path.join(baseDir, "../../input_files/cube_overset_ffd.xyz") + DVGeo = DVGeometry(ffdFile) + DVGeo.addLocalDV("shape", lower=-0.5, upper=0.5, axis="z") + self.CFDSolver.setDVGeo(DVGeo) + + # Apply shape change + dvDict = self.CFDSolver.DVGeo.getValues() + dvDict["shape"][4] = 0.1 + self.CFDSolver.setAeroProblem(self.ap) + self.CFDSolver.DVGeo.setDesignVars(dvDict) + + def test_convergence(self): + + # Solve the flow + self.CFDSolver(self.ap) + + # Check that the flow converged + funcs = {} + self.CFDSolver.checkSolutionFailure(self.ap, funcs) + self.assertFalse(funcs["fail"]) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/reg_tests/test_zipper.py b/tests/reg_tests/test_zipper.py deleted file mode 100644 index 8d7704898..000000000 --- a/tests/reg_tests/test_zipper.py +++ /dev/null @@ -1,63 +0,0 @@ -import unittest -import os -import copy -from adflow import ADFLOW -from reg_default_options import adflowDefOpts -from reg_aeroproblems import ap_simple_cart_cube - - -baseDir = os.path.dirname(os.path.abspath(__file__)) - -# Tests for zipper meshes using a mesh with two overlapping cubes - -gridFile = os.path.join(baseDir, "../../input_files/cube_overset.cgns") -commonTestOptions = { - "gridFile": gridFile, - "equationType": "RANS", - "writeTecplotSurfaceSolution": True, - "monitorVariables": ["cpu", "resrho", "resturb", "cd"], - "volumeVariables": ["resrho", "resturb", "cp", "mach", "blank"], - "mgcycle": "sg", - "L2Convergence": 1e-13, - "nCycles": 500, - "useANKSolver": True, - "useNKSolver": True, - "nearWallDist": 1.0, -} - - -class TestCubeOverset(unittest.TestCase): - N_PROCS = 2 - - def setUp(self): - - super().setUp() - - # Start with the default testing options dictionary - options = copy.copy(adflowDefOpts) - - # Set the output directory - options["outputDirectory"] = os.path.join(baseDir, options["outputDirectory"]) - - # These are the modified options common to these tests - options.update(commonTestOptions) - - # Define the AeroProblem - self.ap = copy.deepcopy(ap_simple_cart_cube) - - # Create the solver - self.CFDSolver = ADFLOW(options=options, debug=False) - - def test_convergence(self): - - # Solve the flow - self.CFDSolver(self.ap) - - # Check that the flow converged - funcs = {} - self.CFDSolver.checkSolutionFailure(self.ap, funcs) - self.assertFalse(funcs["fail"]) - - -if __name__ == "__main__": - unittest.main()