MTHINC interface compression method

Author: sbryngelson

docs/documentation/case.md

@@ -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 |
 | `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) |
@@ -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`).
 
 - `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).

docs/module_categories.json

@@ -7,7 +7,8 @@
             "m_weno",
             "m_riemann_solvers",
             "m_muscl",
-            "m_variables_conversion"
+            "m_variables_conversion",
+            "m_thinc"
         ]
     },
     {

examples/1D_sodshocktube_muscl/case.py

@@ -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,

examples/2D_advection_muscl/case.py

@@ -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,

examples/2D_riemann_test_muscl/case.py

@@ -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,

examples/2D_shockdroplet_muscl/case.py

@@ -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,

examples/3D_rayleigh_taylor_muscl/case.py

@@ -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,

examples/3D_shockdroplet_muscl/case.py

@@ -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,

src/simulation/m_global_parameters.fpp

@@ -14,6 +14,7 @@ module m_global_parameters
 
     use m_derived_types
     use m_helper_basic
+
     ! $:USE_GPU_MODULE()
 
     implicit none
@@ -26,8 +27,10 @@ module m_global_parameters
     character(LEN=path_len) :: case_dir       !< Case folder location
     logical                 :: run_time_info  !< Run-time output flag
     integer                 :: t_step_old     !< Existing IC/grid folder
+
     ! Computational Domain Parameters
     integer :: proc_rank  !< Rank of the local processor
+
     !> @name Number of cells in the x-, y- and z-directions, respectively
     !> @{
     integer :: m, n, p
@@ -57,13 +60,15 @@ module m_global_parameters
     !> @{
     real(wp), target, allocatable, dimension(:) :: x_cc, y_cc, z_cc
     !> @}
-    ! type(bounds_info) :: x_domain, y_domain, z_domain !< Locations of the domain bounds in the x-, y- and z-coordinate directions
+    ! type(bounds_info) :: x_domain, y_domain, z_domain !<
+    !! Locations of the domain bounds in the x-, y- and z-coordinate directions
     !> @name Cell-width distributions in the x-, y- and z-directions, respectively
     !> @{
     real(wp), target, allocatable, dimension(:) :: dx, dy, dz
     !> @}
 
     real(wp) :: dt  !< Size of the time-step
+
     $:GPU_DECLARE(create='[x_cb, y_cb, z_cb, x_cc, y_cc, z_cc, dx, dy, dz, dt, m, n, p]')
 
     !> @name Starting time-step iteration, stopping time-step iteration and the number of time-step iterations between successive
@@ -81,6 +86,7 @@ module m_global_parameters
 
     logical :: cfl_adap_dt, cfl_const_dt, cfl_dt
     integer :: t_step_print  !< Number of time-steps between printouts
+
     ! Simulation Algorithm Parameters
     integer :: model_eqns  !< Multicomponent flow model
     #:if MFC_CASE_OPTIMIZATION
@@ -162,7 +168,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
@@ -206,6 +212,7 @@ module m_global_parameters
     integer  :: relax_model   !< Relaxation model
     real(wp) :: palpha_eps    !< trigger parameter for the p relaxation procedure, phase change model
     real(wp) :: ptgalpha_eps  !< trigger parameter for the pTg relaxation procedure, phase change model
+
     $:GPU_DECLARE(create='[relax, relax_model, palpha_eps, ptgalpha_eps]')
 
     integer :: num_bc_patches
@@ -293,7 +300,9 @@ module m_global_parameters
 
     $:GPU_DECLARE(create='[Re_size, Re_size_max, Re_idx]')
 
-    ! WENO averaging flag: use arithmetic mean or unaltered WENO-reconstructed cell-boundary values
+    ! The WENO average (WA) flag regulates whether the calculation of any cell- average spatial derivatives is carried out in each
+    ! cell by utilizing the arithmetic mean of the left and right, WENO-reconstructed, cell-boundary values or simply, the unaltered
+    ! left and right, WENO-reconstructed, cell- boundary values.
     !> @{
     real(wp) :: wa_flg
     !> @}
@@ -311,24 +320,39 @@ module m_global_parameters
 
     $:GPU_DECLARE(create='[dir_idx, dir_flg, dir_idx_tau]')
 
+    !> The number of cells that are necessary to be able to store enough boundary conditions data to march the solution in the
+    !! physical computational domain to the next time-step.
     integer :: buff_size  !< Number of ghost cells for boundary condition storage
+
     $:GPU_DECLARE(create='[buff_size]')
 
-    integer                  :: shear_num              !< Number of shear stress components
-    integer, dimension(3)    :: shear_indices          !< Indices of the stress components that represent shear stress
-    integer                  :: shear_BC_flip_num      !< Number of shear stress components to reflect for boundary conditions
+    integer               :: shear_num !! Number of shear stress components
+    integer, dimension(3) :: shear_indices      !< Indices of the stress components that represent shear stress
+    integer               :: shear_BC_flip_num  !< Number of shear stress components to reflect for boundary conditions
+    !> Indices of shear stress components to reflect for boundary conditions. Size: (1:3, 1:shear_BC_flip_num) for (x/y/z,
+    !! [indices])
     integer, dimension(3, 2) :: shear_BC_flip_indices  !< Shear stress BC reflection indices (1:3, 1:shear_BC_flip_num)
+
     $:GPU_DECLARE(create='[shear_num, shear_indices, shear_BC_flip_num, shear_BC_flip_indices]')
 
     ! END: Simulation Algorithm Parameters
 
     ! Fluids Physical Parameters
 
+    !> Database of the physical parameters of each of the fluids that is present in the flow. These include the stiffened gas
+    !! equation of state parameters, and the Reynolds numbers.
     type(physical_parameters), dimension(num_fluids_max) :: fluid_pp  !< Stiffened gas EOS parameters and Reynolds numbers per fluid
+
     ! Subgrid Bubble Parameters
     type(subgrid_bubble_physical_parameters) :: bub_pp
-    integer                                  :: fd_order   !< Finite-difference order for CoM and flow probe derivatives
-    integer                                  :: fd_number  !< Finite-difference half-stencil size: MAX(1, fd_order/2)
+
+    !> The order of the finite-difference (fd) approximations of the first-order derivatives that need to be evaluated when the CoM
+    !! or flow probe data files are to be written at each time step
+    integer :: fd_order
+
+    !> The finite-difference number is given by MAX(1, fd_order/2). Essentially, it is a measure of the half-size of the
+    !! finite-difference stencil for the selected order of accuracy.
+    integer :: fd_number
     $:GPU_DECLARE(create='[fd_order, fd_number]')
 
     logical                                              :: probe_wrt
@@ -352,6 +376,9 @@ module m_global_parameters
     type(ib_patch_parameters), dimension(num_patches_max) :: patch_ib  !< Immersed boundary patch parameters
     type(vec3_dt), allocatable, dimension(:)              :: airfoil_grid_u, airfoil_grid_l
     integer                                               :: Np
+    !! Database of the immersed boundary patch parameters for each of the patches employed in the configuration of the initial
+    !! condition. Note that the maximum allowable number of patches, num_patches_max, may be changed in the module
+    !! m_derived_types.f90.
 
     $:GPU_DECLARE(create='[ib, num_ibs, patch_ib, Np, airfoil_grid_u, airfoil_grid_l]')
     !> @}
@@ -479,7 +506,10 @@ module m_global_parameters
     $:GPU_DECLARE(create='[Bx0]')
 
     logical :: fft_wrt
+    !> AMDFlang workaround: keep a dummy logical to avoid a compiler case-optimization bug when a parameter+GPU-kernel conditional
+    !! is false
     logical :: dummy  !< AMDFlang workaround for case-optimization + GPU-kernel bug
+
     !> @name Continuum damage model parameters
     !> @{!
     real(wp) :: tau_star       !< Stress threshold for continuum damage modeling
@@ -502,6 +532,7 @@ contains
     impure subroutine s_assign_default_values_to_user_inputs
 
         integer :: i, j  !< Generic loop iterator
+
         ! Logistics
 
         case_dir = '.'
@@ -561,7 +592,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.
@@ -837,15 +868,15 @@ contains
 
     end subroutine s_assign_default_values_to_user_inputs
 
-    !> Initialize the global parameters module
+    !> The computation of parameters, the allocation of memory, the association of pointers and/or the execution of any other
+    !! procedures that are necessary to setup the module.
     impure subroutine s_initialize_global_parameters_module
 
         integer :: i, j, k
         integer :: fac
 
         #: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
@@ -863,7 +894,9 @@ contains
             $:GPU_UPDATE(device='[igr, igr_order, igr_iter_solver]')
         #:endif
 
-        ! Initialize counts: viscous fluids, surface-tension interfaces, curvature interfaces
+        ! Initializing the number of fluids for which viscous effects will be non-negligible, the number of distinctive material
+        ! interfaces for which surface tension will be important and also, the number of fluids for which the physical and geometric
+        ! curvatures of the interfaces will be computed
         Re_size = 0
         Re_size_max = 0
 
@@ -894,8 +927,10 @@ 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
+                    ! Volume fractions are stored in the indices immediately following the energy equation. IGR tracks a total of
+                    ! (N-1) volume fractions for N fluids, hence the "-1" in adv_idx%end. If num_fluids = 1 then adv_idx%end <
+                    ! adv_idx%beg, which skips all loops over the volume fractions since there is no volume fraction to track
+                    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
@@ -1029,7 +1064,8 @@ contains
                 end if
             end if
 
-            ! Count fluids with non-negligible viscous effects (Re > 0)
+            ! Determining the number of fluids for which the shear and the volume Reynolds numbers, e.g. viscous effects, are
+            ! important
             do i = 1, num_fluids
                 if (fluid_pp(i)%Re(1) > 0) Re_size(1) = Re_size(1) + 1
                 if (fluid_pp(i)%Re(2) > 0) Re_size(2) = Re_size(2) + 1
@@ -1157,7 +1193,9 @@ contains
             end do
         end if
 
-        ! Configure WENO averaging flag (arithmetic mean vs. unaltered values)
+        ! Configuring the WENO average flag that will be used to regulate whether any spatial derivatives are to computed in each
+        ! cell by using the arithmetic mean of left and right, WENO-reconstructed, cell-boundary values or otherwise, the unaltered
+        ! left and right, WENO-reconstructed, cell-boundary values
         wa_flg = 0._wp; if (weno_avg) wa_flg = 1._wp
         $:GPU_UPDATE(device='[wa_flg]')
 

src/simulation/m_mpi_proxy.fpp

@@ -22,7 +22,12 @@ module m_mpi_proxy
 
     implicit none
 
+    !> This variable is utilized to pack and send the buffer of the immersed boundary markers, for a single computational domain
+    !! boundary at the time, to the relevant neighboring processor.
     integer, private, allocatable, dimension(:) :: ib_buff_send  !< IB marker send buffer for halo exchange
+
+    !> q_cons_buff_recv is utilized to receive and unpack the buffer of the immersed boundary markers, for a single computational
+    !! domain boundary at the time, from the relevant neighboring processor.
     integer, private, allocatable, dimension(:) :: ib_buff_recv  !< IB marker receive buffer for halo exchange
     integer                                     :: i_halo_size
     $:GPU_DECLARE(create='[i_halo_size]')
@@ -76,7 +81,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
 
@@ -92,7 +97,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
@@ -245,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