working precision

sbryngelson · sbryngelson · commit 186428543a7f · 2024-12-16T12:46:25.000-05:00
diff --git a/src/common/m_helper.fpp b/src/common/m_helper.fpp
@@ -478,11 +478,11 @@ contains
     recursive function unassociated_legendre(x, l) result(P)
 
         integer, intent(in) :: l
-        real(kind(0d0)), intent(in) :: x
-        real(kind(0d0)) :: P
+        real(wp), intent(in) :: x
+        real(wp) :: P
 
         if (l == 0) then
-            P = 1d0
+            P = 1_wp
         else if (l == 1) then
             P = x
         else
@@ -500,15 +500,15 @@ contains
     recursive function spherical_harmonic_func(x, phi, l, m) result(Y)
 
         integer, intent(in) :: l, m
-        real(kind(0d0)), intent(in) :: x, phi
-        real(kind(0d0)) :: Y, prefactor, pi
+        real(wp), intent(in) :: x, phi
+        real(wp) :: Y, prefactor, pi
 
-        pi = acos(-1d0)
+        pi = acos(-1_wp)
         prefactor = sqrt((2*l + 1)/(4*pi)*factorial(l - m)/factorial(l + m)); 
         if (m == 0) then
             Y = prefactor*associated_legendre(x, l, m); 
         elseif (m > 0) then
-            Y = (-1d0)**m*sqrt(2d0)*prefactor*associated_legendre(x, l, m)*cos(m*phi); 
+            Y = (-1_wp)**m*sqrt(2_wp)*prefactor*associated_legendre(x, l, m)*cos(m*phi); 
         end if
 
     end function spherical_harmonic_func
@@ -522,8 +522,8 @@ contains
     recursive function associated_legendre(x, l, m) result(P)
 
         integer, intent(in) :: l, m
-        real(kind(0d0)), intent(in) :: x
-        real(kind(0d0)) :: P
+        real(wp), intent(in) :: x
+        real(wp) :: P
 
         if (m <= 0 .and. l <= 0) then
             P = 1; 
diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp
@@ -221,9 +221,9 @@ contains
 
         implicit none
         integer, intent(in) :: counts          ! Array of vector lengths for each process
-        real(kind(0d0)), intent(in), dimension(counts) :: my_vector   ! Input vector on each process
+        real(wp), intent(in), dimension(counts) :: my_vector   ! Input vector on each process
         integer, intent(in) :: root               ! Rank of the root process
-        real(kind(0d0)), allocatable, intent(out) :: gathered_vector(:) ! Gathered vector on the root process
+        real(wp), allocatable, intent(out) :: gathered_vector(:) ! Gathered vector on the root process
 
         integer :: i, offset, ierr
         integer, allocatable :: recounts(:), displs(:)
diff --git a/src/common/m_phase_change.fpp b/src/common/m_phase_change.fpp
@@ -85,11 +85,11 @@ contains
     subroutine s_infinite_relaxation_k(q_cons_vf)
 
         type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
-        real(kind(0.0_wp)) :: pS, pSOV, pSSL !< equilibrium pressure for mixture, overheated vapor, and subcooled liquid
-        real(kind(0.0_wp)) :: TS, TSOV, TSSL, TSatOV, TSatSL !< equilibrium temperature for mixture, overheated vapor, and subcooled liquid. Saturation Temperatures at overheated vapor and subcooled liquid
-        real(kind(0.0_wp)) :: rhoe, dynE, rhos !< total internal energy, kinetic energy, and total entropy
-        real(kind(0.0_wp)) :: rho, rM, m1, m2, MCT !< total density, total reacting mass, individual reacting masses
-        real(kind(0.0_wp)) :: TvF !< total volume fraction
+        real(wp) :: pS, pSOV, pSSL !< equilibrium pressure for mixture, overheated vapor, and subcooled liquid
+        real(wp) :: TS, TSOV, TSSL, TSatOV, TSatSL !< equilibrium temperature for mixture, overheated vapor, and subcooled liquid. Saturation Temperatures at overheated vapor and subcooled liquid
+        real(wp) :: rhoe, dynE, rhos !< total internal energy, kinetic energy, and total entropy
+        real(wp) :: rho, rM, m1, m2, MCT !< total density, total reacting mass, individual reacting masses
+        real(wp) :: TvF !< total volume fraction
 
         !$acc declare create(pS, pSOV, pSSL, TS, TSOV, TSatOV, TSatSL, TSSL, rhoe, dynE, rhos, rho, rM, m1, m2, MCT, TvF)
 
@@ -295,13 +295,13 @@ contains
 
         ! initializing variables
         integer, intent(in) :: j, k, l, MFL
-        real(kind(0.0_wp)), intent(out) :: pS
-        real(kind(0.0_wp)), dimension(num_fluids), intent(out) :: p_infpT
-        real(kind(0.0_wp)), intent(in) :: rM
+        real(wp), intent(out) :: pS
+        real(wp), dimension(num_fluids), intent(out) :: p_infpT
+        real(wp), intent(in) :: rM
         type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
-        real(kind(0.0_wp)), intent(in) :: rhoe
-        real(kind(0.0_wp)), intent(out) :: TS
-        real(kind(0.0_wp)) :: gp, gpp, hp, pO, mCP, mQ !< variables for the Newton Solver
+        real(wp), intent(in) :: rhoe
+        real(wp), intent(out) :: TS
+        real(wp) :: gp, gpp, hp, pO, mCP, mQ !< variables for the Newton Solver
 
         integer :: i, ns !< generic loop iterators
 
@@ -399,17 +399,17 @@ contains
 #endif
 
         integer, intent(in) :: j, k, l
-        real(kind(0.0_wp)), intent(inout) :: pS
-        real(kind(0.0_wp)), dimension(num_fluids), intent(in) :: p_infpT
-        real(kind(0.0_wp)), intent(in) :: rhoe
+        real(wp), intent(inout) :: pS
+        real(wp), dimension(num_fluids), intent(in) :: p_infpT
+        real(wp), intent(in) :: rhoe
         type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
-        real(kind(0.0_wp)), intent(inout) :: TS
+        real(wp), intent(inout) :: TS
 
-        real(kind(0.0_wp)), dimension(num_fluids) :: p_infpTg !< stiffness for the participating fluids for pTg-equilibrium
-        real(kind(0.0_wp)), dimension(2, 2) :: Jac, InvJac, TJac !< matrices for the Newton Solver
-        real(kind(0.0_wp)), dimension(2) :: R2D, DeltamP !< residual and correction array
-        real(kind(0.0_wp)) :: Om ! underrelaxation factor
-        real(kind(0.0_wp)) :: mCP, mCPD, mCVGP, mCVGP2, mQ, mQD ! auxiliary variables for the pTg-solver
+        real(wp), dimension(num_fluids) :: p_infpTg !< stiffness for the participating fluids for pTg-equilibrium
+        real(wp), dimension(2, 2) :: Jac, InvJac, TJac !< matrices for the Newton Solver
+        real(wp), dimension(2) :: R2D, DeltamP !< residual and correction array
+        real(wp) :: Om ! underrelaxation factor
+        real(wp) :: mCP, mCPD, mCVGP, mCVGP2, mQ, mQD ! auxiliary variables for the pTg-solver
 
         !< Generic loop iterators
         integer :: i, ns
@@ -525,9 +525,9 @@ contains
 
         !> @name variables for the correction of the reacting partial densities
         !> @{
-        real(kind(0.0_wp)), intent(out) :: MCT
+        real(wp), intent(out) :: MCT
         type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
-        real(kind(0.0_wp)), intent(inout) :: rM
+        real(wp), intent(inout) :: rM
         integer, intent(in) :: j, k, l
         !> @}
         if (rM < 0.0_wp) then
@@ -586,15 +586,15 @@ contains
         !$acc routine seq
 #endif
 
-        real(kind(0.0_wp)), dimension(2, 2), intent(out) :: InvJac
+        real(wp), dimension(2, 2), intent(out) :: InvJac
         integer, intent(in) :: j
-        real(kind(0.0_wp)), dimension(2, 2), intent(out) :: Jac
+        real(wp), dimension(2, 2), intent(out) :: Jac
         integer, intent(in) :: k, l
-        real(kind(0.0_wp)), intent(in) :: mCPD, mCVGP, mCVGP2, pS
+        real(wp), intent(in) :: mCPD, mCVGP, mCVGP2, pS
         type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
-        real(kind(0.0_wp)), dimension(2, 2), intent(out) :: TJac
+        real(wp), dimension(2, 2), intent(out) :: TJac
 
-        real(kind(0.0_wp)) :: ml, mT, TS, dFdT, dTdm, dTdp ! mass of the reacting fluid, total reacting mass, and auxiliary variables
+        real(wp) :: ml, mT, TS, dFdT, dTdm, dTdp ! mass of the reacting fluid, total reacting mass, and auxiliary variables
 
         ! mass of the reacting liquid
         ml = q_cons_vf(lp + contxb - 1)%sf(j, k, l)
@@ -694,12 +694,12 @@ contains
 #endif
 
         integer, intent(in) :: j, k, l
-        real(kind(0.0_wp)), intent(in) :: mCPD, mCVGP, mQD
+        real(wp), intent(in) :: mCPD, mCVGP, mQD
         type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
-        real(kind(0.0_wp)), intent(in) :: pS, rhoe
-        real(kind(0.0_wp)), dimension(2), intent(out) :: R2D
+        real(wp), intent(in) :: pS, rhoe
+        real(wp), dimension(2), intent(out) :: R2D
 
-        real(kind(0.0_wp)) :: ml, mT, TS !< mass of the reacting liquid, total reacting mass, equilibrium temperature
+        real(wp) :: ml, mT, TS !< mass of the reacting liquid, total reacting mass, equilibrium temperature
 
         ! mass of the reacting liquid
         ml = q_cons_vf(lp + contxb - 1)%sf(j, k, l)
@@ -744,11 +744,11 @@ contains
         !$acc routine seq
 #endif
 
-        real(kind(0.0_wp)), intent(in) :: pSat
-        real(kind(0.0_wp)), intent(out) :: TSat
-        real(kind(0.0_wp)), intent(in) :: TSIn
+        real(wp), intent(in) :: pSat
+        real(wp), intent(out) :: TSat
+        real(wp), intent(in) :: TSIn
 
-        real(kind(0.0_wp)) :: dFdT, FT, Om !< auxiliary variables
+        real(wp) :: dFdT, FT, Om !< auxiliary variables
 
         ! Generic loop iterators
         integer :: ns
diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp
@@ -1038,7 +1038,7 @@ contains
         allocate (x_d1(m*n))
         allocate (y_d1(m*n))
         counter = 0
-        maxalph_loc = 0d0
+        maxalph_loc = 0_wp
         do k = 0, p
             do j = 0, n
                 do i = 0, m
@@ -1060,7 +1060,7 @@ contains
             cent = 0
         end if
 
-        thres = 0.9d0*maxalph_glb
+        thres = 0.9_wp*maxalph_glb
         do k = 0, n
             do j = 0, m
                 axp = q_prim_vf(E_idx + 2)%sf(j + 1, k, cent)
@@ -1110,7 +1110,7 @@ contains
                         x_td(i), y_td(i), size(x_td)
                 else
                     write (211, '(F12.9,1X,F12.9,1X,F3.1)') &
-                        x_td(i), y_td(i), 0d0
+                        x_td(i), y_td(i), 0_wp
                 end if
             end do
         end if
@@ -1119,40 +1119,40 @@ contains
 
     subroutine s_write_energy_data_file(q_prim_vf, q_cons_vf)
         type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf, q_cons_vf
-        real(kind(0d0)) :: Elk, Egk, Elp, Egint, Vb, Vl, pres_av, Et
-        real(kind(0d0)) :: rho, pres, dV, tmp, gamma, pi_inf, MaxMa, MaxMa_glb, maxvel, c, Ma, H
-        real(kind(0d0)), dimension(num_dims) :: vel
-        real(kind(0d0)), dimension(num_fluids) :: gammas, pi_infs, adv
+        real(wp) :: Elk, Egk, Elp, Egint, Vb, Vl, pres_av, Et
+        real(wp) :: rho, pres, dV, tmp, gamma, pi_inf, MaxMa, MaxMa_glb, maxvel, c, Ma, H
+        real(wp), dimension(num_dims) :: vel
+        real(wp), dimension(num_fluids) :: gammas, pi_infs, adv
         integer :: i, j, k, l, s !looping indices
         integer :: ierr, counter, root !< number of data points extracted to fit shape to SH perturbations
 
-        Egk = 0d0
-        Elp = 0d0
-        Egint = 0d0
-        Vb = 0d0
-        maxvel = 0d0
-        MaxMa = 0d0
-        Vl = 0d0
-        Elk = 0d0
-        Et = 0d0
-        Vb = 0d0
-        dV = 0d0
-        pres_av = 0d0
-        pres = 0d0
+        Egk = 0_wp
+        Elp = 0_wp
+        Egint = 0_wp
+        Vb = 0_wp
+        maxvel = 0_wp
+        MaxMa = 0_wp
+        Vl = 0_wp
+        Elk = 0_wp
+        Et = 0_wp
+        Vb = 0_wp
+        dV = 0_wp
+        pres_av = 0_wp
+        pres = 0_wp
         do k = 0, p
             do j = 0, n
                 do i = 0, m
-                    pres = 0d0
+                    pres = 0_wp
                     dV = dx(i)*dy(j)*dz(k)
-                    rho = 0d0
-                    gamma = 0d0
-                    pi_inf = 0d0
+                    rho = 0_wp
+                    gamma = 0_wp
+                    pi_inf = 0_wp
                     pres = q_prim_vf(E_idx)%sf(i, j, k)
                     Egint = Egint + q_prim_vf(E_idx + 2)%sf(i, j, k)*(fluid_pp(2)%gamma*pres)*dV
                     do s = 1, num_dims
                         vel(s) = q_prim_vf(num_fluids + s)%sf(i, j, k)
-                        Egk = Egk + 0.5d0*q_prim_vf(E_idx + 2)%sf(i, j, k)*q_prim_vf(2)%sf(i, j, k)*vel(s)*vel(s)*dV
-                        Elk = Elk + 0.5d0*q_prim_vf(E_idx + 1)%sf(i, j, k)*q_prim_vf(1)%sf(i, j, k)*vel(s)*vel(s)*dV
+                        Egk = Egk + 0.5_wp*q_prim_vf(E_idx + 2)%sf(i, j, k)*q_prim_vf(2)%sf(i, j, k)*vel(s)*vel(s)*dV
+                        Elk = Elk + 0.5_wp*q_prim_vf(E_idx + 1)%sf(i, j, k)*q_prim_vf(1)%sf(i, j, k)*vel(s)*vel(s)*dV
                         if (abs(vel(s)) > maxvel) then
                             maxvel = abs(vel(s))
                         end if
@@ -1164,14 +1164,14 @@ contains
                         rho = rho + adv(l)*q_prim_vf(l)%sf(i, j, k)
                     end do
 
-                    H = ((gamma + 1d0)*pres + pi_inf)/rho
+                    H = ((gamma + 1_wp)*pres + pi_inf)/rho
 
                     call s_compute_speed_of_sound(pres, rho, &
                                                   gamma, pi_inf, &
-                                                  H, adv, 0d0, 0d0, c)
+                                                  H, adv, 0_wp, 0_wp, c)
 
                     Ma = maxvel/c
-                    if (Ma > MaxMa .and. adv(1) > 1.0d0 - 1.0d-10) then
+                    if (Ma > MaxMa .and. adv(1) > 1.0_wp - 1.0d-10) then
                         MaxMa = Ma
                     end if
                     Vl = Vl + adv(1)*dV
diff --git a/src/pre_process/m_assign_variables.fpp b/src/pre_process/m_assign_variables.fpp
@@ -507,7 +507,7 @@ contains
                 theta = atan2(y_cc(k), x_cc(j))
                 phi = atan2(sqrt(x_cc(j)**2 + y_cc(k)**2), z_cc(l))
                 !spherical coord, assuming Rmax=1
-                xi_sph = (rcoord**3 - R0ref**3 + 1d0)**(1d0/3d0)
+                xi_sph = (rcoord**3 - R0ref**3 + 1_wp)**(1_wp/3_wp)
                 xi_cart(1) = xi_sph*sin(phi)*cos(theta)
                 xi_cart(2) = xi_sph*sin(phi)*sin(theta)
                 xi_cart(3) = xi_sph*cos(phi)
@@ -520,7 +520,7 @@ contains
             ! assigning the reference map to the q_prim vector field
             do i = 1, num_dims
                 q_prim_vf(i + xibeg - 1)%sf(j, k, l) = eta*xi_cart(i) + &
-                                                       (1d0 - eta)*orig_prim_vf(i + xibeg - 1)
+                                                       (1_wp - eta)*orig_prim_vf(i + xibeg - 1)
             end do
 
         end if
diff --git a/src/pre_process/m_patches.fpp b/src/pre_process/m_patches.fpp
@@ -1412,9 +1412,9 @@ contains
         integer, intent(INOUT), dimension(0:m, 0:n, 0:p) :: patch_id_fp
         type(scalar_field), dimension(1:sys_size) :: q_prim_vf
 
-        real(kind(0d0)) :: r, x_p, eps, phi
-        real(kind(0d0)), dimension(2:9) :: as, Ps
-        real(kind(0d0)) :: radius, x_centroid, y_centroid, z_centroid, eta, smooth_coeff
+        real(wp) :: r, x_p, eps, phi
+        real(wp), dimension(2:9) :: as, Ps
+        real(wp) :: radius, x_centroid, y_centroid, z_centroid, eta, smooth_coeff
         logical :: non_axis_sym
 
         integer :: i, j, k !< generic loop iterators
@@ -1485,7 +1485,7 @@ contains
                             if (patch_icpp(patch_id)%smoothen) then
                                 eta = tanh(smooth_coeff/min(dx, dy, dz)* &
                                            ((r - as(2)*Ps(2) - as(3)*Ps(3) - as(4)*Ps(4) - as(5)*Ps(5) - as(6)*Ps(6) - as(7)*Ps(7)) &
-                                            - radius))*(-0.5d0) + 0.5d0
+                                            - radius))*(-0.5_wp) + 0.5_wp
                             end if
 
                             call s_assign_patch_primitive_variables(patch_id, i, j, k, &
@@ -1502,7 +1502,7 @@ contains
 
                     if (non_axis_sym) then
                         phi = atan(((y_cc(j) - y_centroid) + eps)/((x_cc(i) - x_centroid) + eps))
-                        r = dsqrt((x_cc(i) - x_centroid)**2d0 + (y_cc(j) - y_centroid)**2d0) + eps
+                        r = dsqrt((x_cc(i) - x_centroid)**2_wp + (y_cc(j) - y_centroid)**2_wp) + eps
                         x_p = (eps)/r
                         Ps(2) = spherical_harmonic_func(x_p, phi, 2, 2)
                         Ps(3) = spherical_harmonic_func(x_p, phi, 3, 3)
@@ -1513,7 +1513,7 @@ contains
                         Ps(8) = spherical_harmonic_func(x_p, phi, 8, 8)
                         Ps(9) = spherical_harmonic_func(x_p, phi, 9, 9)
                     else
-                        r = dsqrt((x_cc(i) - x_centroid)**2d0 + (y_cc(j) - y_centroid)**2d0) + eps
+                        r = dsqrt((x_cc(i) - x_centroid)**2_wp + (y_cc(j) - y_centroid)**2_wp) + eps
                         x_p = abs(x_cc(i) - x_centroid + eps)/r
                         Ps(2) = unassociated_legendre(x_p, 2)
                         Ps(3) = unassociated_legendre(x_p, 3)
diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp
@@ -932,15 +932,15 @@ contains
     subroutine s_write_com_files(t_step, c_mass) ! -------------------
 
         integer, intent(in) :: t_step
-        real(kind(0d0)), dimension(num_fluids, 5), intent(in) :: c_mass
+        real(wp), dimension(num_fluids, 5), intent(in) :: c_mass
         integer :: i, j !< Generic loop iterator
-        real(kind(0d0)) :: nondim_time !< Non-dimensional time
+        real(wp) :: nondim_time !< Non-dimensional time
 
         ! Non-dimensional time calculation
         if (t_step_old /= dflt_int) then
-            nondim_time = real(t_step + t_step_old, kind(0d0))*dt
+            nondim_time = real(t_step + t_step_old, wp)*dt
         else
-            nondim_time = real(t_step, kind(0d0))*dt
+            nondim_time = real(t_step, wp)*dt
         end if
 
         if (proc_rank == 0) then
diff --git a/src/simulation/m_derived_variables.f90 b/src/simulation/m_derived_variables.f90
@@ -326,14 +326,14 @@ end subroutine s_derive_acceleration_component
     !!  @param c_m Mass,x-location,y-location,z-location
     subroutine s_derive_center_of_mass(q_vf, c_m)
         type(scalar_field), dimension(sys_size), intent(IN) :: q_vf
-        real(kind(0d0)), dimension(1:num_fluids, 1:5), intent(INOUT) :: c_m
+        real(wp), dimension(1:num_fluids, 1:5), intent(INOUT) :: c_m
         integer :: i, j, k, l !< Generic loop iterators
-        real(kind(0d0)) :: tmp, tmp_out !< Temporary variable to store quantity for mpi_allreduce
-        real(kind(0d0)) :: dV !< Discrete cell volume
+        real(wp) :: tmp, tmp_out !< Temporary variable to store quantity for mpi_allreduce
+        real(wp) :: dV !< Discrete cell volume
 
         do i = 1, num_fluids
             do j = 1, 5
-                c_m(i, j) = 0.0d0
+                c_m(i, j) = 0.0_wp
             end do
         end do
 
diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp
diff --git a/src/simulation/m_hyperelastic.fpp b/src/simulation/m_hyperelastic.fpp
diff --git a/src/simulation/m_hypoelastic.fpp b/src/simulation/m_hypoelastic.fpp
diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp
