MHD & RMHD

docs/documentation/case.md

@@ -372,7 +372,7 @@ 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 |
-| `riemann_solver`       | Integer | Riemann solver algorithm: [1] HLL*; [2] HLLC; [3] Exact*	 |
+| `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  |
 | `wave_speeds`          | Integer | Wave-speed estimation: [1] Direct (Batten et al. 1997); [2] Pressure-velocity* (Toro 1999)	 |
@@ -450,6 +450,7 @@ It is recommended to set `weno_eps` to $10^{-6}$ for WENO-JS, and to $10^{-40}$
 
 - `riemann_solver` specifies the choice of the Riemann solver that is used in simulation by an integer from 1 through 3.
 `riemann_solver = 1`, `2`, and `3` correspond to HLL, HLLC, and Exact Riemann solver, respectively ([Toro, 2013](references.md#Toro13)).
+`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#Miyoshi05)).
 
 - `low_Mach` specifies the choice of the low Mach number correction scheme for the HLLC Riemann solver. `low_Mach = 0` is default value and does not apply any correction scheme. `low_Mach = 1` and `2` apply the anti-dissipation pressure correction method ([Chen et al., 2022](references.md#Chen22)) and the improved velocity reconstruction method ([Thornber et al., 2008](references.md#Thornber08)). This feature requires `riemann_solver = 2` and `model_eqns = 2`.
 
@@ -848,7 +849,34 @@ $$ a_{x[y,z]} = g_{x[y,z]} + k_{x[y,z]}\sin\left(w_{x[y,z]}t + p_{x[y,z]}\right)
 
 By convention, positive accelerations in the `x[y,z]` direction are in the positive `x[y,z]` direction.
 
-### 14. Cylindrical Coordinates
+### 14. Magnetohydrodynamics (MHD)
+
+| Parameter         | Type    | Description                                         |
+| ---:              | :---:   | :---                                                |
+| `mhd`             | Logical | Enable ideal MHD simulation                         |
+| `relativity`      | Logical | Enable relativistic MHD simulation                  |
+| `powell`          | Logical | Enable Powell's method for solenoidal constraint    |
+| `fd_order`        | Integer | Finite difference order for Powell's method         |
+| `Bx[y,z]`         | Real    | Initial magnetic field in the x[y,z] direction      |
+| `Bx0`             | Real    | Constant magnetic field in the x direction (1D only)|
+
+- `mhd` is currently only available for single-component flows and 5-equation model. Its compatibility with most other features is work in progress.
+
+- `relativity` only works for `mhd` enabled and activates relativistic MHD (RMHD) simulation.
+
+- `powell` activates Powell's eight-wave method to impose the solenoidal constraint in the MHD simulation [Powell (1994)](references.md#Powell94). It should not be used in conjunction with HLLD (`riemann_solver = 4`).
+
+- `fd_order` specifies the finite difference order for computing the RHS of the Powell's method. `fd_order = 1`, `2`, and `4` are allowed.
+
+- `Bx0` is only used in 1D simulations to specify the constant magnetic field in the x direction. It must be specified in 1D simulations. `Bx` must not be used in 1D simulations.
+
+- `Bx`, `By`, and `Bz` are used to specify the initial magnetic field in the x, y, and z directions, respectively. They must be specified in all 1D/2D/3D MHD simulations, with the exception of `Bx` in 1D simulations.
+
+Note: In 1D/2D/3D simulations, all three velocity components are treated as state variables and must be specified in the case file.
+
+Note: For relativistic flow, the conservative and primitive densities are different. `rho_wrt` outputs the primitive (rest mass) density.
+
+### 15. Cylindrical Coordinates
 
 When ``cyl_coord = 'T'`` is set in 3D the following constraints must be met:
 

docs/documentation/references.md

@@ -57,3 +57,7 @@
 - <a id="Tiwari13">Tiwari, A., Freund, J. B., and Pantano, C. (2013). A diffuse interface model with immiscibility preservation. Journal of computational physics, 252:290–309.</a>
 
 - <a id="Toro13">Toro, E. F. (2013). Riemann solvers and numerical methods for fluid dynamics: a practical introduction. Springer Science & Business Media.</a>
+
+- <a id="Miyoshi05">Miyoishi, T., and Kusano, K. (2005). A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics. Journal of Computational Physics, 208(1), 315-344.</a>
+
+- <a id="Powell94">Powell, K. G. (1994). An approximate Riemann solver for magnetohydrodynamics: (That works in more than one dimension). In Upwind and high-resolution schemes (pp. 570-583). Springer.</a>

examples/1D_brio_wu/case.py

@@ -0,0 +1,79 @@
+#!/usr/bin/env python3
+import json
+
+# An upwind differencing scheme for the equations of ideal magnetohydrodynamics
+# M. Brio and C. C. Wu
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": 0,
+            "x_domain%end": 1.0,
+            "m": 800,
+            "n": 0,
+            "p": 0,
+            "dt": 0.0002,
+            "t_step_start": 0,
+            "t_step_stop": 1000,
+            "t_step_save": 100,
+            # Simulation Algorithm Parameters
+            "num_patches": 2,
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 1,
+            "mpp_lim": "F",
+            "mixture_err": "F",
+            "time_stepper": 3,
+            "weno_order": 5,
+            "mapped_weno": "T",
+            "weno_eps": 1.0e-16,
+            "null_weights": "F",
+            "mp_weno": "F",
+            "riemann_solver": 1,
+            "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",
+            "rho_wrt": "T",
+            "parallel_io": "T",
+            # MHD
+            "mhd": "T",
+            "Bx0": 0.75,
+            # Patch 1 Left
+            "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)%vel(2)": 0.0,
+            "patch_icpp(1)%vel(3)": 0.0,
+            "patch_icpp(1)%pres": 1.0,
+            "patch_icpp(1)%By": 1.0,
+            "patch_icpp(1)%Bz": 0.0,
+            "patch_icpp(1)%alpha_rho(1)": 1.0,
+            "patch_icpp(1)%alpha(1)": 1.0,
+            # Patch 2 Right
+            "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)%vel(2)": 0.0,
+            "patch_icpp(2)%vel(3)": 0.0,
+            "patch_icpp(2)%pres": 0.1,
+            "patch_icpp(2)%By": -1.0,
+            "patch_icpp(2)%Bz": 0.0,
+            "patch_icpp(2)%alpha_rho(1)": 0.125,
+            "patch_icpp(2)%alpha(1)": 1.0,
+            # Fluids Physical Parameters
+            "fluid_pp(1)%gamma": 1.0e00 / (2.0e00 - 1.0e00),
+            "fluid_pp(1)%pi_inf": 0.0,
+        }
+    )
+)

examples/1D_brio_wu_hlld/case.py

@@ -0,0 +1,80 @@
+#!/usr/bin/env python3
+import json
+
+# An upwind differencing scheme for the equations of ideal magnetohydrodynamics
+# M. Brio and C. C. Wu
+# Note: HLLD is not used in the paper
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": 0,
+            "x_domain%end": 1.0,
+            "m": 800,
+            "n": 0,
+            "p": 0,
+            "dt": 0.0002,
+            "t_step_start": 0,
+            "t_step_stop": 1000,
+            "t_step_save": 100,
+            # Simulation Algorithm Parameters
+            "num_patches": 2,
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 1,
+            "mpp_lim": "F",
+            "mixture_err": "F",
+            "time_stepper": 3,
+            "weno_order": 5,
+            "mapped_weno": "T",
+            "weno_eps": 1.0e-16,
+            "null_weights": "F",
+            "mp_weno": "F",
+            "riemann_solver": 4,
+            "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",
+            "rho_wrt": "T",
+            "parallel_io": "T",
+            # MHD
+            "mhd": "T",
+            "Bx0": 0.75,
+            # Patch 1 Left
+            "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)%vel(2)": 0.0,
+            "patch_icpp(1)%vel(3)": 0.0,
+            "patch_icpp(1)%pres": 1.0,
+            "patch_icpp(1)%By": 1.0,
+            "patch_icpp(1)%Bz": 0.0,
+            "patch_icpp(1)%alpha_rho(1)": 1.0,
+            "patch_icpp(1)%alpha(1)": 1.0,
+            # Patch 2 Right
+            "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)%vel(2)": 0.0,
+            "patch_icpp(2)%vel(3)": 0.0,
+            "patch_icpp(2)%pres": 0.1,
+            "patch_icpp(2)%By": -1.0,
+            "patch_icpp(2)%Bz": 0.0,
+            "patch_icpp(2)%alpha_rho(1)": 0.125,
+            "patch_icpp(2)%alpha(1)": 1.0,
+            # Fluids Physical Parameters
+            "fluid_pp(1)%gamma": 1.0e00 / (2.0e00 - 1.0e00),
+            "fluid_pp(1)%pi_inf": 0.0,
+        }
+    )
+)

examples/1D_brio_wu_rmhd/case.py

@@ -0,0 +1,81 @@
+#!/usr/bin/env python3
+import json
+
+# An HLLC Riemann solver for relativistic flows – II. Magnetohydrodynamics
+# A. Mignone and G. Bodo
+# Problem 1 Figure 3
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": 0,
+            "x_domain%end": 1.0,
+            "m": 1600,
+            "n": 0,
+            "p": 0,
+            "dt": 0.0002,
+            "t_step_start": 0,
+            "t_step_stop": 2000,
+            "t_step_save": 100,
+            # Simulation Algorithm Parameters
+            "num_patches": 2,
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 1,
+            "mpp_lim": "F",
+            "mixture_err": "F",
+            "time_stepper": 3,
+            "weno_order": 5,
+            "mapped_weno": "T",
+            "weno_eps": 1.0e-16,
+            "null_weights": "F",
+            "mp_weno": "F",
+            "riemann_solver": 1,
+            "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",
+            "rho_wrt": "T",
+            "parallel_io": "T",
+            # RMHD
+            "mhd": "T",
+            "relativity": "T",
+            "Bx0": 0.5,
+            # Patch 1 Left
+            "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)%vel(2)": 0.0,
+            "patch_icpp(1)%vel(3)": 0.0,
+            "patch_icpp(1)%pres": 1.0,
+            "patch_icpp(1)%By": 1.0,
+            "patch_icpp(1)%Bz": 0.0,
+            "patch_icpp(1)%alpha_rho(1)": 1.0,
+            "patch_icpp(1)%alpha(1)": 1.0,
+            # Patch 2 Right
+            "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)%vel(2)": 0.0,
+            "patch_icpp(2)%vel(3)": 0.0,
+            "patch_icpp(2)%pres": 0.1,
+            "patch_icpp(2)%By": -1.0,
+            "patch_icpp(2)%Bz": 0.0,
+            "patch_icpp(2)%alpha_rho(1)": 0.125,
+            "patch_icpp(2)%alpha(1)": 1.0,
+            # Fluids Physical Parameters
+            "fluid_pp(1)%gamma": 1.0e00 / (2.0e00 - 1.0e00),
+            "fluid_pp(1)%pi_inf": 0.0,
+        }
+    )
+)

examples/1D_dai_woodward/case.py

@@ -0,0 +1,88 @@
+#!/usr/bin/env python3
+import json
+import math
+
+# A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics
+# T. Miyoshi and K. Kusano
+# Section 6.1 Figure 1
+
+# Normalize magnetic field
+sqrt_4pi = math.sqrt(4 * math.pi)
+Bx0 = 2.0 / sqrt_4pi  # Typo in Miyoshi paper
+By_left = 3.6 / sqrt_4pi
+Bz_left = 2.0 / sqrt_4pi
+By_right = 4.0 / sqrt_4pi
+Bz_right = 2.0 / sqrt_4pi
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": 0,
+            "x_domain%end": 1.0,
+            "m": 800,
+            "n": 0,
+            "p": 0,
+            "dt": 0.0002,
+            "t_step_start": 0,
+            "t_step_stop": 1000,
+            "t_step_save": 100,
+            # Simulation Algorithm Parameters
+            "num_patches": 2,
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 1,
+            "mpp_lim": "F",
+            "mixture_err": "F",
+            "time_stepper": 3,
+            "weno_order": 1,
+            "weno_eps": 1.0e-16,
+            "null_weights": "F",
+            "mp_weno": "F",
+            "riemann_solver": 1,
+            "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",
+            "rho_wrt": "T",
+            "parallel_io": "T",
+            # MHD Settings
+            "mhd": "T",
+            "Bx0": Bx0,
+            # Patch 1 (Left State)
+            "patch_icpp(1)%geometry": 1,
+            "patch_icpp(1)%x_centroid": 0.25,
+            "patch_icpp(1)%length_x": 0.5,
+            "patch_icpp(1)%vel(1)": 1.2,
+            "patch_icpp(1)%vel(2)": 0.01,
+            "patch_icpp(1)%vel(3)": 0.5,
+            "patch_icpp(1)%pres": 0.95,
+            "patch_icpp(1)%By": By_left,
+            "patch_icpp(1)%Bz": Bz_left,
+            "patch_icpp(1)%alpha_rho(1)": 1.08,
+            "patch_icpp(1)%alpha(1)": 1.0,
+            # Patch 2 (Right State)
+            "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)%vel(2)": 0.0,
+            "patch_icpp(2)%vel(3)": 0.0,
+            "patch_icpp(2)%pres": 1.0,
+            "patch_icpp(2)%By": By_right,
+            "patch_icpp(2)%Bz": Bz_right,
+            "patch_icpp(2)%alpha_rho(1)": 1.0,
+            "patch_icpp(2)%alpha(1)": 1.0,
+            # Fluids Physical Parameters
+            "fluid_pp(1)%gamma": 1.0e00 / ((5.0 / 3.0) - 1.0e00),
+            "fluid_pp(1)%pi_inf": 0.0,
+        },
+    )
+)