Skip to content

Commit

Permalink
Several features on tecplot io, slices, MG PC, cavitation constraints…
Browse files Browse the repository at this point in the history 
…, and overset hole cutting. (#231)

* hack to get cpmin with ks

* differentiated the hacky version

* fixed adjoint time print in Fortran

* use modified Gram-Schmidt

* less janky cavitation

* implemented orthogonalization option

* minor fix to get mg pc working with coupled ank

* new config file to do the avx2 stuff by default

* added moment and related outputs to the lift distribution file

* rearranged and updated the slice integration to use the local quarter chord for the moment computations

* added the missing deallocate call to mg pc

* adding cylindrical slicing. WIP

* added the missing mgpc destroy routine

* properly implemented the cavitation constraint

* changed the behavior of cavitation number option

* updates to the slice business

* fixes for the bad merge

* added optional displacement vector when adding integration surfaces

* added cutoff tolerance for the pc update algorithm based on rel convergence. also fixed the mgpc setup for ank turb ksp

* reverted the agmg changes for turbksp. needs a proper implementation to handle two different agmg pcs kicking around at the same time

* added option to always use approx sa with ank solver variants

* added the same update to the turbksp solver

* added python option to disable the overset debug print of connectivity errors

* removed dumb typo

* reran tapenade

* added arbitrary slices

* added the option to disable lift and slices in writesolution routine

* added function to explicitly flag overset cells that are inside provided surface geometries

* added a few more options to the explicit surface cutting

* re-implementing the ks-based cpmin. added the original cavitation sensor back

* reran tapenade

* more fixes to recover the og cavitation sensor

* remove high quality debug print

* variable cleanup

* bug fixes

* complexify fixes

* fix black

* add cavitation test trained with the og cavitation sensor

* bugfixes and reran tapenade

* updated cavitation test to include the ks-based cpmin

* formatting

* update cavitation tests

* further increase the weird test tolerance

* added the explicit surface blanking doc in options

* added the coord transfer changes

* added missing option docs

* more docs

* minor modification to options docs

* added note on plot3d file creation

* bug fix for explicit hole cutting

* addressing sabet's review part 1

* addressing sabet's review part 2: fix the import comment in tests

* part 3: overset API typos

* part 4: applying the same fixes to the actuator region code, which is where the overset api changes mostly came from

* part 5: comments and typos on the slice changes

* part 5: trying to fix the cylindrical slice edge case

* part 5.1: divide the strings up

* part 6: surface integration modifications

* clarified computecavitation option behavior

* changes for the cpmin computation so it returns zero when computeCavitation is set to False.

* rerun tapenade

* fixed variable declaration order. added test output to gitignore.

* black formatting

* updated avx2 config file to be consistent with the newer config files

* minor update in docstrings

* cleaned up slice names

* black formatting

* black formatter is gaslighting me

* first pass at eirikur's comments

* reduced some code duplication between slice addition methods

* camelcase updates

* rename cpmin_exact to cpmin_family

* black fix and reran tapenade

* fixed cpmin computation for TS cases

* missed a capitalization

* changed the number of default cylindrical slices

Co-authored-by: sseraj <[email protected]>
  • Loading branch information
@anilyil @sseraj
anilyil and sseraj authored Dec 16, 2022
1 parent b4ceb9a commit d5fe462
Show file tree
Hide file tree
Showing 31 changed files with 1,997 additions and 291 deletions.
1 change: 1 addition & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ tests/*out

# for the .cgnsTimestep... files
tests/output_files/*.cgns*
tests/output_files/*.plt

*.cgns
*.txt
Expand Down
361 changes: 327 additions & 34 deletions adflow/pyADflow.py
View file

Large diffs are not rendered by default.

55 changes: 55 additions & 0 deletions config/defaults/config.LINUX_INTEL_AVX2.mk
View file
Original file line number Diff line number Diff line change
@@ -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
3 changes: 2 additions & 1 deletion doc/conf.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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),
}
139 changes: 138 additions & 1 deletion doc/options.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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.
"surfFile": "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.
"blockIDs": fuselage_blocks,
# Optional entry: coordXfer
# 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.
"coordXfer": coordXfer,
},
# Similar entries for the fuselage.
"fuselage": {
"surfFile": "fuselage_surface.x",
"blockIDs": 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.
"surfFile": "wing_surface.x",
"blockIDs": 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 ``coordXfer`` 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.
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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: >
Expand Down
43 changes: 41 additions & 2 deletions src/NKSolver/NKSolvers.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
4 changes: 4 additions & 0 deletions src/adjoint/adjointUtils.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -1632,6 +1633,9 @@ subroutine destroyPETScVars

call KSPDestroy(adjointKSP, ierr)
call EChk(ierr,__FILE__,__LINE__)

call destroyAGMG()

adjointPETScVarsAllocated = .False.
end if

Expand Down
Loading

0 comments on commit d5fe462

Please sign in to comment.