Merge branch 'main' into igr-integration

StevenXue98 · StevenXue98 · commit 20033cc059d6 · 2026-06-23T00:20:49.000-04:00
diff --git a/src/common/include/1dHardcodedIC.fpp b/src/common/include/1dHardcodedIC.fpp
@@ -51,9 +51,7 @@
         molar_mass_inv = y1/31.998_wp + y2/18.01508_wp + y3/16.04256_wp + y4/28.0134_wp
 
         q_prim_vf(eqn_idx%cont%beg)%sf(i, 0, 0) = 1.01325_wp*(10.0_wp)**5/(temp*8.3144626_wp*1000.0_wp*molar_mass_inv)
-
     case(191)  ! 1D Dual Isothermal case
-
         q_prim_vf(eqn_idx%E)%sf(i, 0, 0) = 101325.0_wp
         q_prim_vf(eqn_idx%mom%beg)%sf(i, 0, 0) = 0.0_wp
         q_prim_vf(eqn_idx%species%beg)%sf(i, 0, 0) = 1.0_wp
diff --git a/src/common/include/2dHardcodedIC.fpp b/src/common/include/2dHardcodedIC.fpp
@@ -251,7 +251,7 @@
             ! q_prim_vf(eqn_idx%B%beg)%sf(i,j,0) = 1._wp/(4._wp*pi) * (alpha**8 - 2._wp*alpha**4 + 1._wp)
             ! q_prim_vf(eqn_idx%E)%sf(i,j,0) = 6._wp - q_prim_vf(eqn_idx%B%beg)%sf(i,j,0)**2/2._wp
         end if
-    case (262)  ! Tilted 2D MHD shock‐tube at α = arctan2 (≈63.4°)
+    case (262)  ! Tilted 2D MHD shock-tube at \alpha = arctan2 (\approx63.4 deg)
         ! rotate by \alpha = atan(2)
         alpha = atan(2._wp)
         cosA = cos(alpha)
diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp
@@ -278,9 +278,10 @@ module m_derived_types
 
         ! Geometry 13 (2D modal Fourier): fourier_cos(n), fourier_sin(n) for mode n
         real(wp), dimension(1:max_2d_fourier_modes) :: fourier_cos, fourier_sin
-        logical :: modal_clip_r_to_min  !< When true, clip boundary radius: R(theta) = max(R(theta), modal_r_min) (Non-exp form only)
-        real(wp) :: modal_r_min  !< Minimum boundary radius when modal_clip_r_to_min is true (Non-exp form only)
-        logical :: modal_use_exp_form  !< When true, boundary = radius*exp(Fourier series)
+        !> When true, clip boundary radius: R(theta) = max(R(theta), modal_r_min) (Non-exp form only)
+        logical  :: modal_clip_r_to_min
+        real(wp) :: modal_r_min         !< Minimum boundary radius when modal_clip_r_to_min is true (Non-exp form only)
+        logical  :: modal_use_exp_form  !< When true, boundary = radius*exp(Fourier series)
         ! Geometry 14 (3D spherical harmonic): sph_har_coeff(l,m) for real Y_lm
         real(wp), dimension(0:max_sph_harm_degree,-max_sph_harm_degree:max_sph_harm_degree) :: sph_har_coeff
         real(wp), dimension(3) :: normal  !< Patch orientation normal vector (x, y, z)
diff --git a/src/common/m_helper_basic.fpp b/src/common/m_helper_basic.fpp
@@ -111,7 +111,6 @@ contains
 
     !> Compute ghost-cell buffer size and set interior/buffered coordinate index bounds.
     subroutine s_configure_coordinate_bounds(recon_type, weno_polyn, muscl_polyn, igr_order, buff_size, idwint, idwbuff, viscous, &
-
         & bubbles_lagrange, m, n, p, num_dims, igr, ib)
 
         integer, intent(in)                                :: recon_type, weno_polyn, muscl_polyn
diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp
@@ -291,7 +291,6 @@ contains
     !! single process, within its assigned section of the computational domain. Finally, note that the global extrema values are
     !! only bookkeept on the rank 0 processor.
     impure subroutine s_mpi_reduce_stability_criteria_extrema(icfl_max_loc, vcfl_max_loc, Rc_min_loc, icfl_max_glb, vcfl_max_glb, &
-
         & Rc_min_glb)
 
         real(wp), intent(in)  :: icfl_max_loc
@@ -367,13 +366,13 @@ contains
     !> Reduce a local integer value to its global sum across all MPI ranks.
     impure subroutine s_mpi_allreduce_integer_sum(var_loc, var_glb)
 
-        integer, intent(in)  :: var_loc
-        integer, intent(out) :: var_glb
+        integer(kind=8), intent(in)  :: var_loc
+        integer(kind=8), intent(out) :: var_glb
 
 #ifdef MFC_MPI
         integer :: ierr  !< Generic flag used to identify and report MPI errors
 
-        call MPI_ALLREDUCE(var_loc, var_glb, 1, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierr)
+        call MPI_ALLREDUCE(var_loc, var_glb, 1, MPI_INTEGER8, MPI_SUM, MPI_COMM_WORLD, ierr)
 #else
         var_glb = var_loc
 #endif
@@ -1444,7 +1443,6 @@ contains
 
         if (mpi_dir == 1) then
             if (pbc_loc == -1) then  ! PBC at the beginning
-
                 if (bc_x%end >= 0) then  ! PBC at the beginning and end
                     call MPI_SENDRECV(dx(m - buff_size + 1), buff_size, mpi_p, bc_x%end, 0, dx(-buff_size), buff_size, mpi_p, &
                                       & bc_x%beg, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr)
@@ -1463,7 +1461,6 @@ contains
             end if
         else if (mpi_dir == 2) then
             if (pbc_loc == -1) then  ! PBC at the beginning
-
                 if (bc_y%end >= 0) then  ! PBC at the beginning and end
                     call MPI_SENDRECV(dy(n - buff_size + 1), buff_size, mpi_p, bc_y%end, 0, dy(-buff_size), buff_size, mpi_p, &
                                       & bc_y%beg, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr)
@@ -1482,7 +1479,6 @@ contains
             end if
         else
             if (pbc_loc == -1) then  ! PBC at the beginning
-
                 if (bc_z%end >= 0) then  ! PBC at the beginning and end
                     call MPI_SENDRECV(dz(p - buff_size + 1), buff_size, mpi_p, bc_z%end, 0, dz(-buff_size), buff_size, mpi_p, &
                                       & bc_z%beg, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr)
diff --git a/src/simulation/m_body_forces.fpp b/src/simulation/m_body_forces.fpp
@@ -99,7 +99,6 @@ contains
         $:END_GPU_PARALLEL_LOOP()
 
         if (bf_x) then  ! x-direction body forces
-
             $:GPU_PARALLEL_LOOP(private='[j, k, l]', collapse=3)
             do l = 0, p
                 do k = 0, n
@@ -114,7 +113,6 @@ contains
         end if
 
         if (bf_y) then  ! y-direction body forces
-
             $:GPU_PARALLEL_LOOP(private='[j, k, l]', collapse=3)
             do l = 0, p
                 do k = 0, n
@@ -130,7 +128,6 @@ contains
         end if
 
         if (bf_z) then  ! z-direction body forces
-
             $:GPU_PARALLEL_LOOP(private='[j, k, l]', collapse=3)
             do l = 0, p
                 do k = 0, n
diff --git a/src/simulation/m_bubbles.fpp b/src/simulation/m_bubbles.fpp
@@ -333,8 +333,7 @@ contains
     !> Adaptive time stepping routine for subgrid bubbles (See Heirer, E. Hairer S.P.Norsett G. Wanner, Solving Ordinary
     !! Differential Equations I, Chapter II.4)
     subroutine s_advance_step(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, fntait, fBtait, f_bub_adv_src, f_divu, bub_id, fmass_v, &
-
-        & fmass_g, fbeta_c, fbeta_t, fCson, adap_dt_stop)
+                              & fmass_g, fbeta_c, fbeta_t, fCson, adap_dt_stop)
         $:GPU_ROUTINE(function_name='s_advance_step',parallelism='[seq]', cray_inline=True)
 
         real(wp), intent(inout) :: fR, fV, fpb, fmass_v
@@ -488,8 +487,7 @@ contains
 
     !> Integrate bubble variables over the given time step size, h, using a third-order accurate embedded Runge-Kutta scheme.
     subroutine s_advance_substep(err, fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, fntait, fBtait, f_bub_adv_src, f_divu, bub_id, &
-
-        & fmass_v, fmass_g, fbeta_c, fbeta_t, fCson, h, myR_tmp, myV_tmp, myPb_tmp, myMv_tmp)
+                                 & fmass_v, fmass_g, fbeta_c, fbeta_t, fCson, h, myR_tmp, myV_tmp, myPb_tmp, myMv_tmp)
         $:GPU_ROUTINE(function_name='s_advance_substep',parallelism='[seq]', cray_inline=True)
 
         real(wp), intent(out)               :: err
diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp
@@ -72,8 +72,8 @@ contains
     !> Initializes the values of various IBM variables, such as ghost points and image points.
     impure subroutine s_ibm_setup()
 
-        integer :: i, j, k
-        integer :: max_num_gps
+        integer         :: i, j, k
+        integer(kind=8) :: max_num_gps
 
         call nvtxStartRange("SETUP-IBM-MODULE")
 
@@ -119,10 +119,10 @@ contains
         ! find the number of ghost points and set them to be the maximum total across ranks
         call s_find_num_ghost_points(num_gps)
         if (moving_immersed_boundary_flag) then
-            call s_mpi_allreduce_integer_sum(num_gps, max_num_gps)
-            max_num_gps = min(max_num_gps*2, (m + 1)*(n + 1)*(p + 1))
+            call s_mpi_allreduce_integer_sum(int(num_gps, 8), max_num_gps)
+            max_num_gps = min(max_num_gps*2_8, int(m + 1, 8)*int(n + 1, 8)*int(p + 1, 8))
         else
-            max_num_gps = num_gps
+            max_num_gps = int(num_gps, 8)
         end if
 
         ! set the size of the ghost point arrays to be the amount of points total, plus a factor of 2 buffer
@@ -1071,7 +1071,8 @@ contains
     subroutine s_compute_centroid_offset(ib_marker)
 
         integer, intent(in)      :: ib_marker
-        integer                  :: i, j, k, num_cells, num_cells_local, decoded_gbl_id
+        integer                  :: i, j, k, num_cells_local, decoded_gbl_id
+        integer(kind=8)          :: num_cells
         real(wp), dimension(1:3) :: center_of_mass, center_of_mass_local
 
         ! Offset only needs to be computes for specific geometries
@@ -1098,7 +1099,7 @@ contains
             end do
 
             ! reduce the mass contribution over all MPI ranks and compute COM
-            call s_mpi_allreduce_integer_sum(num_cells_local, num_cells)
+            call s_mpi_allreduce_integer_sum(int(num_cells_local, 8), num_cells)
             if (num_cells /= 0) then
                 call s_mpi_allreduce_sum(center_of_mass_local(1), center_of_mass(1))
                 call s_mpi_allreduce_sum(center_of_mass_local(2), center_of_mass(2))
diff --git a/src/simulation/m_muscl.fpp b/src/simulation/m_muscl.fpp
@@ -84,9 +84,7 @@ contains
     end subroutine s_initialize_muscl_module
 
     !> Perform MUSCL reconstruction of left and right cell-boundary values from cell-averaged variables
-    subroutine s_muscl(v_vf, vL_rs_vf_x, vR_rs_vf_x, muscl_dir, is1_muscl_d, &
-
-        & is2_muscl_d, is3_muscl_d)
+    subroutine s_muscl(v_vf, vL_rs_vf_x, vR_rs_vf_x, muscl_dir, is1_muscl_d, is2_muscl_d, is3_muscl_d)
 
         type(scalar_field), dimension(1:), intent(in) :: v_vf
         real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:), intent(inout) :: vL_rs_vf_x, vR_rs_vf_x
diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp
@@ -458,8 +458,7 @@ contains
 
     !> Compute the right-hand side of the semi-discrete governing equations for a single time stage
     impure subroutine s_compute_rhs(q_cons_vf, q_T_sf, q_prim_vf, bc_type, rhs_vf, pb_in, rhs_pb, mv_in, rhs_mv, t_step, &
-
-        & time_avg, stage)
+                                    & time_avg, stage)
 
         type(scalar_field), dimension(sys_size), intent(inout)                                     :: q_cons_vf
         type(scalar_field), intent(inout)                                                          :: q_T_sf
@@ -1375,7 +1374,6 @@ contains
         integer                                                :: i, j, k, l
 
         if (idir == 1) then  ! x-direction
-
             if (surface_tension) then
                 $:GPU_PARALLEL_LOOP(private='[j, k, l]', collapse=3)
                 do l = 0, p
diff --git a/src/simulation/m_riemann_solver_hll.fpp b/src/simulation/m_riemann_solver_hll.fpp
@@ -25,10 +25,9 @@ module m_riemann_solver_hll
 contains
 
     !> HLL approximate Riemann solver, Harten et al. SIAM Review (1983)
-    subroutine s_hll_riemann_solver(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, &
-
-        & dqL_prim_dz_vf, qL_prim_vf, qR_prim_rsx_vf, dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, qR_prim_vf, q_prim_vf, &
-            & flux_vf, flux_src_vf, flux_gsrc_vf, norm_dir, ix, iy, iz)
+    subroutine s_hll_riemann_solver(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, qL_prim_vf, qR_prim_rsx_vf, &
+                                    & dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, qR_prim_vf, q_prim_vf, flux_vf, &
+                                    & flux_src_vf, flux_gsrc_vf, norm_dir, ix, iy, iz)
 
         real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:), intent(inout) :: qL_prim_rsx_vf, qR_prim_rsx_vf
         type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
diff --git a/src/simulation/m_riemann_solver_hllc.fpp b/src/simulation/m_riemann_solver_hllc.fpp
@@ -28,10 +28,9 @@ module m_riemann_solver_hllc
 contains
 
     !> HLLC Riemann solver with contact restoration, Toro et al. Shock Waves (1994)
-    subroutine s_hllc_riemann_solver(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, &
-
-        & dqL_prim_dz_vf, qL_prim_vf, qR_prim_rsx_vf, dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, qR_prim_vf, q_prim_vf, &
-            & flux_vf, flux_src_vf, flux_gsrc_vf, norm_dir, ix, iy, iz)
+    subroutine s_hllc_riemann_solver(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, qL_prim_vf, qR_prim_rsx_vf, &
+                                     & dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, qR_prim_vf, q_prim_vf, flux_vf, &
+                                     & flux_src_vf, flux_gsrc_vf, norm_dir, ix, iy, iz)
 
         real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:), intent(inout) :: qL_prim_rsx_vf, qR_prim_rsx_vf
         type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
diff --git a/src/simulation/m_riemann_solver_hlld.fpp b/src/simulation/m_riemann_solver_hlld.fpp
@@ -18,10 +18,9 @@ module m_riemann_solver_hlld
 contains
 
     !> HLLD Riemann solver for MHD, Miyoshi & Kusano JCP (2005)
-    subroutine s_hlld_riemann_solver(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, &
-
-        & dqL_prim_dz_vf, qL_prim_vf, qR_prim_rsx_vf, dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, qR_prim_vf, q_prim_vf, &
-            & flux_vf, flux_src_vf, flux_gsrc_vf, norm_dir, ix, iy, iz)
+    subroutine s_hlld_riemann_solver(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, qL_prim_vf, qR_prim_rsx_vf, &
+                                     & dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, qR_prim_vf, q_prim_vf, flux_vf, &
+                                     & flux_src_vf, flux_gsrc_vf, norm_dir, ix, iy, iz)
 
         real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:), intent(inout) :: qL_prim_rsx_vf, qR_prim_rsx_vf
         type(scalar_field), allocatable, dimension(:), intent(inout) :: dqL_prim_dx_vf, dqR_prim_dx_vf, dqL_prim_dy_vf, &
diff --git a/src/simulation/m_riemann_solver_lf.fpp b/src/simulation/m_riemann_solver_lf.fpp
@@ -22,10 +22,9 @@ module m_riemann_solver_lf
 contains
 
     !> Lax-Friedrichs (Rusanov) approximate Riemann solver
-    subroutine s_lf_riemann_solver(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, &
-
-        & dqL_prim_dz_vf, qL_prim_vf, qR_prim_rsx_vf, dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, qR_prim_vf, q_prim_vf, &
-            & flux_vf, flux_src_vf, flux_gsrc_vf, norm_dir, ix, iy, iz)
+    subroutine s_lf_riemann_solver(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, qL_prim_vf, qR_prim_rsx_vf, &
+                                   & dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, qR_prim_vf, q_prim_vf, flux_vf, flux_src_vf, &
+                                   & flux_gsrc_vf, norm_dir, ix, iy, iz)
 
         real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:), intent(inout) :: qL_prim_rsx_vf, qR_prim_rsx_vf
         type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp
@@ -26,10 +26,9 @@ contains
 
     !> Dispatch to the subroutines that are utilized to compute the Riemann problem solution. For additional information please
     !! reference: 1) s_hll_riemann_solver 2) s_hllc_riemann_solver 3) s_lf_riemann_solver 4) s_hlld_riemann_solver
-    subroutine s_riemann_solver(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, &
-
-        & qL_prim_vf, qR_prim_rsx_vf, dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, qR_prim_vf, q_prim_vf, flux_vf, &
-            & flux_src_vf, flux_gsrc_vf, norm_dir, ix, iy, iz)
+    subroutine s_riemann_solver(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, qL_prim_vf, qR_prim_rsx_vf, &
+                                & dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, qR_prim_vf, q_prim_vf, flux_vf, flux_src_vf, &
+                                & flux_gsrc_vf, norm_dir, ix, iy, iz)
 
         real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:), intent(inout) :: qL_prim_rsx_vf, qR_prim_rsx_vf
         type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
diff --git a/src/simulation/m_riemann_state.fpp b/src/simulation/m_riemann_state.fpp
@@ -62,7 +62,6 @@ contains
     !! geometries. For more information please refer to: 1) s_compute_cartesian_viscous_source_flux 2)
     !! s_compute_cylindrical_viscous_source_flux
     subroutine s_compute_viscous_source_flux(velL_vf, dvelL_dx_vf, dvelL_dy_vf, dvelL_dz_vf, velR_vf, dvelR_dx_vf, dvelR_dy_vf, &
-
         & dvelR_dz_vf, flux_src_vf, q_prim_vf, norm_dir, ix, iy, iz)
 
         type(scalar_field), dimension(num_vels), intent(in) :: velL_vf, velR_vf, dvelL_dx_vf, dvelR_dx_vf, dvelL_dy_vf, &
@@ -84,9 +83,8 @@ contains
     end subroutine s_compute_viscous_source_flux
 
     !> Populate the left and right Riemann state variable buffers based on boundary conditions
-    subroutine s_populate_riemann_states_variables_buffers(qL_prim_rsx_vf, dqL_prim_dx_vf, &
-
-        & dqL_prim_dy_vf, dqL_prim_dz_vf, qR_prim_rsx_vf, dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, norm_dir, ix, iy, iz)
+    subroutine s_populate_riemann_states_variables_buffers(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, &
+        & qR_prim_rsx_vf, dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, norm_dir, ix, iy, iz)
 
         real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:), intent(inout) :: qL_prim_rsx_vf, qR_prim_rsx_vf
         type(scalar_field), allocatable, dimension(:), intent(inout) :: dqL_prim_dx_vf, dqR_prim_dx_vf, dqL_prim_dy_vf, &
@@ -176,7 +174,6 @@ contains
             end if
 
             if (bc_x%end == BC_RIEMANN_EXTRAP) then  ! Riemann state extrap. BC at end
-
                 $:GPU_PARALLEL_LOOP(collapse=3)
                 do i = 1, sys_size
                     do l = is3%beg, is3%end
@@ -274,7 +271,6 @@ contains
             end if
 
             if (bc_y%end == BC_RIEMANN_EXTRAP) then  ! Riemann state extrap. BC at end
-
                 $:GPU_PARALLEL_LOOP(collapse=3)
                 do i = 1, sys_size
                     do l = is3%beg, is3%end
@@ -366,7 +362,6 @@ contains
             end if
 
             if (bc_z%end == BC_RIEMANN_EXTRAP) then  ! Riemann state extrap. BC at end
-
                 $:GPU_PARALLEL_LOOP(collapse=3)
                 do i = 1, sys_size
                     do k = is2%beg, is2%end
@@ -565,7 +560,6 @@ contains
 
     !> Compute cylindrical viscous source flux contributions for momentum and energy
     subroutine s_compute_cylindrical_viscous_source_flux(velL_vf, dvelL_dx_vf, dvelL_dy_vf, dvelL_dz_vf, velR_vf, dvelR_dx_vf, &
-
         & dvelR_dy_vf, dvelR_dz_vf, flux_src_vf, q_prim_vf, norm_dir, ix, iy, iz)
 
         type(scalar_field), dimension(num_dims), intent(in)    :: velL_vf, velR_vf
@@ -788,7 +782,6 @@ contains
 
     !> Compute Cartesian viscous source flux contributions for momentum and energy
     subroutine s_compute_cartesian_viscous_source_flux(dvelL_dx_vf, dvelL_dy_vf, dvelL_dz_vf, dvelR_dx_vf, dvelR_dy_vf, &
-
         & dvelR_dz_vf, flux_src_vf, q_prim_vf, norm_dir)
 
         ! Arguments
diff --git a/src/simulation/m_viscous.fpp b/src/simulation/m_viscous.fpp
diff --git a/src/simulation/m_weno.fpp b/src/simulation/m_weno.fpp
diff --git a/tests/coverage_map.json.gz b/tests/coverage_map.json.gz
diff --git a/toolchain/pyproject.toml b/toolchain/pyproject.toml
