changes to several fils

so they can become more readable

Author: JRChreim

src/common/m_checker_common.fpp

@@ -253,7 +253,7 @@ contains
         @: PROHIBIT(cyl_coord .and. n == 0, "n must be positive (2D or 3D) for cylindrical coordinates")
         @: PROHIBIT(sph_coord .and. n == 0, "n must be positive (1D) for axysimmetric spherical coordinates")
         @: PROHIBIT(cyl_coord .and. p == 0 .and. bc_y%beg /= -2, "bc_y%beg must be -2 for 2D cylindrical coordinates (p = 0)")
-!        @: PROHIBIT(sph_coord .and. p == 0 .and. bc_y%beg /= -2, "bc_y%beg must be -2 axysimmetric spherical coordinates (p = 0)")
+        @: PROHIBIT(sph_coord .and. p == 0 .and. bc_y%beg /= -2, "bc_y%beg must be -2 axysimmetric spherical coordinates (p = 0)")
         @: PROHIBIT(cyl_coord .and. p > 0 .and. bc_y%beg /= -14, "bc_y%beg must be -14 for 3D cylindrical coordinates (p > 0)")
         @: PROHIBIT(cyl_coord .and. (bc_y%end > -1 .or. bc_y%end < -16), "bc_y%end must be between -1 and -16")
         @: PROHIBIT(sph_coord .and. (bc_y%end > -1 .or. bc_y%end < -16), "bc_y%end must be between -1 and -16")

src/common/m_phase_change.fpp

@@ -110,17 +110,16 @@ contains
                         ! original volume fractions, before any relaxation
                         alpha0k(i) = q_cons_vf(i + advxb - 1)%sf(j, k, l)
 
-                        ! original mass fractions, before any relaxation
+                        ! original partial density, before any relaxation
                         m0k(i) = q_cons_vf(i + contxb - 1)%sf(j, k, l)
 
                         ! Mixture density
-                        rho = rho + q_cons_vf(i + contxb - 1)%sf(j, k, l)
+                        rho = rho + m0k(i)
 
                         ! calculating the total internal energy such that the energy-fraction for each of the
                         ! fluids can be proportionally distributed when the sum of the internal energies differs from the
-
-                        if ( ieee_is_nan( q_cons_vf(i + contxb - 1)%sf(j, k, l) ) ) then
-                            print *, 'partial densities', q_cons_vf(i + contxb - 1)%sf(j, k, l)
+                        if ( ieee_is_nan( m0k(i) ) ) then
+                            print *, 'partial densities', m0k(i)
                         end if
 
                         if (model_eqns == 3) then
@@ -177,7 +176,7 @@ contains
                             ! if ( ( pS < pCr ) .and. ( pS > 0 ) .and. &
                             if ( ( pS < pCr ) .and. &
                             ! 2.1 Homogeneous pTg-equilibrium criterium
-                            ( ( ( pS < 0 ) .and. ( pS + minval(p_infpT) > 0.0 ) ) &
+                            ( ( ( pS < 0 ) .and. ( pS + minval(p_infpT) > 0.0_wp ) ) &
                             .or. &
                             ! 2.2. Heterogeneous pTg-equilibrium.
                             ( (alphak(lp) > palpha_eps) .and. (alphak(vp) > palpha_eps) ) ) &
@@ -278,7 +277,6 @@ contains
 
                         ! updating conservative variables through either p- or pT-equilibrium
                         call update_conservative_vars( j, k, l, pS, q_cons_vf, Tk )
-
                     else
                         !$acc loop seq
                         do i = 1, num_fluids

src/simulation/m_rhs.fpp

@@ -624,7 +624,7 @@ contains
         real(wp), intent(inout) :: time_avg
 
         real(wp), dimension(0:m, 0:n, 0:p) :: nbub
-        real(wp) :: t_start, t_finish
+        real(wp) :: t_start, t_finish, multip
         integer :: i, j, k, l, id !< Generic loop iterators
 
         call nvtxStartRange("COMPUTE-RHS")
@@ -644,7 +644,14 @@ contains
             end do
         end do
 
-        ! Converting Conservative to Primitive Variables
+        ! determining what the multiplying factor will be for the source terms in the cylindrical/spherical coordinates
+        if (cyl_coord) then
+            multip = 5e-1_wp
+        else if (sph_coord) then
+            multip = 1.0_wp
+        else 
+            multip = 0.0_wp
+        end if
 
         if (mpp_lim .and. bubbles_euler) then
             !$acc parallel loop collapse(3) gang vector default(present)
@@ -817,7 +824,8 @@ contains
                                                  iden_cons_qp, &
                                                  iden_prim_qp, &
                                                  nill_cons_qp, &
-                                                 nill_prim_qp)
+                                                 nill_prim_qp, &
+                                                 multip)
             call nvtxEndRange
 
             ! RHS additions for hypoelasticity
@@ -836,7 +844,8 @@ contains
                                                       flux_src_n(id)%vf, &
                                                       dq_prim_dx_qp(1)%vf, &
                                                       dq_prim_dy_qp(1)%vf, &
-                                                      dq_prim_dz_qp(1)%vf)
+                                                      dq_prim_dz_qp(1)%vf, &
+                                                      multip)
                 call nvtxEndRange
             end if
 
@@ -937,29 +946,20 @@ contains
 
     end subroutine s_compute_rhs
 
-    subroutine s_compute_advection_source_term(idir, rhs_vf, q_cons_vf, q_prim_vf, flux_src_n_vf, iden_cons_vf, iden_prim_vf, nill_cons_vf, nill_prim_vf)
+    subroutine s_compute_advection_source_term(idir, rhs_vf, q_cons_vf, q_prim_vf, flux_src_n_vf, iden_cons_vf, iden_prim_vf, nill_cons_vf, nill_prim_vf, multip)
 
         integer, intent(in) :: idir
         type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf
         type(vector_field), intent(inout) :: q_cons_vf
         type(vector_field), intent(inout) :: q_prim_vf
         type(vector_field), intent(inout) :: flux_src_n_vf
         type(vector_field), intent(in) :: iden_cons_vf, iden_prim_vf, nill_cons_vf, nill_prim_vf
-        real(wp) :: multip
+        real(wp), intent(in) :: multip
 
         real(wp), dimension(0:m, 0:n, 0:p) :: KNeg, KNill
 
         integer :: i, j, k, l, q
 
-        ! determining what the multiplying factor will be for the source terms in the cylindrical/spherical coordinates
-        if (cyl_coord) then
-            multip = 5e-1_wp
-        else if (sph_coord) then
-            multip = 1.0_wp
-        else 
-            multip = 0.0_wp
-        end if        
-
         if (alt_soundspeed) then
             !$acc parallel loop collapse(3) gang vector default(present)
             do l = 0, p
@@ -979,15 +979,14 @@ contains
 
                         Kterm(j, k, l) = alpha1(j, k, l)*alpha2(j, k, l)*(blkmod2(j, k, l) - blkmod1(j, k, l))/ &
                                          (alpha1(j, k, l)*blkmod2(j, k, l) + alpha2(j, k, l)*blkmod1(j, k, l))
-        
-                        KNill(j, k, l) = 0.0_wp
-
-                        KNeg(j, k, l) = -1.0_wp * Kterm(j, k, l)
                     end do
                 end do
             end do
+            KNeg(0:m, 0:n, 0:p) = -1.0_wp * Kterm(0:m, 0:n, 0:p)
         end if
-        
+
+        KNill(0:m, 0:n, 0:p) = 0.0_wp
+
         if (idir == 1) then
 
             if (bc_x%beg <= -5 .and. bc_x%beg >= -13) then
@@ -1034,7 +1033,6 @@ contains
             end if
 
         elseif (idir == 2) then
-
             ! RHS Contribution in y-direction
             ! Applying the Riemann fluxes
 
@@ -1084,7 +1082,6 @@ contains
             end if
 
             if (grid_geometry == 3) then ! Cylindrical Coordinates
-
                 !$acc parallel loop collapse(4) gang vector default(present)
                 do j = 1, sys_size
                     do k = 0, p
@@ -1150,6 +1147,7 @@ contains
                     ! call s_calculate_rhs_operator(idir, 1, sys_size, 0, m, 0, n, 0, p, multip, rhs_vf, flux_n, flux_gsrc_n, &
                     !                                 KNill, iden_cons_vf, iden_prim_vf, 0, 1, KNill, iden_cons_vf, iden_prim_vf, 0, 1, 'L')
                 else
+
                     if (alt_soundspeed) then
                         do j = advxb, advxe
                             if ((j == advxe) .and. (bubbles_euler .neqv. .true.)) then
@@ -1244,24 +1242,17 @@ contains
     end subroutine s_compute_advection_source_term
 
     subroutine s_compute_additional_physics_rhs(idir, q_prim_vf, rhs_vf, flux_src_n, &
-                                                dq_prim_dx_vf, dq_prim_dy_vf, dq_prim_dz_vf)
+                                                dq_prim_dx_vf, dq_prim_dy_vf, dq_prim_dz_vf, multip)
 
         integer, intent(in) :: idir
         type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
         type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf
         type(scalar_field), dimension(sys_size), intent(in) :: flux_src_n
         type(scalar_field), dimension(sys_size), intent(in) :: dq_prim_dx_vf, dq_prim_dy_vf, dq_prim_dz_vf
-        real(wp) :: multip
+        real(wp), intent(in) :: multip
 
         integer :: i, j, k, l
 
-        ! determining what the multiplying factor will be for the source terms in the cylindrical/spherical coordinates
-        if (cyl_coord) then
-            multip = 5e-1_wp
-        else if (sph_coord) then
-            multip = 1.0_wp
-        end if
-
         if (idir == 1) then ! x-direction
 
             if (surface_tension) then

src/simulation/m_time_steppers.fpp

@@ -394,7 +394,7 @@ contains
         if (grid_geometry == 3) call s_apply_fourier_filter(q_cons_ts(1)%vf)
 
         if (model_eqns == 3 .and. (.not. relax)) then
-            call s_pressure_relaxation_procedure(q_cons_ts(2)%vf)
+            call s_pressure_relaxation_procedure(q_cons_ts(1)%vf)
         end if
 
         if (adv_n) call s_comp_alpha_from_n(q_cons_ts(1)%vf)