Adding MUSCL Reconstruction

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 |
 | `weno_order`	         | Integer | WENO order [1,3,5] |
 | `weno_eps`	         | Real    | WENO perturbation (avoid division by zero) |
 | `mapped_weno`	         | Logical | WENO-M (WENO with mapping of nonlinear weights) |
@@ -388,6 +389,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  |
@@ -474,7 +480,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,
+        }
+    )
+)