MHD Hyperbolic Divergence Cleaning (#1086)

Co-authored-by: Daniel J. Vickers <dnlvickers5@gmail.com>
Co-authored-by: Spencer Bryngelson <sbryngelson@gmail.com>

Author: 3 people

.typos.toml

@@ -21,6 +21,7 @@ TKE = "TKE"
 HSA = "HSA"
 infp = "infp"
 Sur = "Sur"
+iz = "iz"
 
 [files]
 extend-exclude = ["docs/documentation/references*", "tests/", "toolchain/cce_simulation_workgroup_256.sh"]

docs/documentation/case.md

@@ -939,22 +939,25 @@ By convention, positive accelerations in the `x[y,z]` direction are in the posit
 
 ### 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)|
+| Parameter                | Type    | Description                                              |
+| ---:                     | :---:   | :---                                                     |
+| `mhd`                    | Logical | Enable ideal MHD simulation                              |
+| `relativity`             | Logical | Enable relativistic MHD simulation                       |
+| `hyper_cleaning`         | Logical | Enable hyperbolic (GLM) divergence cleaning for `div B`  |
+| `hyper_cleaning_speed`   | Real    | Cleaning wave speed `c_h`                                |
+| `hyper_cleaning_tau`     | Real    | Cleaning damping timescale `tau`                         |
+| `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). It should not be used in conjunction with HLLD (`riemann_solver = 4`).
+- `hyper_cleaning` [Dedner et al., 2002](references.md) only works with `mhd` in 2D/3D and reduces numerical `div B` errors by propagation and damping. Currently not compatible with HLLD (`riemann_solver = 4`).
+
+- `hyper_cleaning_speed` sets the propagation speed of divergence-cleaning waves.
 
-- `fd_order` specifies the finite difference order for computing the RHS of the Powell's method. `fd_order = 1`, `2`, and `4` are allowed.
+- `hyper_cleaning_tau` sets the decay timescale for divergence-cleaning.
 
 - `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.
 

docs/documentation/references.md

@@ -66,7 +66,7 @@
 
 - <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>
+- <a id="Dedner02">Dedner, A., Kemm, F., Kröner, D., Munz, C. D., Schnitzer, T., & Wesenberg, M. (2002). Hyperbolic divergence cleaning for the MHD equations. Journal of Computational Physics, 175(2), 645-673.</a>
 
 - <a id="Cao19">Cao, S., Zhang, Y., Liao, D., Zhong, P., and Wang, K. G. (2019). Shock-induced damage and dynamic fracture in cylindrical bodies submerged in liquid. International Journal of Solids and Structures, 169:55–71. Elsevier.</a>
 

examples/2D_mhd_magnetic_vortex/case.py

@@ -54,7 +54,7 @@
             "parallel_io": "T",
             # MHD
             "mhd": "T",
-            "patch_icpp(1)%hcid": 252,
+            "patch_icpp(1)%hcid": 253,
             "patch_icpp(1)%geometry": 3,
             "patch_icpp(1)%x_centroid": 0.0,
             "patch_icpp(1)%y_centroid": 0.0,

examples/2D_mhd_rotor/case.py

@@ -0,0 +1,84 @@
+#!/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": 399,
+            "n": 399,
+            "p": 0,
+            "dt": 0.0004,
+            "t_step_start": 0,
+            "t_step_stop": 375,
+            "t_step_save": 5,
+            # Simulation Algorithm Parameters
+            "num_patches": 1,
+            "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-6,
+            "null_weights": "F",
+            "mp_weno": "F",
+            "riemann_solver": 1,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -1,
+            "bc_x%end": -1,
+            "bc_y%beg": -1,
+            "bc_y%end": -1,
+            # Formatted Database Files Structure Parameters
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "rho_wrt": "T",
+            "parallel_io": "T",
+            # MHD
+            "mhd": "T",
+            "hyper_cleaning": "T",
+            "hyper_cleaning_speed": 2.5,
+            "hyper_cleaning_tau": 0.004,
+            # Patch 1 - 2D MHD Rotor Problem
+            # gamma = 1.4
+            # Ambient medium (r > 0.1):
+            #   rho = 1
+            #   p = 1
+            #   v = (0, 0, 0)
+            #   B = (1, 0, 0)
+            # Rotor (r <= 0.1):
+            #   rho = 10
+            #   v has angular velocity of 20
+            "patch_icpp(1)%hcid": 252,
+            "patch_icpp(1)%geometry": 3,
+            "patch_icpp(1)%x_centroid": 0.5,
+            "patch_icpp(1)%y_centroid": 0.5,
+            "patch_icpp(1)%length_x": 1.0,
+            "patch_icpp(1)%length_y": 1.0,
+            "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)%Bx": 5.0 / math.sqrt(math.pi * 4.0),
+            "patch_icpp(1)%By": 0.0,
+            "patch_icpp(1)%Bz": 0.0,
+            "patch_icpp(1)%alpha_rho(1)": 1.0,
+            "patch_icpp(1)%alpha(1)": 1.0,
+            # Fluids Physical Parameters
+            "fluid_pp(1)%gamma": 1.0e00 / (1.4e00 - 1.0e00),
+            "fluid_pp(1)%pi_inf": 0.0,
+        }
+    )
+)

examples/2D_orszag_tang_powell/case.py …es/2D_orszag_tang_hyper_cleaning/case.pyexamples/2D_orszag_tang_powell/case.py renamed to examples/2D_orszag_tang_hyper_cleaning/case.py examples/2D_orszag_tang_powell/case.py renamed to examples/2D_orszag_tang_hyper_cleaning/case.py

@@ -16,9 +16,9 @@
             "m": 512,
             "n": 512,
             "p": 0,
-            "dt": 0.0005,
+            "dt": 0.0004,
             "t_step_start": 0,
-            "t_step_stop": 2000,
+            "t_step_stop": 2500,
             "t_step_save": 100,
             # Simulation Algorithm Parameters
             "num_patches": 1,
@@ -28,7 +28,8 @@
             "mpp_lim": "F",
             "mixture_err": "F",
             "time_stepper": 3,
-            "weno_order": 1,
+            "weno_order": 5,
+            "mapped_weno": "T",
             "weno_eps": 1.0e-6,
             "null_weights": "F",
             "mp_weno": "F",
@@ -47,8 +48,9 @@
             "parallel_io": "T",
             # MHD
             "mhd": "T",
-            "powell": "T",
-            "fd_order": 2,
+            "hyper_cleaning": "T",
+            "hyper_cleaning_speed": 2.5,
+            "hyper_cleaning_tau": 0.004,
             # Patch 1 - Analytical for v and B
             # gamma = 5/3
             #   rho = 25/(36*pi)

src/common/include/2dHardcodedIC.fpp

@@ -1,9 +1,14 @@
 #:def Hardcoded2DVariables()
     ! Place any declaration of intermediate variables here
-    real(wp) :: eps
+    real(wp) :: eps, eps_mhd, C_mhd
     real(wp) :: r, rmax, gam, umax, p0
     real(wp) :: rhoH, rhoL, pRef, pInt, h, lam, wl, amp, intH, intL, alph
     real(wp) :: factor
+    real(wp) :: r0, alpha, r2
+    real(wp) :: sinA, cosA
+
+    real(wp) :: r_sq
+
     ! # 207
     real(wp) :: sigma, gauss1, gauss2
     ! # 208
@@ -183,7 +188,43 @@
             q_prim_vf(E_idx)%sf(i, j, 0) = 3.e-5_wp
         end if
 
-    case (252) ! MHD Smooth Magnetic Vortex
+        ! case 252 is for the 2D MHD Rotor problem
+    case (252) ! 2D MHD Rotor Problem
+        ! Ambient conditions are set in the JSON file.
+        ! This case imposes the dense, rotating cylinder.
+        !
+        ! gamma = 1.4
+        ! Ambient medium (r > 0.1):
+        !   rho = 1, p = 1, v = 0, B = (1,0,0)
+        ! Rotor (r <= 0.1):
+        !   rho = 10, p = 1
+        !   v has angular velocity w=20, giving v_tan=2 at r=0.1
+
+        ! Calculate distance squared from the center
+        r_sq = (x_cc(i) - 0.5_wp)**2 + (y_cc(j) - 0.5_wp)**2
+
+        ! inner radius of 0.1
+        if (r_sq <= 0.1**2) then
+            ! -- Inside the rotor --
+            ! Set density uniformly to 10
+            q_prim_vf(contxb)%sf(i, j, 0) = 10._wp
+
+            ! Set vup constant rotation of rate v=2
+            ! v_x = -omega * (y - y_c)
+            ! v_y =  omega * (x - x_c)
+            q_prim_vf(momxb)%sf(i, j, 0) = -20._wp*(y_cc(j) - 0.5_wp)
+            q_prim_vf(momxb + 1)%sf(i, j, 0) = 20._wp*(x_cc(i) - 0.5_wp)
+
+            ! taper width of 0.015
+        else if (r_sq <= 0.115**2) then
+            ! linearly smooth the function between r = 0.1 and 0.115
+            q_prim_vf(contxb)%sf(i, j, 0) = 1._wp + 9._wp*(0.115_wp - sqrt(r_sq))/(0.015_wp)
+
+            q_prim_vf(momxb)%sf(i, j, 0) = -(2._wp/sqrt(r_sq))*(y_cc(j) - 0.5_wp)*(0.115_wp - sqrt(r_sq))/(0.015_wp)
+            q_prim_vf(momxb + 1)%sf(i, j, 0) = (2._wp/sqrt(r_sq))*(x_cc(i) - 0.5_wp)*(0.115_wp - sqrt(r_sq))/(0.015_wp)
+        end if
+
+    case (253) ! MHD Smooth Magnetic Vortex
         ! Section 5.2 of
         ! Implicit hybridized discontinuous Galerkin methods for compressible magnetohydrodynamics
         ! C. Ciuca, P. Fernandez, A. Christophe, N.C. Nguyen, J. Peraire
@@ -199,6 +240,62 @@
         ! pressure
         q_prim_vf(E_idx)%sf(i, j, 0) = 1._wp + (1 - 2._wp*(x_cc(i)**2 + y_cc(j)**2))*exp(1 - (x_cc(i)**2 + y_cc(j)**2))/((2._wp*pi)**3)
 
+    case (260)  ! Gaussian Divergence Pulse
+        !  Bx(x) = 1 + C * erf((x-0.5)/σ)
+        !  ⇒  ∂Bx/∂x = C * (2/√π) * exp[-((x-0.5)/σ)**2] * (1/σ)
+        !  Choose C = ε * σ * √π / 2  ⇒ ∂Bx/∂x = ε * exp[-((x-0.5)/σ)**2]
+        !  ψ is initialized to zero everywhere.
+
+        eps_mhd = patch_icpp(patch_id)%a(2)
+        sigma = patch_icpp(patch_id)%a(3)
+        C_mhd = eps_mhd*sigma*sqrt(pi)*0.5_wp
+
+        ! B-field
+        q_prim_vf(B_idx%beg)%sf(i, j, 0) = 1._wp + C_mhd*erf((x_cc(i) - 0.5_wp)/sigma)
+
+    case (261)  ! Blob
+        r0 = 1._wp/sqrt(8._wp)
+        r2 = x_cc(i)**2 + y_cc(j)**2
+        r = sqrt(r2)
+        alpha = r/r0
+        if (alpha < 1) then
+            q_prim_vf(B_idx%beg)%sf(i, j, 0) = 1._wp/sqrt(4._wp*pi)*(alpha**8 - 2._wp*alpha**4 + 1._wp)
+            ! q_prim_vf(B_idx%beg)%sf(i,j,0) = 1._wp/sqrt(4000._wp*pi) * (4096._wp*r2**4 - 128._wp*r2**2 + 1._wp)
+            ! q_prim_vf(B_idx%beg)%sf(i,j,0) = 1._wp/(4._wp*pi) * (alpha**8 - 2._wp*alpha**4 + 1._wp)
+            ! q_prim_vf(E_idx)%sf(i,j,0) = 6._wp - q_prim_vf(B_idx%beg)%sf(i,j,0)**2/2._wp
+        end if
+
+    case (262)  ! Tilted 2D MHD shock‐tube at α = arctan2 (≈63.4°)
+        ! rotate by α = atan(2)
+        alpha = atan(2._wp)
+        cosA = cos(alpha)
+        sinA = sin(alpha)
+        ! projection along shock normal
+        r = x_cc(i)*cosA + y_cc(j)*sinA
+
+        if (r <= 0.5_wp) then
+            ! LEFT state: ρ=1, v∥=+10, v⊥=0, p=20, B∥=B⊥=5/√(4π)
+            q_prim_vf(contxb)%sf(i, j, 0) = 1._wp
+            q_prim_vf(momxb)%sf(i, j, 0) = 10._wp*cosA
+            q_prim_vf(momxb + 1)%sf(i, j, 0) = 10._wp*sinA
+            q_prim_vf(E_idx)%sf(i, j, 0) = 20._wp
+            q_prim_vf(B_idx%beg)%sf(i, j, 0) = (5._wp/sqrt(4._wp*pi))*cosA &
+                                               - (5._wp/sqrt(4._wp*pi))*sinA
+            q_prim_vf(B_idx%beg + 1)%sf(i, j, 0) = (5._wp/sqrt(4._wp*pi))*sinA &
+                                                   + (5._wp/sqrt(4._wp*pi))*cosA
+        else
+            ! RIGHT state: ρ=1, v∥=−10, v⊥=0, p=1, B∥=B⊥=5/√(4π)
+            q_prim_vf(contxb)%sf(i, j, 0) = 1._wp
+            q_prim_vf(momxb)%sf(i, j, 0) = -10._wp*cosA
+            q_prim_vf(momxb + 1)%sf(i, j, 0) = -10._wp*sinA
+            q_prim_vf(E_idx)%sf(i, j, 0) = 1._wp
+            q_prim_vf(B_idx%beg)%sf(i, j, 0) = (5._wp/sqrt(4._wp*pi))*cosA &
+                                               - (5._wp/sqrt(4._wp*pi))*sinA
+            q_prim_vf(B_idx%beg + 1)%sf(i, j, 0) = (5._wp/sqrt(4._wp*pi))*sinA &
+                                                   + (5._wp/sqrt(4._wp*pi))*cosA
+        end if
+        ! v^z and B^z remain zero by default
+
     case (270)
         ! This hardcoded case extrudes a 1D profile to initialize a 2D simulation domain
         @: HardcodedReadValues()

src/common/m_variables_conversion.fpp

@@ -882,6 +882,7 @@ contains
 
                     if (cont_damage) qK_prim_vf(damage_idx)%sf(j, k, l) = qK_cons_vf(damage_idx)%sf(j, k, l)
 
+                    if (hyper_cleaning) qK_prim_vf(psi_idx)%sf(j, k, l) = qK_cons_vf(psi_idx)%sf(j, k, l)
 #ifdef MFC_POST_PROCESS
                     if (bubbles_lagrange) qK_prim_vf(beta_idx)%sf(j, k, l) = qK_cons_vf(beta_idx)%sf(j, k, l)
 #endif
@@ -1152,6 +1153,8 @@ contains
 
                     if (cont_damage) q_cons_vf(damage_idx)%sf(j, k, l) = q_prim_vf(damage_idx)%sf(j, k, l)
 
+                    if (hyper_cleaning) q_cons_vf(psi_idx)%sf(j, k, l) = q_prim_vf(psi_idx)%sf(j, k, l)
+
                 end do
             end do
         end do

src/post_process/m_data_output.fpp

@@ -363,6 +363,9 @@ contains
             ! Damage state variable
             if (cont_damage) dbvars = dbvars + 1
 
+            ! Hyperbolic cleaning for MHD
+            if (hyper_cleaning) dbvars = dbvars + 1
+
             ! Magnetic field
             if (mhd) then
                 if (n == 0) then

src/post_process/m_global_parameters.fpp

@@ -117,6 +117,7 @@ module m_global_parameters
     integer :: b_size          !< Number of components in the b tensor
     integer :: tensor_size     !< Number of components in the nonsymmetric tensor
     logical :: cont_damage     !< Continuum damage modeling
+    logical :: hyper_cleaning  !< Hyperbolic cleaning for MHD
     logical :: igr             !< enable IGR
     integer :: igr_order       !< IGR reconstruction order
     logical, parameter :: chemistry = .${chemistry}$. !< Chemistry modeling
@@ -143,6 +144,7 @@ module m_global_parameters
     integer :: c_idx                               !< Index of color function
     type(int_bounds_info) :: species_idx           !< Indexes of first & last concentration eqns.
     integer :: damage_idx                          !< Index of damage state variable (D) for continuum damage model
+    integer :: psi_idx                                 !< Index of hyperbolic cleaning state variable for MHD
     !> @}
 
     ! Cell Indices for the (local) interior points (O-m, O-n, 0-p).
@@ -407,6 +409,7 @@ contains
         b_size = dflt_int
         tensor_size = dflt_int
         cont_damage = .false.
+        hyper_cleaning = .false.
         igr = .false.
 
         bc_x%beg = dflt_int; bc_x%end = dflt_int
@@ -812,6 +815,13 @@ contains
                 damage_idx = dflt_int
             end if
 
+            if (hyper_cleaning) then
+                psi_idx = sys_size + 1
+                sys_size = psi_idx
+            else
+                psi_idx = dflt_int
+            end if
+
         end if
 
         if (chemistry) then