Skip to content

Commit

Permalink
Merge branch 'main' into multi_bc_fix
Browse files Browse the repository at this point in the history
  • Loading branch information
@eirikurj
eirikurj committed Jun 28, 2024
2 parents efbda44 + 6706920 commit a973beb
Show file tree
Hide file tree
Showing 10 changed files with 343 additions and 125 deletions.
252 changes: 194 additions & 58 deletions adflow/pyADflow.py
View file

Large diffs are not rendered by default.

16 changes: 15 additions & 1 deletion doc/options.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,7 @@ isosurface:
desc: >
Dictionary specifying the type and values to be used for isosurfaces.
Any of the ``volumeVariables`` may be used.
An example dictionary is : ``{"Vx":-0.001, "shock":1.0}``.
An example dictionary is : ``{"vx":-0.001, "shock":1.0}``.
This will place an isosurface at (essentially) 0 x-velocity and an isosurface at the shock sensor value of 1 (used to visualize the shock region).
isoVariables:
Expand Down Expand Up @@ -255,6 +255,12 @@ useApproxWallDistance:
After the geometry deforms (such as during an optimization) the spatial search algorithm is not run, but the distance between the (new) parametric location and the (new) grid cell center is computed and taken as the wall distance.
This is substantially faster and permits efficient wall-distance updates for use in aerostructural analysis.
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.
eulerWallTreatment:
desc: >
Specifies how wall boundary conditions are implemented for inviscid simulations.
Expand Down Expand Up @@ -712,6 +718,14 @@ oversetPriority:
A lower factor will encourage the usage of that block mesh.
This option may be required to get the flooding algorithm working properly.
recomputeOverlapMatrix:
desc: >
Flag to determine if the domain overlap matrix is re-computed when a full overset update is done.
This overlap matrix determines the connections between domains for the donor search algorithm; for the code to search a domain for potential donors, that domain must overlap with the current domain.
This matrix is re-computed when a full overset update is done because the domain overlaps can change.
If the overlap matrix is outdated, the code will fail to find donors in cases where there are many potential donors available from overlapping domains that are not represented in the matrix.
If the user knows that the overlap matrix will not change, and they want to improve the efficiency of the full overset update, they can set this option to False.
oversetDebugPrint:
desc: >
Flag to enable or disable the debug printout from the overset algorithm when a hole cutting process fails.
Expand Down
2 changes: 2 additions & 0 deletions src/f2py/adflow.pyf
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1009,6 +1009,7 @@ python module libadflow
logical :: lumpeddiss
real(kind=realtype) :: sigma
logical :: useapproxwalldistance
logical :: updatewallassociations
logical :: lowspeedpreconditioner
logical :: useblockettes
end module inputdiscretization
Expand Down Expand Up @@ -1261,6 +1262,7 @@ python module libadflow
logical :: debugzipper
logical :: usezippermesh
logical :: useoversetwallscaling
logical :: recomputeoverlapmatrix
logical :: oversetdebugprint
real(kind=realtype) :: selfzipcutoff
end module inputoverset
Expand Down
2 changes: 2 additions & 0 deletions src/inputParam/inputParamRoutines.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4042,6 +4042,8 @@ subroutine setDefaultValues
lumpedDiss = .False.
approxSA = .False.
useApproxWallDistance = .False.
updateWallAssociations = .False.
recomputeOverlapMatrix = .True.
cflLimit = 3.0
adjointPETScVarsAllocated = .False.
adjointPETScPreProcVarsAllocated = .False.
Expand Down
4 changes: 3 additions & 1 deletion src/modules/inputParam.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,7 @@ module inputDiscretization
! computations. Typically only used for
! repeated calls when the wall distance would
! not have changed significantly
! updateWallAssociation : Logical to determine if the full wall distance
! updateWallAssociations: Logical to determine if the full wall distance
! assocation is to be performed on the next
! wall distance calculation. This is only
! significant when useApproxWallDistance is
Expand Down Expand Up @@ -92,6 +92,7 @@ module inputDiscretization
logical :: radiiNeededFine, radiiNeededCoarse

logical :: useApproxWallDistance
logical :: updateWallAssociations
logical :: lowSpeedPreconditioner
end module inputDiscretization

Expand Down Expand Up @@ -880,5 +881,6 @@ module inputOverset
integer(kind=intType) :: nFloodIter
logical :: useZipperMesh
logical :: useOversetWallScaling
logical :: recomputeOverlapMatrix
logical :: oversetDebugPrint
end module inputOverset
2 changes: 1 addition & 1 deletion src/modules/wallDistanceData.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ module wallDistanceData
real(kind=realType), dimension(:), pointer :: xSurf

logical, dimension(:), allocatable :: wallDistanceDataAllocated
logical, dimension(:), allocatable :: updateWallAssociation
logical, dimension(:), allocatable :: updateLevelWallAssociation

#ifndef USE_TAPENADE
real(kind=realType), dimension(:), pointer :: xSurfd
Expand Down
149 changes: 103 additions & 46 deletions src/overset/oversetAPI.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2105,7 +2105,7 @@ subroutine flagCellsInSurface(pts, npts, conn, nconn, flag, ncell, blockids, nbl
use ADTData
use blockPointers, only: x, il, jl, kl, nDom, iBlank, vol, nbkGlobal, kBegOr
use adjointVars, only: nCellsLocal
use utils, only: setPointers, EChk
use utils, only: setPointers, EChk, mynorm2
use sorting, only: famInList
implicit none

Expand All @@ -2122,7 +2122,7 @@ subroutine flagCellsInSurface(pts, npts, conn, nconn, flag, ncell, blockids, nbl
integer(kind=intType) :: iblock, cell_counter, k_cgns
real(kind=realType) :: dStar, frac, volLocal
real(kind=realType), dimension(3) :: minX, maxX, v1, v2, v3, axisVec
real(kind=realType), dimension(3) :: diag1, diag2, diag3, diag4
real(kind=realType), dimension(3) :: diag1, diag2, diag3, diag4, vertexPt
real(kind=realType) :: dd1, dd2, dd3, dd4, diag_max
type(adtType) :: ADT
real(kind=realType), dimension(:, :), allocatable :: norm
Expand Down Expand Up @@ -2237,51 +2237,66 @@ subroutine flagCellsInSurface(pts, npts, conn, nconn, flag, ncell, blockids, nbl
! we can work with squared distances for the sake of efficiency. the projection routine
! will also return the squared distance for the same reason.

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

! reset the "closest point already found" variable.
coor(4) = dStar
intInfo(3) = 0
call minDistancetreeSearchSinglePoint(ADT, coor, intInfo, &
uvw, dummy, 0, BB, &
frontLeaves, frontLeavesNew)
cellID = intInfo(3)
if (cellID > 0) then
! Now check if this was successful or not:
if (checkInside()) then
! Whoohoo! We are inside the region. Flag this cell
flag(cell_counter) = 1
exit
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
! First test the cell center
coor(1:3) = eighth * ( &
x(i - 1, j - 1, k - 1, :) + &
x(i, j - 1, k - 1, :) + &
x(i - 1, j, k - 1, :) + &
x(i, j, k - 1, :) + &
x(i - 1, j - 1, k, :) + &
x(i, j - 1, k, :) + &
x(i - 1, j, k, :) + &
x(i, j, k, :))

! reset the "closest point already found" variable.
coor(4) = dStar
intInfo(3) = 0
call minDistancetreeSearchSinglePoint(ADT, coor, intInfo, &
uvw, dummy, 0, BB, &
frontLeaves, frontLeavesNew)
cellID = intInfo(3)

if (cellID > 0) then

! Now check if this was successful or not:
if (checkInside()) then
! We are inside the region. Flag this cell
flag(cell_counter) = 1
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

! Now check if this was successful or not:
if (checkNearby()) then
! 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

cell_counter = cell_counter + 1
Expand Down Expand Up @@ -2332,6 +2347,48 @@ function checkInside()
end if
end function checkInside

function checkNearby()

implicit none
logical :: checkNearby
integer(kind=intType) :: jj
real(kind=realType) :: shp(4), xp(3), normal(3), v1(3), projLen
real(kind=realType) :: vertexVec(3), projVec(3), projDist

! bi-linear shape functions (CCW ordering)
shp(1) = (one - uvw(1)) * (one - uvw(2))
shp(2) = (uvw(1)) * (one - uvw(2))
shp(3) = (uvw(1)) * (uvw(2))
shp(4) = (one - uvw(1)) * (uvw(2))

xp = zero
do jj = 1, 4
xp = xp + shp(jj) * pts(:, conn(jj, cellID))
end do

! this is the vector that points from the cell center to the vertex
vertexVec = vertexPt - coor(1:3)

! this is the vector from the cell center to the projection
projVec = xp - coor(1:3)

! get the projection distance
projDist = mynorm2(projVec)
! normalize it
projVec = projVec / projDist

! compute the projected length of the vertexVec in the direction of projVec
projLen = vertexVec(1) * projVec(1) + vertexVec(2) * projVec(2) + vertexVec(3) * projVec(3)

! we flag this cell if this dot product is greater than the projection length
if (projLen .gt. projDist) then
checkNearby = .True.
else
checkNearby = .False.
end if

end function

end subroutine flagCellsInSurface

subroutine updateOverset(flag, n, closedFamList, nFam)
Expand All @@ -2352,7 +2409,7 @@ subroutine updateOverset(flag, n, closedFamList, nFam)
! oversetComm routine.

use constants
use inputOverset, only: oversetUpdateMode
use inputOverset, only: oversetUpdateMode, recomputeOverlapMatrix
use block, only: flowDOms
use inputTimeSpectral, only: nTimeIntervalsSpectral
use oversetCommUtilities, onlY: updateOversetConnectivity
Expand Down Expand Up @@ -2408,7 +2465,7 @@ subroutine updateOverset(flag, n, closedFamList, nFam)
do level = 1, nLevels
if (level == 1) then
call setExplicitHoleCut(flag)
call oversetComm(level, .False., .false., closedFamList)
call oversetComm(level, recomputeOverlapMatrix, .false., closedFamList)
else
call oversetComm(level, .False., .True., closedFamList)
end if
Expand Down
23 changes: 12 additions & 11 deletions src/overset/oversetUtilities.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1541,7 +1541,7 @@ subroutine qsortPocketEdgeType(arr, nn)

end subroutine qsortPocketEdgeType

subroutine checkOverset(level, sps, totalOrphans, printBadCells)
subroutine checkOverset(level, sps, totalOrphans, lastCall)

!
! CheckOverset checks the integrity of the overset connectivity
Expand All @@ -1560,7 +1560,7 @@ subroutine checkOverset(level, sps, totalOrphans, printBadCells)
! Input/Output
integer(kind=intType), intent(in) :: level, sps
integer(kind=intType), intent(out) :: totalOrphans
logical, intent(in) :: printBadCells
logical, intent(in) :: lastCall

! Working
integer(kind=intType) :: i, j, k, nn, ii, jj, kk, n, ierr
Expand Down Expand Up @@ -1588,13 +1588,8 @@ subroutine checkOverset(level, sps, totalOrphans, printBadCells)
badCell = .True.
end if
end do stencilLoop
if (badCell .and. printBadCells) then
if (oversetDebugPrint) &
print *, 'Error in connectivity at :', nbkglobal, i + iBegOr, j + jBegOr, k + kBegOr
! we can modify iBlankLast because this is the last checkOverset call.
! we set iBlankLast to -5 to mark orphan cells, this value will then
! be moved to iBlank after we are done with other loops.
flowDoms(nn, level, sps)%iBlankLast(i, j, k) = -5
if (badCell .and. lastCall .and. oversetDebugPrint) then
print *, 'Error in connectivity at :', nbkglobal, i + iBegOr, j + jBegOr, k + kBegOr
end if
end if
end do
Expand Down Expand Up @@ -1670,6 +1665,12 @@ subroutine checkOverset(level, sps, totalOrphans, printBadCells)
! This cell is an orphan:
n = n + 1
orphans(:, n) = (/ii, jj, kk/)
if (lastCall) then
! we can modify iBlankLast because this is the last checkOverset call.
! we set iBlankLast to -5 to mark orphan cells, this value will then
! be moved to iBlank after we are done with other loops.
flowDoms(nn, level, sps)%iBlankLast(ii, jj, kk) = -5
end if
end if
end if
end do stencilLoop3
Expand All @@ -1688,10 +1689,10 @@ subroutine checkOverset(level, sps, totalOrphans, printBadCells)
print *, 'Total number of orphans:', totalOrphans
end if

! if this is the last checkOverset call with printBadCells = .True., then
! if this is the last checkOverset call with lastCall = .True., then
! we can move the orphan information to iBlank because we have done all
! the checking required.
if (printBadCells .and. (totalOrphans .gt. 0)) then
if (lastCall .and. (totalOrphans .gt. 0)) then
! loop over the cells and write the orphan information to iBlank
do nn = 1, nDom
call setPointers(nn, level, sps)
Expand Down
6 changes: 3 additions & 3 deletions src/preprocessing/preprocessingAPI.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ subroutine preprocessing
use inputTimeSpectral
use section
use wallDistance, only: xVolumeVec, xSurfVec, wallScatter, &
wallDistanceDataAllocated, updateWallAssociation, &
wallDistanceDataAllocated, updateLevelWallAssociation, &
computeWallDistance
use oversetData, only: cumDomProc, nDomProc, wallFringes, nDomTotal, &
overlapMatrix, oversetPresent, localWallFringes
Expand Down Expand Up @@ -159,9 +159,9 @@ subroutine preprocessing

! Allocate some data of size nLevels for the fast wall distance calc
allocate (xVolumeVec(nLevels), xSurfVec(nLevels, mm), wallScatter(nLevels, mm), &
wallDistanceDataAllocated(nLevels), updateWallAssociation(nLevels))
wallDistanceDataAllocated(nLevels), updateLevelWallAssociation(nLevels))
wallDistanceDataAllocated = .False.
updateWallAssociation = .True.
updateLevelWallAssociation = .True.

! Nullify the wallFringe poiter as initialization
nullify (wallFringes, localWallFringes)
Expand Down
Loading

0 comments on commit a973beb

Please sign in to comment.