implement fixes and expand to cover more cases of double precision

Author: arciyer123

src/common/m_constants.fpp

@@ -8,11 +8,11 @@ module m_constants
 
     character, parameter :: dflt_char = ' ' !< Default string value
 
-    real(wp), parameter :: dflt_real = -1d6                !< Default real value
-    real(wp), parameter :: sgm_eps = 1e-16               !< Segmentation tolerance
-    real(wp), parameter :: small_alf = 1e-11                !< Small alf tolerance
+    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 :: pi = 3.141592653589793_wp !< Pi
-    real(wp), parameter :: verysmall = 1.e-12              !< Very small number
+    real(wp), parameter :: verysmall = 1.e-12_wp              !< Very small number
 
     integer, parameter :: num_stcls_min = 5    !< Minimum # of stencils
     integer, parameter :: path_len = 400  !< Maximum path length
@@ -24,7 +24,7 @@ module m_constants
     integer, parameter :: num_patches_max = 10
     integer, parameter :: pathlen_max = 400
     integer, parameter :: nnode = 4    !< Number of QBMM nodes
-    real(wp), parameter :: capillary_cutoff = 1e-6 !< color function gradient magnitude at which to apply the surface tension fluxes
+    real(wp), parameter :: capillary_cutoff = 1e-6_wp !< 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
 

src/common/m_helper.fpp

@@ -153,7 +153,7 @@ contains
         if (thermal == 2) gamma_m = 1._wp
 
         temp = 293.15_wp
-        D_m = 0.242e-4
+        D_m = 0.242e-4_wp
         uu = sqrt(pl0/rhol0)
 
         omega_ref = 3._wp*k_poly*Ca + 2._wp*(3._wp*k_poly - 1._wp)/Web
@@ -260,17 +260,17 @@ contains
         real(wp), dimension(nb) :: phi
 
         ! nondiml. min. & max. initial radii for numerical quadrature
-        !sd   = 0.05D0
-        !R0mn = 0.75D0
-        !R0mx = 1.3D0
+        !sd   = 0.05e0_wp
+        !R0mn = 0.75e0_wp
+        !R0mx = 1.3e0_wp
 
-        !sd   = 0.3D0
-        !R0mn = 0.3D0
-        !R0mx = 6.D0
+        !sd   = 0.3e0_wp
+        !R0mn = 0.3e0_wp
+        !R0mx = 6.e0_wp
 
-        !sd   = 0.7D0
-        !R0mn = 0.12D0
-        !R0mx = 150.D0
+        !sd   = 0.7e0_wp
+        !R0mn = 0.12e0_wp
+        !R0mx = 150.e0_wp
 
         sd = poly_sigma
         R0mn = 0.8_wp*exp(-2.8_wp*sd)

src/common/m_helper_basic.f90

@@ -22,7 +22,7 @@ module m_helper_basic
     !> This procedure checks if two floating point numbers of wp are within tolerance.
     !! @param a First number.
     !! @param b Second number.
-    !! @param tol_input Relative error (default = 1e-6).
+    !! @param tol_input Relative error (default = 1e-6_wp).
     !! @return Result of the comparison.
     logical function f_approx_equal(a, b, tol_input) result(res)
         !$acc routine seq
@@ -35,7 +35,7 @@ logical function f_approx_equal(a, b, tol_input) result(res)
         if (present(tol_input)) then
             tol = tol_input
         else
-            tol = 1e-6
+            tol = 1e-6_wp
         end if
 
         if (a == b) then

src/common/m_phase_change.fpp

@@ -51,10 +51,10 @@ module m_phase_change
 
     !> @name Parameters for the first order transition phase change
     !> @{
-    integer, parameter :: max_iter = 1e8        !< max # of iterations
-    real(wp), parameter :: pCr = 4.94d7   !< Critical water pressure
+    integer, parameter :: max_iter = 1e8_wp        !< max # of iterations
+    real(wp), parameter :: pCr = 4.94e7_wp   !< Critical water pressure
     real(wp), parameter :: TCr = 385.05 + 273.15  !< Critical water temperature
-    real(wp), parameter :: mixM = 1.0e-8 !< threshold for 'mixture cell'. If Y < mixM, phase change does not happen
+    real(wp), parameter :: mixM = 1.0e-8_wp !< threshold for 'mixture cell'. If Y < mixM, phase change does not happen
     integer, parameter :: lp = 1    !< index for the liquid phase of the reacting fluid
     integer, parameter :: vp = 2    !< index for the vapor phase of the reacting fluid
     !> @}
@@ -148,7 +148,7 @@ contains
                     !$acc loop seq
                     do i = momxb, momxe
 
-                        dynE = dynE + 5.0e-1*q_cons_vf(i)%sf(j, k, l)**2/rho
+                        dynE = dynE + 5.0e-1_wp*q_cons_vf(i)%sf(j, k, l)**2/rho
 
                     end do
 
@@ -352,12 +352,12 @@ contains
 
         ! Maybe improve this condition afterwards. As long as the initial guess is in between -min(ps_inf)
         ! and infinity, a solution should be able to be found.
-        pS = 1.0d4
+        pS = 1.0e4_wp
 
         ! Newton Solver for the pT-equilibrium
         ns = 0
-        ! change this relative error metric. 1E4 is just arbitrary
-        do while ((abs(pS - pO) > palpha_eps) .and. (abs((pS - pO)/pO) > palpha_eps/1e4) .or. (ns == 0))
+        ! change this relative error metric. 1e4_wp is just arbitrary
+        do while ((abs(pS - pO) > palpha_eps) .and. (abs((pS - pO)/pO) > palpha_eps/1e4_wp) .or. (ns == 0))
 
             ! increasing counter
             ns = ns + 1
@@ -425,17 +425,17 @@ contains
         ns = 0
 
         ! Relaxation factor
-        Om = 1.0e-3
+        Om = 1.0e-3_wp
 
         p_infpTg = p_infpT
 
         if (((pS < 0.0_wp) .and. ((q_cons_vf(lp + contxb - 1)%sf(j, k, l) &
                                    + q_cons_vf(vp + contxb - 1)%sf(j, k, l)) > ((rhoe &
                                                                                  - gs_min(lp)*ps_inf(lp)/(gs_min(lp) - 1))/qvs(lp)))) .or. &
-            ((pS >= 0.0_wp) .and. (pS < 1.0e-1))) then
+            ((pS >= 0.0_wp) .and. (pS < 1.0e-1_wp))) then
 
             ! improve this initial condition
-            pS = 1.0d4
+            pS = 1.0e4_wp
 
         end if
 
@@ -447,7 +447,7 @@ contains
         R2D(1) = 0.0_wp; R2D(2) = 0.0_wp
         DeltamP(1) = 0.0_wp; DeltamP(2) = 0.0_wp
         do while (((sqrt(R2D(1)**2 + R2D(2)**2) > ptgalpha_eps) &
-                   .and. ((sqrt(R2D(1)**2 + R2D(2)**2)/rhoe) > (ptgalpha_eps/1d6))) &
+                   .and. ((sqrt(R2D(1)**2 + R2D(2)**2)/rhoe) > (ptgalpha_eps/1e6_wp))) &
                   .or. (ns == 0))
 
             ! Updating counter for the iterative procedure
@@ -754,7 +754,7 @@ contains
             ns = 0
 
             ! underrelaxation factor
-            Om = 1.0e-3
+            Om = 1.0e-3_wp
             do while ((abs(FT) > ptgalpha_eps) .or. (ns == 0))
                 ! increasing counter
                 ns = ns + 1

src/common/m_precision_select.f90

@@ -14,7 +14,7 @@ module m_precision_select
     integer, parameter :: double_precision = selected_real_kind(15, 307)
 
     ! Set the working precision (wp) to single or double precision
-    integer, parameter :: wp = single_precision  ! Change this to single_precision if needed
+    integer, parameter :: wp = double_precision  ! Change this to single_precision if needed
 
 #ifdef MFC_MPI
     ! Declare mpi_p as a module variable

src/common/m_variables_conversion.fpp

@@ -189,7 +189,7 @@ contains
                 Y_rs(i) = rhoYks(i)/rho
             end do
 
-            if (sum(Y_rs) > 1e-16) then
+            if (sum(Y_rs) > 1e-16_wp) then
                 call get_temperature(.true., energy - dyn_p, 1200._wp, Y_rs, T)
                 call get_pressure(rho, T, Y_rs, pres)
             else
@@ -295,7 +295,7 @@ contains
                 alpha_K(i) = min(max(0._wp, alpha_K(i)), 1._wp)
             end do
 
-            alpha_K = alpha_K/max(sum(alpha_K), 1e-16)
+            alpha_K = alpha_K/max(sum(alpha_K), 1e-16_wp)
 
         end if
 
@@ -420,7 +420,7 @@ contains
                 alpha_K(i) = min(max(0._wp, alpha_K(i)), 1._wp)
             end do
 
-            alpha_K = alpha_K/max(sum(alpha_K), 1e-16)
+            alpha_K = alpha_K/max(sum(alpha_K), 1e-16_wp)
 
         end if
 
@@ -987,7 +987,7 @@ contains
                         if (model_eqns /= 4) then
                             qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l) &
                                                         /rho_K
-                            dyn_pres_K = dyn_pres_K + 5e-1*qK_cons_vf(i)%sf(j, k, l) &
+                            dyn_pres_K = dyn_pres_K + 5e-1_wp*qK_cons_vf(i)%sf(j, k, l) &
                                          *qK_prim_vf(i)%sf(j, k, l)
                         else
                             qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l) &
@@ -1349,7 +1349,7 @@ contains
 
                     ! Computing the energy from the pressure
                     E_K = gamma_K*pres_K + pi_inf_K &
-                          + 5e-1*rho_K*vel_K_sum + qv_K
+                          + 5e-1_wp*rho_K*vel_K_sum + qv_K
 
                     ! mass flux, this should be \alpha_i \rho_i u_i
                     !$acc loop seq
@@ -1468,7 +1468,7 @@ contains
                     (rho*(1._wp - adv(num_fluids)))
             end if
         else
-            c = ((H - 5e-1*vel_sum)/gamma)
+            c = ((H - 5e-1_wp*vel_sum)/gamma)
         end if
 
         if (mixture_err .and. c < 0._wp) then

src/post_process/m_derived_variables.fpp

@@ -212,7 +212,7 @@ contains
                     end if
 
                     if (mixture_err .and. q_sf(i, j, k) < 0._wp) then
-                        q_sf(i, j, k) = 1e-16
+                        q_sf(i, j, k) = 1e-16_wp
                     else
                         q_sf(i, j, k) = sqrt(q_sf(i, j, k))
                     end if
@@ -285,8 +285,8 @@ contains
                         end if
                     end if
 
-                    if (abs(top) < 1e-8) top = 0._wp
-                    if (abs(bottom) < 1e-8) bottom = 0._wp
+                    if (abs(top) < 1e-8_wp) top = 0._wp
+                    if (abs(bottom) < 1e-8_wp) bottom = 0._wp
 
                     if (top == bottom) then
                         slope = 1._wp
@@ -295,20 +295,20 @@ contains
                         !           (bottom == 0._wp .AND. top /= 0._wp)) THEN
                         !           slope = 0._wp
                     else
-                        slope = (top*bottom)/(bottom**2._wp + 1e-16)
+                        slope = (top*bottom)/(bottom**2._wp + 1e-16_wp)
                     end if
 
                     ! Flux limiter function
                     if (flux_lim == 1) then ! MINMOD (MM)
                         q_sf(j, k, l) = max(0._wp, min(1._wp, slope))
                     elseif (flux_lim == 2) then ! MUSCL (MC)
-                        q_sf(j, k, l) = max(0._wp, min(2._wp*slope, 5e-1*(1._wp + slope), 2._wp))
+                        q_sf(j, k, l) = max(0._wp, min(2._wp*slope, 5e-1_wp*(1._wp + slope), 2._wp))
                     elseif (flux_lim == 3) then ! OSPRE (OP)
-                        q_sf(j, k, l) = (15e-1*(slope**2._wp + slope))/(slope**2._wp + slope + 1._wp)
+                        q_sf(j, k, l) = (15e-1_wp*(slope**2._wp + slope))/(slope**2._wp + slope + 1._wp)
                     elseif (flux_lim == 4) then ! SUPERBEE (SB)
                         q_sf(j, k, l) = max(0._wp, min(1._wp, 2._wp*slope), min(slope, 2._wp))
                     elseif (flux_lim == 5) then ! SWEBY (SW) (beta = 1.5)
-                        q_sf(j, k, l) = max(0._wp, min(15e-1*slope, 1._wp), min(slope, 15e-1))
+                        q_sf(j, k, l) = max(0._wp, min(15e-1_wp*slope, 1._wp), min(slope, 15e-1_wp))
                     elseif (flux_lim == 6) then ! VAN ALBADA (VA)
                         q_sf(j, k, l) = (slope**2._wp + slope)/(slope**2._wp + 1._wp)
                     elseif (flux_lim == 7) then ! VAN LEER (VL)

src/pre_process/include/2dHardcodedIC.fpp

@@ -5,7 +5,7 @@
 
     real(wp) :: rhoH, rhoL, pRef, pInt, h, lam, wl, amp, intH, alph
 
-    eps = 1e-9
+    eps = 1e-9_wp
 
 #:enddef
 
@@ -71,7 +71,7 @@
     case (204) ! Rayleigh-Taylor instability
         rhoH = 3
         rhoL = 1
-        pRef = 1e5
+        pRef = 1e5_wp
         pInt = pRef
         h = 0.7
         lam = 0.2
@@ -80,7 +80,7 @@
 
         intH = amp*sin(2*pi*x_cc(i)/lam - pi/2) + h
 
-        alph = 5e-1*(1 + tanh((y_cc(j) - intH)/2.5e-3))
+        alph = 5e-1_wp*(1 + tanh((y_cc(j) - intH)/2.5e-3_wp))
 
         if (alph < eps) alph = eps
         if (alph > 1 - eps) alph = 1 - eps

src/pre_process/include/3dHardcodedIC.fpp

@@ -5,7 +5,7 @@
 
     real(wp) :: eps
 
-    eps = 1e-9
+    eps = 1e-9_wp
 #:enddef
 
 #:def Hardcoded3D()
@@ -14,7 +14,7 @@
     case (300) ! Rayleigh-Taylor instability
         rhoH = 3
         rhoL = 1
-        pRef = 1e5
+        pRef = 1e5_wp
         pInt = pRef
         h = 0.7
         lam = 0.2
@@ -23,7 +23,7 @@
 
         intH = amp*(sin(2*pi*x_cc(i)/lam - pi/2) + sin(2*pi*z_cc(k)/lam - pi/2)) + h
 
-        alph = 5e-1*(1 + tanh((y_cc(j) - intH)/2.5e-3))
+        alph = 5e-1_wp*(1 + tanh((y_cc(j) - intH)/2.5e-3_wp))
 
         if (alph < eps) alph = eps
         if (alph > 1 - eps) alph = 1 - eps

src/pre_process/m_assign_variables.fpp

@@ -197,7 +197,7 @@ contains
         #:endif
 
         ! Updating the patch identities bookkeeping variable
-        if (1._wp - eta < 1e-16) patch_id_fp(j, k, l) = patch_id
+        if (1._wp - eta < 1e-16_wp) patch_id_fp(j, k, l) = patch_id
 
     end subroutine s_assign_patch_mixture_primitive_variables
 
@@ -216,7 +216,7 @@ contains
         real(wp) :: R3bar, n0, ratio, nH, vfH, velH, rhoH, deno
 
         p0 = 101325
-        pres_mag = 1e-1
+        pres_mag = 1e-1_wp
         loc = x_cc(177)
         n_tait = fluid_pp(1)%gamma
         B_tait = fluid_pp(1)%pi_inf
@@ -669,7 +669,7 @@ contains
         end if
 
         ! Updating the patch identities bookkeeping variable
-        if (1._wp - eta < 1e-16) patch_id_fp(j, k, l) = patch_id
+        if (1._wp - eta < 1e-16_wp) patch_id_fp(j, k, l) = patch_id
 
     end subroutine s_assign_patch_species_primitive_variables