Skip to content

Commit

Permalink
address Eirikur's comments
Browse files Browse the repository at this point in the history
  • Loading branch information
@anilyil
anilyil committed Jun 20, 2024
1 parent c88f8c8 commit f8e1c0b
Show file tree
Hide file tree
Showing 5 changed files with 120 additions and 114 deletions.
126 changes: 74 additions & 52 deletions adflow/pyADflow.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -302,11 +302,14 @@ def __init__(self, comm=None, options=None, debug=False, dtype="d"):
flag = numpy.zeros(n)

# Call the user supplied callback if necessary
self.oversetCutCallback(flag)
self._oversetCutCallback(flag)
cutCallBackTime = time.time()

# also run through the surface callback routine
self.oversetExplicitSurfaceCallback(flag, firstCall=True)
# the initialization is done separately in case we want to do
# full overset updates later on.
self._initializeExplicitSurfaceCallback()
self._oversetExplicitSurfaceCallback(flag)
explicitSurfaceCutTime = time.time()

# Need to reset the oversetPriority option since the CGNSGrid
Expand Down Expand Up @@ -4336,10 +4339,10 @@ def updateGeometryInfo(self, warpMesh=True):
# Only need to call the cutCallBack and regenerate the zipper mesh
# if we're doing a full update.
if self.getOption("oversetUpdateMode") == "full":
self.oversetCutCallback(flag)
self._oversetCutCallback(flag)

# also run the explicit surface blanking
self.oversetExplicitSurfaceCallback(flag, firstCall=False)
self._oversetExplicitSurfaceCallback(flag)

# Verify previous mesh failures
self.adflow.killsignals.routinefailed = self.comm.allreduce(
Expand Down Expand Up @@ -4367,38 +4370,82 @@ def updateGeometryInfo(self, warpMesh=True):
)
self.adflow.killsignals.fatalfail = self.adflow.killsignals.routinefailed

def oversetCutCallback(self, flag):
# this routine goes through the explicit callback routine provided by the user and
# modifies the flag array in place
def _oversetCutCallback(self, flag):
"""This routine goes through the explicit callback routine provided by the user
and modifies the flag array in place. The most common use case for this is to
blank out the cells on the wrong side of the symmetry plane.
Parameters
----------
flag : ndarray
Array that is used as a mask to select the explicitly blanked cells.
This is modified in place.
"""
cutCallBack = self.getOption("cutCallBack")
if cutCallBack is not None:
n = len(flag)
xCen = self.adflow.utils.getcellcenters(1, n).T
cellIDs = self.adflow.utils.getcellcgnsblockids(1, n)
cutCallBack(xCen, self.CGNSZoneNameIDs, cellIDs, flag)

def oversetExplicitSurfaceCallback(self, flag, firstCall=True):
# exclude the cells inside closed surfaces if we are provided with them
explicitSurfaceCallback = self.getOption("explicitSurfaceCallback")
if explicitSurfaceCallback is not None:
# the user wants to exclude cells that lie within a list of surfaces.
def _initializeExplicitSurfaceCallback(self):
"""Rotuine that loads the external surfaces provided by the user for explicit blanking.
We can do this just once because there may be subsequent calls with the same surfaces.
"""

# get the number of cells
n = len(flag)
self.explicitSurfaceCallback = self.getOption("explicitSurfaceCallback")

if firstCall:
# first, call the callback function with cgns zone name IDs.
# this need to return us a dictionary with the surface mesh information,
# as well as which blocks in the cgns mesh to include in the search.
# we dont need to update this dictionary on subsequent calls.
self.blankingSurfDict = explicitSurfaceCallback(self.CGNSZoneNameIDs)
blankingSurfData = {}
if self.explicitSurfaceCallback is not None:
# first, call the callback function with cgns zone name IDs.
# this need to return us a dictionary with the surface mesh information,
# as well as which blocks in the cgns mesh to include in the search.
# we dont need to update this dictionary on subsequent calls.
self.blankingSurfDict = self.explicitSurfaceCallback(self.CGNSZoneNameIDs)

Check warning on line 4403 in adflow/pyADflow.py

View check run for this annotation

Codecov / codecov/patch

adflow/pyADflow.py#L4403 

Added line #L4403 was not covered by tests

# also keep track of a second dictionary that saves the surface info. We do this
# separately because the surface files can be shared across different blanking calls.
# in this case, we dont want to process the surfaces multiple times.
blankingSurfData = {}

Check warning on line 4408 in adflow/pyADflow.py

View check run for this annotation

Codecov / codecov/patch

adflow/pyADflow.py#L4408 

Added line #L4408 was not covered by tests

for surf in self.blankingSurfDict:

Check warning on line 4410 in adflow/pyADflow.py

View check run for this annotation

Codecov / codecov/patch

adflow/pyADflow.py#L4410 

Added line #L4410 was not covered by tests
# this is the plot3d surface that defines the closed volume
surfFile = self.blankingSurfDict[surf]["surfFile"]

Check warning on line 4412 in adflow/pyADflow.py

View check run for this annotation

Codecov / codecov/patch

adflow/pyADflow.py#L4412 

Added line #L4412 was not covered by tests

# if this is the first call, we need to load the surface meshes.
# we might have duplicate mesh files, so no need to load them again
# if we have already loadded one copy
if surfFile not in blankingSurfData:

Check warning on line 4417 in adflow/pyADflow.py

View check run for this annotation

Codecov / codecov/patch

adflow/pyADflow.py#L4417 

Added line #L4417 was not covered by tests
# optional coordinate transformation to do general manipulation of the coordinates
if "coordXfer" in self.blankingSurfDict[surf]:
coordXfer = self.blankingSurfDict[surf]["coordXfer"]

Check warning on line 4420 in adflow/pyADflow.py

View check run for this annotation

Codecov / codecov/patch

adflow/pyADflow.py#L4419-L4420 

Added lines #L4419 - L4420 were not covered by tests
else:
coordXfer = None

Check warning on line 4422 in adflow/pyADflow.py

View check run for this annotation

Codecov / codecov/patch

adflow/pyADflow.py#L4422 

Added line #L4422 was not covered by tests

# read the plot3d surface
pts, conn = self._readPlot3DSurfFile(surfFile, convertToTris=False, coordXfer=coordXfer)

Check warning on line 4425 in adflow/pyADflow.py

View check run for this annotation

Codecov / codecov/patch

adflow/pyADflow.py#L4425 

Added line #L4425 was not covered by tests

blankingSurfData[surfFile] = {"pts": pts, "conn": conn}

Check warning on line 4427 in adflow/pyADflow.py

View check run for this annotation

Codecov / codecov/patch

adflow/pyADflow.py#L4427 

Added line #L4427 was not covered by tests

self.blankingSurfData = blankingSurfData

Check warning on line 4429 in adflow/pyADflow.py

View check run for this annotation

Codecov / codecov/patch

adflow/pyADflow.py#L4429 

Added line #L4429 was not covered by tests

def _oversetExplicitSurfaceCallback(self, flag):
"""This routine runs the explicit surface callback algorithm if user adds
surfaces that define the boundaries of the compute domain. This approach
will work more robustly than the automated flooding algorithm for complex
grids.
Parameters
----------
flag : ndarray
Array that is used as a mask to select the explicitly blanked cells.
This is modified in place.
"""

if self.explicitSurfaceCallback is not None:

# get the number of cells
n = len(flag)

Check warning on line 4447 in adflow/pyADflow.py

View check run for this annotation

Codecov / codecov/patch

adflow/pyADflow.py#L4447 

Added line #L4447 was not covered by tests

# loop over the surfaces
for surf in self.blankingSurfDict:
if self.comm.rank == 0:
Expand All @@ -4417,31 +4464,9 @@ def oversetExplicitSurfaceCallback(self, flag, firstCall=True):
else:
kMin = -1

Check warning on line 4465 in adflow/pyADflow.py

View check run for this annotation

Codecov / codecov/patch

adflow/pyADflow.py#L4465 

Added line #L4465 was not covered by tests

# if this is the first call, we need to load the surface meshes.
# we might have duplicate mesh files, so no need to load them again
# if we have already loadded one copy
if firstCall:
if surfFile not in blankingSurfData:
# optional coordinate transformation to do general manipulation of the coordinates
if "coordXfer" in self.blankingSurfDict[surf]:
coordXfer = self.blankingSurfDict[surf]["coordXfer"]
else:
coordXfer = None

# read the plot3d surface
pts, conn = self._readPlot3DSurfFile(surfFile, convertToTris=False, coordXfer=coordXfer)

blankingSurfData[surfFile] = {"pts": pts, "conn": conn}
else:
# just load from the local dict
pts = blankingSurfData[surfFile]["pts"]
conn = blankingSurfData[surfFile]["conn"]

else:
# we can just load the pts and conn from the dictionary
# these are only saved if we have a DVGeo, a Mesh, and we are doing a full overset update.
pts = self.blankingSurfData[surfFile]["pts"]
conn = self.blankingSurfData[surfFile]["conn"]
# the surf file is loaded in initialization
pts = self.blankingSurfData[surfFile]["pts"]
conn = self.blankingSurfData[surfFile]["conn"]

Check warning on line 4469 in adflow/pyADflow.py

View check run for this annotation

Codecov / codecov/patch

adflow/pyADflow.py#L4468-L4469 

Added lines #L4468 - L4469 were not covered by tests

# get a new flag array
surfFlag = numpy.zeros(n, "intc")

Check warning on line 4472 in adflow/pyADflow.py

View check run for this annotation

Codecov / codecov/patch

adflow/pyADflow.py#L4472 

Added line #L4472 was not covered by tests
Expand All @@ -4453,12 +4478,9 @@ def oversetExplicitSurfaceCallback(self, flag, firstCall=True):
# update the flag array with the new info
flag[:] = numpy.any([flag, surfFlag], axis=0)

Check warning on line 4479 in adflow/pyADflow.py

View check run for this annotation

Codecov / codecov/patch

adflow/pyADflow.py#L4479 

Added line #L4479 was not covered by tests

# finally before we quit, check if we need to save the blanking surface data
if firstCall and (self.getOption("oversetUpdateMode") == "full"):
# save the points and conn if we are doing full overset updates.
# we will add the points to the potential DVGeo and re-use the conn info
# during hole cutting update
self.blankingSurfData = blankingSurfData
# we can delete the surface info if we are not running in full overset update mode
if self.getOption("oversetUpdateMode") != "full":
del self.blankingSurfData

Check warning on line 4483 in adflow/pyADflow.py

View check run for this annotation

Codecov / codecov/patch

adflow/pyADflow.py#L4482-L4483 

Added lines #L4482 - L4483 were not covered by tests

def getAdjointResNorms(self):
"""
Expand Down
3 changes: 1 addition & 2 deletions doc/options.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -259,8 +259,7 @@ updateWallAssociations:
desc: >
Flag to update wall associations, even if the approximate wall distance routines are used.
By default, we don't update the associations because the update itself cannot be differentiated.
However, for large mesh changes, users might still want to update the associations for
more accurate wall distance values.
However, for large mesh changes, users might still want to update the associations for more accurate wall distance values.
eulerWallTreatment:
desc: >
Expand Down
1 change: 1 addition & 0 deletions src/inputParam/inputParamRoutines.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4043,6 +4043,7 @@ subroutine setDefaultValues
approxSA = .False.
useApproxWallDistance = .False.
updateWallAssociations = .False.
recomputeOverlapMatrix = .True.
cflLimit = 3.0
adjointPETScVarsAllocated = .False.
adjointPETScPreProcVarsAllocated = .False.
Expand Down
75 changes: 30 additions & 45 deletions src/overset/oversetAPI.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2128,7 +2128,6 @@ subroutine flagCellsInSurface(pts, npts, conn, nconn, flag, ncell, blockids, nbl
real(kind=realType), dimension(:, :), allocatable :: norm
integer(kind=intType), dimension(:), allocatable :: normCount
integer(kind=intType), dimension(:, :), pointer :: tmp
logical :: keepChecking

! ADT Type required data
integer(kind=intType), dimension(:), pointer :: frontLeaves, frontLeavesNew
Expand Down Expand Up @@ -2257,61 +2256,47 @@ subroutine flagCellsInSurface(pts, npts, conn, nconn, flag, ncell, blockids, nbl
frontLeaves, frontLeavesNew)
cellID = intInfo(3)

keepChecking = .True.

if (cellID > 0) then

! Now check if this was successful or not:
if (checkInside()) then
! Whoohoo! We are inside the region. Flag this cell
! We are inside the region. Flag this cell
flag(cell_counter) = 1
keepChecking = .False.
else
! we are outside. now check if the projection distance is larger than
! the max diagonal. if so, we can quit early here.
if (uvw(4) .gt. diag_max) then
! projection is larger than our biggest diagonal.
! other nodes wont be in the surface, so we can exit the cell early here
keepChecking = .False.
end if
end if
end if
else if (uvw(4) .lt. diag_max) then
! We are outside, but the projection is smaller than the biggest diagonal.
! This means we are outside, but not sufficiently far.
! Check vectors between the center and each vertex

do l = 1, 8
select case (l)
case (1)
vertexPt = x(i - 1, j - 1, k - 1, :)
case (2)
vertexPt = x(i, j - 1, k - 1, :)
case (3)
vertexPt = x(i, j, k - 1, :)
case (4)
vertexPt = x(i - 1, j, k - 1, :)
case (5)
vertexPt = x(i - 1, j - 1, k, :)
case (6)
vertexPt = x(i, j - 1, k, :)
case (7)
vertexPt = x(i, j, k, :)
case (8)
vertexPt = x(i - 1, j, k, :)
end select

if (keepChecking) then
! we are outside, but not sufficiently far.
! Check vectors between the center and each vertex

! actually test each node
do l = 1, 8
select case (l)
case (1)
vertexPt = x(i - 1, j - 1, k - 1, :)
case (2)
vertexPt = x(i, j - 1, k - 1, :)
case (3)
vertexPt = x(i, j, k - 1, :)
case (4)
vertexPt = x(i - 1, j, k - 1, :)
case (5)
vertexPt = x(i - 1, j - 1, k, :)
case (6)
vertexPt = x(i, j - 1, k, :)
case (7)
vertexPt = x(i, j, k, :)
case (8)
vertexPt = x(i - 1, j, k, :)
end select

if (cellID > 0) then
! Now check if this was successful or not:
if (checkNearby()) then
! Whoohoo! We are inside the region. Flag this cell
! We are inside the region. Flag this cell.
flag(cell_counter) = 1
exit
end if
end if
end do
end if
end do

end if
end if
end if

cell_counter = cell_counter + 1
Expand Down
Loading

0 comments on commit f8e1c0b

Please sign in to comment.