diff --git a/adflow/pyADflow.py b/adflow/pyADflow.py index 790a2ed2c..2b0f9ff06 100644 --- a/adflow/pyADflow.py +++ b/adflow/pyADflow.py @@ -5521,6 +5521,8 @@ def _getDefaultOptions(): ], "turbulenceOrder": [str, ["first order", "second order"]], "turbResScale": [(float, list, type(None)), None], + "meshMaxSkewness": [float, 1.0], + "useSkewnessCheck": [bool, False], "turbulenceProduction": [str, ["strain", "vorticity", "Kato-Launder"]], "useQCR": [bool, False], "useRotationSA": [bool, False], @@ -5696,6 +5698,7 @@ def _getDefaultOptions(): "setMonitor": [bool, True], "printWarnings": [bool, True], "printNegativeVolumes": [bool, False], + "printBadlySkewedCells": [bool, False], "monitorVariables": [list, ["cpu", "resrho", "resturb", "cl", "cd"]], "surfaceVariables": [list, ["cp", "vx", "vy", "vz", "mach"]], "volumeVariables": [list, ["resrho"]], @@ -5783,6 +5786,7 @@ def _getImmutableOptions(self): "zippersurfacefamily", "cutcallback", "explicitsurfacecallback", + "useskewnesscheck", ) def _getOptionMap(self): @@ -5894,6 +5898,8 @@ def _getOptionMap(self): }, "turbulenceorder": {"first order": 1, "second order": 2, "location": ["discr", "orderturb"]}, "turbresscale": ["iter", "turbresscale"], + "meshmaxskewness": ["iter", "meshmaxskewness"], + "useskewnesscheck": ["iter", "useskewnesscheck"], "turbulenceproduction": { "strain": self.adflow.constants.strain, "vorticity": self.adflow.constants.vorticity, @@ -6095,6 +6101,7 @@ def _getOptionMap(self): "printiterations": ["iter", "printiterations"], "printwarnings": ["iter", "printwarnings"], "printnegativevolumes": ["iter", "printnegativevolumes"], + "printbadlyskewedcells": ["iter", "printbadlyskewedcells"], "printtiming": ["adjoint", "printtiming"], "setmonitor": ["adjoint", "setmonitor"], "storeconvhist": ["io", "storeconvinneriter"], diff --git a/doc/options.yaml b/doc/options.yaml index 35bef4e48..94d49a9e5 100644 --- a/doc/options.yaml +++ b/doc/options.yaml @@ -1248,6 +1248,10 @@ printNegativeVolumes: desc: > Flag to print the block indices, cell center coordinates, and volume for each negative volume cell in the mesh. +printBadlySkewedCells: + desc: > + Flag to print the block indices, cell center coordinates, and skewness for each cell whose skewness is above the value defined in ``meshMaxSkewness``. Only used when ``useSkewnessCheck`` is active. + monitorVariables: desc: > List of the variables whose convergence should be monitored. @@ -1545,4 +1549,13 @@ cavSensorSharpness: cavExponent: desc: > - The exponent for the numerator term (- Cp - cavitationnumber) of the cavitation sensor. \ No newline at end of file + The exponent for the numerator term (- Cp - cavitationnumber) of the cavitation sensor. + +meshMaxSkewness: + desc: > + Adflow throws an error and fails if the `skewness` of the mesh is above this value. `Skewness` is defined as described `here `__. Only used when ``useSkewnessCheck`` is active. + +useSkewnessCheck: + desc: > + When set to true, ADflow computes the `skewness` of each cell and throws an error if it is above ``meshMaxSkewness``. See also ``printBadlySkewedCells``. + diff --git a/src/f2py/adflow.pyf b/src/f2py/adflow.pyf index 1c569a5bd..5cc9678d0 100644 --- a/src/f2py/adflow.pyf +++ b/src/f2py/adflow.pyf @@ -1082,10 +1082,13 @@ python module libadflow real(kind=realtype) allocatable,dimension(:) :: etark real(kind=realtype) allocatable,dimension(:) :: cdisrk real(kind=realtype), dimension(4) :: turbresscale + real(kind=realtype) :: meshmaxskewness + logical :: useskewnesscheck logical :: freezeturbsource logical :: printiterations logical :: printwarnings logical :: printnegativevolumes + logical :: printbadlyskewedcells real(kind=realtype) ::maxl2deviationfactor logical :: uselinresmonitor logical :: usedisscontinuation @@ -1501,4 +1504,4 @@ python module libadflow end python module libadflow_cs #else end python module libadflow -#endif \ No newline at end of file +#endif diff --git a/src/modules/block.F90 b/src/modules/block.F90 index ecf098902..1fdb9b3fe 100644 --- a/src/modules/block.F90 +++ b/src/modules/block.F90 @@ -382,6 +382,7 @@ module block ! needed for unsteady problems on ! deforming grids. Only allocated on ! the finest grid level. + ! skew(0:ib,0:jb,0:kb) - Skewness of Volume ! uv(2,2:il,2:jl,2:kl) - Parametric location on elemID for each cell. ! Only used for fast wall distance calcs. ! porI(1:il,2:jl,2:kl) - Porosity in the i direction. @@ -433,6 +434,7 @@ module block real(kind=realType), dimension(:, :, :), pointer :: vol real(kind=realType), dimension(:, :, :, :), pointer :: volOld real(kind=realType), dimension(:, :, :), pointer :: volref + real(kind=realType), dimension(:, :, :), pointer :: skew real(kind=realType), dimension(:, :, :, :), pointer :: uv integer(kind=intType), dimension(:, :, :, :), pointer :: surfNodeIndices diff --git a/src/modules/blockPointers.F90 b/src/modules/blockPointers.F90 index e2e5e9682..4705f9fd7 100644 --- a/src/modules/blockPointers.F90 +++ b/src/modules/blockPointers.F90 @@ -97,6 +97,7 @@ module blockPointers real(kind=realType), dimension(:, :, :), pointer :: vol real(kind=realType), dimension(:, :, :), pointer :: volref real(kind=realType), dimension(:, :, :, :), pointer :: volOld + real(kind=realType), dimension(:, :, :), pointer :: skew real(kind=realType), dimension(:, :, :, :), pointer :: dadidata integer(kind=porType), dimension(:, :, :), pointer :: porI, porJ, porK diff --git a/src/modules/constants.F90 b/src/modules/constants.F90 index 4dc4a23c9..4653d56ec 100644 --- a/src/modules/constants.F90 +++ b/src/modules/constants.F90 @@ -511,4 +511,8 @@ module constants integer(kind=intType), parameter :: iZippWallY = 8 integer(kind=intType), parameter :: iZippWallZ = 9 + ! Metric types + integer(kind=intType), parameter :: MetricVolume = 1 + integer(kind=intType), parameter :: MetricSkewness = 2 + end module constants diff --git a/src/modules/inputParam.F90 b/src/modules/inputParam.F90 index 5f7e41444..be8a4a729 100644 --- a/src/modules/inputParam.F90 +++ b/src/modules/inputParam.F90 @@ -247,9 +247,11 @@ module inputIteration ! cdisRk: Dissipative coefficients in the runge kutta ! scheme. The values depend on the number of ! stages specified. - ! printIterations: If True, iterations are printed to stdout - ! turbresscale: Scaling factor for turbulent residual. Necessary for - ! NKsolver with RANS. Only tested on SA. + ! printIterations: If True, iterations are printed to stdout + ! turbresscale: Scaling factor for turbulent residual. Necessary for + ! NKsolver with RANS. Only tested on SA. + ! meshMaxSkewness If one cell has a highe skewness than this, the Solver + ! errors out. ! iterType : String used for specifying which type of iteration was taken ! ! Definition of the string, which stores the multigrid cycling @@ -285,7 +287,10 @@ module inputIteration logical :: printIterations logical :: printWarnings logical :: printNegativeVolumes + logical :: printBadlySkewedCells real(kind=realType), dimension(4) :: turbResScale + real(kind=realType) :: meshMaxSkewness + logical :: useSkewnessCheck logical :: useDissContinuation real(kind=realType) :: dissContMagnitude, dissContMidpoint, dissContSharpness diff --git a/src/preprocessing/preprocessingAPI.F90 b/src/preprocessing/preprocessingAPI.F90 index e1ad7289d..8275ed3c2 100644 --- a/src/preprocessing/preprocessingAPI.F90 +++ b/src/preprocessing/preprocessingAPI.F90 @@ -2593,6 +2593,7 @@ subroutine allocateMetric(level) use inputTimeSpectral use iteration use inputUnsteady + use inputIteration, only: useSkewnessCheck use utils, only: terminate, setPointers implicit none ! @@ -2640,6 +2641,16 @@ subroutine allocateMetric(level) "Memory allocation failure for & &normals and volumes") + ! Allocate the memory for the skewness if desired + if (useSkewnessCheck) then + allocate (flowDoms(nn, level, sps)%skew(0:ib, 0:jb, 0:kb), & + stat=ierr) + if (ierr /= 0) & + call terminate("allocateMetric", & + "Memory allocation failure for & + &skewness") + end if + ! Added by HDN ! Added s[I,J,K]ALE if (equationMode == unSteady .and. useALE) then @@ -2721,7 +2732,8 @@ subroutine metric(level) use communication use inputTimeSpectral use checkVolBlock - use inputIteration, only: printWarnings, printNegativeVolumes + use inputIteration, only: printWarnings, printNegativeVolumes, & + useSkewnessCheck, meshMaxSkewness, printBadlySkewedCells use utils, only: setPointers, terminate, returnFail use commonFormats, only: stringSpace, stringInt1 implicit none @@ -2743,11 +2755,14 @@ subroutine metric(level) integer(kind=intType) :: nn, mm, sps integer(kind=intType) :: nVolNeg, nVolPos integer(kind=intType) :: nVolBad, nVolBadGlobal + integer(kind=intType) :: nBadSkew, nBadSkewGlobal integer(kind=intType) :: nBlockBad, nBlockBadGlobal real(kind=realType) :: fact, mult real(kind=realType) :: xp, yp, zp, vp1, vp2, vp3, vp4, vp5, vp6 + real(kind=realType), dimension(24) :: angles + real(kind=realType), dimension(3) :: v1, v2 real(kind=realType), dimension(:, :, :), pointer :: ss @@ -2758,6 +2773,7 @@ subroutine metric(level) logical :: badVolume, iBlankAllocated logical, dimension(:, :, :), pointer :: volumeIsNeg + logical, dimension(:, :, :), pointer :: volumeIsSkewed type(checkVolBlockType), & dimension(nDom, nTimeIntervalsSpectral) :: checkVolDoms @@ -2766,6 +2782,7 @@ subroutine metric(level) nVolBad = 0 nBlockBad = 0 + nBadSkew = 0 ! Loop over the number of spectral solutions and local blocks. @@ -2789,6 +2806,16 @@ subroutine metric(level) call terminate("metric", & "Memory allocation failure for volumeIsNeg") volumeIsNeg => checkVolDoms(nn, sps)%volumeIsNeg + + if (useSkewnessCheck) then + allocate (checkVolDoms(nn, sps)%volumeIsSkewed(2:il, 2:jl, 2:kl), & + stat=ierr) + if (ierr /= 0) & + call terminate("metric", & + "Memory allocation failure for volumeIsSkewed") + volumeIsSkewed => checkVolDoms(nn, sps)%volumeIsSkewed + end if + ! ! Volume and block orientation computation. ! @@ -2898,8 +2925,33 @@ subroutine metric(level) vol(i, j, k) = sixth * (vp1 + vp2 + vp3 + vp4 + vp5 + vp6) + ! Calculate the min and max angles of each face + ! This is necessairy for the skewness calculation + ! If we are not interested in skewness, skip it + + if (useSkewnessCheck) then + call minMaxAngle(angles(1:4), & + x(i, m, k, :), x(i, j, k, :), x(l, j, k, :), x(l, m, k, :)) + + call minMaxAngle(angles(5:8), & + x(i, m, n, :), x(i, j, n, :), x(l, j, n, :), x(l, m, n, :)) + + call minMaxAngle(angles(9:12), & + x(i, m, k, :), x(i, m, n, :), x(i, j, n, :), x(i, j, k, :)) + + call minMaxAngle(angles(13:16), & + x(l, m, k, :), x(l, m, n, :), x(l, j, n, :), x(l, j, k, :)) + + call minMaxAngle(angles(17:20), & + x(i, j, k, :), x(i, j, n, :), x(l, j, n, :), x(l, j, k, :)) + + call minMaxAngle(angles(21:24), & + x(i, m, k, :), x(i, m, n, :), x(l, m, n, :), x(l, m, k, :)) + end if + ! Check the volume and update the number of positive ! and negative volumes if needed. + ! Also calculate the skewness of the cell if (checkAll) then @@ -2940,6 +2992,24 @@ subroutine metric(level) if (badVolume) nVolBad = nVolBad + 1 + ! Calculate the skewness and update nBadSkew if it + ! is above the threshold + ! Skip it if we don't care about skewness + + if (useSkewnessCheck) then + skew(i, j, k) = max((maxval(angles) - (pi / two)) / & + (pi / two), & + ((pi / two) - minval(angles)) / & + (pi / two)) + + if (skew(i, j, k) > meshMaxSkewness) then + volumeIsSkewed(i, j, k) = .true. + nBadSkew = nBadSkew + 1 + else + volumeIsSkewed(i, j, k) = .false. + end if + end if + end if ! Set the volume to the absolute value. @@ -3025,6 +3095,18 @@ subroutine metric(level) else checkVolDoms(nn, sps)%blockHasNegVol = .false. end if + + ! Check if the block has badly skewed volumes. Only do this if we are + ! actually interested + + if (useSkewnessCheck) then + if (nBadSkew > 0) then + checkVolDoms(nn, sps)%blockHasSkewedVol = .true. + else + checkVolDoms(nn, sps)%blockHasSkewedVol = .false. + end if + end if + ! ! Computation of the face normals in i-, j- and k-direction. ! Formula's are valid for a right handed block; for a left @@ -3269,7 +3351,7 @@ subroutine metric(level) ! list of the bad volumes and terminate executation. if (printNegativeVolumes) then - call writeNegVolumes(checkVolDoms) + call writeBadVolumes(checkVolDoms, MetricVolume) end if call returnFail("metric", & @@ -3292,6 +3374,47 @@ subroutine metric(level) end if end if + ! Determine the global number of badly skewed blocks. The result must be + ! known on all processors and thus an allreduce is needed. + + if (useSkewnessCheck) then + call mpi_allreduce(nBadSkew, nBadSkewGlobal, 1, adflow_integer, & + mpi_sum, ADflow_comm_world, ierr) + + ! Test if badly skewed cells are present in the grid. If so, the action + ! taken depends on the grid level. + + if (nBadSkewGlobal > 0) then + if (level == 1) then + + ! Badly skewed volumes present on the fine grid level. Print a + ! list of the bad volumes and terminate executation. + + if (printBadlySkewedCells) then + call writeBadVolumes(checkVolDoms, MetricSkewness) + end if + + call returnFail("metric", & + "Badly skewed cells present in grid.") + call mpi_barrier(ADflow_comm_world, ierr) + + else + + ! Coarser grid level. The fine grid is okay, but due to the + ! coarsening badly skewed volumes are introduced. Print a warning. + + if (myID == 0) then + print "(a)", "#" + print "(a)", "# Warning" + print stringInt1, "#* Badly skewed cells present on coarse grid level ", level, "." + print "(a)", "#* Computation continues, but be aware of this" + print "(a)", "#" + end if + + end if + end if + end if + ! Determine the global number of bad volumes. The result will ! only be known on processor 0. The quality volume check will ! only be done for the finest grid level. @@ -3326,6 +3449,18 @@ subroutine metric(level) end do end do + ! Release memory for skewness again + if (useSkewnessCheck) then + do sps = 1, nTimeIntervalsSpectral + do nn = 1, nDom + deallocate (checkVolDoms(nn, sps)%volumeIsSkewed, stat=ierr) + if (ierr /= 0) & + call terminate("metric", & + "Deallocation failure for volumeIsSkewed") + end do + end do + end if + contains ! ================================================================ @@ -3363,13 +3498,53 @@ end function volpym end subroutine metric + subroutine minMaxAngle(angles, p1, p2, p3, p4) + ! + ! Calculates the min and max angles inside a quadrilateral + ! + use constants + use utils, only: norm2cplx + implicit none + ! + ! Subroutine arguments. + ! + real(kind=realType), dimension(3), intent(in) :: p1, p2, p3, p4 + real(kind=realType), dimension(4), intent(out) :: angles + ! + ! Local variables. + ! + integer(kind=intType) :: i + real(kind=realType), dimension(4, 3) :: v + + ! initialize variables + v(1, :) = p2 - p1 + v(2, :) = p3 - p2 + v(3, :) = p4 - p3 + v(4, :) = p1 - p4 + + ! loop through vectors and calculate the first 3 angles + do i = 1, 3 + angles(i) = acos(dot_product(-v(i, :), v(i + 1, :)) / & + (norm2cplx(-v(i, :)) * norm2cplx(v(i + 1, :)))) + end do + + ! since it is a quadrilateral, the last angle must be: + ! 360° - sum(first 3 angles) + angles(4) = two * pi - sum(angles(:3)) + + end subroutine minMaxAngle + ! ================================================================== - subroutine writeNegVolumes(checkVolDoms) + subroutine writeBadVolumes(checkVolDoms, mode) ! - ! writeNegVolumes writes the negative volumes of a block to - ! stdout. If a block is flagged to have negative volumes it is - ! assumed that the block is intended to be a right handed block. + ! Writes the bad volumes of a block to stdout. + ! Depending on the mode, it writes those metrics: + ! 1: negative volumes - If a block is flagged + ! to have negative volumes it is assumed + ! that the block is intended to be a + ! right handed block. + ! 2: skewed volumes ! use constants use blockPointers @@ -3384,6 +3559,7 @@ subroutine writeNegVolumes(checkVolDoms) ! ! Subroutine arguments. ! + integer, intent(in) :: mode type(checkVolBlockType), & dimension(nDom, nTimeIntervalsSpectral), intent(in) :: checkVolDoms ! @@ -3392,7 +3568,22 @@ subroutine writeNegVolumes(checkVolDoms) integer :: proc, ierr integer(kind=intType) :: nn, sps, i, j, k real(kind=realType), dimension(3) :: xc + real(kind=realType) :: quantity + character(len=10) :: intString1, intString2, intString3 + character(len=20) :: modeString, descString + + logical :: blockIsBad, volumeIsBad + + ! prepare the string to use according to the mode it is operating in + select case (mode) + case (MetricVolume) + modeString = "Negative volumes" + descString = "Volume" + case (MetricSkewness) + modeString = "Badly skewed cells" + descString = "Skewness" + end select ! Processor 0 prints a message that negative volumes are present ! in the grid. @@ -3400,8 +3591,8 @@ subroutine writeNegVolumes(checkVolDoms) if (myID == 0) then print "(a)", "#" print "(a)", "# Error" - print "(a)", "# Negative volumes found in the grid." - print "(a)", "# A list of the negative volumes is printed below" + print "(3(a))", "# ", trim(modeString), " found in the grid." + print "(a)", "# A list of those volumes is printed below" print "(a)", "#" end if @@ -3421,8 +3612,17 @@ subroutine writeNegVolumes(checkVolDoms) domains: do nn = 1, nDom ! Test for a bad block. + blockIsBad = .false. + select case (mode) + case (MetricVolume) + if (checkVolDoms(nn, sps)%blockHasNegVol) & + blockIsBad = .true. + case (MetricSkewness) + if (checkVolDoms(nn, sps)%blockHasSkewedVol) & + blockIsBad = .true. + end select - testBad: if (checkVolDoms(nn, sps)%blockHasNegVol) then + testBad: if (blockIsBad) then ! Set the pointers for this block. @@ -3438,8 +3638,9 @@ subroutine writeNegVolumes(checkVolDoms) case (steady, unsteady) print "(a)", "#" + print stringSpace, "# Block", trim(cgnsDoms(nbkGlobal)%zoneName), & - "contains the following negative volumes" + "contains the following ", trim(modeString) print "(a)", "#--------------------------------------------------------------------" print "(a)", "#" @@ -3452,7 +3653,7 @@ subroutine writeNegVolumes(checkVolDoms) print "(a)", "#" print stringSpace, "# Spectral solution", trim(intString1), "block", & - trim(cgnsDoms(nbkGlobal)%zoneName), "contains the following negative volumes" + trim(cgnsDoms(nbkGlobal)%zoneName), "contains the following ", trim(modeString) print "(a)", "#--------------------------------------------------------------------" print "(a)", "#" @@ -3464,8 +3665,20 @@ subroutine writeNegVolumes(checkVolDoms) do k = 2, kl do j = 2, jl do i = 2, il - if (checkVolDoms(nn, sps)%volumeIsNeg(i, j, k)) then + ! Test for a bad volume. + + volumeIsBad = .false. + select case (mode) + case (MetricVolume) + if (checkVolDoms(nn, sps)%volumeIsNeg(i, j, k)) & + volumeIsBad = .true. + case (MetricSkewness) + if (checkVolDoms(nn, sps)%volumeIsSkewed(i, j, k)) & + volumeIsBad = .true. + end select + + if (volumeIsBad) then xc(1:3) = eighth * (x(i - 1, j - 1, k - 1, 1:3) & + x(i, j - 1, k - 1, 1:3) & + x(i - 1, j, k - 1, 1:3) & @@ -3483,10 +3696,16 @@ subroutine writeNegVolumes(checkVolDoms) intString2 = adjustl(intString2) intString3 = adjustl(intString3) - print "(7(A), 4(ES10.3, A))", "# Indices (", trim(intString1), & - ",", trim(intString2), ",", trim(intString3), "), coordinates (", & - xc(1), ",", xc(2), ",", xc(3), "), Volume: ", -vol(i, j, k) + select case (mode) + case (MetricVolume) + quantity = -vol(i, j, k) + case (MetricSkewness) + quantity = skew(i, j, k) + end select + print "(7(A), 3(ES10.3, A), 2(A), F6.3)", "# Indices (", trim(intString1), & + ",", trim(intString2), ",", trim(intString3), "), coordinates (", & + xc(1), ",", xc(2), ",", xc(3), "), ", trim(descString), ": ", quantity end if end do end do @@ -3504,7 +3723,7 @@ subroutine writeNegVolumes(checkVolDoms) end do procLoop - end subroutine writeNegVolumes + end subroutine writeBadVolumes subroutine faceRotationMatrices(level, allocMem) ! ! faceRotationMatrices computes the rotation matrices on the diff --git a/src/preprocessing/preprocessingModules.F90 b/src/preprocessing/preprocessingModules.F90 index 4269176ba..2b569d1e3 100644 --- a/src/preprocessing/preprocessingModules.F90 +++ b/src/preprocessing/preprocessingModules.F90 @@ -775,6 +775,8 @@ module checkVolBlock logical :: blockHasNegVol logical, dimension(:, :, :), pointer :: volumeIsNeg + logical :: blockHasSkewedVol + logical, dimension(:, :, :), pointer :: volumeIsSkewed end type checkVolBlockType diff --git a/src/utils/utils.F90 b/src/utils/utils.F90 index 1652107c9..b50cc8c17 100644 --- a/src/utils/utils.F90 +++ b/src/utils/utils.F90 @@ -3379,6 +3379,8 @@ subroutine setPointers(nn, mm, ll) volRef => flowDoms(nn, mm, ll)%volRef volOld => flowDoms(nn, 1, ll)%volOld + skew => flowDoms(nn, mm, ll)%skew + porI => flowDoms(nn, mm, 1)%porI porJ => flowDoms(nn, mm, 1)%porJ porK => flowDoms(nn, mm, 1)%porK @@ -4823,6 +4825,7 @@ subroutine deallocateBlock(nn, level, sps) use inputUnsteady use inputPhysics use iteration + use inputIteration, only: useSkewnessCheck use block, only: viscSubfaceType, BCDataType, flowDoms implicit none ! @@ -5194,6 +5197,12 @@ subroutine deallocateBlock(nn, level, sps) deallocate (flowDoms(nn, level, sps)%vol, stat=ierr) if (ierr /= 0) deallocationFailure = .true. + if (useSkewnessCheck) then + if (associated(flowDoms(nn, level, sps)%skew)) & + deallocate (flowDoms(nn, level, sps)%skew, stat=ierr) + if (ierr /= 0) deallocationFailure = .true. + end if + if (associated(flowDoms(nn, level, sps)%volRef)) & deallocate (flowDoms(nn, level, sps)%volRef, stat=ierr) if (ierr /= 0) deallocationFailure = .true. @@ -6645,5 +6654,28 @@ subroutine getRotMatrix(vec1, vec2, rotMat) end subroutine getRotMatrix + real(kind=realType) function norm2cplx(v) + ! if input is real: + ! returns the norm of the input array + ! + ! if input is complex: + ! returns a complex number where the real part represents the norm of + ! all the real input parts and the complex part represents the norm of + ! all the complex input parts. + + use constants + + implicit none + + real(kind=realType), dimension(:), intent(in) :: v + +#ifdef USE_COMPLEX + norm2cplx = cmplx(NORM2(real(v)), NORM2(aimag(v))) +#else + norm2cplx = NORM2(v) +#endif + + end function norm2cplx + #endif end module utils diff --git a/tests/unit_tests/test_mesh_skewness.py b/tests/unit_tests/test_mesh_skewness.py new file mode 100644 index 000000000..6cd22180d --- /dev/null +++ b/tests/unit_tests/test_mesh_skewness.py @@ -0,0 +1,94 @@ +import unittest +import os +from adflow import ADFLOW +from idwarp import USMesh +import numpy as np +import sys + +baseDir = os.path.dirname(os.path.abspath(__file__)) +sys.path.append(os.path.join(baseDir, "../reg_tests")) + +from reg_aeroproblems import ap_tutorial_wing # noqa E402 + + +class BasicTests(unittest.TestCase): + N_PROCS = 1 + + def setUp(self): + gridFile = "input_files/cube_1x1x1.cgns" + self.options = { + "gridfile": os.path.join(baseDir, "../../", gridFile), + "useSkewnessCheck": True, + "meshMaxSkewness": 0.5, + "printBadlySkewedCells": True, + } + self.idwarp_options = { + "gridfile": os.path.join(baseDir, "../../", gridFile), + } + eps = 1e-12 + self.warp_no_fail = 1 - eps + self.warp_fail = 1 + eps + + # only here for debugging + self.write_mesh = False + + def test_skew_x(self): + hold, move = 1, 0 + + # warp the mesh slightly + CFDSolver = self.create_solver() + self.warp_mesh(CFDSolver, hold, move, self.warp_no_fail) + self.assertFalse(CFDSolver.adflow.killsignals.fatalfail) + + # warp it more + CFDSolver = self.create_solver() + self.warp_mesh(CFDSolver, hold, move, self.warp_fail) + self.assertTrue(CFDSolver.adflow.killsignals.fatalfail) + + def test_skew_y(self): + hold, move = 0, 1 + + # warp the mesh slightly + CFDSolver = self.create_solver() + self.warp_mesh(CFDSolver, hold, move, self.warp_no_fail) + self.assertFalse(CFDSolver.adflow.killsignals.fatalfail) + + # warp it more + CFDSolver = self.create_solver() + self.warp_mesh(CFDSolver, hold, move, self.warp_fail) + self.assertTrue(CFDSolver.adflow.killsignals.fatalfail) + + def test_skew_z(self): + hold, move = 0, 2 + + # warp the mesh slightly + CFDSolver = self.create_solver() + self.warp_mesh(CFDSolver, hold, move, self.warp_no_fail) + self.assertFalse(CFDSolver.adflow.killsignals.fatalfail) + + # warp it more + CFDSolver = self.create_solver() + self.warp_mesh(CFDSolver, hold, move, self.warp_fail) + self.assertTrue(CFDSolver.adflow.killsignals.fatalfail) + + def create_solver(self): + CFDSolver = ADFLOW(options=self.options, debug=False) + mesh = USMesh(options=self.idwarp_options) + CFDSolver.setMesh(mesh) + return CFDSolver + + def warp_mesh(self, CFDSolver, hold, move, amount): + # setup warping + coords = CFDSolver.mesh.getSurfaceCoordinates().copy() + index = np.argwhere(coords[:, hold] == 0) + coords[index, move] += amount + + # warp + CFDSolver.mesh.setSurfaceCoordinates(coords) + CFDSolver.setAeroProblem(ap_tutorial_wing) # dummy AP + if self.write_mesh: + CFDSolver.writeMeshFile(f"skewed_{hold}_{move}_{amount}.cgns") + + +if __name__ == "__main__": + unittest.main()