Overset update fixes

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(

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: >

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

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

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

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.

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

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()