fix: correct MTHINC normals on non-uniform grids (#1401)

Co-authored-by: Ben Wilfong <bwilfong3@gatech.edu>
Co-authored-by: Ben Wilfong <48168887+wilfonba@users.noreply.github.com>

Author: sbryngelson

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_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)