add normal snapping for near grid aligned normals in MTHINC

Author: wilfonba

src/common/m_constants.fpp

@@ -46,9 +46,10 @@ module m_constants
     integer, parameter :: WENO_TYPE = 1   !< Using WENO for reconstruction type
     integer, parameter :: MUSCL_TYPE = 2  !< Using MUSCL for reconstruction type
     ! Interface Compression
-    real(wp), parameter :: dflt_ic_eps = 1e-4_wp    !< Ensure compression is only applied to surface cells in THINC
-    real(wp), parameter :: dflt_ic_beta = 1.6_wp    !< Sharpness parameter's default value used in THINC
-    real(wp), parameter :: moncon_cutoff = 1e-8_wp  !< Monotonicity constraint's limiter to prevent extremas in THINC
+    real(wp), parameter :: dflt_ic_eps = 1e-4_wp        !< Ensure compression is only applied to surface cells in THINC
+    real(wp), parameter :: dflt_ic_beta = 1.6_wp        !< Sharpness parameter's default value used in THINC
+    real(wp), parameter :: moncon_cutoff = 1e-8_wp      !< Monotonicity constraint's limiter to prevent extremas in THINC
+    real(wp), parameter :: mthinc_align_tol = 1.e-8_wp  !< Grid aligned snap tolerance for MTHINC normals
     ! Chemistry
     real(wp), parameter :: dflt_T_guess = 1200._wp  !< Default guess for temperature (when a previous value is not available)
 

src/simulation/m_muscl.fpp

@@ -211,7 +211,7 @@ contains
             #:endfor
         end if
 
-        if (int_comp > 0 .and. v_size >= advxe) then
+        if (int_comp > 0 .and. v_size >= eqn_idx%adv%end) then
             call nvtxStartRange("WENO-INTCOMP")
             #:for MUSCL_DIR, XYZ in [(1, 'x'), (2, 'y'), (3, 'z')]
                 if (muscl_dir == ${MUSCL_DIR}$) then

src/simulation/m_thinc.fpp

@@ -212,7 +212,7 @@ contains
 
         type(scalar_field), dimension(:), intent(in) :: v_vf
         integer                                      :: j, k, l
-        real(wp)                                     :: nr_x, nr_y, nr_z, nmag, ac
+        real(wp)                                     :: nr_x, nr_y, nr_z, nmag, nmax, ac
         type(int_bounds_info), dimension(3)          :: id_norm
 
         $:GPU_PARALLEL_LOOP(collapse=3, private='[j, k, l]')
@@ -239,7 +239,7 @@ contains
         end if
 
         ! Compute unit normal and solve for d at interior
-        $:GPU_PARALLEL_LOOP(collapse=3, private='[j, k, l, nr_x, nr_y, nr_z, nmag, ac]', copyin='[id_norm]')
+        $:GPU_PARALLEL_LOOP(collapse=3, private='[j, k, l, nr_x, nr_y, nr_z, nmag, nmax, ac]', copyin='[id_norm]')
         do l = id_norm(3)%beg, id_norm(3)%end
             do k = id_norm(2)%beg, id_norm(2)%end
                 do j = id_norm(1)%beg, id_norm(1)%end
@@ -265,6 +265,16 @@ contains
                             nr_y = nr_y/nmag
                             nr_z = nr_z/nmag
 
+                            ! Snap near-grid-aligned normals to exact alignment
+                            nmax = max(abs(nr_x), abs(nr_y), abs(nr_z))
+                            if (abs(nr_x) < mthinc_align_tol*nmax) nr_x = 0._wp
+                            if (abs(nr_y) < mthinc_align_tol*nmax) nr_y = 0._wp
+                            if (abs(nr_z) < mthinc_align_tol*nmax) nr_z = 0._wp
+                            nmag = sqrt(nr_x*nr_x + nr_y*nr_y + nr_z*nr_z)
+                            nr_x = nr_x/nmag
+                            nr_y = nr_y/nmag
+                            nr_z = nr_z/nmag
+
                             mthinc_nhat(1, j, k, l) = nr_x
                             mthinc_nhat(2, j, k, l) = nr_y
                             mthinc_nhat(3, j, k, l) = nr_z

src/simulation/m_weno.fpp

@@ -1372,7 +1372,7 @@ contains
             #:endif
         end if
 
-        if (int_comp > 0 .and. v_size >= advxe) then
+        if (int_comp > 0 .and. v_size >= eqn_idx%adv%end) then
             call nvtxStartRange("WENO-INTCOMP")
             #:for WENO_DIR, XYZ in [(1, 'x'), (2, 'y'), (3, 'z')]
                 if (weno_dir == ${WENO_DIR}$) then