Updates to several overset update routines

doc/options.yaml

@@ -255,6 +255,13 @@ 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.
@@ -712,6 +719,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.

src/f2py/adflow.pyf

@@ -1006,6 +1006,7 @@ python module libadflow
          logical :: lumpeddiss
          real(kind=realtype) :: sigma
          logical :: useapproxwalldistance
+         logical :: updatewallassociations
          logical :: lowspeedpreconditioner
          logical :: useblockettes
        end module inputdiscretization
@@ -1258,6 +1259,7 @@ python module libadflow
          logical :: debugzipper
          logical :: usezippermesh
          logical :: useoversetwallscaling
+         logical :: recomputeoverlapmatrix
          logical :: oversetdebugprint
          real(kind=realtype) :: selfzipcutoff
        end module inputoverset

src/inputParam/inputParamRoutines.F90

@@ -4042,6 +4042,7 @@ subroutine setDefaultValues
         lumpedDiss = .False.
         approxSA = .False.
         useApproxWallDistance = .False.
+        updateWallAssociations = .False.
         cflLimit = 3.0
         adjointPETScVarsAllocated = .False.
         adjointPETScPreProcVarsAllocated = .False.

src/modules/inputParam.F90

@@ -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
@@ -92,6 +92,7 @@ module inputDiscretization
     logical :: radiiNeededFine, radiiNeededCoarse
 
     logical :: useApproxWallDistance
+    logical :: updateWallAssociations
     logical :: lowSpeedPreconditioner
 end module inputDiscretization
 
@@ -880,5 +881,6 @@ module inputOverset
     integer(kind=intType) :: nFloodIter
     logical :: useZipperMesh
     logical :: useOversetWallScaling
+    logical :: recomputeOverlapMatrix
     logical :: oversetDebugPrint
 end module inputOverset

src/modules/wallDistanceData.F90

@@ -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

src/overset/oversetAPI.F90

@@ -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
 
@@ -2122,12 +2122,13 @@ 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
         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
@@ -2237,51 +2238,82 @@ 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)
+
+                                ! increment the cell counter before the checks in case we exit early
+
+                                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
+                                        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
+
+                                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
+                                                flag(cell_counter) = 1
                                                 exit
                                             end if
                                         end if
-                                    end if
-                                end do
+                                    end do
+                                end if
+
                             end if
 
                             cell_counter = cell_counter + 1
@@ -2332,6 +2364,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)
@@ -2352,7 +2426,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
@@ -2408,7 +2482,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

src/overset/oversetUtilities.F90

@@ -1591,10 +1591,6 @@ subroutine checkOverset(level, sps, totalOrphans, printBadCells)
                             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
                             end if
                         end if
                     end do
@@ -1670,6 +1666,12 @@ subroutine checkOverset(level, sps, totalOrphans, printBadCells)
                                         ! This cell is an orphan:
                                         n = n + 1
                                         orphans(:, n) = (/ii, jj, kk/)
+                                        if (printBadCells) 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

src/preprocessing/preprocessingAPI.F90

@@ -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
@@ -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)