Several features on tecplot io, slices, MG PC, cavitation constraints, and overset hole cutting.

.gitignore

@@ -42,6 +42,7 @@ tests/*out
 
 # for the .cgnsTimestep... files
 tests/output_files/*.cgns*
+tests/output_files/*.plt
 
 *.cgns
 *.txt

config/defaults/config.LINUX_INTEL_AVX2.mk

@@ -0,0 +1,55 @@
+# ----------------------------------------------------------------------
+# Config file for Intel ifort
+# ----------------------------------------------------------------------
+
+# ------- Define a possible parallel make (use PMAKE = make otherwise)--
+PMAKE = make -j 4
+
+# ------- Define the MPI Compilers--------------------------------------
+FF90 = mpiifort
+CC   = mpiicc
+
+# ------- Define Precision Flags ---------------------------------------
+# Options for Integer precision flags: -DUSE_LONG_INT
+# Options for Real precision flags: -DUSE_SINGLE_PRECISION, -DUSE_QUADRUPLE_PRECISION
+# Default (nothing specified) is 4 byte integers and 8 byte reals
+
+FF90_INTEGER_PRECISION_FLAG =
+FF90_REAL_PRECISION_FLAG    =
+CC_INTEGER_PRECISION_FLAG   =
+CC_REAL_PRECISION_FLAG      =
+
+# ------- Define CGNS Inlcude and linker flags -------------------------
+# Define the CGNS include directory and linking flags for the CGNS library.
+# We are assuming that HDF5 came from PETSc so it is included in ${PETSC_LIB}.
+# Otherwise you will have to specify the HDF5 library.
+CGNS_INCLUDE_FLAGS=-I$(CGNS_HOME)/include
+CGNS_LINKER_FLAGS=-L$(CGNS_HOME)/lib -lcgns
+
+# ------- Define Compiler Flags ----------------------------------------
+FF77_FLAGS   = -fPIC -r8
+FF90_FLAGS = $(FF77_FLAGS) -std08
+FFXX_OPT_FLAGS = -O2 -xCORE-AVX2
+C_FLAGS      = -fPIC -O -xCORE-AVX2
+
+# ------- Define Archiver  and Flags -----------------------------------
+AR       = ar
+AR_FLAGS = -rvs
+
+# ------- Define Linker Flags ------------------------------------------
+LINKER       = $(FF90)
+LINKER_FLAGS = -nofor-main
+
+# ------- Define Petsc Info ---
+include ${PETSC_DIR}/lib/petsc/conf/variables
+PETSC_INCLUDE_FLAGS=${PETSC_CC_INCLUDES} -I$(PETSC_DIR)
+PETSC_LINKER_FLAGS=${PETSC_LIB}
+
+# Combine flags from above -- don't modify here
+FF90_PRECISION_FLAGS = $(FF90_INTEGER_PRECISION_FLAG)$(FF90_REAL_PRECISION_FLAG)
+CC_PRECISION_FLAGS   = $(CC_INTEGER_PRECISION_FLAG) $(CC_REAL_PRECISION_FLAG)
+
+# Define potentially different python, python-config and f2py executables:
+PYTHON = python
+PYTHON-CONFIG = python3-config # use python-config for python 2
+F2PY = f2py

doc/conf.py

@@ -33,5 +33,6 @@
 
 # intersphinx
 intersphinx_mapping = {
-    "mach-aero": (f"https://mdolab-mach-aero.readthedocs-hosted.com/en/latest", None),
+    "mach-aero": ("https://mdolab-mach-aero.readthedocs-hosted.com/en/latest", None),
+    "pygeo": ("https://mdolab-pygeo.readthedocs-hosted.com/en/latest", None),
 }

doc/options.yaml

@@ -332,6 +332,12 @@ cavitationNumber:
   desc: >
     The -Cp value used to trigger the cavitation sensor.
 
+cpMinRho:
+  desc: >
+    The rho parameter used for the KS aggregation used for computing minimum Cp.
+    KS aggregation is used to compute the differentiable min Cp within the given family group.
+    A higher rho value will approach the actual min Cp more accurately, at the cost of a more nonlinear function.
+
 nCycles:
   desc: >
     The maximum number of "iterations" to run.
@@ -483,6 +489,116 @@ cutCallback:
     Therefore, we are flagging every cell where the y-component of the center coordinate (index 0:x, 1:y, 2:z) has a negative value.
     The flags array is initialized to zero internally by pyADflow, and a value of 1 for a cell tags it for flooding.
 
+explicitSurfaceCallback:
+  desc: >
+    A user-provided Python function that is used as a callback method to pass the information required for the explicit surface blanking method.
+    To be able to understand what this option does, you should really understand how the overset flooding method works.
+
+    With some configurations, the automatic flooding algorithm makes it extremely challenging to prevent the entire mesh from flooding.
+    The overlap of volume blocks can even be such that setting the ``"oversetPriority"`` option for each block is not enough.
+    To handle these extreme cases, we have a surface based explicit blanking method.
+    The idea is very similar to the ``cutCallback`` option above; this option is used to tag cells that are *inside* a closed geometry provided by the user.
+    This way, we try to eliminate as much of the flood seeds and flooded cells as possible.
+    The approach is not perfect; we only tag the cells that have at least one vertex inside the provided surface geometry as explicitly blanked (-4).
+    This results in two layers of fringe cells that are automatically added, which prevent flooding under most cases.
+    However, the flooding algorithm is actually more involved than this; normally, potential wall donors are tagged as flood seeds, not just the cells that have a vertex inside a wall BC.
+    Regardless, the method works well enough that it can give you a valid hole cutting with cases that are practically impossible to get a valid hole cutting with just the default approach.
+
+    The way user interacts with this method is through this callback function.
+    The reason we have the callback function is that we need the CGNS block names for each "computational block".
+    These names are used to determine which volume blocks will be explicitly blanked by which surface; we do not want every surface to act on every computational domain because that would result in most near wall cells getting tagged and would be very error prone.
+
+    To understand how this option is used, consider an example where we have a wing and a fuselage, in an overset mesh.
+    The wing and fuselage component meshes are joined together with a collar mesh and a background cartesian mesh.
+    Some of the wing mesh is inside the fuselage surface, and some of the fuselage mesh is overlapping with the wing surface.
+    Normally, this case would be handled by the automatic flooding, which works fairly well for most applications including this example, but here, we are just demonstrating the alternative approach.
+    With the regular flooding algorithm, some cells in these overlap regions would get tagged as flood seeds on both the wing and the fuselage.
+    Then, these flood seeds are used to flood the cells that are located inside these geometries; e.g. all of the wing volume cells that are inside the fuselage surface.
+
+    Instead of relying on the flooding, here, we want to explicitly blank cells that are inside the other component's surface.
+    To do this, we want to explicitly tag the wing mesh volume cells that have at least one vertex inside the fuselage surface, and similarly, tag the fuselage volume cells that have at least one vertex inside the wing surface.
+    The ``explicitSurfaceCallback`` option that would achieve this would look something like this:
+
+      .. code-block:: python
+
+        def explicitSurfaceCallback(CGNSZoneNameIDs):
+           # arrays to save the integer CGNS block IDs for each component
+           fuselage_blocks = []
+           wing_blocks = []
+
+           for name, id in CGNSZoneNameIDs.items():
+               # We can only get the blocks we want by checking the CGNS blocks' name.
+               # In this example, we don't want the explicit blanking algorithm to work
+               # on the cartesian background or the collar meshes.
+               if "wing" in name:
+                   wing_blocks.append(id)
+               elif "fuselage" in name:
+                   fuselage_blocks.append(id)
+
+           surf_dict = {
+               "wing": {
+                   # The surface file for each entry must be provided.
+                   # This is a closed plot3d surface mesh. It does not
+                   # necessarily need to be closed; we project the cell
+                   # vertices on this mesh and if the projection direction
+                   # is coming from the inside of the surface, the cell is
+                   # tagged. So a CFD mesh with a symmetry plane can
+                   # have a plot3d surface thats open at the symmetry plane.
+                   # The normals must always point outside. An easy way to
+                   # obtain this surface file is to convert the surface
+                   # CGNS grid used to generate the mesh into the plot3d
+                   # format using the cgns2plot3d method in cgnsutils.
+                   "surf_file": "wing_surface.x",
+
+                   # The block IDs we will search and blank using this surface.
+                   # Here, we only blank the fuselage cells that are inside the wing
+                   # and don't want to consider the background or collar meshes.
+                   "block_ids": fuselage_blocks,
+
+                   # Optional entry: coord_xfer
+                   # This is another callback function that performs a coordinate
+                   # transformation on the loaded plot3d mesh nodes. See the
+                   # addPointSet method in DVGeo.py at
+                   # https://github.com/mdolab/pygeo for the call signature.
+                   # We use this approach so that the users can have complete
+                   # control over the surface mesh nodes after they are loaded.
+                   "coord_xfer": coord_xfer,
+               },
+
+               # Similar entries for the fuselage.
+               "fuselage": {
+                   "surf_file": "fuselage_surface.x",
+                   "block_ids": wing_blocks,
+               },
+
+               # We can have multiple entries with the same surface to blank different
+               # subsets of CGNS blocks. The entries of this top level dictionary
+               # is just for bookkeeping, and is not important. The code prints
+               # which entry it is blanking as it works through the dictionary.
+               # With large meshes (millions of cells) on few number of processors (<10)
+               # this process can take up to a minute or so. The process should run quicker
+               # with more processors. The surface mesh is loaded on all processors
+               # so there is a memory limit on the surface mesh refinement.
+               "wing_self": {
+                   # Here, assume we also want to blank some of the wing cells because
+                   # they intersect the wing geometry itself somehow. E.g. very large splay.
+                   # However, we only want to blank cells that have a k-index that is
+                   # larger than 32 to avoid tagging cells near the wall. The
+                   # optional "kmin" entry is used to achieve this.
+                   "surf_file": "wing_surface.x",
+                   "block_ids": wing_blocks,
+                   "kmin": 32,
+               },
+           }
+
+           # pyADflow will loop over the entries of this dictionary and perform the
+           # explicit blanking with the information provided.
+           return surf_dict
+
+    The call signature of the ``coord_xfer`` option is documented in the DVGeometry's addPointset method:
+    :meth:`addPointset <pygeo:pygeo.parameterization.DVGeo.DVGeometry.addPointSet>`
+
+
 oversetUpdateMode:
   desc: >
     How to update the overset connectivity after mesh warping.
@@ -536,6 +652,11 @@ oversetPriority:
     A lower factor will encourage the usage of that block mesh.
     This option may be required to get the flooding algorithm working properly.
 
+oversetDebugPrint:
+  desc: >
+    Flag to enable or disable the debug printout from the overset algorithm when a hole cutting process fails.
+    This information is not useful anymore and users should obtain the volume solution after an unsuccessful hole cutting and debug the mesh by plotting the cells with incomplete connectivities, which receive an iblank value of -5.
+
 timeIntegrationScheme:
   desc: >
     The type of time integration scheme to use for unsteady analysis.
@@ -774,6 +895,12 @@ ANKUseTurbDADI:
     Only applies to decoupled ANK.
     If False, an internal ANK-like solver is used for the turbulence in the decoupled mode.
 
+ANKUseApproxSA:
+  desc: >
+    The ANK solver switches from the approximate to the exact SA implicit formulation when the solver switches from the first order to the second order routines (printed as SANK or CSANK).
+    Setting this flag to false will force the SA implicit treatment to always use the approximate formulation with the ANK solver variants.
+    This should have no effect on the final solution, but can help with convergence with challenging cases.
+
 ANKSwitchTol:
   desc: >
     Relative convergence in the residual norm before we switch to the ANK solver.
@@ -913,6 +1040,15 @@ ANKPCUpdateTol:
   desc: >
     If the ANK solver converges by this amount relative to the last nonlinear iteration where the preconditioner was updated, the preconditioner is updated again.
 
+ANKPCUpdateCutoff:
+  desc: >
+    The cutoff tolerance below which we disable the PC update based on nonlinear convergence introduced in the option above.
+    E.g. the default value of 1e-6 will trigger PC updates with the ANK if the solver converges by a certain amount, but past a relative L2 convergence of 1e-6, this algorithm is disabled and the preconditioner is only updated with the mechanisms other than relative convergence.
+    This is useful when CSANK is used, which can take steps that converge the L2 norm by multiple orders of magnitude in each iteration towards the end of the convergence.
+    Here, we do not need to update the PC at each solver iteration, so we disable it with this option.
+    This option only changes the PC update behavior with the coupled ANK variants.
+    Decoupled ANK variants don't get affected by this.
+
 ANKADPC:
   desc: >
     Whether or not to use the AD-based preconditioner for the ANK solver.
@@ -1263,7 +1399,8 @@ sepSensorSharpness:
 
 computeCavitation:
   desc: >
-    Whether or not to compute cavitation.
+    Whether or not to compute cavitation related cost functions, which are `cavitation` and `cpMin`.
+    If this option is not set to `True`, the code will return zero for these two cost functions.
 
 cavSensorOffset:
   desc: >

src/NKSolver/NKSolvers.F90

@@ -1672,6 +1672,7 @@ module ANKSolver
   real(kind=realType)   :: ANK_switchTol
   real(kind=realType)   :: ANK_divTol = 10
   logical :: ANK_useTurbDADI
+  logical :: ANK_useApproxSA
   real(kind=realType) :: ANK_turbcflscale
   logical :: ANK_useFullVisc
   logical :: ANK_ADPC
@@ -1686,6 +1687,7 @@ module ANKSolver
   real(kind=realType) :: ANK_secondOrdSwitchTol, ANK_coupledSwitchTol
   real(kind=realType) :: ANK_physLSTol, ANK_unstdyLSTol
   real(kind=realType) :: ANK_pcUpdateTol
+  real(kind=realType) :: ANK_pcUpdateCutoff
   real(kind=realType) :: lambda
   logical :: ANK_solverSetup=.False.
   integer(kind=intTYpe) :: ANK_iter
@@ -2093,6 +2095,12 @@ subroutine FormJacobianANK
                       ! get the global cell index
                       irow = globalCell(i, j, k)
 
+                      if (useCoarseMats) then
+                         do lvl=1, agmgLevels-1
+                            coarseRows(lvl+1) = coarseIndices(nn, lvl)%arr(i, j, k)
+                         end do
+                      end if
+
                       ! Add the contribution to the matrix in PETSc
                       call setBlock()
                    end do
@@ -3470,6 +3478,10 @@ subroutine ANKTurbSolveKSP
             rtol = min(ANK_rtol, rtol)
         end if
 
+        ! also check if we are using approxSA always
+        if (ANK_useApproxSA) &
+           approxSA = .True.
+
         ! Record the total residual and relative convergence for next iteration
         totalR_old = totalR
         rtolLast = rtol
@@ -3529,6 +3541,10 @@ subroutine ANKTurbSolveKSP
             approxSA = .False.
         end if
 
+        ! put back the approxsa flag if we were using it
+        if (ANK_useApproxSA) &
+          approxSA = .False.
+
         ! Compute the maximum step that will limit the change
         ! in SA variable to some user defined fraction.
         call physicalityCheckANKTurb(lambdaTurb)
@@ -3700,7 +3716,7 @@ subroutine ANKStep(firstCall)
     real(kind=realType) :: atol, val, v2, factK, gm1
     real(kind=alwaysRealType) :: rtol, totalR_dummy, linearRes, norm
     real(kind=alwaysRealType) :: resHist(ANK_maxIter+1)
-    real(kind=alwaysRealType) :: unsteadyNorm, unsteadyNorm_old
+    real(kind=alwaysRealType) :: unsteadyNorm, unsteadyNorm_old, rel_pcUpdateTol
     logical :: correctForK, LSFailed
 
     ! Enter this check if this is the first ANK step OR we are switching to the coupled ANK solver
@@ -3766,13 +3782,28 @@ subroutine ANKStep(firstCall)
        ANK_iter = ANK_iter + 1
     end if
 
+    ! figure out if we want to scale the ANKPCUpdateTol
+    if (.not. ANK_coupled) then
+      rel_pcUpdateTol = ANK_pcUpdateTol
+    else
+      ! for coupled ANK, we dont want to update the PC as frequently,
+      ! so we reduce the relative tol by 4 orders of magnitude,
+      ! *if* we are converged past pc update cutoff wrt free stream already
+      if (totalR / totalR0 .lt. ANK_pcUpdateCutoff) then
+         rel_pcUpdateTol = ANK_pcUpdateTol * 1e-4_realType
+      else
+         ! if we are not that far down converged, use the option directly
+         rel_pcUpdateTol = ANK_pcUpdateTol
+      end if
+    end if
+
     ! Compute the norm of rVec, which is identical to the
     ! norm of the unsteady residual vector.
     call VecNorm(rVec, NORM_2, unsteadyNorm_old, ierr)
     call EChk(ierr, __FILE__, __LINE__)
 
     ! Determine if if we need to form the Preconditioner
-    if (mod(ANK_iter, ANK_jacobianLag) == 0 .or. totalR/totalR_pcUpdate < ANK_pcUpdateTol) then
+    if (mod(ANK_iter, ANK_jacobianLag) == 0 .or. totalR/totalR_pcUpdate < rel_pcUpdateTol) then
 
        ! First of all, update the minimum cfl wrt the overall convergence
        ANK_CFLMin = min(ANK_CFLLimit, ANK_CFLMinBase*(totalR0/totalR)**ANK_CFLExponent)
@@ -3887,6 +3918,10 @@ subroutine ANKStep(firstCall)
        rtol = min(ANK_rtol, rtol)
     end if
 
+    ! also check if we are using approxSA always
+    if (ANK_useApproxSA) &
+      approxSA = .True.
+
     ! Record the total residual and relative convergence for next iteration
     totalR_old = totalR
     rtolLast = rtol
@@ -3951,6 +3986,10 @@ subroutine ANKStep(firstCall)
 
     end if
 
+   ! put back the approxsa flag if we were using it
+    if (ANK_useApproxSA) &
+      approxSA = .False.
+
     ! Compute the maximum step that will limit the change in pressure
     ! and energy to some user defined fraction.
     call physicalityCheckANK(lambda)

src/adjoint/adjointUtils.F90

@@ -1614,6 +1614,7 @@ subroutine destroyPETScVars
     use constants
     use ADjointPETSc, only : dRdWT, dRdwPreT, adjointKSP, adjointPETScVarsAllocated
     use inputAdjoint, only : approxPC
+    use agmg, only : destroyAGMG
     use utils, only : EChk
     implicit none
 
@@ -1632,6 +1633,9 @@ subroutine destroyPETScVars
 
        call KSPDestroy(adjointKSP, ierr)
        call EChk(ierr,__FILE__,__LINE__)
+
+       call destroyAGMG()
+
        adjointPETScVarsAllocated = .False.
     end if