A short note on reconstruction variables in shock capturing schemes for magnetohydrodynamics
AA short note on reconstruction variables in shock capturingschemes for magnetohydrodynamics
Takahiro Miyoshi a, ∗ , Takashi Minoshima b a Graduate School of Advanced Science and Engineering, Hiroshima University, Higashihiroshima 739-8526, Japan b Center for Mathematical Science and Advanced Technology, Japan Agency for Marine-Earth Science and Technology,Yokohama 236-0001, Japan
Abstract
We propose a set of quick and easy approximate characteristic variables for higher-order recon-structions of shock capturing schemes for magnetohydrodynamics (MHD). Numerical experi-ments suggest that the reconstructions using the approximate characteristic variables are morerobust than those using the conservative or primitive variables, while their computational e ffi -ciencies are comparable. The approximate characteristic variables are simple compared to thefull characteristic variables for MHD, and can be a practical choice of reconstruction variables. Keywords: magnetohydrodynamics, shock capturing scheme, reconstruction
1. Introduction
High-order shock capturing schemes such as monotonic upstream-centered scheme for con-servation laws (MUSCL) [23], essentially non-oscillatory scheme (ENO) [9], weighted essen-tially non-oscillatory scheme (WENO) [12], and monotonicity-preserving scheme (MP) [21] aredesigned based on the theory of hyperbolic conservation laws. The reconstruction techniquesused in the aforementioned schemes extend a monotone first-order scheme to higher orders,while suppressing numerical oscillations near discontinuities. For a scalar conservation law, thevariable to be reconstructed is the conservative variable, which is identical to the characteristicvariable. By contrast, reconstruction variables in a system of nonlinear hyperbolic conservationlaws are not unique because the conservative variables and characteristic ones are di ff erent ingeneral.In order to avoid spurious numerical oscillations, the reconstruction should be performed inthe characteristic fields [9]. However, in the magnetohydrodynamic (MHD) system, the char-acteristic variables are much more complicated than those in the hydrodynamic system [5, 19],and as a result, the reconstruction of those demands high computational cost. For example, inan experiment using an open source MHD simulation code CANS + [13], where the fifth-orderMP of the characteristic variables is adopted as the reconstruction, the computational time ofthe reconstruction steps is more than 60% of the total time. To reduce the computational time,a particular high-order multidimensional scheme has been proposed that reduces the number of ∗ Corresponding author
Email address: [email protected] (Takahiro Miyoshi)
Preprint submitted to Elsevier September 3, 2020 a r X i v : . [ phy s i c s . c o m p - ph ] S e p ppearances of the characteristic variables while maintaining robustness [4]. Meanwhile, theconservative variables or the primitive variables have often been used instead of the character-istic variables due to computational e ffi ciency in many practical applications. In this paper, wepropose a set of quick and easy reconstruction variables for high-order MHD schemes, where thereconstructions using the proposed variables are more robust than those using the conservativeor primitive variables, while their e ffi ciencies are comparable.
2. Set of variables of the MHD equations
Consider the one-dimensional conservation laws: ∂ U ∂ t + ∂ F ∂ x = , (1)where U and F are the state vector of the conservative variables and corresponding flux vector.The set of the conservative variables of MHD is given by U = (cid:16) ρ, m x , m y , m z , B y , B z , e (cid:17) T , (2)where ρ , m , B , and e are the density, momentum, magnetic field, and total energy density, re-spectively. The subscripts x , y , and z denote the x -, y -, and z -components of vector fields. Inone-dimension, B x is constant in space and time due to the solenoidal condition of the magneticfield. The flux vector of the ideal MHD equations is F = m xm x ρ + p + B − B xm x ρ m y − B x B ym x ρ m z − B x B zm x ρ B y − B x m y ρ m x ρ B z − B x m z ρ m x ρ ( e + p + B / − B x ( m x B x + m y B y + m z B z ) ρ . (3)When an ideal gas equation of state is considered, the pressure p is determined from p = ( γ − (cid:32) e − m ρ − B (cid:33) , (4)where γ denotes the ratio of specific heats. The conservation laws of MHD can be rewritten as a quasilinear form for primitive variables V [18]: ∂ V ∂ t + A P ∂ V ∂ x = , (5)where A P is a coe ffi cient matrix. An often used set of the primitive variables is given by V = (cid:16) ρ, v x , v y , v z , B y , B z , p (cid:17) T , (6)2here v = m /ρ is the velocity. The temperature or the entropy density can become an alterna-tive to the pressure. Using the transformation matrix (see (19) and (20) in [18]), infinitesimalvariations of the conservative and primitive variables are related by d V = M − d U , M = ∂ U ∂ V . (7) Since the ideal MHD equations are hyperbolic, the equations are decomposed into a systemof nonlinear advection equations as R − P ∂ V ∂ t + R − P A P R P R − P ∂ V ∂ x ≡ ∂ W ∂ t + Λ ∂ W ∂ x = , (8)where R P and Λ are right eigenvectors of A P and diagonal matrix of the eigenvalues of A P ,respectively. The characteristic variables W are transformed into the others by d W = R − P d V = R − d U , (9)where R = MR P . The eigenvectors must be appropriately normalized so as to be well-behaved[5, 19] because the eigenvectors can be singular when the eigenvalues are degenerate. The trans-formation to and from the characteristic variables of the MHD equations is time-consumingcompared with that of the Euler equations and can be a computational bottleneck. Let us consider approximate subsystems reduced from the MHD equations. The compress-ible subsystem is obtained by taking the limit of B x → B y , B z , and p , we find ∂∂ t ρ v x p (cid:48) T + v x ρ v x ρ ρ c (cid:48) f v x ∂∂ x ρ v x p (cid:48) T = , (10)where p (cid:48) T = p + (cid:16) B y + B z (cid:17) , ρ c (cid:48) f = γ p + B x + B y + B z . (11)The second term on the right side of the last equation has been added to include a contributionof B x . The incompressible subsystems, on the other hand, are derived in the limit of ∂ v x /∂ x → ∂∂ t (cid:32) v y , z B y , z (cid:33) + (cid:32) v x − ρ B x − B x v x (cid:33) ∂∂ x (cid:32) v y , z B y , z (cid:33) = . (12)The subsystems (10) and (12) readily lead to the relations between the characteristic variables ofthe subsystems, W (cid:48) m , and the primitive variables: dW (cid:48) dW (cid:48) dW (cid:48) = c (cid:48) f − ρ c (cid:48) f − ρ c (cid:48) f d ρ dv x d p (cid:48) T , (cid:32) dW (cid:48) , dW (cid:48) , (cid:33) = (cid:32) √ ρ √ ρ − (cid:33) (cid:32) dv y , z dB y , z (cid:33) . (13)3eanwhile, the inverse transformations are given by d ρ dv x d p (cid:48) T = c (cid:48) f c (cid:48) f c (cid:48) f ρ c (cid:48) f − ρ c (cid:48) f
12 12 dW (cid:48) dW (cid:48) dW (cid:48) , (cid:32) dv y , z dB y , z (cid:33) = (cid:32) √ ρ √ ρ − (cid:33) (cid:32) dW (cid:48) , dW (cid:48) , (cid:33) . (14)Using the relation d p (cid:48) T = d p + B y dB y + B z dB z , d p = (cid:16) dW (cid:48) + dW (cid:48) − B y dW (cid:48) + B y dW (cid:48) − B z dW (cid:48) + B z dW (cid:48) (cid:17) . (15)Here W (cid:48) m is referred to as approximate characteristic variables. In vector form, the approximatecharacteristic variables are related to the conservative and primitive variables by d W (cid:48) = R (cid:48) P − d V = R (cid:48)− d U (16)where R (cid:48) = MR (cid:48) P . The transformation matrices R (cid:48) P − consisting of (13) and R (cid:48) P of (14) areconsiderably simplified compared to the eigenvectors of the complete MHD system.
3. Reconstruction procedures
In the finite volume approach, higher-order spatial accuracy can be achieved by reconstruct-ing the cell-averaged value and evaluating the values on the left and right faces of the cell. Sim-ilarly, in the finite di ff erence approach, higher-order accuracy is realized by reconstructing thepoint value and evaluating the values on the left and right sides of the midpoint.Consider first the scalar conservation law for u . Let u Lj + / denote u on the left of the cell faceor the midpoint, j + /
2. In general, introducing a well-controlled interpolation function I tosuppress numerical oscillations, u Lj + / = I (cid:16) · · · , u j − , u j , u j + , · · · (cid:17) , (17)where u j + s is the cell-averaged or the point value at j + s . Similarly, the right value u Rj − / isobtained by reversing the order of u j + s in I .Correspondingly, we evaluate the variables of the MHD on j + / I as the following subsections. Compute the conservative variables on the left U Lj + / using U j + s : U Lj + / = I (cid:16) · · · , U j − , U j , U j + , · · · (cid:17) , (18)where I acts on each component of the vector. Hereinafter, (18) is referred to as U - reconstruction . Compute the primitive variables V Lj + / using V j + s ≡ V ( U j + s ) as V Lj + / = I (cid:16) · · · , V j − , V j , V j + , · · · (cid:17) . (19)Hereinafter referred to as V - reconstruction . 4 .3. Characteristic variable reconstruction Compute V Lj + / by interpolating the characteristic variables applied with the eigenvectors forthe primitive system at j as [1, 2] V Lj + / = R P j I (cid:16) · · · , ¯ W j − , ¯ W j , ¯ W j + , · · · (cid:17) , ¯ W j + s = R P − j V j + s . (20)Hereinafter referred to as W - reconstruction . Note that, for example in the finite-di ff erenceWENO schemes [10, 3], the characteristic decomposition is performed using averaged eigen-vectors at the midpoint j + /
2. This paper, however, does not focus on that. Also, U Lj + / can bedirectly computed by applying R : U Lj + / = R j I (cid:16) · · · , ˜ W j − , ˜ W j , ˜ W j + , · · · (cid:17) , ˜ W j + s = R − j U j + s . (21)Although U Lj + / (cid:44) U ( V Lj + / ) in general, both lead to similar results in later tests (but not shown). Compute V Lj + / applying the transformation matrix R (cid:48) P to the approximate characteristicvariables W (cid:48) (referred to as W (cid:48) - reconstruction ): V Lj + / = R (cid:48) P j I (cid:16) · · · , ¯ W (cid:48) j − , ¯ W (cid:48) j , ¯ W (cid:48) j + , · · · (cid:17) , ¯ W (cid:48) j + s = R (cid:48) P j − V (cid:48) j + s . (22)Alternatively, compute U Lj + / applying the transformation matrix R (cid:48) to the approximate charac-teristic variables W (cid:48) (results not shown): U Lj + / = R (cid:48) j I (cid:16) · · · , ˜ W (cid:48) j − , ˜ W (cid:48) j , ˜ W (cid:48) j + , · · · (cid:17) , ˜ W (cid:48) j + s = R (cid:48) j − U (cid:48) j + s . (23)
4. Numerical experiments
Typical numerical experiments are presented. The Harten-Lax-van Leer Discontinuities(HLLD) approximate Riemann solver [16] is employed as a base scheme for MHD. The second-order MUSCL and the fifth-order weighted compact nonlinear scheme (WCNS), which is con-structed by combining the fourth-order interpolation and fourth-order finite di ff erence [14], areadopted as high-order schemes. The third-order TVD Runge-Kutta method [20] is used for thetime integration in all the experiments.In Fig. 1, the accuracy of the high-order schemes with the di ff erent reconstruction variablesis examined by performing the one-dimensional circularly polarized Alfv´en wave test [22]. Thetime step size is fixed small enough so as not to be a ff ected by errors of the time integral. TheCFL number at the highest resolution is on the order of 10 − . The results indicate that W - reconstruction in the MUSCL scheme with the minmod limiter is more accurate than U -, V -,and W (cid:48) - reconstructions because the limiter is not activated in the distribution of the exact char-acteristic variables. On the other hand, the higher-order reconstructions, including the MUSCLscheme with the Koren limiter [11], do not depend on the reconstruction variables and achievethe expected order of accuracy. Note that although the Koren limiter can reconstruct a smoothdistribution with third-order accuracy, it is limited to second-order in this problem.We present a typical one-dimensional shock tube problem [6], which contains two fast shocks,two slow shocks, two rotational discontinuities, and one contact discontinuity. The left and rightstates are initialized as V L = (1 . , . , . , . , . / √ π, / √ π, . T and V R = (1 , , , , / √ π, / √ π, T B x = / √ π , respectively. The CFL number is 0 . reconstructions produce the similar solutions without nu-merical oscillations. The results obtained by the WCNS scheme, however, reveal that consid-erable overshoots are observed in U - and V - reconstructions as shown in Fig. 2. By contrast,the WCNS scheme with W (cid:48) - reconstruction is as robust as that with W - reconstruction thoughsmall bumps are observed. Another typical shock tube problem [5] is also demonstrated. Thisproblem contains a number of di ff erent waves: two fast rarefaction waves, one slow shock,one slow compound wave, and a contact discontinuity. The left and right states are initiallygiven by V L = (1 , , , , , , T and V R = (0 . , , , , − , , . T with B x = .
75. TheMUSCL-minmod scheme with each reconstruction gives similar reasonable results. On theother hand, the WCNS scheme with each reconstruction yields slightly di ff erent solutions asshown in Fig. 3. Large numerical oscillations are observed in U - and V - reconstructions . In W (cid:48) - reconstruction , however, numerical oscillations are reduced though those are still larger than in W - reconstruction . Note that a numerical oscillation with a relatively long wavelength attachedto the right-moving fast rarefaction wave is observed even in W - reconstruction . Furthermore,we performed several shock tube problems consisting of pure waves such as isolated shocks andrarefaction waves [8], and confirmed that the results for W (cid:48) - reconstruction are very similar tothose for W - reconstruction (not shown).Finally, let us demonstrate the Orszag-Tang vortex problem [17], which is a standard test formultidimensional MHD. Considering the trade-o ff between computational cost and accuracy, thehigh-order schemes are necessary, especially for practical multidimensional problems. To main-tain the solenoidal condition of the magnetic field, the hyperbolic divergence cleaning method[7] is used. The number of the grid is 400 for x , y ∈ [0 , π ], and the CFL number is set to 0 . W - and W (cid:48) - reconstructions yields a similar result to the refer-ence solution. In contrast to that, visible numerical oscillations attached to shocks are observedin high-compression regions when U - (not shown) and V - reconstructions are adopted.
5. Conclusions
We have proposed a set of quick and easy reconstruction variables for high-order shock cap-turing schemes for MHD. The numerical experiments suggest that W (cid:48) - reconstruction is capableof the higher-order reconstructions almost without numerical oscillations even when B x (cid:44) dv x / dx (cid:44)
0. In the present WCNS code, the computational time of W - reconstruction is more than3 . . U - and V - reconstructions ,while the time of W (cid:48) - reconstruction is about 1 . . W (cid:48) - reconstruction are more robust than those with U -and V - reconstructions , while their computational e ffi ciencies are comparable. Thus we concludethat W (cid:48) - reconstruction can be a practical choice for the higher-order reconstructions in shockcapturing schemes for MHD.In terms of robustness or accuracy, W (cid:48) - reconstruction can be switched to W - reconstruction orothers adaptively in a very problematic situation. Indeed, another set of reconstruction variablesthat improves the accuracy in low Mach number MHD flows has been proposed [15]. We expectthat the present approach of introducing reasonably approximated characteristic fields can beextended to more complex systems. The system of relativistic MHD may be an example sincethe corresponding eigensystem is too complicated to be used for the higher-order reconstructions.6 cknowledgement This work was supported by MEXT / JSPS KAKENHI Grant Numbers JP20K11851, JP20H00156,JP19H01928 (T. Miyoshi).
References [1] D. Balsara, Total variation diminishing scheme for adiabatic and isothermal magnetohydrodynamics, Astrophys. J.Suppl. 116 (1998) 133-153.[2] D. Balsara, Second-order-accurate schemes for magnetohydrodynamics with divergence-free reconstruction, As-trophys. J. Suppl. 151 (2004) 149-184.[3] D. Balsara, C. W. Shu, Monotonicity preserving weighted essentially non-oscillatory schemes with increasinglyhigh order of accuracy, J. Comput. Phys. 160 (2000) 405-452.[4] D. Balsara, T. Rumpf, M. Dumbser, C. D. Munz, E ffi cient, high accuracy ADER-WENO schemes for hydrody-namics and divergence-free magnetohydrodynamics, J. Comput. Phys. 228 (2009) 2480-2516.[5] M. Brio, C. C. Wu, An upwind di ff erencing scheme for the equations of ideal magnetohydrodynamics, J. Comput.Phys. 75 (1988) 400-422.[6] W. Dai, P. R. Woodward, An approximate Riemann solver for ideal magnetohydrodynamics, J. Comput. Phys. 111(1994) 354-372.[7] A. Dedner, F. Kemm, D. Kr¨oner, C. D. Munz, T. Schnitzer, M. Wesenberg, Hyperbolic divergence cleaning for theMHD equations, J. Comput. Phys. 175 (2002) 645-673.[8] S. A. E. G. Falle, S. S. Komissarov, P. Joarder, A multidimensional upwind scheme for magnetohydrodynamics,Mon. Not. R. Astron. Soc. 297 (1998) 265-277.[9] A. Harten, B. Engquist, S. Osher, S. R. Chakravarthy, Uniformly high order accurate essentially non-oscillatoryschemes, III, J. Comput. Phys. 71 (1987) 231-303.[10] G. S. Jiang, C. C. Wu, A high-order WENO finite di ff erence scheme for the equations of ideal magnetohydrody-namics, J. Comput. Phys. 150 (1999) 561-594.[11] B. Koren, A robust upwind discretization method for advection, di ff usion and source terms, in: C. B. Vreugdenhil,B. Koren (Eds.), Numerical Methods for AdvectionDi ff usion Problems, Vieweg, Braunschweig, Germany, 1993,pp. 117-138.[12] X. D. Liu, S. Osher, T. Chan, Weighted essentially non-oscillatory schemes, J. Comput. Phys. 115 (1994) 200-212.[13] Y. Matsumoto, Y. Asahina, Y. Kudoh, T. Kawashima, J. Matsumoto, H. R. Takahashi, T. Minoshima, S. Zenitani,T. Miyoshi, R. Matsumoto, Magnetohydrodynamic simulation code CANS + : Assessments and applications, Publ.Astron. Soc. Japan 71 (2019) 83(1-26).[14] T. Minoshima, T. Miyoshi, Y. Matsumoto, A high-order weighted finite di ff erence scheme with a multistate approx-imate Riemann solver for divergence-free magnetohydrodynamic simulations, Astrophys. J. 242 (2019) 14(1-29).[15] T. Minoshima, K. Kitamura, T. Miyoshi, A multistate low-dissipation advection upstream splitting method for idealmagnetohydrodynamics, Astrophys. J. 248 (2020) 12(1-21).[16] T. Miyoshi, K. Kusano, A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics, J. Com-put. Phys. 208 (2005) 315-344.[17] A. Orszag, C. M. Tang, Small-scale structure of two-dimensional magnetohydrodynamic turbulence, J. Fluid.Mech. 90 (1979) 129-143.[18] K. G. Powell, P. L. Roe, T. J. Linde, T. I. Gombosi, D. L. De Zeeuw, A solution-adaptive upwind scheme for idealmagnetohydrodynamics, J. Comput. Phys. 154 (1999) 284-309.[19] P. L. Roe, D. S. Balsara, Notes on the eigensystem of magnetohydrodynamics, SIAM J. Appl. Math. 56 (1996)57-67.[20] C. W. Shu, S. Osher, E ffi cient implementation of essentially non-oscillatory shock-capturing schemes, J. Comput.Phys. 77 (1988) 439-471.[21] A. Suresha, H. T. Huynh, Accurate monotonicity-preserving schemes with RungeKutta time stepping, J. Comput.Phys. 136 (1997) 83-99.[22] G. T´oth, The ∇ · B constraint in shock-capturing magnetohydrodynamic codes, J. Comput. Phys. 161 (2000) 605-652.[23] B. van Leer, Towards the ultimate conservative di ff erence scheme, V. A second order sequel to Godunov’s method,J. Comput. Phys. 32 (1979) 101-136. -8 -7 -6 -5 -4 -3 -2 -1 Δ xL UVWW' -8 -7 -6 -5 -4 -3 -2 -1 Δ xL UVWW'
Figure 1: The L -norm errors of the magnetic field strength for di ff erent grid resolutions. Each color corresponds toa di ff erent reconstruction . (Left) The base scheme (solid line), MUSCL scheme with the minmod limiter (dash-dottedlines), and WCNS scheme (dashed lines). Orange open circles show the second-order linear reconstruction. (Right) TheMUSCL scheme with the Koren limiter (dash-dotted lines). The base and WCNS schemes are also included because ofvisibility. Orange open circles show the third-order linear reconstruction. .01.5-1.0 -0.5 0.0 0.5 1.0 B y UVWW' B y UVWW'
Figure 2: Dai-Woodward shock tube problem computed by the WCNS scheme. B y is plotted with the color correspond-ing to each set of reconstruction variables. The right is an enlarged view. .00.5-1.0 -0.5 0.0 0.5 1.0 v x UVWW' v x UVWW'
Figure 3: Brio-Wu shock tube problem computed by the WCNS scheme. V x is plotted with the color corresponding toeach set of reconstruction variables. The right is an enlarged view. a) (b) (c) (d) 0.00.51.0 Figure 4: Schlieren-like images calculated as exp ( − |∇ ρ | / |∇ ρ | max ) for the Orszag-Tang vortex obtained by the WCNSscheme with (a) V - reconstruction , (b) W - reconstruction , and (c) W (cid:48) - reconstruction . (d) The reference solution obtainedby the MUSCL-minmod scheme with W - reconstruction on 4000 grid. Since the resolution of the grid is 10 times higher,exp ( − |∇ ρ | / |∇ ρ | max ) is shown for comparison.) is shown for comparison.