Merge branch 'master' into weno-beta-sos

Author: sbryngelson

.github/workflows/test.yml

@@ -172,7 +172,7 @@ jobs:
 
         include:
           - os:    ubuntu
-            mpi:   no-mpi
+            mpi:   mpi
             precision: single
             debug: no-debug
             intel: false
@@ -386,11 +386,12 @@ jobs:
       - name: Test
         if:   '!matrix.nvhpc'
         run:  |
-          /bin/bash mfc.sh test -v --max-attempts 3 -j $(nproc) $ONLY_CHANGES $TEST_ALL $TEST_PCT
+          /bin/bash mfc.sh test -v --max-attempts 3 -j $(nproc) $ONLY_CHANGES $TEST_ALL $TEST_PCT $PRECISION
         env:
           TEST_ALL: ${{ matrix.mpi == 'mpi' && '--test-all' || '' }}
           TEST_PCT: ${{ matrix.debug == 'reldebug' && '-% 20' || '' }}
           ONLY_CHANGES: ${{ github.event_name == 'pull_request' && '--only-changes' || '' }}
+          PRECISION: ${{ matrix.precision != '' && format('--{0}', matrix.precision) || '' }}
 
       # ── NVHPC build + test (via docker exec into long-lived container) ──
       - name: Build (NVHPC)

examples/2D_advection_muscl/case.py

@@ -21,7 +21,7 @@
             "t_step_save": 100,
             # Simulation Algorithm Parameters
             "num_patches": 2,
-            "model_eqns": 3,
+            "model_eqns": 2,
             "alt_soundspeed": "F",
             "num_fluids": 2,
             "mpp_lim": "T",

src/simulation/m_ib_patches.fpp

@@ -21,12 +21,7 @@ module m_ib_patches
 
     implicit none
 
-    private; public :: s_apply_ib_patches, s_update_ib_rotation_matrix, f_convert_cyl_to_cart, s_instantiate_STL_models, &
-        & s_decode_patch_periodicity
-
-    real(wp) :: cart_y, cart_z
-    $:GPU_DECLARE(create='[cart_y, cart_z]')
-    ! Variables to be used to hold cell locations in Cartesian coordinates if 3D simulation is using cylindrical coordinates
+    private; public :: s_apply_ib_patches, s_update_ib_rotation_matrix, s_instantiate_STL_models, s_decode_patch_periodicity
 
 contains
 
@@ -542,18 +537,12 @@ contains
 
         ! Checking whether the sphere covers a particular cell in the domain and verifying whether the current patch has permission
         ! to write to that cell. If both queries check out, the primitive variables of the current patch are assigned to this cell.
-        $:GPU_PARALLEL_LOOP(private='[i, j, k, cart_y, cart_z]', copyin='[encoded_patch_id, center, radius]', collapse=3)
+        $:GPU_PARALLEL_LOOP(private='[i, j, k]', copyin='[encoded_patch_id, center, radius]', collapse=3)
         do k = kl, kr
             do j = jl, jr
                 do i = il, ir
-                    if (grid_geometry == 3) then
-                        call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k))
-                    else
-                        cart_y = y_cc(j)
-                        cart_z = z_cc(k)
-                    end if
                     ! Updating the patch identities bookkeeping variable
-                    if (((x_cc(i) - center(1))**2 + (cart_y - center(2))**2 + (cart_z - center(3))**2 <= radius**2)) then
+                    if (((x_cc(i) - center(1))**2 + (y_cc(j) - center(2))**2 + (z_cc(k) - center(3))**2 <= radius**2)) then
                         ib_markers%sf(i, j, k) = encoded_patch_id
                     end if
                 end do
@@ -603,19 +592,12 @@ contains
         ! Checking whether the cuboid covers a particular cell in the domain and verifying whether the current patch has permission
         ! to write to to that cell. If both queries check out, the primitive variables of the current patch are assigned to this
         ! cell.
-        $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local, cart_y, cart_z]', copyin='[encoded_patch_id, center, length, &
-                            & inverse_rotation]', collapse=3)
+        $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local]', copyin='[encoded_patch_id, center, length, inverse_rotation]', &
+                            & collapse=3)
         do k = kl, kr
             do j = jl, jr
                 do i = il, ir
-                    if (grid_geometry == 3) then
-                        ! TODO :: This does not work and is not covered by any tests. This should be fixed
-                        call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k))
-                    else
-                        cart_y = y_cc(j)
-                        cart_z = z_cc(k)
-                    end if
-                    xyz_local = [x_cc(i), cart_y, cart_z] - center  ! get coordinate frame centered on IB
+                    xyz_local = [x_cc(i), y_cc(j), z_cc(k)] - center  ! get coordinate frame centered on IB
                     xyz_local = matmul(inverse_rotation, xyz_local)  ! rotate the frame into the IB's coordinates
 
                     if (-0.5*length(1) <= xyz_local(1) .and. 0.5*length(1) >= xyz_local(1) .and. -0.5*length(2) <= xyz_local(2) &
@@ -672,18 +654,12 @@ contains
         ! Checking whether the cylinder covers a particular cell in the domain and verifying whether the current patch has the
         ! permission to write to that cell. If both queries check out, the primitive variables of the current patch are assigned to
         ! this cell.
-        $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local, cart_y, cart_z]', copyin='[encoded_patch_id, center, length, radius, &
+        $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local]', copyin='[encoded_patch_id, center, length, radius, &
                             & inverse_rotation]', collapse=3)
         do k = kl, kr
             do j = jl, jr
                 do i = il, ir
-                    if (grid_geometry == 3) then
-                        call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k))
-                    else
-                        cart_y = y_cc(j)
-                        cart_z = z_cc(k)
-                    end if
-                    xyz_local = [x_cc(i), cart_y, cart_z] - center  ! get coordinate frame centered on IB
+                    xyz_local = [x_cc(i), y_cc(j), z_cc(k)] - center  ! get coordinate frame centered on IB
                     xyz_local = matmul(inverse_rotation, xyz_local)  ! rotate the frame into the IB's coordinates
 
                     if (((.not. f_is_default(length(1)) .and. xyz_local(2)**2 + xyz_local(3)**2 <= radius**2 .and. &
@@ -963,30 +939,6 @@ contains
 
     end subroutine s_update_ib_rotation_matrix
 
-    !> Convert cylindrical (r, theta) coordinates to Cartesian (y, z)
-    subroutine s_convert_cylindrical_to_cartesian_coord(cyl_y, cyl_z)
-
-        $:GPU_ROUTINE(parallelism='[seq]')
-
-        real(wp), intent(in) :: cyl_y, cyl_z
-
-        cart_y = cyl_y*sin(cyl_z)
-        cart_z = cyl_y*cos(cyl_z)
-
-    end subroutine s_convert_cylindrical_to_cartesian_coord
-
-    !> Convert a 3D cylindrical coordinate vector (x, r, theta) to Cartesian (x, y, z)
-    pure function f_convert_cyl_to_cart(cyl) result(cart)
-
-        $:GPU_ROUTINE(parallelism='[seq]')
-
-        real(wp), dimension(1:3), intent(in) :: cyl
-        real(wp), dimension(1:3)             :: cart
-
-        cart = (/cyl(1), cyl(2)*sin(cyl(3)), cyl(2)*cos(cyl(3))/)
-
-    end function f_convert_cyl_to_cart
-
     subroutine get_bounding_indices(left_bound, right_bound, cell_centers, left_index, right_index)
 
         real(wp), intent(in)                         :: left_bound, right_bound

src/simulation/m_muscl.fpp

@@ -17,6 +17,7 @@ module m_muscl
     use m_mpi_proxy
     use m_helper
     use m_thinc
+    use m_nvtx
 
     private; public :: s_initialize_muscl_module, s_muscl, s_finalize_muscl_module
 

src/simulation/m_rhs.fpp

@@ -670,7 +670,16 @@ contains
                                                                 & qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, id)
                     end if
                 else
-                    if (all(Re_size == 0)) then
+                    if (int_comp > 0) then
+                        ! THINC reads cont and adv from v_rs_ws; must reconstruct full sys_size range to populate both
+                        iv%beg = 1; iv%end = sys_size
+                        call s_reconstruct_cell_boundary_values(q_prim_qp%vf(1:sys_size), qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, &
+                                                                & qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, id)
+                        ! Surface tension requires first-order energy; overwrite the higher-order result from the full pass above
+                        iv%beg = eqn_idx%E; iv%end = eqn_idx%E
+                        call s_reconstruct_cell_boundary_values_first_order(q_prim_qp%vf(eqn_idx%E), qL_rsx_vf, qL_rsy_vf, &
+                            & qL_rsz_vf, qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, id)
+                    else if (all(Re_size == 0)) then
                         iv%beg = 1; iv%end = eqn_idx%E - 1
                         call s_reconstruct_cell_boundary_values(q_prim_qp%vf(iv%beg:iv%end), qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, &
                                                                 & qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, id)

src/simulation/m_thinc.fpp

@@ -30,32 +30,14 @@ module m_thinc
     real(wp) :: gq3_pts(3) = [-5e-1_wp*0.7745966692414834_wp, 0._wp, 5e-1_wp*0.7745966692414834_wp]
     !> Weights: 5/18, 8/18, 5/18
     real(wp) :: gq3_wts(3) = [5._wp/18._wp, 8._wp/18._wp, 5._wp/18._wp]
-    !> ln(2)
-    real(wp) :: ln2 = 0.6931471805599453_wp
-    $:GPU_DECLARE(create='[gq3_pts, gq3_wts, ln2]')
+    $:GPU_DECLARE(copyin='[gq3_pts, gq3_wts]')
 
     real(wp), allocatable, dimension(:,:,:,:) :: mthinc_nhat !> Unit normal vector
     real(wp), allocatable, dimension(:,:,:)   :: mthinc_d    !> Interface position parameter
     $:GPU_DECLARE(create='[mthinc_nhat, mthinc_d]')
 
 contains
 
-    !> @brief Stable computation of ln(cosh(x))
-    function f_log_cosh(x) result(res)
-
-        $:GPU_ROUTINE(parallelism='[seq]')
-        real(wp), intent(in) :: x
-        real(wp)             :: res, ax
-
-        ax = abs(x)
-        if (ax > 20._wp) then
-            res = ax - ln2
-        else
-            res = ax + log(1._wp + exp(-2._wp*ax)) - ln2
-        end if
-
-    end function f_log_cosh
-
     !> @brief Stable difference: log_cosh(a+h) - log_cosh(a-h) = 2*atanh(tanh(a)*tanh(h)). Avoids catastrophic cancellation when h
     !! is small relative to a.
     function f_log_cosh_diff(a, h) result(res)
@@ -259,16 +241,19 @@ contains
                     ac = v_vf(eqn_idx%adv%beg)%sf(j, k, l)
 
                     if (ac >= ic_eps .and. ac <= 1._wp - ic_eps) then
-                        nr_x = (v_vf(eqn_idx%adv%beg)%sf(j + 1, k, l) - v_vf(eqn_idx%adv%beg)%sf(j - 1, k, l))*5e-1_wp
+                        nr_x = (v_vf(eqn_idx%adv%beg)%sf(j + 1, k, l) - v_vf(eqn_idx%adv%beg)%sf(j - 1, k, &
+                                & l))*(x_cb(j) - x_cb(j - 1))/(x_cc(j + 1) - x_cc(j - 1))
 
                         nr_y = 0._wp
                         if (n > 0) then
-                            nr_y = (v_vf(eqn_idx%adv%beg)%sf(j, k + 1, l) - v_vf(eqn_idx%adv%beg)%sf(j, k - 1, l))*5e-1_wp
+                            nr_y = (v_vf(eqn_idx%adv%beg)%sf(j, k + 1, l) - v_vf(eqn_idx%adv%beg)%sf(j, k - 1, &
+                                    & l))*(y_cb(k) - y_cb(k - 1))/(y_cc(k + 1) - y_cc(k - 1))
                         end if
 
                         nr_z = 0._wp
                         if (p > 0) then
-                            nr_z = (v_vf(eqn_idx%adv%beg)%sf(j, k, l + 1) - v_vf(eqn_idx%adv%beg)%sf(j, k, l - 1))*5e-1_wp
+                            nr_z = (v_vf(eqn_idx%adv%beg)%sf(j, k, l + 1) - v_vf(eqn_idx%adv%beg)%sf(j, k, &
+                                    & l - 1))*(z_cb(l) - z_cb(l - 1))/(z_cc(l + 1) - z_cc(l - 1))
                         end if
 
                         nmag = sqrt(nr_x*nr_x + nr_y*nr_y + nr_z*nr_z)