Skip to content

Commit

Permalink
Overset update fixes (#225)
Browse files Browse the repository at this point in the history
* 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
  • Loading branch information
sseraj authored Jul 28, 2022
1 parent ea9b0ed commit b3ed666
Show file tree
Hide file tree
Showing 9 changed files with 164 additions and 88 deletions.
3 changes: 2 additions & 1 deletion adflow/pyADflow.py
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down
6 changes: 5 additions & 1 deletion doc/options.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -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: >
Expand Down
4 changes: 1 addition & 3 deletions src/NKSolver/blockette.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
32 changes: 18 additions & 14 deletions src/adjoint/masterRoutines.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down
4 changes: 3 additions & 1 deletion src/overset/finalOversetCommStructures.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
8 changes: 4 additions & 4 deletions src/overset/oversetAPI.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
17 changes: 16 additions & 1 deletion src/overset/oversetCommUtilites.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
115 changes: 115 additions & 0 deletions tests/reg_tests/test_overset.py
Original file line number Diff line number Diff line change
@@ -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()
63 changes: 0 additions & 63 deletions tests/reg_tests/test_zipper.py

This file was deleted.

0 comments on commit b3ed666

Please sign in to comment.