MFlowCode · sbryngelson · Mar 26, 2026 · Mar 26, 2026 · Mar 27, 2026 · Mar 29, 2026
diff --git a/docs/documentation/case.md b/docs/documentation/case.md
@@ -446,8 +446,8 @@ See @ref equations "Equations" for the mathematical models these parameters cont
 | `muscl_order`              | Integer | MUSCL order [1,2] |
 | `muscl_lim`                | Integer | MUSCL Slope Limiter: [1] minmod; [2] monotonized central; [3] Van Albada; [4] Van Leer; [5] SUPERBEE |
 | `muscl_eps`                | Real    | MUSCL limiter slope-product threshold (default: hard-coded thresholds; set to 0 for textbook behavior) |
+| `int_comp`                 | Integer | Interface Compression [0] Off [1] THINC [2] MTHINC (default 0) |
 | `flux_lim`                 | Integer | Flux limiter for post-process: [1] minmod; [2] MUSCL; [3] OSPRE; [4] SUPERBEE |
-| `int_comp`                 | Logical | THINC Interface Compression |
 | `ic_eps`                   | Real    | Interface compression threshold (default: 1e-4) |
 | `ic_beta`                  | Real    | Interface compression sharpness parameter (default: 1.6) |
 | `riemann_solver`           | Integer | Riemann solver algorithm: [1] HLL*; [2] HLLC; [3] Exact*; [4] HLLD	(only for MHD) |
@@ -547,7 +547,7 @@ It is recommended to set `weno_eps` to $10^{-6}$ for WENO-JS, and to $10^{-40}$
 When not set (default), the threshold is 1e-9 for minmod/MC, and 1e-6 for others.
 Setting `muscl_eps = 0` gives textbook limiter behavior where limiters activate whenever both slopes have the same sign.
 
-- `int_comp` activates interface compression using THINC used in MUSCL Reconstruction, with control parameters (`ic_eps`, and `ic_beta`).
+- `int_comp` activates interface compression using [1] THINC or [2] MTHINC (default off) used in variable reconstruction, with control parameters (`ic_eps`, and `ic_beta`).
 
 - `riemann_solver` specifies the choice of the Riemann solver that is used in simulation by an integer from 1 through 4.
 `riemann_solver = 1`, `2`, and `3` correspond to HLL, HLLC, and Exact Riemann solver, respectively (\cite Toro09).

diff --git a/docs/module_categories.json b/docs/module_categories.json
@@ -7,7 +7,8 @@
             "m_weno",
             "m_riemann_solvers",
             "m_muscl",
-            "m_variables_conversion"
+            "m_variables_conversion",
+            "m_thinc"
         ]
     },
     {

@@ -36,7 +36,6 @@
             "recon_type": 2,
             "muscl_order": 2,
             "muscl_lim": 2,
-            "int_comp": "T",
             "riemann_solver": 2,
             "wave_speeds": 1,
             "avg_state": 2,

@@ -30,7 +30,7 @@
             "recon_type": 2,
             "muscl_order": 2,
             "muscl_lim": 2,
-            "int_comp": "T",
+            "int_comp": 1,
             "null_weights": "F",
             "riemann_solver": 2,
             "wave_speeds": 1,

@@ -34,7 +34,6 @@
             # "weno_eps": 1e-16,
             "muscl_order": 2,
             "muscl_lim": 1,
-            "int_comp": "T",
             "riemann_solver": 2,
             "wave_speeds": 1,
             "avg_state": 2,

@@ -59,7 +59,7 @@
             "recon_type": 2,
             "muscl_order": 2,
             "muscl_lim": 4,
-            "int_comp": "T",
+            "int_comp": 1,
             "null_weights": "F",
             "riemann_solver": 2,
             "wave_speeds": 1,

@@ -57,7 +57,7 @@
     "recon_type": 2,
     "muscl_order": 2,
     "muscl_lim": 4,
-    "int_comp": "T",
+    "int_comp": 1,
     "avg_state": 2,
     "riemann_solver": 2,
     "wave_speeds": 1,

@@ -202,7 +202,7 @@
             "recon_type": 2,
             "muscl_order": 2,
             "muscl_lim": 3,
-            "int_comp": "T",
+            "int_comp": 1,
             "riemann_solver": 2,
             "wave_speeds": 1,
             "avg_state": 2,

diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp
@@ -47,9 +47,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)
 

diff --git a/src/simulation/m_checker.fpp b/src/simulation/m_checker.fpp
@@ -80,6 +80,8 @@ contains
         @:PROHIBIT(p + 1 < min(1, p)*num_stcls_min*muscl_order, &
                    & "For 3D simulation, p must be greater than or equal to (num_stcls_min*muscl_order - 1), whose value is " &
                    & // trim(numStr))
+        @:PROHIBIT(muscl_order == 1 .and. int_comp > 0, &
+                   & "int_comp requires muscl_order >= 2 (muscl_order=1 leaves the reconstruction workspace uninitialised)")
 
     end subroutine s_check_inputs_muscl
 

diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp
@@ -148,24 +148,25 @@ module m_global_parameters
     logical :: nv_uvm_pref_gpu  !< Enable explicit gpu memory hints (default FALSE)
     !> @}
 
-    real(wp)           :: muscl_eps                    !< MUSCL limiter slope-product threshold
-    real(wp)           :: weno_eps                     !< Binding for the WENO nonlinear weights
-    real(wp)           :: teno_CT                      !< Smoothness threshold for TENO
-    logical            :: mp_weno                      !< Monotonicity preserving (MP) WENO
-    logical            :: weno_avg                     !< Average left/right cell-boundary states
-    logical            :: weno_Re_flux                 !< WENO reconstruct velocity gradients for viscous stress tensor
-    integer            :: riemann_solver               !< Riemann solver algorithm
-    integer            :: low_Mach                     !< Low Mach number fix to HLLC Riemann solver
-    integer            :: wave_speeds                  !< Wave speeds estimation method
-    integer            :: avg_state                    !< Average state evaluation method
-    logical            :: alt_soundspeed               !< Alternate mixture sound speed
-    logical            :: null_weights                 !< Null undesired WENO weights
-    logical            :: mixture_err                  !< Mixture properties correction
-    logical            :: hypoelasticity               !< hypoelasticity modeling
-    logical            :: hyperelasticity              !< hyperelasticity modeling
-    logical            :: int_comp                     !< THINC interface compression
-    real(wp)           :: ic_eps                       !< THINC Epsilon to compress on surface cells
-    real(wp)           :: ic_beta                      !< THINC Sharpness Parameter
+    real(wp) :: muscl_eps        !< MUSCL limiter slope-product threshold
+    real(wp) :: weno_eps         !< Binding for the WENO nonlinear weights
+    real(wp) :: teno_CT          !< Smoothness threshold for TENO
+    logical  :: mp_weno          !< Monotonicity preserving (MP) WENO
+    logical  :: weno_avg         !< Average left/right cell-boundary states
+    logical  :: weno_Re_flux     !< WENO reconstruct velocity gradients for viscous stress tensor
+    integer  :: riemann_solver   !< Riemann solver algorithm
+    integer  :: low_Mach         !< Low Mach number fix to HLLC Riemann solver
+    integer  :: wave_speeds      !< Wave speeds estimation method
+    integer  :: avg_state        !< Average state evaluation method
+    logical  :: alt_soundspeed   !< Alternate mixture sound speed
+    logical  :: null_weights     !< Null undesired WENO weights
+    logical  :: mixture_err      !< Mixture properties correction
+    logical  :: hypoelasticity   !< hypoelasticity modeling
+    logical  :: hyperelasticity  !< hyperelasticity modeling
+    integer  :: int_comp         !< Interface compression: 0=off, 1=THINC, 2=MTHINC
+    real(wp) :: ic_eps           !< THINC Epsilon to compress on surface cells
+    real(wp) :: ic_beta          !< THINC Sharpness Parameter
+    $:GPU_DECLARE(create='[int_comp, ic_eps, ic_beta]')
     integer            :: hyper_model                  !< hyperelasticity solver algorithm
     logical            :: elasticity                   !< elasticity modeling, true for hyper or hypo
     logical, parameter :: chemistry = .${chemistry}$.  !< Chemistry modeling
@@ -546,7 +547,7 @@ contains
         ptgalpha_eps = dflt_real
         hypoelasticity = .false.
         hyperelasticity = .false.
-        int_comp = .false.
+        int_comp = 0
         ic_eps = dflt_ic_eps
         ic_beta = dflt_ic_beta
         elasticity = .false.
@@ -841,7 +842,6 @@ contains
 
         #:if not MFC_CASE_OPTIMIZATION
             ! Determining the degree of the WENO polynomials
-
             if (recon_type == WENO_TYPE) then
                 weno_polyn = (weno_order - 1)/2
                 if (teno) then
@@ -1178,7 +1178,7 @@ contains
             fd_number = max(1, fd_order/2)
         end if
 
-        if (mhd) then  ! TODO merge with above; waiting for hyperelasticity PR
+        if (mhd) then
             fd_number = max(1, fd_order/2)
         end if
 
@@ -1201,7 +1201,7 @@ contains
             grid_geometry = 1
         else if (cyl_coord .and. p == 0) then  ! Axisymmetric cylindrical grid
             grid_geometry = 2
-        else  ! Fully 3D cylindrical grid
+        else
             grid_geometry = 3
         end if
 
@@ -1232,6 +1232,7 @@ contains
             $:GPU_UPDATE(device='[num_fluids, num_dims, viscous, num_vels, nb, muscl_lim]')
         #:endif
 
+        $:GPU_UPDATE(device='[int_comp, ic_eps, ic_beta]')
         $:GPU_UPDATE(device='[muscl_eps]')
         $:GPU_UPDATE(device='[dir_idx, dir_flg, dir_idx_tau]')
 

diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp
@@ -76,7 +76,7 @@ contains
             & 'num_probes', 'num_integrals', 'bubble_model', 'thermal',        &
             & 'num_source', 'relax_model', 'num_ibs', 'n_start',    &
             & 'num_bc_patches', 'num_igr_iters', 'num_igr_warm_start_iters', &
-            & 'adap_dt_max_iters', 'collision_model' ]
+            & 'adap_dt_max_iters', 'int_comp', 'collision_model' ]
             call MPI_BCAST(${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
         #:endfor
 
@@ -94,7 +94,7 @@ contains
             & 'bc_x%isothermal_out', 'bc_y%isothermal_out', 'bc_z%isothermal_out', &
             & 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt', 'surface_tension',       &
             & 'shear_stress', 'bulk_stress', 'bubbles_lagrange',                &
-            & 'hyperelasticity', 'down_sample', 'int_comp','fft_wrt', &
+            & 'hyperelasticity', 'down_sample', 'fft_wrt', &
             & 'hyper_cleaning', 'ib_state_wrt']
             call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
         #:endfor
@@ -250,7 +250,6 @@ contains
 
 #ifdef MFC_MPI
         integer :: ierr  !< Generic flag used to identify and report MPI errors
-
         call MPI_BCAST(phi_rn, num_freq, mpi_p, 0, MPI_COMM_WORLD, ierr)
 #endif
 

diff --git a/src/simulation/m_muscl.fpp b/src/simulation/m_muscl.fpp
@@ -16,8 +16,10 @@ 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, s_interface_compression
+    private; public :: s_initialize_muscl_module, s_muscl, s_finalize_muscl_module
 
     integer :: v_size
     $:GPU_DECLARE(create='[v_size]')
@@ -209,85 +211,19 @@ contains
             #:endfor
         end if
 
-        if (int_comp) then
-            call s_interface_compression(vL_rs_vf_x, vL_rs_vf_y, vL_rs_vf_z, vR_rs_vf_x, vR_rs_vf_y, vR_rs_vf_z, muscl_dir, &
-                                         & is1_muscl_d, is2_muscl_d, is3_muscl_d)
+        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
+                    call s_thinc_compression(v_rs_ws_${XYZ}$_muscl, vL_rs_vf_x, vL_rs_vf_y, vL_rs_vf_z, vR_rs_vf_x, vR_rs_vf_y, &
+                                             & vR_rs_vf_z, muscl_dir, is1_muscl, is2_muscl, is3_muscl)
+                end if
+            #:endfor
+            call nvtxEndRange()
         end if
 
     end subroutine s_muscl
 
-    !> Apply THINC interface-compression to sharpen volume-fraction reconstructions at material interfaces
-    subroutine s_interface_compression(vL_rs_vf_x, vL_rs_vf_y, vL_rs_vf_z, vR_rs_vf_x, vR_rs_vf_y, vR_rs_vf_z, muscl_dir, &
-
-        & is1_muscl_d, is2_muscl_d, is3_muscl_d)
-
-        real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:), intent(inout) :: vL_rs_vf_x, vL_rs_vf_y, &
-             & vL_rs_vf_z, vR_rs_vf_x, vR_rs_vf_y, vR_rs_vf_z
-        integer, intent(in)               :: muscl_dir
-        type(int_bounds_info), intent(in) :: is1_muscl_d, is2_muscl_d, is3_muscl_d
-        integer                           :: j, k, l
-        real(wp)                          :: aCL, aCR, aC, aTHINC, qmin, qmax, A, B, C, sign, moncon
-
-        #:for MUSCL_DIR, XYZ in [(1, 'x'), (2, 'y'), (3, 'z')]
-            if (muscl_dir == ${MUSCL_DIR}$) then
-                $:GPU_PARALLEL_LOOP(collapse=3,private='[j, k, l, aCL, aC, aCR, aTHINC, moncon, sign, qmin, qmax]')
-                do l = is3_muscl%beg, is3_muscl%end
-                    do k = is2_muscl%beg, is2_muscl%end
-                        do j = is1_muscl%beg, is1_muscl%end
-                            aCL = v_rs_ws_${XYZ}$_muscl(j - 1, k, l, eqn_idx%adv%beg)
-                            aC = v_rs_ws_${XYZ}$_muscl(j, k, l, eqn_idx%adv%beg)
-                            aCR = v_rs_ws_${XYZ}$_muscl(j + 1, k, l, eqn_idx%adv%beg)
-
-                            moncon = (aCR - aC)*(aC - aCL)
-
-                            if (aC >= ic_eps .and. aC <= 1._wp - ic_eps .and. moncon > moncon_cutoff) then  ! Interface cell
-
-                                if (aCR - aCL > 0._wp) then
-                                    sign = 1._wp
-                                else
-                                    sign = -1._wp
-                                end if
-
-                                qmin = min(aCR, aCL)
-                                qmax = max(aCR, aCL) - qmin
-
-                                C = (aC - qmin + sgm_eps)/(qmax + sgm_eps)
-                                B = exp(sign*ic_beta*(2._wp*C - 1._wp))
-                                A = (B/cosh(ic_beta) - 1._wp)/tanh(ic_beta)
-
-                                ! Left reconstruction
-                                aTHINC = qmin + 5e-1_wp*qmax*(1._wp + sign*A)
-                                if (aTHINC < ic_eps) aTHINC = ic_eps
-                                if (aTHINC > 1 - ic_eps) aTHINC = 1 - ic_eps
-                                vL_rs_vf_${XYZ}$ (j, k, l, eqn_idx%cont%beg) = vL_rs_vf_${XYZ}$ (j, k, l, &
-                                                  & eqn_idx%cont%beg)/vL_rs_vf_${XYZ}$ (j, k, l, eqn_idx%adv%beg)*aTHINC
-                                vL_rs_vf_${XYZ}$ (j, k, l, eqn_idx%cont%end) = vL_rs_vf_${XYZ}$ (j, k, l, &
-                                                  & eqn_idx%cont%end)/(1._wp - vL_rs_vf_${XYZ}$ (j, k, l, &
-                                                  & eqn_idx%adv%beg))*(1._wp - aTHINC)
-                                vL_rs_vf_${XYZ}$ (j, k, l, eqn_idx%adv%beg) = aTHINC
-                                vL_rs_vf_${XYZ}$ (j, k, l, eqn_idx%adv%end) = 1 - aTHINC
-
-                                ! Right reconstruction
-                                aTHINC = qmin + 5e-1_wp*qmax*(1._wp + sign*(tanh(ic_beta) + A)/(1._wp + A*tanh(ic_beta)))
-                                if (aTHINC < ic_eps) aTHINC = ic_eps
-                                if (aTHINC > 1 - ic_eps) aTHINC = 1 - ic_eps
-                                vR_rs_vf_${XYZ}$ (j, k, l, eqn_idx%cont%beg) = vL_rs_vf_${XYZ}$ (j, k, l, &
-                                                  & eqn_idx%cont%beg)/vL_rs_vf_${XYZ}$ (j, k, l, eqn_idx%adv%beg)*aTHINC
-                                vR_rs_vf_${XYZ}$ (j, k, l, eqn_idx%cont%end) = vL_rs_vf_${XYZ}$ (j, k, l, &
-                                                  & eqn_idx%cont%end)/(1._wp - vL_rs_vf_${XYZ}$ (j, k, l, &
-                                                  & eqn_idx%adv%beg))*(1._wp - aTHINC)
-                                vR_rs_vf_${XYZ}$ (j, k, l, eqn_idx%adv%beg) = aTHINC
-                                vR_rs_vf_${XYZ}$ (j, k, l, eqn_idx%adv%end) = 1 - aTHINC
-                            end if
-                        end do
-                    end do
-                end do
-                $:END_GPU_PARALLEL_LOOP()
-            end if
-        #:endfor
-
-    end subroutine s_interface_compression
-
     !> Reshape cell-averaged variable data into direction-local work arrays for MUSCL reconstruction
     subroutine s_initialize_muscl(v_vf, muscl_dir)