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 1 commit into
base: master
Choose a base branch
from
+2,353 −2,018
Draft

MTHINC #1303

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
6 changes: 3 additions & 3 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,9 +162,9 @@ 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
real(wp) :: ic_beta !< THINC/MTHINC Sharpness Parameter
integer :: hyper_model !< hyperelasticity solver algorithm
logical :: elasticity !< elasticity modeling, true for hyper or hypo
logical, parameter :: chemistry = .${chemistry}$. !< Chemistry modeling
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
4 changes: 2 additions & 2 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
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
16 changes: 15 additions & 1 deletion src/simulation/m_rhs.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@ module m_rhs
use m_body_forces
use m_chemistry
use m_igr
use m_thinc
use m_pressure_relaxation

implicit none
Expand Down Expand Up @@ -641,6 +642,12 @@ contains
call nvtxEndRange
end if

if (int_comp == 2 .and. n > 0) then
call nvtxStartRange("RHS-COMPRESSION-NORMALS")
call s_compute_mthinc_normals(q_prim_qp%vf)
call nvtxEndRange
end if
Comment on lines +645 to +649
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 | 🟠 Major

Use the same THINC enablement guard here as in startup.

src/simulation/m_start_up.fpp:1114-1122 only initializes THINC when .not. igr .or. dummy, but this block runs for any int_comp == 2. In a pure IGR run, q_prim_qp%vf is not prepared for THINC and the MTHINC workspace was never allocated.

Minimal fix 
-        if (int_comp == 2 .and. n > 0) then
+        if (int_comp == 2 .and. n > 0 .and. (.not. igr .or. dummy)) then
             call nvtxStartRange("RHS-COMPRESSION-NORMALS")
             call s_compute_mthinc_normals(q_prim_qp%vf)
             call nvtxEndRange
         end if
📝 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
if (int_comp == 2 .and. n > 0) then
call nvtxStartRange("RHS-COMPRESSION-NORMALS")
call s_compute_mthinc_normals(q_prim_qp%vf)
call nvtxEndRange
end if
if (int_comp == 2 .and. n > 0 .and. (.not. igr .or. dummy)) then
call nvtxStartRange("RHS-COMPRESSION-NORMALS")
call s_compute_mthinc_normals(q_prim_qp%vf)
call nvtxEndRange
end if


! Loop over coordinate directions for dimensional splitting
do id = 1, num_dims
if (igr .or. dummy) then
Expand Down Expand Up @@ -675,7 +682,7 @@ contains
if ((.not. igr) .or. dummy) then ! Finite volume solve

! Reconstructing Primitive/Conservative Variables
call nvtxStartRange("RHS-WENO")
call nvtxStartRange("RHS-RECONSTRUCTION")

if (.not. surface_tension) then
if (all(Re_size == 0)) then
Expand Down Expand Up @@ -1661,6 +1668,13 @@ contains
call s_${SCHEME}$ (v_vf(iv%beg:iv%end), vL_x(:,:,:,iv%beg:iv%end), vL_y(:,:,:,:), vL_z(:,:,:,:), vR_x(:,:,:, &
& iv%beg:iv%end), vR_y(:,:,:,:), vR_z(:,:,:,:), recon_dir, is1, is2, is3)
end if

if (int_comp > 0 .and. iv%beg <= advxb .and. iv%end >= advxe) then
! Only run interface compression when the volume fractions are in the reconstructed variables
call nvtxStartRange("RHS-RECONSTRUCTION-COMPRESSION")
call s_thinc_compression(q_prim_qp%vf, vL_x, vL_y, vL_z, vR_x, vR_y, vR_z, recon_dir, is1, is2, is3)
call nvtxEndRange()
end if
end if
#:endfor

Expand Down
3 changes: 3 additions & 0 deletions src/simulation/m_start_up.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ module m_start_up
use m_variables_conversion
use m_weno
use m_muscl
use m_thinc
use m_riemann_solvers
use m_cbc
use m_boundary_common
Expand Down Expand Up @@ -926,6 +927,7 @@ contains
else if (recon_type == MUSCL_TYPE) then
call s_initialize_muscl_module()
end if
if (int_comp > 0) call s_initialize_thinc_module()
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

Mirror THINC finalization with the initialization predicate.

Initialization happens inside .not. igr .or. dummy, but finalization only runs in the igr == .false. branch. In igr=.true., dummy=.true. runs, the new THINC/MTHINC state can be allocated and never released.

As per coding guidelines, Always pair @:ALLOCATE calls with matching @:DEALLOCATE calls in the corresponding finalization subroutine.

Also applies to: 1304-1304

call s_initialize_cbc_module()
call s_initialize_riemann_solvers_module()
end if
Expand Down Expand Up @@ -1097,6 +1099,7 @@ contains
else if (recon_type == MUSCL_TYPE) then
call s_finalize_muscl_module()
end if
if (int_comp > 0) call s_finalize_thinc_module()
end if
call s_finalize_variables_conversion_module()
if (grid_geometry == 3) call s_finalize_fftw_module
Expand Down
Loading
Loading