raise Error if a Mesh Cell is skewed above threshold (#235)

* raise Error if a Mesh Cell is skewed above threshold

* Fix calculation of skewness, correct typos

* Add skewness test

* Run flacke8 and black

* Add documentation

* fix print typo

* Update error printing to new format

* Run fprettify

* Fix bug where skewness code is executed when deactivated

* Add constants for choosen metric

* run fprettify

* deallocate skew variable

* Make complex to compile

* run fprettify

* Make error prints nice again

* Tidy up skewness branch

* Run fprettify

* Add back lost space

* inflating my commit counts by pushing newline changes

---------

Co-authored-by: Eirikur Jonsson <36180221+eirikurj@users.noreply.github.com>
Co-authored-by: Anil Yildirim <anily@umich.edu>

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"],

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.
+    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 <https://www.simscale.com/docs/simulation-setup/meshing/mesh-quality/#skewness>`__. 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``.
+

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
+#endif

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
 

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

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

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