Skip to content

MTHINC #1303

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wilfonba wants to merge 2 commits into
base: master
Choose a base branch
from
Draft

MTHINC #1303

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
5a8ca06
MTHINC interface compression method
sbryngelson Mar 26, 2026
b7f85c6
fix: remove stale THINC comment from m_muscl.fpp
sbryngelson Mar 26, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions docs/documentation/case.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -433,8 +433,8 @@ See @ref equations "Equations" for the mathematical models these parameters cont
| `mp_weno` | Logical | Monotonicity preserving WENO |
| `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 |
| `int_comp` | Integer | Interface Compression [1] THINC [2] MTHINC |
Copy link

Copilot AI Mar 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The parameter table entry for int_comp was updated to an integer but no longer documents the 0=off value (it currently only lists [1] THINC [2] MTHINC). This can confuse users and contradicts other docs/tooling that treat 0 as the default/off value. Please update the table text to include 0 explicitly.

Suggested change
| `int_comp` | Integer | Interface Compression [1] THINC [2] MTHINC |
| `int_comp` | Integer | Interface Compression: [0] off; [1] THINC; [2] MTHINC |

Copilot uses AI. Check for mistakes.
Copy link
Contributor

@coderabbitai coderabbitai bot Mar 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟡 Minor

Document all three valid states for int_comp parameter.

The table entry only lists values [1] and [2], but according to the PR objectives, int_comp is a three-state enumeration: 0 (off/disabled), 1 (THINC), and 2 (MTHINC). Users need to know that setting int_comp=0 disables interface compression.

📝 Suggested documentation fix 
-| `int_comp`                 | Integer | Interface Compression [1] THINC [2] MTHINC |
+| `int_comp`                 | Integer | Interface Compression: [0] Off; [1] THINC; [2] MTHINC |
📝 Committable suggestion

‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.

Suggested change
| `int_comp` | Integer | Interface Compression [1] THINC [2] MTHINC |
| `int_comp` | Integer | Interface Compression: [0] Off; [1] THINC; [2] MTHINC |

| `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) |
Expand Down Expand Up @@ -530,7 +530,7 @@ It is recommended to set `weno_eps` to $10^{-6}$ for WENO-JS, and to $10^{-40}$
- `muscl_lim` specifies the slope limiter that is used in 2nd order MUSCL Reconstruction by an integer from 1 through 5.
`muscl_lim = 1`, `2`, `3`, `4`, and `5` correspond to minmod, monotonized central, Van Albada, Van Leer, and SUPERBEE, respectively.

- `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 used in variable reconstruction, with control parameters (`ic_eps`, and `ic_beta`).
Copy link
Contributor

@coderabbitai coderabbitai bot Mar 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟡 Minor

Clarify all three states and the behavior when int_comp=0.

The description states "int_comp activates interface compression using [1] THINC or [2] MTHINC," but doesn't document the disabled state (0) or explain what happens when int_comp=0. Since this parameter changed from a boolean to a three-state integer enumeration, all valid values should be clearly documented to help users understand how to disable the feature.

📝 Suggested documentation improvement 
-- `int_comp` activates interface compression using [1] THINC or [2] MTHINC used in variable reconstruction, with control parameters (`ic_eps`, and `ic_beta`).
+- `int_comp` controls interface compression in variable reconstruction with control parameters (`ic_eps`, and `ic_beta`). `int_comp = 0` disables interface compression; `int_comp = 1` enables THINC compression; `int_comp = 2` enables MTHINC compression.

Based on learnings, modifications to public parameters affect downstream users. This parameter changed from boolean to integer semantics, so complete documentation of all states is essential for users migrating from the previous boolean interface.


- `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).
Expand Down
3 changes: 2 additions & 1 deletion docs/module_categories.json
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,8 @@
"m_weno",
"m_riemann_solvers",
"m_muscl",
"m_variables_conversion"
"m_variables_conversion",
"m_thinc"
]
},
{
Expand Down
2 changes: 1 addition & 1 deletion examples/1D_sodshocktube_muscl/case.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@
"recon_type": 2,
"muscl_order": 2,
"muscl_lim": 2,
"int_comp": "T",
"int_comp": 1,
"riemann_solver": 2,
"wave_speeds": 1,
"avg_state": 2,
Expand Down
2 changes: 1 addition & 1 deletion examples/2D_advection_muscl/case.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
2 changes: 1 addition & 1 deletion examples/2D_riemann_test_muscl/case.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@
# "weno_eps": 1e-16,
"muscl_order": 2,
"muscl_lim": 1,
"int_comp": "T",
"int_comp": 1,
"riemann_solver": 2,
"wave_speeds": 1,
"avg_state": 2,
Expand Down
2 changes: 1 addition & 1 deletion examples/2D_shockdroplet_muscl/case.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
2 changes: 1 addition & 1 deletion examples/3D_rayleigh_taylor_muscl/case.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
2 changes: 1 addition & 1 deletion examples/3D_shockdroplet_muscl/case.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
18 changes: 8 additions & 10 deletions src/simulation/m_global_parameters.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -162,7 +162,7 @@ module m_global_parameters
logical :: mixture_err !< Mixture properties correction
logical :: hypoelasticity !< hypoelasticity modeling
logical :: hyperelasticity !< hyperelasticity modeling
logical :: int_comp !< THINC interface compression
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
integer :: hyper_model !< hyperelasticity solver algorithm
Expand Down Expand Up @@ -561,7 +561,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.
Expand Down Expand Up @@ -845,7 +845,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
Expand Down Expand Up @@ -894,8 +893,7 @@ contains
E_idx = mom_idx%end + 1

if (igr) then
! IGR: volume fractions after energy (N-1 for N fluids; skipped when num_fluids=1)
adv_idx%beg = E_idx + 1 ! Alpha for fluid 1
adv_idx%beg = E_idx + 1
adv_idx%end = E_idx + num_fluids - 1
else
! Volume fractions are stored in the indices immediately following the energy equation. WENO/MUSCL + Riemann
Expand Down Expand Up @@ -990,11 +988,11 @@ contains
internalEnergies_idx%end = adv_idx%end + num_fluids
sys_size = internalEnergies_idx%end
else if (model_eqns == 4) then
cont_idx%beg = 1 ! one continuity equation
cont_idx%beg = 1
cont_idx%end = 1 ! num_fluids
mom_idx%beg = cont_idx%end + 1 ! one momentum equation in each direction
mom_idx%beg = cont_idx%end + 1
mom_idx%end = cont_idx%end + num_vels
E_idx = mom_idx%end + 1 ! one energy equation
E_idx = mom_idx%end + 1
adv_idx%beg = E_idx + 1
adv_idx%end = adv_idx%beg ! one volume advection equation
alf_idx = adv_idx%end
Expand Down Expand Up @@ -1173,7 +1171,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

Expand All @@ -1196,7 +1194,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

Expand Down
5 changes: 2 additions & 3 deletions src/simulation/m_mpi_proxy.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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' ]
& 'adap_dt_max_iters', 'int_comp' ]
call MPI_BCAST(${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
#:endfor

Expand All @@ -92,7 +92,7 @@ contains
& 'bc_z%grcbc_in', 'bc_z%grcbc_out', 'bc_z%grcbc_vel_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
Expand Down Expand Up @@ -245,7 +245,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

Expand Down
77 changes: 1 addition & 76 deletions src/simulation/m_muscl.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ module m_muscl
use m_mpi_proxy
use m_helper

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]')
Expand Down Expand Up @@ -210,83 +210,8 @@ 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)
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, advxb)
aC = v_rs_ws_${XYZ}$_muscl(j, k, l, advxb)
aCR = v_rs_ws_${XYZ}$_muscl(j + 1, k, l, advxb)

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, contxb) = vL_rs_vf_${XYZ}$ (j, k, l, contxb)/vL_rs_vf_${XYZ}$ (j, k, &
& l, advxb)*aTHINC
vL_rs_vf_${XYZ}$ (j, k, l, contxe) = vL_rs_vf_${XYZ}$ (j, k, l, &
& contxe)/(1._wp - vL_rs_vf_${XYZ}$ (j, k, l, advxb))*(1._wp - aTHINC)
vL_rs_vf_${XYZ}$ (j, k, l, advxb) = aTHINC
vL_rs_vf_${XYZ}$ (j, k, l, advxe) = 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, contxb) = vL_rs_vf_${XYZ}$ (j, k, l, contxb)/vL_rs_vf_${XYZ}$ (j, k, &
& l, advxb)*aTHINC
vR_rs_vf_${XYZ}$ (j, k, l, contxe) = vL_rs_vf_${XYZ}$ (j, k, l, &
& contxe)/(1._wp - vL_rs_vf_${XYZ}$ (j, k, l, advxb))*(1._wp - aTHINC)
vR_rs_vf_${XYZ}$ (j, k, l, advxb) = aTHINC
vR_rs_vf_${XYZ}$ (j, k, l, advxe) = 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)

Expand Down
Loading
Loading