Adding MUSCL Reconstruction

README.md

@@ -165,7 +165,9 @@ They are organized below.
 
 * Shock and interface capturing schemes
 	* First-order upwinding
- 	* WENO reconstructions of order 3, 5, and 7
+ 	* MUSCL (order 2)
+  		* Slope limiters: minmod, monotonized central, Van Albada, Van Leer, superbee
+ 	* WENO reconstructions (orders 3, 5, and 7)
   	* WENO variants: WENO-JS, WENO-M, WENO-Z, TENO
    	* Monotonicity-preserving reconstructions
 	* Reliable handling of large density ratios
@@ -187,7 +189,7 @@ They are organized below.
 * Ideal weak scaling to 100% of the largest GPU and superchip supercomputers
  	* \>36K AMD APUs (MI300A) on [LLNL El Capitan](https://hpc.llnl.gov/hardware/compute-platforms/el-capitan)
    	* \>3K AMD APUs (MI300A) on [LLNL Tuolumne](https://hpc.llnl.gov/hardware/compute-platforms/tuolumne)
-	* \>33K AMD GPUs (MI250X) on the first exascale computer, [OLCF Frontier](https://www.olcf.ornl.gov/frontier/) 
+	* \>33K AMD GPUs (MI250X) on [OLCF Frontier](https://www.olcf.ornl.gov/frontier/) 
 	* \>10K NVIDIA GPUs (V100) on [OLCF Summit](https://www.olcf.ornl.gov/summit/) 
 * Near compute roofline behavior
 * RDMA (remote data memory access; GPU-GPU direct communication) via GPU-aware MPI on NVIDIA (CUDA-aware MPI) and AMD GPU systems
@@ -197,7 +199,7 @@ They are organized below.
 
 * [Fypp](https://fypp.readthedocs.io/en/stable/fypp.html) metaprogramming for code readability, performance, and portability
 * Continuous Integration (CI)
-	* \>300 Regression tests with each PR.
+	* Approx. 500 Regression tests with each PR.
  		* Performed with GNU (GCC), Intel (oneAPI), Cray (CCE), and NVIDIA (NVHPC) compilers on NVIDIA and AMD GPUs.
 		* Line-level test coverage reports via [Codecov](https://app.codecov.io/gh/MFlowCode/MFC) and `gcov`
 	* Benchmarking to avoid performance regressions and identify speed-ups

docs/documentation/case.md

@@ -379,6 +379,7 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](r
 | `mixture_err`          | Logical | Mixture properties correction |
 | `time_stepper`         | Integer | Runge--Kutta order [1-3] |
 | `adap_dt`              | Logical | Strang splitting scheme with adaptive time stepping |
+| `recon_type`           | Integer | Reconstruction Type: [1] WENO; [2] MUSCL |
 | `adap_dt_tol`          | Real    | Tolerance for adaptive time stepping in Strang splitting scheme|
 | `adap_dt_max_iters`    | Integer | Max iteration for adaptive time stepping in Strang splitting scheme |
 | `weno_order`	         | Integer | WENO order [1,3,5] |
@@ -390,6 +391,11 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](r
 | `teno_CT`              | Real    | TENO threshold for smoothness detection |
 | `null_weights`         | Logical | Null WENO weights at boundaries |
 | `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`             | 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) |
 | `low_Mach`             | Integer | Low Mach number correction for HLLC Riemann solver: [0] None; [1] Pressure (Chen et al. 2022); [2] Velocity (Thornber et al. 2008)	 |
 | `avg_state`	         | Integer | Averaged state evaluation method: [1] Roe average*; [2] Arithmetic mean  |
@@ -476,7 +482,14 @@ It is recommended to set `weno_eps` to $10^{-6}$ for WENO-JS, and to $10^{-40}$
 
 - `mp_weno` activates monotonicity preservation in the WENO reconstruction (MPWENO) such that the values of reconstructed variables do not reside outside the range spanned by WENO stencil ([Balsara and Shu, 2000](references.md); [Suresh and Huynh, 1997](references.md)).
 
-- `riemann_solver` specifies the choice of the Riemann solver that is used in simulation by an integer from 1 through 3.
+- `muscl_order` specifies the order of the MUSCL scheme that is used for spatial reconstruction of variables by an integer of 1, or 2, that corresponds to the 1st, and 2nd order respectively. When using `muscl_order = 2`, `muscl_lim` must be defined. 
+
+- `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`).
+
+- `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 ([Toro, 2013](references.md)).
 `riemann_solver = 4` is only for MHD simulations. It resolves 5 of the full seven-wave structure of the MHD equations ([Miyoshi and Kusano, 2005](references.md)).
 

examples/1D_sodshocktube_muscl/case.py

@@ -0,0 +1,71 @@
+#!/usr/bin/env python3
+import math
+import json
+
+# Numerical setup
+Nx = 399
+dx = 1.0 / (1.0 * (Nx + 1))
+
+Tend, Nt = 0.1, 1000
+mydt = Tend / (1.0 * Nt)
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": 0.0e00,
+            "x_domain%end": 1.0e00,
+            "m": Nx,
+            "n": 0,
+            "p": 0,
+            "dt": mydt,
+            "t_step_start": 0,
+            "t_step_stop": int(Nt),
+            "t_step_save": int(math.ceil(Nt / 10.0)),
+            # Simulation Algorithm Parameters
+            "num_patches": 2,
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 1,
+            "mpp_lim": "F",
+            "mixture_err": "F",
+            "time_stepper": 3,
+            "recon_type": 2,
+            "muscl_order": 2,
+            "muscl_lim": 2,
+            "int_comp": "T",
+            "riemann_solver": 2,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -3,
+            "bc_x%end": -3,
+            # Formatted Database Files Structure Parameters
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "parallel_io": "T",
+            # Patch 1 L
+            "patch_icpp(1)%geometry": 1,
+            "patch_icpp(1)%x_centroid": 0.25,
+            "patch_icpp(1)%length_x": 0.5,
+            "patch_icpp(1)%vel(1)": 0.0,
+            "patch_icpp(1)%pres": 1.0,
+            "patch_icpp(1)%alpha_rho(1)": 1.0e00,
+            "patch_icpp(1)%alpha(1)": 1.0,
+            # Patch 2 R
+            "patch_icpp(2)%geometry": 1,
+            "patch_icpp(2)%x_centroid": 0.75,
+            "patch_icpp(2)%length_x": 0.5,
+            "patch_icpp(2)%vel(1)": 0.0,
+            "patch_icpp(2)%pres": 0.1,
+            "patch_icpp(2)%alpha_rho(1)": 0.125e00,
+            "patch_icpp(2)%alpha(1)": 1.0,
+            # Fluids Physical Parameters
+            "fluid_pp(1)%gamma": 1.0e00 / (1.4 - 1.0e00),
+            "fluid_pp(1)%pi_inf": 0.0,
+        }
+    )
+)

examples/2D_advection_muscl/case.py

@@ -0,0 +1,83 @@
+#!/usr/bin/env python3
+import json
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": 0.0e00,
+            "x_domain%end": 1.0e00,
+            "y_domain%beg": 0.0e00,
+            "y_domain%end": 1.0e00,
+            "m": 99,
+            "n": 99,
+            "p": 0,
+            "dt": 5.0e-07,
+            "t_step_start": 0,
+            "t_step_stop": 1000,
+            "t_step_save": 100,
+            # Simulation Algorithm Parameters
+            "num_patches": 2,
+            "model_eqns": 3,
+            "alt_soundspeed": "F",
+            "num_fluids": 2,
+            "mpp_lim": "T",
+            "mixture_err": "T",
+            "time_stepper": 3,
+            "recon_type": 2,
+            "muscl_order": 2,
+            "muscl_lim": 2,
+            "int_comp": "T",
+            "null_weights": "F",
+            "riemann_solver": 2,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -3,
+            "bc_x%end": -3,
+            "bc_y%beg": -3,
+            "bc_y%end": -3,
+            # Formatted Database Files Structure Parameters
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "parallel_io": "T",
+            # Patch 1: Base
+            "patch_icpp(1)%geometry": 3,
+            "patch_icpp(1)%x_centroid": 0.5e00,
+            "patch_icpp(1)%y_centroid": 0.5e00,
+            "patch_icpp(1)%length_x": 1.0e00,
+            "patch_icpp(1)%length_y": 1.0e00,
+            "patch_icpp(1)%vel(1)": 100.0e00,
+            "patch_icpp(1)%vel(2)": 100.0e00,
+            "patch_icpp(1)%pres": 1.0e05,
+            "patch_icpp(1)%alpha_rho(1)": 1000.0e00,
+            "patch_icpp(1)%alpha_rho(2)": 1.0,
+            "patch_icpp(1)%alpha(1)": 1.0e-12,
+            "patch_icpp(1)%alpha(2)": 1.0 - 1.0e-12,
+            # Patch 2: Density to transport
+            "patch_icpp(2)%geometry": 2,
+            "patch_icpp(2)%smoothen": "T",
+            "patch_icpp(2)%smooth_patch_id": 1,
+            "patch_icpp(2)%smooth_coeff": 0.5e00,
+            "patch_icpp(2)%x_centroid": 0.1e00,
+            "patch_icpp(2)%y_centroid": 0.1e00,
+            "patch_icpp(2)%radius": 0.1e00,
+            "patch_icpp(2)%alter_patch(1)": "T",
+            "patch_icpp(2)%vel(1)": 100.0e00,
+            "patch_icpp(2)%vel(2)": 100.0e00,
+            "patch_icpp(2)%pres": 1.0e05,
+            "patch_icpp(2)%alpha_rho(1)": 1.0,
+            "patch_icpp(2)%alpha_rho(2)": 1.0,
+            "patch_icpp(2)%alpha(1)": 0,
+            "patch_icpp(2)%alpha(2)": 1.0,
+            # Fluids Physical Parameters
+            "fluid_pp(1)%gamma": 1.0e00 / (2.35e00 - 1.0e00),
+            "fluid_pp(1)%pi_inf": 2.35e00 * 1.0e09 / (2.35e00 - 1.0e00),
+            "fluid_pp(2)%gamma": 1.0e00 / (1.4e00 - 1.0e00),
+            "fluid_pp(2)%pi_inf": 0.0e00,
+        }
+    )
+)

examples/2D_riemann_test_muscl/case.py

@@ -0,0 +1,99 @@
+#!/usr/bin/env python3
+import json
+import math
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": 0.0,
+            "x_domain%end": 1.0,
+            "y_domain%beg": 0.0,
+            "y_domain%end": 1.0,
+            "m": 499,
+            "n": 499,
+            "p": 0,
+            "dt": 8e-05,
+            "t_step_start": 0,
+            "t_step_stop": 1000,
+            "t_step_save": 100,
+            # Simulation Algorithm Parameters
+            "num_patches": 4,
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 1,
+            "mpp_lim": "F",
+            "mixture_err": "T",
+            "time_stepper": 3,
+            # "mp_weno": "F",
+            "recon_type": 2,
+            # "weno_order": 5,
+            # "weno_eps": 1e-16,
+            "muscl_order": 2,
+            "muscl_lim": 1,
+            "int_comp": "T",
+            "riemann_solver": 2,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -3,
+            "bc_x%end": -3,
+            "bc_y%beg": -3,
+            "bc_y%end": -3,
+            # Formatted Database Files Structure Parameters
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "parallel_io": "T",
+            # Patch 1: Base
+            "patch_icpp(1)%geometry": 3,
+            "patch_icpp(1)%x_centroid": 0.4,
+            "patch_icpp(1)%y_centroid": 0.4,
+            "patch_icpp(1)%length_x": 0.8,
+            "patch_icpp(1)%length_y": 0.8,
+            "patch_icpp(1)%vel(1)": 4 / math.sqrt(11),
+            "patch_icpp(1)%vel(2)": 4 / math.sqrt(11),
+            "patch_icpp(1)%pres": 9 / 310,
+            "patch_icpp(1)%alpha_rho(1)": 77 / 558,
+            "patch_icpp(1)%alpha(1)": 1,
+            # Patch 1: Base
+            "patch_icpp(2)%geometry": 3,
+            "patch_icpp(2)%x_centroid": 0.4,
+            "patch_icpp(2)%y_centroid": 0.9,
+            "patch_icpp(2)%length_x": 0.8,
+            "patch_icpp(2)%length_y": 0.2,
+            "patch_icpp(2)%vel(1)": 4 / math.sqrt(11),
+            "patch_icpp(2)%vel(2)": 0,
+            "patch_icpp(2)%pres": 0.3,
+            "patch_icpp(2)%alpha_rho(1)": 33 / 62,
+            "patch_icpp(2)%alpha(1)": 1,
+            # Patch 1: Base
+            "patch_icpp(3)%geometry": 3,
+            "patch_icpp(3)%x_centroid": 0.9,
+            "patch_icpp(3)%y_centroid": 0.4,
+            "patch_icpp(3)%length_x": 0.2,
+            "patch_icpp(3)%length_y": 0.8,
+            "patch_icpp(3)%vel(1)": 0,
+            "patch_icpp(3)%vel(2)": 4 / math.sqrt(11),
+            "patch_icpp(3)%pres": 0.3,
+            "patch_icpp(3)%alpha_rho(1)": 33 / 62,
+            "patch_icpp(3)%alpha(1)": 1,
+            # Patch 1: Base
+            "patch_icpp(4)%geometry": 3,
+            "patch_icpp(4)%x_centroid": 0.9,
+            "patch_icpp(4)%y_centroid": 0.9,
+            "patch_icpp(4)%length_x": 0.2,
+            "patch_icpp(4)%length_y": 0.2,
+            "patch_icpp(4)%vel(1)": 0,
+            "patch_icpp(4)%vel(2)": 0,
+            "patch_icpp(4)%pres": 1.5,
+            "patch_icpp(4)%alpha_rho(1)": 1.5,
+            "patch_icpp(4)%alpha(1)": 1.0,
+            # Fluids Physical Parameters
+            "fluid_pp(1)%gamma": 1.0e00 / (1.4e00 - 1.0e00),
+            "fluid_pp(1)%pi_inf": 0.0e00,
+        }
+    )
+)