Merge branch 'master' into mixlayer

Author: sbryngelson

.github/workflows/frontier/submit.sh

@@ -33,10 +33,9 @@ sbatch <<EOT
 #SBATCH -A CFD154                  # charge account
 #SBATCH -N 1                       # Number of nodes required
 $sbatch_device_opts
-#SBATCH -t 01:59:00                # Duration of the job (Ex: 15 mins)
+#SBATCH -t 03:59:00                # Duration of the job (Ex: 15 mins)
 #SBATCH -o$job_slug.out            # Combined output and error messages file
 #SBATCH -p extended                # Extended partition for shorter queues
-#SBATCH -q debug                   # Use debug QOS - only one job per user allowed in queue!
 #SBATCH -W                         # Do not exit until the submitted job terminates.
 
 set -e

.github/workflows/lint-source.yml

@@ -42,6 +42,10 @@ jobs:
       run: |
         ! grep -iR -e '\.\.\.' -e '\-\-\-' -e '===' ./src/*
 
+    - name: Looking for false integers
+      run: |
+        ! grep -onRP '(?<![0-9.eE\-])\b[0-9]*_wp\b' src/
+
     - name: Looking for junk comments in examples
       run: |
         ! grep -R '# ===' ./benchmarks **/*.py

.github/workflows/pmd.yml

@@ -30,7 +30,102 @@ jobs:
               unzip -q pmd.zip
               PMD_HOME="pmd-bin-${PMD_VERSION}"
           
+              SOURCE_DIR="${1:-src}"
+              total_files=$(find "$SOURCE_DIR" -type f \( -name "*.f" -o -name "*.f90" -o -name "*.for" -o -name "*.fpp" -o -name "*.F" -o -name "*.F90" \) | wc -l)
+              processed=0
+
+              find "$SOURCE_DIR" -type f \( -name "*.f" -o -name "*.f90" -o -name "*.for" -o -name "*.fpp" -o -name "*.F" -o -name "*.F90" \) -print0 | 
+              while IFS= read -r -d $'\0' file; do
+                  processed=$((processed + 1))
+                  
+                  # Create a temporary file with same permissions as original
+                  TMP_FILE=$(mktemp)
+                  if [ $? -ne 0 ]; then
+                      echo -e "Failed to create temporary file for $file, skipping"
+                      continue
+                  fi
+                  
+                  # Copy permissions from original file
+                  chmod --reference="$file" "$TMP_FILE"
+                  
+                  # More comprehensive sed command to handle different Fortran comment styles:
+                  # 1. Replace lines that are entirely comments with an empty line:
+                  #    - Lines starting with '!' (free form comments)
+                  #    - Lines starting with 'c', 'C', '*', 'd', 'D' in column 1 (fixed form comments)
+                  # 2. Remove end-of-line comments (anything after '!' that isn't in a string)
+                  # 3. Preserve strings containing '!' characters                  
+                  sed -E '
+                    # First handle & continuation style (modern Fortran)
+                    :ampersand_loop
+                    /&[[:space:]]*$/ {
+                      N
+                      s/&[[:space:]]*\n[[:space:]]*(&)?/ /g
+                      tampersand_loop
+                    }
+                    
+                    # Handle fixed-form continuation (column 6 indicator)
+                    :fixed_form_loop
+                    /^[[:space:]]{0,5}[^[:space:]!&]/ {
+                      N
+                      s/\n[[:space:]]{5}[^[:space:]]/ /g
+                      tfixed_form_loop
+                    }
+                    
+                    # Remove any remaining continuation markers
+                    s/&//g
+                    
+                    # Normalize spacing - replace multiple spaces with single space
+                    s/[[:space:]]{2,}/ /g
+                    
+                    # Remove spaces around mathematical operators
+                    s/[[:space:]]*\*[[:space:]]*/*/g
+                    s/[[:space:]]*\+[[:space:]]*/+/g
+                    s/[[:space:]]*-[[:space:]]*/-/g
+                    s/[[:space:]]*\/[[:space:]]*/\//g
+                    s/[[:space:]]*\*\*[[:space:]]*/\*\*/g
+                    
+                    # Remove spaces in common Fortran constructs (array indexing, function calls)
+                    s/\([[:space:]]*([^,)[:space:]]+)[[:space:]]*,/(\1,/g      # First argument
+                    s/,[[:space:]]*([^,)[:space:]]+)[[:space:]]*,/,\1,/g       # Middle arguments
+                    s/,[[:space:]]*([^,)[:space:]]+)[[:space:]]*\)/,\1)/g      # Last argument
+                    s/\([[:space:]]*([^,)[:space:]]+)[[:space:]]*\)/(\1)/g     # Single argument
+                    
+                    # Remove spaces around brackets and parentheses
+                    s/\[[[:space:]]*/</g
+                    s/\[[[:space:]]*/>/g
+                    s/\[[[:space:]]*/</g
+                    s/[[:space:]]*\]/]/g
+                    s/\([[:space:]]*/(/g
+                    s/[[:space:]]*\)/)/g
+
+                    # Remove spaces around comparison operators
+                    s/[[:space:]]*<=[[:space:]]*/</g
+                    s/[[:space:]]*>=[[:space:]]*/>/g
+                    s/[[:space:]]*<[[:space:]]*/</g
+                    s/[[:space:]]*>[[:space:]]*/>/g
+                    s/[[:space:]]*==[[:space:]]*/==/g
+
+                    # Remove full-line comments
+                    /^\s*!/d
+                    /^[cC*dD]/d
+                    /^[ \t]*[cC*dD]/d
+                    /^[[:space:]]*$/d
+
+                    # Remove end-of-line comments, preserving quoted strings
+                    s/([^"'\''\\]*("[^"]*")?('\''[^'\'']*'\''?)?[^"'\''\\]*)[!].*$/\1/
+                  ' "$file" > "$TMP_FILE"
+
+                  if cmp -s "$file" "$TMP_FILE"; then
+                      rm "$TMP_FILE"
+                  else
+                      # Overwrite the original file with the processed content
+                      mv "$TMP_FILE" "$file"
+                  fi
+              done
+
               "${PMD_HOME}/bin/pmd" cpd \
                 --dir src \
                 --language fortran \
-                --minimum-tokens=40
+                --minimum-tokens=20 \
+                --no-fail-on-violation \
+                --no-fail-on-error

src/common/m_checker_common.fpp

@@ -322,7 +322,7 @@ contains
         @:PROHIBIT(surface_tension .and. sigma < 0._wp, &
             "sigma must be greater than or equal to zero")
 
-        @:PROHIBIT(surface_tension .and. sigma == dflt_real, &
+        @:PROHIBIT(surface_tension .and. f_approx_equal(sigma, dflt_real), &
             "sigma must be set if surface_tension is enabled")
 
         @:PROHIBIT(.not. f_is_default(sigma) .and. .not. surface_tension, &
@@ -347,9 +347,12 @@ contains
         !! Called by s_check_inputs_common for all three stages
     impure subroutine s_check_inputs_moving_bc
         #:for X, VB2, VB3 in [('x', 'vb2', 'vb3'), ('y', 'vb3', 'vb1'), ('z', 'vb1', 'vb2')]
-            if (any((/bc_${X}$%vb1, bc_${X}$%vb2, bc_${X}$%vb3/) /= 0._wp)) then
+            if (.not. (f_approx_equal(bc_${X}$%vb1, 0._wp) .and. &
+                       f_approx_equal(bc_${X}$%vb2, 0._wp) .and. &
+                       f_approx_equal(bc_${X}$%vb3, 0._wp))) then
                 if (bc_${X}$%beg == BC_SLIP_WALL) then
-                    if (any((/bc_${X}$%${VB2}$, bc_${X}$%${VB3}$/) /= 0._wp)) then
+                    if (.not. (f_approx_equal(bc_${X}$%${VB2}$, 0._wp) .and. &
+                               f_approx_equal(bc_${X}$%${VB3}$, 0._wp))) then
                         call s_mpi_abort("bc_${X}$%beg must be -15 if "// &
                                          "bc_${X}$%${VB2}$ or bc_${X}$%${VB3}$ "// &
                                          "is set. Exiting.", CASE_FILE_ERROR_CODE)
@@ -362,9 +365,12 @@ contains
         #:endfor
 
         #:for X, VE2, VE3 in [('x', 've2', 've3'), ('y', 've3', 've1'), ('z', 've1', 've2')]
-            if (any((/bc_${X}$%ve1, bc_${X}$%ve2, bc_${X}$%ve3/) /= 0._wp)) then
+            if (.not. (f_approx_equal(bc_${X}$%ve1, 0._wp) .and. &
+                       f_approx_equal(bc_${X}$%ve2, 0._wp) .and. &
+                       f_approx_equal(bc_${X}$%ve3, 0._wp))) then
                 if (bc_${X}$%end == BC_SLIP_WALL) then
-                    if (any((/bc_${X}$%${VE2}$, bc_${X}$%${VE3}$/) /= 0._wp)) then
+                    if (.not. (f_approx_equal(bc_${X}$%${VE2}$, 0._wp) .and. &
+                               f_approx_equal(bc_${X}$%${VE3}$, 0._wp))) then
                         call s_mpi_abort("bc_${X}$%end must be -15 if "// &
                                          "bc_${X}$%${VE2}$ or bc_${X}$%${VE3}$ "// &
                                          "is set. Exiting.", CASE_FILE_ERROR_CODE)

src/common/m_constants.fpp

@@ -8,9 +8,9 @@ module m_constants
 
     character, parameter :: dflt_char = ' ' !< Default string value
 
-    real(wp), parameter :: dflt_real = -1e6_wp                !< Default real value
-    real(wp), parameter :: sgm_eps = 1e-16_wp               !< Segmentation tolerance
-    real(wp), parameter :: small_alf = 1e-11_wp                !< Small alf tolerance
+    real(wp), parameter :: dflt_real = -1.e6_wp                !< Default real value
+    real(wp), parameter :: sgm_eps = 1.e-16_wp               !< Segmentation tolerance
+    real(wp), parameter :: small_alf = 1.e-11_wp                !< Small alf tolerance
     real(wp), parameter :: pi = 3.141592653589793_wp !< Pi
     real(wp), parameter :: verysmall = 1.e-12_wp              !< Very small number
 
@@ -26,7 +26,7 @@ module m_constants
     integer, parameter :: pathlen_max = 400
     integer, parameter :: nnode = 4    !< Number of QBMM nodes
     integer, parameter :: gp_layers = 3 !< Number of ghost point layers for IBM
-    real(wp), parameter :: capillary_cutoff = 1e-6 !< color function gradient magnitude at which to apply the surface tension fluxes
+    real(wp), parameter :: capillary_cutoff = 1.e-6 !< color function gradient magnitude at which to apply the surface tension fluxes
     real(wp), parameter :: acoustic_spatial_support_width = 2.5_wp !< Spatial support width of acoustic source, used in s_source_spatial
     real(wp), parameter :: dflt_vcfl_dt = 100._wp !< value of vcfl_dt when viscosity is off for computing adaptive timestep size
     real(wp), parameter :: broadband_spectral_level_constant = 20._wp !< The constant to scale the spectral level at the lower frequency bound
@@ -41,25 +41,25 @@ module m_constants
     integer, parameter :: Ifactor_bary_3D = 20 !< Multiple factor of the ratio (triangle area to cell face area) for interpolation on triangle facets for 3D models
     integer, parameter :: num_ray = 20 !< Default number of rays traced per cell
     real(wp), parameter :: ray_tracing_threshold = 0.9_wp !< Threshold above which the cell is marked as the model patch
-    real(wp), parameter :: threshold_vector_zero = 1e-10 !< Threshold to treat the component of a vector to be zero
-    real(wp), parameter :: threshold_edge_zero = 1e-10 !< Threshold to treat two edges to be overlapped
-    real(wp), parameter :: threshold_bary = 1e-1 !< Threshold to interpolate a barycentric facet
-    real(wp), parameter :: initial_distance_buffer = 1e12_wp !< Initialized levelset distance for the shortest path pair algorithm
+    real(wp), parameter :: threshold_vector_zero = 1.e-10_wp !< Threshold to treat the component of a vector to be zero
+    real(wp), parameter :: threshold_edge_zero = 1.e-10_wp !< Threshold to treat two edges to be overlapped
+    real(wp), parameter :: threshold_bary = 1.e-1_wp !< Threshold to interpolate a barycentric facet
+    real(wp), parameter :: initial_distance_buffer = 1.e12_wp !< Initialized levelset distance for the shortest path pair algorithm
 
     ! Lagrange bubbles constants
     integer, parameter :: mapCells = 3 !< Number of cells around the bubble where the smoothening function will have effect
     real(wp), parameter :: R_uni = 8314._wp ! Universal gas constant - J/kmol/K
 
     ! Strang Splitting constants
-    real(wp), parameter :: dflt_adap_dt_tol = 1e-4_wp !< Default tolerance for adaptive step size
+    real(wp), parameter :: dflt_adap_dt_tol = 1.e-4_wp !< Default tolerance for adaptive step size
     integer, parameter :: adap_dt_max_iters = 100 !< Maximum number of iterations
     ! Constants of the algorithm described by Heirer, E. Hairer S.P.Nørsett G. Wanner, Solving Ordinary Differential Equations I, Chapter II.4
     ! to choose the initial time step size for the adaptive time stepping routine
-    real(wp), parameter :: threshold_first_guess = 1e-5_wp
-    real(wp), parameter :: threshold_second_guess = 1e-15_wp
-    real(wp), parameter :: scale_first_guess = 1e-3_wp
-    real(wp), parameter :: scale_guess = 1e-2_wp
-    real(wp), parameter :: small_guess = 1e-6_wp
+    real(wp), parameter :: threshold_first_guess = 1.e-5_wp
+    real(wp), parameter :: threshold_second_guess = 1.e-15_wp
+    real(wp), parameter :: scale_first_guess = 1.e-3_wp
+    real(wp), parameter :: scale_guess = 1.e-2_wp
+    real(wp), parameter :: small_guess = 1.e-6_wp
 
     ! Relativity
     integer, parameter :: relativity_cons_to_prim_max_iter = 100

src/common/m_eigen_solver.f90

@@ -10,6 +10,8 @@ module m_eigen_solver
 
     use m_precision_select
 
+    use m_helper_basic         !< Functions to compare floating point numbers
+
     implicit none
 
     private; 
@@ -124,7 +126,7 @@ pure subroutine cbal(nm, nl, ar, ai, low, igh, scale)
 
             do 110 i = 1, l
                 if (i == j) go to 110
-                if (ar(j, i) /= 0.0_wp .or. ai(j, i) /= 0.0_wp) go to 120
+                if (.not. f_approx_equal(ar(j, i), 0.0_wp) .or. .not. f_approx_equal(ai(j, i), 0.0_wp)) go to 120
 110         end do
 
             ml = l
@@ -140,7 +142,7 @@ pure subroutine cbal(nm, nl, ar, ai, low, igh, scale)
 
             do 150 i = k, l
                 if (i == j) go to 150
-                if (ar(i, j) /= 0.0_wp .or. ai(i, j) /= 0.0_wp) go to 170
+                if (.not. f_approx_equal(ar(i, j), 0.0_wp) .or. .not. f_approx_equal(ai(i, j), 0.0_wp)) go to 170
 150         end do
 
             ml = k
@@ -164,7 +166,7 @@ pure subroutine cbal(nm, nl, ar, ai, low, igh, scale)
                 r = r + abs(ar(i, j)) + abs(ai(i, j))
 200         end do
 !     guard against zero c or r due to underflow
-            if (c == 0.0_wp .or. r == 0.0_wp) go to 270
+            if (f_approx_equal(c, 0.0_wp) .or. f_approx_equal(r, 0.0_wp)) go to 270
             g = r/radix
             f = 1.0_wp
             s = c + r
@@ -242,7 +244,7 @@ pure subroutine corth(nm, nl, low, igh, ar, ai, ortr, orti)
             do 90 i = ml, igh
                 scale = scale + abs(ar(i, ml - 1)) + abs(ai(i, ml - 1))
 90          end do
-            if (scale == 0._wp) go to 180
+            if (f_approx_equal(scale, 0._wp)) go to 180
             mp = ml + igh
 !     for i=igh step -1 until ml do
             do 100 ii = ml, igh
@@ -254,7 +256,7 @@ pure subroutine corth(nm, nl, low, igh, ar, ai, ortr, orti)
 
             g = sqrt(h)
             call pythag(ortr(ml), orti(ml), f)
-            if (f == 0._wp) go to 103
+            if (f_approx_equal(f, 0._wp)) go to 103
             h = h + f*g
             g = g/f
             ortr(ml) = (1.0_wp + g)*ortr(ml)
@@ -374,8 +376,8 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
 !     for i=igh-1 step -1 until low+1 do
 105     do 140 ii = 1, iend
             i = igh - ii
-            if (abs(ortr(i)) == 0._wp .and. abs(orti(i)) == 0._wp) go to 140
-            if (abs(hr(i, i - 1)) == 0._wp .and. abs(hi(i, i - 1)) == 0._wp) go to 140
+            if (f_approx_equal(abs(ortr(i)), 0._wp) .and. f_approx_equal(abs(orti(i)), 0._wp)) go to 140
+            if (f_approx_equal(abs(hr(i, i - 1)), 0._wp) .and. f_approx_equal(abs(hi(i, i - 1)), 0._wp)) go to 140
 !     norm below is negative of h formed in corth
             norm = hr(i, i - 1)*ortr(i) + hi(i, i - 1)*orti(i)
             ip1 = i + 1
@@ -410,7 +412,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
 
         do 170 i = l, igh
             ll = min0(i + 1, igh)
-            if (abs(hi(i, i - 1)) == 0._wp) go to 170
+            if (f_approx_equal(abs(hi(i, i - 1)), 0._wp)) go to 170
             call pythag(hr(i, i - 1), hi(i, i - 1), norm)
             yr = hr(i, i - 1)/norm
             yi = hi(i, i - 1)/norm
@@ -458,7 +460,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
             tst1 = abs(hr(l - 1, l - 1)) + abs(hi(l - 1, l - 1)) &
                    + abs(hr(l, l)) + abs(hi(l, l))
             tst2 = tst1 + abs(hr(l, l - 1))
-            if (tst2 == tst1) go to 300
+            if (f_approx_equal(tst2, tst1)) go to 300
 260     end do
 !     form shift
 300     if (l == en) go to 660
@@ -468,7 +470,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
         si = hi(en, en)
         xr = hr(enm1, en)*hr(en, enm1)
         xi = hi(enm1, en)*hr(en, enm1)
-        if (xr == 0.0_wp .and. xi == 0.0_wp) go to 340
+        if (f_approx_equal(xr, 0.0_wp) .and. f_approx_equal(xi, 0.0_wp)) go to 340
         yr = (hr(enm1, enm1) - sr)/2.0_wp
         yi = (hi(enm1, enm1) - si)/2.0_wp
         call csroot(yr**2 - yi**2 + xr, 2.0_wp*yr*yi + xi, zzr, zzi)
@@ -522,7 +524,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
 500     end do
 
         si = hi(en, en)
-        if (abs(si) == 0._wp) go to 540
+        if (f_approx_equal(abs(si), 0._wp)) go to 540
         call pythag(hr(en, en), si, norm)
         sr = hr(en, en)/norm
         si = si/norm
@@ -567,7 +569,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
 590         end do
 600     end do
 
-        if (abs(si) == 0._wp) go to 240
+        if (f_approx_equal(abs(si), 0._wp)) go to 240
 
         do 630 i = 1, en
             yr = hr(i, en)
@@ -602,7 +604,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
             end do
         end do
 
-        if (nl == 1 .or. norm == 0._wp) go to 1001
+        if (nl == 1 .or. f_approx_equal(norm, 0._wp)) go to 1001
 !     for en=nl step -1 until 2 do
         do 800 nn = 2, nl
             en = nl + 2 - nn
@@ -625,7 +627,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
 
                 yr = xr - wr(i)
                 yi = xi - wi(i)
-                if (yr /= 0.0_wp .or. yi /= 0.0_wp) go to 765
+                if (.not. f_approx_equal(yr, 0.0_wp) .or. .not. f_approx_equal(yi, 0.0_wp)) go to 765
                 tst1 = norm
                 yr = tst1
 760             yr = 0.01_wp*yr
@@ -635,7 +637,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
                 call cdiv(zzr, zzi, yr, yi, hr(i, en), hi(i, en))
 !     overflow control
                 tr = abs(hr(i, en)) + abs(hi(i, en))
-                if (tr == 0.0_wp) go to 780
+                if (f_approx_equal(tr, 0.0_wp)) go to 780
                 tst1 = tr
                 tst2 = tst1 + 1.0_wp/tst1
                 if (tst2 > tst1) go to 780
@@ -796,11 +798,11 @@ pure elemental subroutine pythag(a, b, c)
 
         real(wp) :: p, r, s, t, u
         p = max(abs(a), abs(b))
-        if (p == 0.0_wp) go to 20
+        if (f_approx_equal(p, 0.0_wp)) go to 20
         r = (min(abs(a), abs(b))/p)**2
 10      continue
         t = 4.0_wp + r
-        if (t == 4.0_wp) go to 20
+        if (f_approx_equal(t, 4.0_wp)) go to 20
         s = r/t
         u = 1.0_wp + 2.0_wp*s
         p = u*p