A Spatial-Temporal asymptotic preserving scheme for radiation magnetohydrodynamics
AA SPATIAL-TEMPORAL ASYMPTOTIC PRESERVING SCHEME FOR RADIATIONMAGNETOHYDRODYNAMICS
SHI JIN, MIN TANG, AND XIAOJIANG ZHANG
Abstract.
The radiation magnetohydrodynamics (RMHD) system couples the ideal magnetohydrodynam-ics equations with a gray radiation transfer equation. The main challenge is that the radiation travels at thespeed of light while the magnetohydrodynamics changes with fluid. The time scales of these two processescan vary dramatically. In order to use mesh sizes and time steps that are independent of the speed of light,asymptotic preserving (AP) schemes in both space and time are desired. In this paper, we develop an APscheme in both space and time for the RMHD system. Two different scalings are considered, one results inan equilibrium diffusion limit system, while the other results in a non-equilibrium system. The main idea isto decompose the radiative intensity into three parts, each part is treated differently. The performances ofthe semi-implicit method are presented, for both optically thin and thick regions, as well as for the radiativeshock problem. Comparisons with the semi-analytic solution are given to verify the accuracy and asymptoticproperties of the method. Introduction
Radiation magnetohydrodynamics (RMHD) is concerned about the dynamical behaviors of magnetizedfluids that have nonnegligible exchange of energy and momentum with radiation, which is important in hightemperature flow systems, solar and space physics and astrophysics. The radiation transfer equation (RTE)for the radiative intensity I and the fluid temperature T in the mixed frame is governed by ∂I∂t + C n · ∇ I = Cσ a (cid:18) a r T π − I (cid:19) + Cσ s ( J − I ) + 3 n · v σ a (cid:18) a r T π − J (cid:19) + n · v ( σ a + σ s )( I + 3 J ) − σ s v · H − ( σ a − σ s ) v · v C J − ( σ a − σ s ) v · ( v · K ) C , (1.1)where σ a and σ s represent absorption and scattering opacities respectively, n ∈ V is the angular variable, a r is the radiation constant and C is the speed of light [7]. The ideal MHD equations with radiation energyand momentum source terms are ∂ρ∂t + ∇ · ( ρ v ) = 0 ,∂ ( ρ v ) ∂t + ∇ · ( ρ vv − BB + P ∗ ) = − S rp ,∂E∂t + ∇ · [( E + P ∗ ) v − B ( B · v )] = − CS re ,∂ B ∂t + ∇ × ( v × B ) = 0 , (1.2)where S re = σ a ( a r T − E r ) + ( σ a − σ s ) v C · [ F r − ( v E r + v · P r )] , S rp = − σ s + σ a C [ F r − ( v E r + v · P r )] + v C σ a ( a r T − E r ) , Date : January 12, 2021. a r X i v : . [ m a t h . NA ] J a n SHI JIN, MIN TANG, AND XIAOJIANG ZHANG with the energy density E r = 4 πJ = 4 π | V | (cid:90) V I d n = 4 π (cid:104) I (cid:105) , the radiation flux F r = 4 πC H = 4 πC | V | (cid:90) V n I d n = 4 πC (cid:104) n I (cid:105) , and the radiation pressure P r = 4 πK = 4 π | V | (cid:90) V n I d n = 4 π (cid:10) n I (cid:11) . In particular, the operators (cid:104)·(cid:105) , (cid:104) n ·(cid:105) and (cid:10) n · (cid:11) are respectively | V | (cid:82) V · d n , | V | (cid:82) V n · d n and | V | (cid:82) V n · d n in context. Moreover, ρ is fluid density, v is fluid velocity, B is magnetic flux density, p is static pressure, P ∗ is full pressure of fluid, P ∗ ≡ ( p + B / E = E g + 12 ρv + B , where E g is the internal energy density, v = v · v and B = B · B . The radiation MHD system is closed bythe perfect gas equation of state: E g = p/ ( γ − , T = p/ ( R ideal ρ ) , where γ is adiabatic index for an ideal gas, and R ideal is the ideal gas constant.Numerical approximations to RTE have been extensively studied in [12,14,15,21,24,30] and for ideal MHDsee for example [2,18,29]. The RMHD coupled system has been investigated lately [1,6,7,13,27,28,31]. Sincethe transport term in the RTE and the advection term in the MHD equations are at different time scales,the RMHD system can be stiff, as can been seen more clearly after nondimensionalization [19]. To solve theRMHD system, classical numerical discretizations require the space and time steps resolve the speed of light,which leads to very expensive computational cost. Asymptotic-Preserving (AP) schemes provide a genericframework for such multiscale problems [8, 9]. AP schemes were first studied for steady neutron transportproblems in diffusive regimes [16, 17] and then [5, 10] for boundary value problems, and unsteady transportproblems [11,14]. In the discrete setting, from the microscopic scale to the macroscopic, AP schemes preservethe asymptotic limit. When the scaling parameter ε in the multiscale system can not be resolved numerically,AP schemes should automatically become a good solver for the macroscopic models.AP schemes for RMHD have been proposed in the literature but since most of them are based on operatorsplitting, they are only AP in space. In [1], Simon et al. have developed a second-order scheme in both spaceand time for the Euler equations coupled with a gray radiation S model. The scheme in [1] uses MUSCL-Hancock method to solve the Euler equations and the lumped linear-discontinuous Galerkin method for theradiation S model. For the time discretization, the TR/BDF2 time integration method is employed. Jianget al. devised an Implicit/Explicit (IMEX) scheme for the splitting system that treats the transport term inthe RTE and convective term in the fluid equations explicitly, and the source term implicitly [7]. Recently,another method by Sun etc. [31] uses the gas-kinetic scheme (GKS) to deal with Euler equations and theunified gas-kinetic scheme (UGKS) for the RTE. However, all the aforementioned methods either require atime step that satisfies a hyperbolic constraint with CFL number proportional to the light speed [1, 7], oruse nonlinear iterations involving (1.1) in order to get rid of the time step constraint [31]. Both approachesare expensive. It is desired to design a scheme that can use a time step that is independent of the lightspeed and employs nonlinear iterations that involve only macroscopic quantities, i.e. the ρ , v , E , B or themoments of radiative intensity.In this paper, we aim at developing a scheme for RMHD that is AP in both space and time. Twodifferent parameter regimes are considered. One is for σ a small and σ s large, which results in the non-equilibrium diffusion limit system. The other is for σ a large, which gives the equilibrium diffusion limit SPATIAL-TEMPORAL ASYMPTOTIC PRESERVING SCHEME FOR RADIATION MAGNETOHYDRODYNAMICS 3 system. The proposed scheme preserves both limits. The main idea is to decompose the intensity into threeparts, two parts correspond to the zeroth and first order moments, while the third part is the residual. Thetwo macroscopic moments are treated implicitly while the residual is explicit. Thanks to the properties ofresidual term, one can update the macroscopic quantities first. For the space discretization, we use Roe’smethod in the Athena code to solve the convective part in the ideal MHD equations and the UGKS for theRTE. AP property of both semi-discretized and fully discretizad are proved.This paper is organized as follows. Section 2 gives the asymptotic limit of the coupled system, and theequilibrium and non-equilibrium diffusion limit systems under different scalings are obtained. The semi-discretization and one dimensional fully discretized scheme for RMHD are illustrated in Sections 3 and 4respectively. Their capability of capturing both diffusion limits are proved. In Section 5, performances ofour AP method are presented in both optically thin and thick regions. The radiative shock problems aretested. By comparing with the semi-analytic solutions in [20], we observe that the scheme is accurate andstable using large time step and meshes that are independent of the light speed. Finally, we conclude inSection 6. 2.
The diffusion limit of RMHD
As presented in [19], the dimensionless form for the coupled system can give a better insight on the relativeimportance of different terms. Consider the following nondimensionalization: x = (cid:98) x(cid:96) ∞ , t = (cid:98) t(cid:96) ∞ /a ∞ , ρ = (cid:98) ρρ ∞ , v = (cid:98) v a ∞ , p = (cid:98) pρ ∞ a ∞ , T = (cid:98) T T ∞ , I = a r T ∞ (cid:98) I,E r = a r T ∞ (cid:99) E r , F r = Ca r T ∞ (cid:99) F r , P r = a r T ∞ (cid:99) P r , σ a = λ a (cid:99) σ a , σ s = λ s (cid:98) σ s , where variables with a hat denote nondimensional quantity and variables with ∞− subscript are the char-acteristic value with unit. More precisely, (cid:96) ∞ , a ∞ , ρ ∞ and T ∞ are respectively the reference length, soundspeed, density and temperature. Moreover, λ a and λ s are respectively the characteristic value of absorbingand scattering coefficients. Then the full dimensionless radiation MHD system becomes ∂I∂t + C n · ∇ I = L a C σ a (cid:18) T π − I (cid:19) + L s C ( J − I ) + 3 L a n · v σ a (cid:18) T π − J (cid:19) + n · v ( L a σ a + L s σ s ) ( I + 3 J ) − L s σ s v · H − ( L a σ a − L s σ s ) v · v C J − ( L a σ a − L s σ s ) v · ( v · K ) C (cid:44) C S,∂ρ∂t + ∇ · ( ρ v ) = 0 ,∂ ( ρ v ) ∂t + ∇ · ( ρ vv − BB + P ∗ ) = −P S rp ,∂E∂t + ∇ · [( E + P ∗ ) v − B ( B · v )] = −CP S re ,∂ B ∂t + ∇ × ( v × B ) = 0 , (2.1)where C = Ca ∞ , P = a r T ∞ ρ ∞ a ∞ , L a = (cid:96) ∞ λ a , L s = (cid:96) ∞ λ s , S re = 4 π (cid:104) S (cid:105) and S rp = 4 π (cid:104) n S (cid:105) . In (2.1), P is a nondimensional constant, which measures the influence of radiation on the flow dynamics. In [4], theauthors considered the equilibrium diffusion regime and non-equilibrium diffusion regime for a simplifiedcoupled system. In the subsequent part, we consider the equilibrium diffusion limit and non-equilibriumdiffusion limit of the RMHD system by using the same scaling assumptions as in [4].2.1. The non-equilibrium diffusion limit.
In the simulation of radiation-hydrodynamics in [1], the speedof light C = 29 . cm ns − , in nonrelativistic flows the sound speed is 0.001386 cm ns − , the radiationconstent a r = 0 . Jk cm − KeV − , the typical computational domain is 0.04 cm , and the units of the SHI JIN, MIN TANG, AND XIAOJIANG ZHANG length, time, temperature and energy are respectively cm, ns, KeV and Jk . Therefore 1 / C = O ( ε ) , P = O (1) . In non-equilibrium regime, we assume L a = ε, L s = 1 /ε , and the radiation intensity I tends to astate which can be characterized by a temperature which is different from the material temperature T . Insummary, we consider the following scaling L a = ε, L s = 1 /ε, P = O (1) , C = c/ε. The RMHD system (2.1) becomes (after dropping the hat): ∂I∂t + cε n · ∇ I = cσ a (cid:18) T π − I (cid:19) + cσ s ε ( J − I ) + 3 ε n · v σ a (cid:18) T π − J (cid:19) + n · v (cid:16) εσ a + σ s ε (cid:17) ( I + 3 J ) − σ s ε v · H − ( ε σ a − σ s ) v · v c J − ( ε σ a − σ s ) v · ( v · K ) c , (2.2a) ∂ρ∂t + ∇ · ( ρ v ) = 0 , (2.2b) ∂ ( ρ v ) ∂t + ∇ · ( ρ vv − BB + P ∗ ) = −P S rp , (2.2c) ∂E∂t + ∇ · [( E + P ∗ ) v − B ( B · v )] = − cε P S re , (2.2d) ∂ B ∂t + ∇ × ( v × B ) = 0 , (2.2e)where S re = εσ a (cid:18) T − π | V | (cid:90) V I d n (cid:19) + ( ε σ a − εσ s ) 4 π v c | V | · (cid:20) cε (cid:90) V n I d n − (cid:18) v (cid:90) V I d n + v (cid:90) V n I d n (cid:19)(cid:21) , S rp = − π ( σ s + ε σ a ) c | V | (cid:20) cε (cid:90) V n I d n − (cid:18) v (cid:90) V I d n + v (cid:90) V n I d n (cid:19)(cid:21) + v c ε σ a (cid:18) T − π | V | (cid:90) V I d n (cid:19) . When ε →
0, we show that the solution to (2.2) can be approximated by the solution to a nonlinear diffusionequation coupled with a MHD system. We assume that the radiation intensity I and temperature T havethe following Chapman-Enskog expansion such that I = I (0) + εI (1) + ε I (2) + · · · ,T = T (0) + εT (1) + ε T (2) + · · · . (2.3)By substituting ansatz (2.3) into equation (2.2a) and collecting the terms of the same order in ε , we have O (cid:18) ε (cid:19) : I (0) = J (0) , (2.4a) O (cid:18) ε (cid:19) : c n · ∇ I (0) = cσ s (cid:16) J (1) − I (1) (cid:17) + n · v σ s (cid:16) I (0) + 3 J (0) (cid:17) . (2.4b)Multiplying both sides of (2.2a) by n , taking its integral with respect to n and combining the obtainedequation with (2.2c), one can get the following momentum conservation equation: ∂ t (cid:18) ρ v P + 4 πεc (cid:104) n I (cid:105) (cid:19) + ∇ · (cid:18) ρ vv − BB + P ∗ P + 4 π (cid:10) n I (cid:11)(cid:19) = 0 . (2.5)By using (2.4a) and (cid:10) n (cid:11) = , when ε →
0, (2.5) gives ∂ t ( ρ (0) v (0) ) + ∇ · (cid:16) ρ (0) v (0) v (0) − B (0) B (0) + ( P ∗ ) (0) (cid:17) = − π P ∇ J (0) . (2.6)By taking the integral with respect to n on both sides of (2.2a), and combining it with (2.2d), one can obtainthe energy conservation equation ∂ t (cid:18) E P + 4 π (cid:104) I (cid:105) (cid:19) + ∇ · (cid:20) ( E + P ∗ ) v − B ( B · v ) P + 4 πc (cid:104) n I (cid:105) ε (cid:21) = 0 . (2.7) SPATIAL-TEMPORAL ASYMPTOTIC PRESERVING SCHEME FOR RADIATION MAGNETOHYDRODYNAMICS 5
By using (2.4a), when ε →
0, (2.7) gives ∂ t (cid:18) E (0) P + 4 π (cid:68) I (0) (cid:69)(cid:19) + ∇ · (cid:34) ( E (0) + ( P ∗ ) (0) ) v (0) − B (0) ( B (0) · v (0) ) P + 4 πc (cid:68) n I (1) (cid:69)(cid:35) = 0 . Noting (2.4b), we find ∂ t (cid:16) E (0) + 4 π P J (0) (cid:17) + ∇ · (cid:20) ( E (0) + ( P ∗ ) (0) ) v (0) − B (0) ( B (0) · v (0) ) + 16 π P v (0) J (0) − πc P σ s ∇ J (0) (cid:21) = 0 . (2.8)(2.2b), (2.6), (2.8) and (2.2e) give four equations for ρ, ρ v , E, J, B , to obtain a closed system, we needone more equation for J . Taking (cid:104)·(cid:105) on both sides of equation (2.2a), yields: ∂ t (cid:104) I (cid:105) + cε ∇ · (cid:104) n I (cid:105) = c πε S re . By using the expansion in (2.3) and (2.4a), the leading order terms are4 π∂ t (cid:68) I (0) (cid:69) +4 πc ∇· (cid:68) n I (1) (cid:69) = cσ a (cid:16) ( T (0) ) − π (cid:68) I (0) (cid:69)(cid:17) − πσ s v (0) c (cid:104) c (cid:68) n I (1) (cid:69) − v (0) (cid:16)(cid:68) I (0) (cid:69) + (cid:68) n I (0) (cid:69)(cid:17)(cid:105) , then from (2.4b), one gets4 π∂ t J (0) + ∇ · (cid:18) π v (0) J (0) − πc σ s ∇ J (0) (cid:19) = cσ a (cid:16) ( T (0) ) − πJ (0) (cid:17) + 4 π v (0) ∇ J (0) . (2.9)Therefore, when ε → ∂ t ρ + ∇ · ( ρ v ) = 0 , (2.10a) ∂ t ( ρ v ) + ∇ · ( ρ vv − BB + P ∗ ) = − π P ∇ J, (2.10b) ∂ t ( E + 4 π P J ) + ∇ · (cid:20) ( E + P ∗ ) v − B ( B · v ) + 16 π P v J (cid:21) = ∇ · (cid:18) πc P σ s ∇ J (cid:19) , (2.10c)4 π∂ t J + ∇ · (cid:18) π v J − πc σ s ∇ J (cid:19) = cσ a (cid:0) T − πJ (cid:1) + 4 π v ∇ J, (2.10d) ∂ t B + ∇ × ( v × B ) = 0 . (2.10e)Here, the non-equilibrium system indicates that J is away from T π , while we will see in section 2.2 that theequilibrium diffusion limit indicates that J ≈ T π .2.2. The equilibrium diffusion limit.
The setting of P and C in this regime are the same with thenon-equilibrium regime. Moreover, we assume L a = 1 /ε, L s = ε , and the radiation intensity I adapt to thematerial temperature. In summary, we consider the following scaling L a = 1 /ε, L s = ε, P = O (1) , C = c/ε. SHI JIN, MIN TANG, AND XIAOJIANG ZHANG
This is called the equilibrium diffusion regime, which will give the classical equilibrium diffusion limit [19,23, 25]. The RMHD system (2.1) becomes (after dropping the hat): ∂I∂t + cε n · ∇ I = cσ a ε (cid:18) T π − I (cid:19) + cσ s ( J − I ) + 3 n · v σ a ε (cid:18) T π − J (cid:19) + n · v (cid:16) σ a ε + εσ s (cid:17) ( I + 3 J ) − εσ s v · H − ( σ a − ε σ s ) v · v c J − ( σ a − ε σ s ) v · ( v · K ) c , (2.11a) ∂ρ∂t + ∇ · ( ρ v ) = 0 , (2.11b) ∂ ( ρ v ) ∂t + ∇ · ( ρ vv − BB + P ∗ ) = −P S rp , (2.11c) ∂E∂t + ∇ · [( E + P ∗ ) v − B ( B · v )] = − cε P S re , (2.11d) ∂ B ∂t + ∇ × ( v × B ) = 0 , (2.11e)where S re = σ a ε (cid:18) T − π | V | (cid:90) V I d n (cid:19) + ( εσ a − ε σ s ) 4 π v c | V | · (cid:20) cε (cid:90) V n I d n − (cid:18) v (cid:90) V I d n + v (cid:90) V n I d n (cid:19)(cid:21) , S rp = − π ( ε σ s + σ a ) c | V | (cid:20) cε (cid:90) V n I d n − (cid:18) v (cid:90) V I d n + v (cid:90) V n I d n (cid:19)(cid:21) + v c σ a (cid:18) T − π | V | (cid:90) V I d n (cid:19) . Using the same expansion as in (2.3), one has T = ( T (0) + εT (1) + ε T (2) + · · · ) = ( T (0) ) + 4 ε ( T (0) ) T (1) + · · · . (2.12)Then substituting (2.12) into (2.11a) and collecting the terms of the same order of ε , one gets O (cid:18) ε (cid:19) : I (0) = ( T (0) ) π , (2.13a) O (cid:18) ε (cid:19) : c n · ∇ I (0) = cσ a (cid:18) T (0) ) T (1) π − I (1) (cid:19) + n · v σ a (cid:18) I (0) + 3( T (0) ) π (cid:19) , (2.13b)by similar calculations in section 2.1, one can get the following equilibrium diffusion limit: ∂ t ρ + ∇ · ( ρ v ) = 0 , (2.14a) ∂ t ( ρ v ) + ∇ · ( ρ vv − BB + P ∗ ) = − P ∇ T , (2.14b) ∂ t (cid:0) E + P T (cid:1) + ∇ · (cid:20) ( E + P ∗ ) v − B ( B · v ) + 4 P v T (cid:21) = ∇ · (cid:18) c P σ a ∇ T (cid:19) , (2.14c) ∂ t B + ∇ × ( v × B ) = 0 . (2.14d)As we can see, when J = T π in (2.10), the non-equilibrium diffusion system (2.10) is the same as theequilibrium diffusion system (2.14), which indicates J is close to T π in the equilibrium system (2.11).3. Time discretization for the RMHD
In this section, based on a decomposition for radiation intensity I , we will present a semi-discretizationfor the RMHD and show its AP property for any dimensions in space and angular variables. SPATIAL-TEMPORAL ASYMPTOTIC PRESERVING SCHEME FOR RADIATION MAGNETOHYDRODYNAMICS 7
A decomposition.
To preserve the asymptotic limit of (2.1) while keeping the computational com-plexity under control, the intensity I will be decomposed into three parts, and each of them is treated in adifferent way. More precisely, we decompose I ( t, x, n ) as follows: I ( t, x, n ) = (cid:104) I (cid:105) + 3 n (cid:104) n I (cid:105) + εQ ( t, x, n ) , (3.1a):= J ( t, x ) + ε n R ( t, x ) + εQ ( t, x, n ) , (3.1b)with 3 (cid:104) n I (cid:105) = εR . Then taking (cid:104) n ·(cid:105) on both sides of (3.1a) and using the condition (cid:10) n (cid:11) = yields (cid:104) n I (cid:105) = (cid:104) n I (cid:105) + ε (cid:104) n Q (cid:105) , which gives (cid:104) n Q (cid:105) = 0 . (3.2)Moreover, taking (cid:104)·(cid:105) on (3.1b), yields (cid:104) I (cid:105) = (cid:104) J + ε n R + εQ (cid:105) = (cid:104) J + εQ (cid:105) = (cid:104) I (cid:105) + ε (cid:104) Q (cid:105) , which implies (cid:104) Q (cid:105) = 0 . (3.3)Substituting the decomposition (3.1) into system (2.2) leads to ∂ t ( J + ε ( n R + Q )) + cε n ∇ · ( J + ε ( n R + Q )) = cσ a (cid:18) T π − J (cid:19) − (cid:16) cεσ a + cσ s ε (cid:17) ( n R + Q ) − σ s v R + nv (cid:16) εσ a + σ s ε (cid:17) (4 J + ε ( n R + Q )) + 3 εσ a nv (cid:18) T π − J (cid:19) − ( ε σ a − σ s ) v · v c (cid:18) J + εK Q (cid:19) , (3.4a) ∂ρ∂t + ∇ · ( ρ v ) = 0 , (3.4b) ∂ ( ρ v ) ∂t + ∇ · ( ρ vv − BB + P ∗ ) = −P (cid:100) S rp , (3.4c) ∂E∂t + ∇ · [( E + P ∗ ) v − B ( B · v )] = − cε P (cid:99) S re , (3.4d) ∂ B ∂t + ∇ × ( v × B ) = 0 , (3.4e)where K Q = (cid:10) n Q (cid:11) and (cid:99) S re = εσ a ( T − πJ ) + 4 π ( ε σ a − εσ s ) v c · (cid:20) c R − v (cid:18) J + εK Q (cid:19)(cid:21) , (cid:100) S rp = − π ( σ s + ε σ a ) c (cid:20) c R − v (cid:18) J + εK Q (cid:19)(cid:21) − ε v c σ a ( T − πJ ) . (3.4a) has three unknown functions, J, R and Q . Thanks to the properties of Q in (3.2),(3.3), the zeroth andfirst moment equations of (3.4a) have no time derivatives with respect to Q . Therefore, instead of solvingthe RMHD system (3.4) directly, we take the zeroth and first moments of (3.4a), and then couple themtogether with (3.4b)–(3.4e) to solve variables ρ, J, R, v , T, B . More precisely, we have ∂ t J + c ∇ · R = c πε (cid:99) S re ,ε ∂ t R + c ε ∇ · J + c ∇ · K Q = c πε (cid:100) S rp , (3.5)coupled with (3.4b)–(3.4e) and by treating K Q explicitly, one can update the macroscopic quantities first.Then Q can be updated implicitly by the first equation in system (3.4a), and the new I can be given by(3.1). The details are illustrated in the next two subsections. SHI JIN, MIN TANG, AND XIAOJIANG ZHANG
The time discretization for the non-equilibrium regime.
Let t s = s ∆ t, s = 0 , , , · · · , and I s ≈ I ( x, t s , n ) , ρ s ≈ ρ ( x, t s , n ) , v s ≈ v ( x, t s , n ) ,T s ≈ T ( x, t s , n ) , B s ≈ B ( x, t s , n ) . We use the following semi-discrete scheme in time for (3.4), which reads: J s +1 − J s ∆ t + ε n R s +1 − R s ∆ t + ε Q s +1 − Q s ∆ t + cε n ∇ · J s +1 + c n ∇ · ( n R s +1 + Q s )= cσ a (cid:18) ( T s +1 ) π − J s +1 (cid:19) − (cid:16) cεσ a + cσ s ε (cid:17) ( n R s +1 + Q s +1 ) − σ s v s +1 R s +1 +3 εσ a nv s +1 (cid:18) ( T s +1 ) π − J s +1 (cid:19) + nv s +1 (cid:16) εσ a + σ s ε (cid:17) (4 J s +1 + ε ( n R s +1 + Q s )) − ( ε σ a − σ s ) v s +1 · v s +1 c (cid:18) J s +1 + εK sQ (cid:19) , (3.6a) ρ s +1 − ρ s ∆ t + ∇ · ( ρ s v s ) = 0 , (3.6b)( ρ v ) s +1 − ( ρ v ) s ∆ t + ∇ · ( ρ s v s v s − B s B s + ( P s ) ∗ ) = −P ( (cid:100) S rp ) s +1 , (3.6c) a ( ρT ) s +1 − ( ρT ) s ∆ t + ρ s +1 ( v s +1 ) − ρ s ( v s ) t + ( B s +1 ) − ( B s ) t + ∇ · [( E s + ( P ∗ ) s ) v s − B s ( B s · v s )] = − cε P ( (cid:99) S re ) s +1 , (3.6d) B s +1 − B s ∆ t + ∇ × ( v s × B s ) = 0 . (3.6e)where K sQ = (cid:10) n Q s (cid:11) and( (cid:99) S re ) s +1 = εσ a (( T s +1 ) − πJ s +1 ) + 4 π ( ε σ a − εσ s ) v s +1 c (cid:20) c R s +1 − v s +1 (cid:18) J s +1 + εK sQ (cid:19)(cid:21) , ( (cid:100) S rp ) s +1 = − π ( σ s + ε σ a ) c (cid:20) c R s +1 − v s +1 (cid:18) J s +1 + εK sQ (cid:19)(cid:21) − ε v s +1 c σ a (( T s +1 ) − πJ s +1 ) . The zeroth and first moment equations of (3.6a) are J s +1 − J s ∆ t + c ∇ · R s +1 = c πε ( (cid:99) S re ) s +1 , (3.7a) ε c R s +1 − R s ∆ t + 13 ∇ · J s +1 + ε ∇ · K sQ = ( (cid:100) S rp ) s +1 π . (3.7b)We first update ρ, B by (3.6b), (3.6e) and then solve (3.6c), (3.6d) and (3.7) to get J, R, v , T . Q s +1 can thenbe updated implicitly by equation (3.6a), and we get the I s +1 by using (3.1). The diffusion limit of (3.7):
Then we will show the AP property of the semi-discretization (3.7). From(3.7b), R s +1 = 4 v s +1 J s +1 c − ∇ · J s +1 σ s + O ( ε ) . (3.8) SPATIAL-TEMPORAL ASYMPTOTIC PRESERVING SCHEME FOR RADIATION MAGNETOHYDRODYNAMICS 9
Multiplying equation (3.7a) by 4 π P , and adding it up with equation (3.6d), the terms on the right handsides cancel. By using (3.8), and sending ε → π P × (3.7a) + (3.6d) yields a ( ρT ) s +1 − ( ρT ) s ∆ t + ρ s +1 ( v s +1 ) − ρ s ( v s ) t + ( B s +1 ) − ( B s ) t + 4 π P J s +1 − J s ∆ t + ∇ · (cid:20) ( E s + ( P s ) ∗ ) v s − B s ( B s · v s ) + 16 π P v s +1 J s +1 (cid:21) = ∇ · (cid:18) πc P σ s ∇ J s +1 (cid:19) . (3.9)This is a semi-discretization for (2.10c). Then multiplying equation (3.7b) by 4 π P , and adding it to equation(3.6c), sending ε →
0, formally one gets( ρ v ) s +1 − ( ρ v ) s ∆ t + ∇ · ( ρ s v s v s − B s B s + ( P ∗ ) s ) = − π P ∇ · J s +1 . (3.10)This is a semi-discretization for (2.10b). Sending ε → π J s +1 − J s ∆ t + 4 πc ∇ · R s +1 = cσ a (( T s +1 ) − πJ s +1 ) − πσ s v s +1 c (cid:18) c R s +1 − v s +1 J s +1 (cid:19) , (3.11)using (3.8) in (3.11), we get4 π J s +1 − J s ∆ t + 16 π ∇ · ( v s +1 J s +1 ) − ∇ · (cid:18) πc σ s ∇ J s +1 (cid:19) = cσ a (( T s +1 ) − πJ s +1 ) + 4 π v s +1 ∇ · J s +1 , (3.12)which is a semi-discretization (2.10d).3.3. The time discretization for the equilibrium regime.
We substitute the decomposition in (3.1)into the RMHD system (2.11), and the time discretization of (2.11) is same with (3.6) except the RTE andsource terms ( (cid:99) S re ) s +1 and ( (cid:100) S rp ) s +1 , which read: J s +1 − J s ∆ t + ε n R s +1 − R s ∆ t + ε Q s +1 − Q s ∆ t + cε n ∇ · J s +1 + c n ∇ · ( n R s +1 + Q s ) = cσ a ε (cid:18) ( T s +1 ) π − J s +1 (cid:19) − (cid:16) cσ a ε + cεσ s (cid:17) ( n R s +1 + Q s +1 ) − ε σ s v s +1 R s +1 − ( σ a − ε σ s ) v s +1 · v s +1 c (cid:18) J s +1 + εK sQ (cid:19) + 3 σ a ε nv s +1 (cid:18) ( T s +1 ) π − J s +1 (cid:19) + nv s +1 (cid:16) σ a ε + εσ s (cid:17) (4 J s +1 + ε ( n R s +1 + Q s )) , (3.13)and ( (cid:99) S re ) s +1 = σ a ε (( T s +1 ) − πJ s +1 ) + 4 π ( εσ a − ε σ s ) v s +1 c (cid:20) c R s +1 − v s +1 (cid:18) J s +1 + εK sQ (cid:19)(cid:21) , (3.14)( (cid:100) S rp ) s +1 = − π ( ε σ s + σ a ) c (cid:20) c R s +1 − v s +1 (cid:18) J s +1 + εK sQ (cid:19)(cid:21) − v s +1 c σ a (( T s +1 ) − πJ s +1 ) . (3.15)Moreover, the procedure of solving the equilibrium case is the same as for (3.6). The diffusion limit of the equilibrium regime:
The derivation of the equilibrium diffusion limit issimilar to the non-equilibrium case, but requires additional analysis for the zeroth moment of equation(3.13). Sending ε → σ a (( T s +1 ) − πJ s +1 ) + O ( ε ) = 0 , which implies 4 πJ s +1 = ( T s +1 ) + O ( ε ) . Therefore, (3.9) and (3.10) in the equilibrium regime become a ( ρT ) s +1 − ( ρT ) s ∆ t + ρ s +1 ( v s +1 ) − ρ s ( v s ) t + ( B s +1 ) − ( B s ) t + P ( T s +1 ) − ( T s ) ∆ t + ∇ · (cid:20) ( E s + ( P s ) ∗ ) v s − B s ( B s · v s ) + 4 P v s +1 ( T s +1 ) (cid:21) = ∇ · (cid:18) πc P σ a ∇ ( T s +1 ) (cid:19) , ( ρ v ) s +1 − ( ρ v ) s ∆ t + ∇ · ( ρ s v s v s − B s B s + ( P ∗ ) s ) = − P ∇ · ( T s +1 ) , which gives a semi-discretization of system (2.14).4. The full discretization for the non-equilibrium regime
For the ease of exposition, we will explain our spatial discretion in 1D. That is, x ∈ [0 , L ], n ∈ [ − , (cid:104) f ( n ) (cid:105) = (cid:82) − f ( n )d n . Recall the RMHD for the non-equilibrium case in the slab geometry: ∂ t I + cε n∂ x I = cσ a (cid:18) T π − I (cid:19) + cσ s ε ( J − I ) + 3 εnv x σ a (cid:18) T π − J (cid:19) + nv x (cid:16) εσ a + σ s ε (cid:17) ( I + 3 J ) − σ s ε v x · H − ( ε σ a − σ s ) v x · v x c J − ( ε σ a − σ s ) v x · ( v x · K ) c , (4.1a) ∂ t ρ + ∂ x ( ρv x ) = 0 , (4.1b) ∂ t ( ρv x ) + ∂ x (cid:0) ρv x + p + B / B x (cid:1) = −P S rp , (4.1c) ∂ t ( ρv y ) + ∂ x ( ρv x v y − B x B y ) = 0 , (4.1d) ∂ t ( ρv z ) + ∂ x ( ρv x v z − B x B z ) = 0 , (4.1e) ∂ t (cid:18) aρT + ρv + B (cid:19) + ∂ x [( E + P ∗ ) v x − ( B · v ) B x ] = − cε P S re , (4.1f) ∂ t B y + ∂ x ( B y v x − B x v y ) = 0 , (4.1g) ∂ t B z + ∂ x ( B z v x − B x v z ) = 0 , (4.1h)where S re = εσ a (cid:18) T − π (cid:90) − I d n (cid:19) + 2 π ( ε σ a − εσ s ) v x c (cid:20) cε (cid:90) − nI d n − (cid:18) v x (cid:90) − I d n + v x (cid:90) − n I d n (cid:19)(cid:21) ,S rp = − π ( σ s + ε σ a ) c (cid:20) cε (cid:90) − nI d n − (cid:18) v x (cid:90) − I d n + v x (cid:90) − n I d n (cid:19)(cid:21) + v x c ε σ a (cid:18) T − π (cid:90) − I d n (cid:19) , and v x , v y and v z are fluid velocity in the x − , y − , and z − directions, respectively. B x , B y and B z are magneticflux in the x − , y − , and z − directions, respectively. The boundary condition becomes I ( t, , n ) = b L ( t, n ) , for n > I ( t, L, n ) = b R ( t, n ) , for n < . (4.2)Higher dimensions can be treated in the dimension by dimension manner. Let [ a, b ] be the computationaldomain, ∆ x = ( b − a ) /N x , and we consider the uniform mesh as follows x i − = a + ( i − x, i = 1 , · · · , N x + 1 , and let x = a < x < x < · · · < x i − < x i < x i + < · · · < x N x < x N x + = b,x i = (cid:16) x i − + x i + (cid:17) / , for i = 1 , · · · , N x . To get a consistent stencil in spatial discretization, we use the unified gas kinetic scheme (UGKS) forspatial discretization [22]. Other space discretization that is AP can be applied as well. The most crucialpoint is that the space discretization of the MHD system has to be consistent with that of RTE.
SPATIAL-TEMPORAL ASYMPTOTIC PRESERVING SCHEME FOR RADIATION MAGNETOHYDRODYNAMICS 11
A finite volume approach.
UGKS is a finite volume method, integrating RMHD system (4.1) over[ t s , t s +1 ] and [ x i − , x i + ] gives J s +1 i − J si ∆ t + εn R s +1 i − R si ∆ t + ε Q s +1 i − Q si ∆ t + 1∆ x ( ζ i + − ζ i − ) = cσ a (cid:18) ( T s +1 i ) π − J s +1 i (cid:19) − (cid:16) cεσ a + cσ s ε (cid:17) ( nR s +1 i + Q s +1 i ) + G i , (4.3a) ρ s +1 i − ρ si ∆ t + 1∆ x ( F s ,i + − F s ,i − ) = 0 , (4.3b)( ρv x ) s +1 i − ( ρv x ) si ∆ t + 1∆ x ( F s ,i + − F s ,i − ) = −P ( (cid:100) S rp ) s +1 i , (4.3c)( ρv y ) s +1 i − ( ρv y ) si ∆ t + 1∆ x ( F s ,i + − F s ,i − ) = 0 , (4.3d)( ρv z ) s +1 i − ( ρv z ) si ∆ t + 1∆ x ( F s ,i + − F s ,i − ) = 0 , (4.3e) a ( ρT ) s +1 i − ( ρT ) si ∆ t + ( ρv + B ) s +1 i − ( ρv + B ) si t + 1∆ x ( F s ,i + − F s ,i − ) = − cε P ( (cid:99) S re ) s +1 i , (4.3f) B s +1 y,i − B sy,i ∆ t + 1∆ x ( F s ,i + − F s ,i − ) = 0 , (4.3g) B s +1 z,i − B sz,i ∆ t + 1∆ x ( F s ,i + − F s ,i − ) = 0 , (4.3h)where G , ( (cid:99) S re ) s +1 i and ( (cid:100) S rp ) s +1 i are G i =3 εσ a,i nv s +1 x,i (cid:18) ( T s +1 i ) π − J s +1 i (cid:19) + nv s +1 x,i (cid:16) εσ a,i + σ s,i ε (cid:17) (cid:0) J s +1 i + ε ( nR s +1 i + Q si ) (cid:1) − σ s,i v s +1 x,i R s +1 i − ( ε σ a,i − σ s,i )( v s +1 x,i ) c (cid:18) J s +1 i + εK sQ,i (cid:19) , ( (cid:99) S re ) s +1 i = 4 πεσ a,i (cid:18) ( T s +1 i ) π − J s +1 i (cid:19) + 4 π ( ε σ a,i − εσ s,i ) v s +1 x,i c (cid:20) c R s +1 i − v s +1 x,i (cid:18) J s +1 i + εK sQ,i (cid:19)(cid:21) , ( (cid:100) S rp ) s +1 i = − π ( σ s,i + ε σ a,i ) c (cid:20) c R s +1 i − v s +1 x,i (cid:18) J s +1 i + εK sQ,i (cid:19)(cid:21) + 4 πε σ a,i v s +1 x,i c (cid:18) ( T s +1 i ) π − J s +1 i (cid:19) . Here the numerical fluxes F s · ,i ± at the interfaces use Roe’s Riemann solver [29], based on the piecewiselinear reconstruction. At last, we discuss the microscopic flux ζ defined at the interface x i + in (4.3a). In[22], all the terms of the microscopic flux are treated explicitly, in this paper all the terms except the initialvalue and K Q are treated implicitly , and the processes of UGKS are detailed in Appendix A.With the expressions (A.3), (A.4) and (A.6), the numerical flux ζ i + is ζ i + = A i + n (cid:16) I s, + i + n> + I s, − i + n< (cid:17) + C i + nJ s +1 i + + C i + π n ( T s +1 i + ) + F i + n (cid:98) G i + + D i + n (cid:16) δ x J s +1 , + i + n> + δ x J s +1 , − i + n< (cid:17) + D i + π n (cid:16) δ x ( T s +1 , + i + ) n> + δ x ( T s +1 , − i + ) n< (cid:17) , (4.4)where n ≶ is the indicator function. Here the two approximations for I , from the left and from the right,at the interface x i + are given by: I s, + i + = I si + ∆ x δ x I si , I s, − i + = I si +1 − ∆ x δ x I si +1 . (4.5) Moreover, the coefficients in the numerical flux are given by: A = c ∆ tεµ (1 − e − µ ∆ t ) ,C = c σ s ∆ tε µ (cid:18) ∆ t − µ (1 − e − µ ∆ t ) (cid:19) , C = c σ a ∆ tεµ (cid:18) ∆ t − µ (1 − e − µ ∆ t ) (cid:19) ,D = − c σ s ∆ tε µ (cid:18) ∆ t (cid:0) e − µ ∆ t (cid:1) − µ (cid:0) − e − µ ∆ t (cid:1)(cid:19) ,D = − c σ a ∆ tε µ (cid:18) ∆ t (cid:0) e − µ ∆ t (cid:1) − µ (cid:0) − e − µ ∆ t (cid:1)(cid:19) ,F = c ∆ tεµ (cid:18) ∆ t − µ (cid:0) − e − µ ∆ t (cid:1)(cid:19) , (4.6)with µ = ε cσ a + cσ s ε . For reader’s convenience, the processes of UGKS are detailed in Appendix A.To obtain a scheme that can update the quantities J , R , T and v x , the zeroth moment and first momentof velocity field n for (4.3a) are needed. Integrating equation (4.3a) for n from -1 to 1, and multiplying theboth sides by 4 π , one can obtain4 π J s +1 i − J si ∆ t + 2 π ∆ x (cid:90) − (cid:98) ζ i + − (cid:98) ζ i − dn = cε ( (cid:99) S re ) s +1 i , (4.7)where (cid:90) − (cid:98) ζ i + dn = (cid:90) − A i + n (cid:16) I s, + i + n> + I s, − i + n< (cid:17) dn + 2 D i + J s +1 i +1 − J s +1 i ∆ x + D i + π ( T s +1 i +1 ) − ( T s +1 i ) ∆ x + F i + εσ a,i + v s +1 x,i + ( T s +1 i + ) π − J s +1 i + + v s +1 x,i + σ s,i + + ε σ a,i + ε (cid:18) J s +1 i + + 2 εK sQ,i + (cid:19) . Multiplying equation (4.3a) by n , then integrating it for n from -1 to 1, and multiplying the both sides by4 π , one can obtain 4 πε R s +1 i − R si ∆ t + 2 π ∆ x (cid:90) − n (cid:16)(cid:98) ζ i + − (cid:98) ζ i − (cid:17) dn = cε ( (cid:100) S rp ) s +1 i , (4.8)where (cid:90) − n (cid:98) ζ i + dn = (cid:90) − A i + n (cid:16) I s, + i + n> + I s, − i + n< (cid:17) dn + C i + J s +1 i +1 + J s +1 i ) + C i + π (( T s +1 i +1 ) + ( T s +1 i ) )+ F i + (cid:20) v s +1 x,i + (cid:16) ε σ a,i + + σ s,i + (cid:17) (cid:18) R s +1 i + + (cid:90) − n Q si + dn (cid:19) − σ s,i + v s +1 x,i + R s +1 i + − ε σ a,i + − σ s,i + )( v s +1 x,i + ) c (cid:18) J s +1 i + + εK sQ,i + (cid:19) . Coupling equations (4.7), (4.8), (4.3c) and (4.3f) as a system and solving the system, quantities J , T , R and v x can be updated. Finally, quantity Q can be updated by equation (4.3a), and the new I can be given byusing (3.1). SPATIAL-TEMPORAL ASYMPTOTIC PRESERVING SCHEME FOR RADIATION MAGNETOHYDRODYNAMICS 13
The diffusion limit of (4.3c) – (4.8) . Firstly, assuming σ a and σ s are positive, as ε →
0, the leadingorder of the coefficients in the numerical flux (4.4) are A (∆ t, ε, c, σ a , σ s ) = 0 + O ( ε ) ,εC (∆ t, ε, c, σ a , σ s ) = c + O ( ε ) , εC = 0 + O ( ε ) ,D = − cσ s + O ( ε ) , D = 0 + O ( ε ) , ε F = 1 σ s + O ( ε ) . This means the zero moment and the first moment of the numerical flux, as ε →
0, have the following limit: (cid:90) − (cid:98) ζ i + dn → v s +1 x,i + J s +1 i + − c σ s,i + J s +1 i +1 − J s +1 i ∆ x , (4.9) ε (cid:90) − n (cid:98) ζ i + dn → c J s +1 i +1 + J s +1 i ) . (4.10)Multiplying equation (4.7) by P , and adding it up with equation (4.3f), the terms on the right hand sidescancel. By sending ε → a ( ρT ) s +1 i − ( ρT ) si ∆ t + ( ρv + B ) s +1 i − ( ρv + B ) si t + 4 π P J s +1 i − J si ∆ t + 1∆ x ( F s ,i + − F s ,i − )+ 16 π P v s +1 x,i + J s +1 i + − v s +1 x,i − J s +1 i − ∆ x = 4 π P c σ s,i + ∆ x J s +1 i +1 − J s +1 i ∆ x − π P c σ s,i − ∆ x J s +1 i − J s +1 i − ∆ x . (4.11)This is a fully discretization for (2.10c). Then multiplying equation (4.8) by εc P , and adding it to equation(4.3c), sending ε → ρv x ) s +1 i − ( ρv x ) si ∆ t + 1∆ x ( F s ,i + − F s ,i − ) + 2 π P J s +1 i +1 − J s +1 i − ∆ x = 0 . (4.12)This is a fully discretization for (2.10b). Sending ε → R s +1 i = 4 c v s +1 x,i J s +1 i − J s +1 i +1 − J s +1 i − xσ s,i , (4.13)using (4.13) in (4.7) and sending ε →
0, one gets4 π J s +1 i − J si ∆ t − π ∆ x (cid:34) c σ s,i + J s +1 i +1 − J s +1 i ∆ x − c σ s,i − J s +1 i − J s +1 i − ∆ x (cid:35) + 16 π v s +1 x,i + J s +1 i + − v s +1 x,i − J s +1 i − ∆ x = cσ a,i (cid:0) ( T s +1 i ) − πJ s +1 i (cid:1) − π v s +1 x,i J s +1 i +1 − J s +1 i − ∆ x , (4.14)which is a fully discretization for (2.10d).4.3. Boundary conditions.
In this section, the numerical fluxes in the system (4.3c)–(4.8) under theboundary condition (4.2) are considered. The construction of the flux on the right boundary is similar to theconstruction on the left boundary, for which we only consider the left boundary case. At the left boundary,the integral representation of the radiative intensity I (A.3) reads: I ( t ) = b L , if n > ,e − µ ( t − t s ) I (cid:16) t s , x + cnε ( t − t s ) (cid:17) + − e − µ
12 ( t − ts ) µ (cid:98) G + t (cid:82) t s e − µ ( t − t s ) (cid:16) cσ s, ε J (cid:16) z, x − cnε ( t − z ) (cid:17) + cσ a, π T (cid:16) z, x − cnε ( t − z ) (cid:17)(cid:17) dz, if n < . (4.15) According to the approximation (A.4) and (A.6) , the reconstruction of I for n < x > x at the leftboundary can be written as: I ( t s , x, n ) = I s , and the reconstruction of J for t in the interval [ t s , t s +1 ]: J ( t, x ) = J s +1 + δ x J s +1 , − ( x − x ) , where J s +1 = ( (cid:104) b L (cid:105) + J s +11 ) and δ x J s +1 , − = J s +11 − J s +112 ∆ x/ . Then the numerical flux at the left boundaryreads: ζ = cnε b L n> + A nI s n< + C nJ s +1 n< + C π n ( T s +1 ) n< + F n (cid:98) G n< + D n δ x J s +1 , − n< + D π n δ x ( T s +1 , − ) n< , which indicates (cid:90) − (cid:98) ζ dn = (cid:90) − (cid:16) cnε b L n> + A nI s n< (cid:17) dn − (cid:32) C (cid:0) J s +11 + (cid:104) b L (cid:105) (cid:1) + C π (( T s +11 ) + ( T s +1 L ) ) (cid:33) − F (cid:98) G x (cid:32) D ( J s +11 − (cid:104) b L (cid:105) ) + D π (( T s +11 ) − ( T s +1 L ) ) (cid:33) , here T s +1 L is the material temperature at the left boundary, and (cid:90) − n (cid:98) ζ dn = (cid:90) − (cid:18) cn ε b L n> + A n I s n< (cid:19) dn + 13 (cid:32) C (cid:0) J s +11 + (cid:104) b L (cid:105) (cid:1) + C π (( T s +11 ) + ( T s +1 L ) ) (cid:33) + F (cid:98) G − x (cid:32) D ( J s +11 − (cid:104) b L (cid:105) ) + D π (( T s +11 ) − ( T s +1 L ) ) (cid:33) . To summarize, we have the following one time step update of the fully discrete version of RMHD.
Algorithm 1: one step of fully discrete update for RMHD
Input: ρ si , v sx,i , v sy,i , v sz,i , T si , B sy,i , B sz,i , I si , J si ( J si = (cid:104) I si (cid:105) ) Output: ρ s +1 i , v s +1 x,i , v s +1 y,i , v s +1 z,i , T s +1 i , B s +1 y,i , B s +1 z,i , I s +1 i , J s +1 i ( J s +1 i = (cid:10) I s +1 i (cid:11) ) obtain ρ s +1 i , v s +1 y,i , v s +1 z,i , B s +1 y,i , B s +1 z,i from (4.3b), (4.3d), (4.3e), (4.3g) and (4.3h); obtain R si , Q si from (3.2) and (3.1); get J s +1 i , T s +1 i , v s +1 x,i and R s +1 i from the system coupled by (4.3c), (4.3f), (4.7) and (4.8); get Q s +1 i from equation (4.3a) and Q si replaced by Q s +1 i and K sQ,i replaced by K s +1 Q,i ; obtain I s +1 i from (3.1). 5. Numerical Examples
In this section, we conduct numerical experiments to test the performance of our proposed method. Theexamples cover both optically thin ε ∼ O (1) and optically thick ε (cid:28) S discrete ordinate method in the angle space. SPATIAL-TEMPORAL ASYMPTOTIC PRESERVING SCHEME FOR RADIATION MAGNETOHYDRODYNAMICS 15
Example 1.
In this example, we will test an ideal MHD shock tube problem as in Section V in [2],where P =0 in (2.2). In this case, the radiation does not affect the fluid, but the fluid provides source forradiation. The computational domain is [ − ,
1] and the initial data is given by( ρ, v x , v y , v z , B y , B z , p )( x,
0) = (cid:40) (1 , , , , , , , x < , (0 . , , , , − , , . , x > . (5.1)The magnetic intensity in the x direction B x = 0 .
75, the ideal gas constant R ideal = 1 and adiabatic indexfor an ideal gas γ = 2, the absorption collision cross-section σ a ( x ) = 1 /
3, the scatter collision cross-section σ s ( x ) = 1 / c = 0 . ε = 10 − in the RMHD system (2.2). The exact solutions at time t > t = 0 . J is solved using the diffusion limit equation (2.10) by the semi-implicit scheme (B.7) and (B.8). For all theschemes, the time step is ∆ t = 0 . x. Example 2.
The uniform accuracy and stability of the AP method at the equilibrium regime are testedin this example. The initial data is the same as in (5.1) in example 1, and the initial temperature and densityfluxes are T (0 , x ) = p/ ( R ideal ρ ) , I (0 , x, n ) = T (0 , x ) π . (5.2)Moreover, we assume that P = 0 . , c = 0 . σ a and the scatter collision cross-section σ s are chosen as σ a ( x ) = σ s ( x ) = 1 / Accuracy . The uniform convergence order of the diffusion limit system (2.14) and the RMHD system (2.11)of our AP scheme are given. In Fig. 3, we ploterror ρ = (cid:107) ρ ∆ x ( · , t max ) − ρ ∆ x/ ( · , t max ) (cid:107) l , error p = (cid:107) p ∆ x ( · , t max ) − p ∆ x/ ( · , t max ) (cid:107) l , error v x = (cid:107) ( v x ) ∆ x ( · , t max ) − ( v x ) ∆ x/ ( · , t max ) (cid:107) l , error v y = (cid:107) ( v y ) ∆ x ( · , t max ) − ( v y ) ∆ x/ ( · , t max ) (cid:107) l , error B y = (cid:107) ( B y ) ∆ x ( · , t max ) − ( B y ) ∆ x/ ( · , t max ) (cid:107) l , error T = (cid:107) T ∆ x ( · , t max ) − T ∆ x/ ( · , t max ) (cid:107) l , (5.3)for the diffusion limit system (2.14) with decreasing ∆ x and t max = 0 .
2. Then in Figs. 4, we plot the sameerrors for RMHD system (2.11) with decreasing ∆ x and t max = 0 . ε = 1,10 − , 10 − ,10 − . Stability condition . To check the stability, we choose ∆ t = C ∆ x with C being constant, and record thelargest possible choice of C in Table 1 for different ε and ∆ x . Uniform hyperbolic stability is observednumerically. Table 1.
Stability test for the RMHD system (2.11). Here we use ∆ t = C ∆ x and thelargest possible C that stabilize the scheme are displayed. ε ∆ x . . . . . . . . . . . . . . . . . . . . -1 0 1 x Reference800 grids (a) Density ρ -1 0 1 x P Reference800 grids (b) Pressure p -1 0 1 x -0.4-0.200.20.40.60.8 v x Reference800 grids (c) Velocity in x-direction v x -1 0 1 x -1.5-1-0.50 v y Reference800 grids (d) Velocity in y-direction v y -1 0 1 x -1-0.500.51 B y Reference800 grids (e) Magnetic field in y-direction B y -1 0 1 x P / Reference800 grids (f) The ratio of the pressure and the density
Figure 1.
The numerical results at time t = 0 . t = 0 . x . Comparison with reference solutions . We consider two cases with ε = 1 and ε = 10 − as in [3]. Tocompute the reference solution with ε = 1, use a fully resolved explicit finite volume solver. For ε (cid:28)
1, wesolve the diffusion limit equation (2.10) by the semi-implicit scheme (B.7) and (B.8). The results are shownin Fig. 5 and Fig. 6. We can see that our results agree well with the reference solutions.5.3.
Example 3.
In the simulation of RMHD, it is a challenging task to produce accurate radiative shocks,especially in the optically thick regime. We test a range of the radiation-hydrodynamic shock problemspresented in [20] in which several shocks are presented in the solution. For each of these shocks, we set P = 10 − , adiabatic index for an ideal gas γ = 5 /
3, absorption collision cross-section σ a = 1 /
3, scatter
SPATIAL-TEMPORAL ASYMPTOTIC PRESERVING SCHEME FOR RADIATION MAGNETOHYDRODYNAMICS 17 -1 0 1 x J Reference800 grids
Figure 2.
Comparison of J of our AP scheme and the reference solution obtained by semi-implicit solver (B.7)–(B.8) for RMHD system (2.2). Here P = 0, ε = 10 − , t = 0 . x = 1 / t = 0 . x for both two schemes. -3 -2 x -3 -2 -1 slope=1error in error in v x error in v y error in perror in B y error in T Figure 3.
Plot of errors (5.3) for diffusion limit system (2.14) with decreasing ∆ x = , , , , . Here ∆ t = 0 . x .collision cross-section σ s = 10 , c = 3 × , ε = 1 in RMHD system (2.2), and the radiation temperature T r ≡ (4 πJ ) . . Moreover, the computational domain is [ − . , .
02] and the results at t = 0 .
04 are displayed.In this example, ∆ x = 1 / , ∆ t = 0 . x . In [20], the authors have obtained some semi-analytic solutionsfor the non-equilibrium diffusion limit system (2.10) at different Mach numbers, we compare the numericalresults with the semi-analytic solutions obtained [20].( ρ, v x , v y , v z , B y , B z , T, T r )( x,
0) = (cid:40) (1 , . , , , , , , , x < , (1 . , . , , , , , . , . , x > , (5.4)( ρ, v x , v y , v z , B y , B z , T, T r )( x,
0) = (cid:40) (1 , , , , , , , , x < , (2 . , . , , , , , . , . , x > . (5.5)First of all, a Mach 1.2 shock is considered, which has no isothermal sonic point (ISP) but a hydrodynamicshock. The initial conditions are shown in (5.4). Fig. 7 compares our numerical results with the semi-analytic -3 -2 x -3 -2 -1 =1 slope=1error in error in v x error in v y error in perror in B y error in T (a) ε = 1 -3 -2 x -3 -2 -1 =0.1 slope=1error in error in v x error in v y error in perror in B y error in T (b) ε = 0 . -3 -2 x -3 -2 -1 =10 -3 slope=1error in error in v x error in v y error in perror in B y error in T (c) ε = 10 − -3 -2 x -3 -2 -1 =10 -5 slope=1error in error in v x error in v y error in perror in B y error in T (d) ε = 10 − Figure 4.
Plot of errors (5.3) for RMHD system (2.11) with decreasing ∆ x = , , , , . Here ∆ t = 0 . x , ε = 1 , . , − and 10 − .solutions, in which we can see good agreement for all quantities, including density, velocity, material andradiation temperature. Due to the hydrodynamic shock, there exhibits discontinuities in the solution profilesof density, velocity and material temperature, and the maximum material temperature is bounded, since thereis no ISP to drive it further.At Mach 2, there are both a hydrodynamic shock and an ISP. The initial conditions are given in (5.5).Fig. 8 compares our numerical results with the semi-analytic solutions. One can see that our numerical resultsare in good agreement with the semi-analytic solutions. In Fig. 8, discontinuities can be seen in materialdensity, velocity and material temperature due to the hydrodynamic shock. Moreover, the Zel’dovich spikecan be observed in material temperature as in Fig. 8(c). The Zel’dovich spike is caused by the ISP embeddedwithin the hydrodynamic shock, which drives up the material temperature at the shock front.6. Conclusion and discussions
In this paper, we have introduced an AP scheme in both space and time for the RMHD system thatcouples the ideal MHD equations with a gray RTE, and two different scalings are considered, one resultsin an equilibrium diffusion limit system, while the other results in nonequilibrium. The main idea is todecompose the intensity into three parts, two parts correspond to the zeroth and first order moments, while
SPATIAL-TEMPORAL ASYMPTOTIC PRESERVING SCHEME FOR RADIATION MAGNETOHYDRODYNAMICS 19 -1 0 1 x explicitAP (a) Density ρ -1 0 1 x P explicitAP (b) Pressure p -1 0 1 x -0.4-0.200.20.40.60.8 v x explicitAP (c) Velocity in x-direction v x -1 0 1 x -1.5-1-0.500.5 v y explicitAP (d) Velocity in y-direction v y -1 0 1 x -1-0.500.51 B y explicitAP (e) Magnetic field in y-direction B y -1 0 1 x P / explicitAP (f) The ratio of the pressure and the density Figure 5.
Comparison of the results of our AP scheme and the reference solution obtainedby fully explicit scheme for RMHD system (2.11). Here ε = 1, t = 0 .
2, ∆ x = 1 / t = 0 . x for the AP scheme, and ∆ x = 1 / t = 0 . x for the explicit solver. -1 0 1 x implicitAP (a) Density ρ -1 0 1 x P implicitAP (b) Pressure p -1 0 1 x -0.4-0.200.20.40.60.8 v x implicitAP (c) Velocity in x-direction v x -1 0 1 x -2-1.5-1-0.500.5 v y implicitAP (d) Velocity in y-direction v y -1 0 1 x -1-0.500.51 B y implicitAP (e) Magnetic field in y-direction B y -1 0 1 x P / implicitAP (f) The ratio of the pressure and the density Figure 6.
Comparison of the results of our AP scheme and the reference solution obtainedby the semi-implicit solver (B.7)–(B.8) for RMHD system (2.11). Here ε = 10 − , t = 0 . x = 1 / t = 0 . x for both two schemes. SPATIAL-TEMPORAL ASYMPTOTIC PRESERVING SCHEME FOR RADIATION MAGNETOHYDRODYNAMICS 21 -0.02 -0.01 0 0.01 0.02 x Semi-analyticnumerical (a) Density ρ -0.02 -0.01 0 0.01 0.02 x v Semi-analyticnumerical (b) Velocity v -0.02 -0.01 0 0.01 0.02 x T Semi-analyticnumerical (c) Temperature T -0.02 -0.01 0 0.01 0.02 x T r Semi-analyticnumerical (d) Radiation temperature T r Figure 7.
Comparison of the results using our AP scheme at time t = 0 .
04 and the semi-analytic solution for Mach number M = 1.2. For AP scheme, we use ∆ x = 1 /
800 and∆ t = 0 . x .the third part is the residual. The two macroscopic moments are treated implicitly while the residual isexplicit. For the space discretization, we use Roe’s method in the Athena code to solve the convective partin the ideal MHD equations and the UGKS for the RTE. Numerical results are presented in the opticallythick region and the radiative shock problem, which are compared to the semi-analytic solution to verifyaccuracy and asymptotic properties of our method.It would be worthwhile to construct a second-order AP method for the radiation magnetohydrodynamicssystem. In the semi-implicit scheme, a simple time-integration method has been used. To obtain a secondscheme, higher order IMEX time-integration method can be considered. Acknowledgement:
The authors would like to thank Dr. Yanfei Jiang from Computational Center forAstrophysics, the Flatiron Institute, for proposing this problem and useful discussions. S. Jin is partiallysupported by NSFC12031013 and the Strategic Priority Research Program of Chinese Academy of Sciences,XDA25010404; M. Tang and X. J. Zhang are partially supported by NSFC11871340 and the Strategic PriorityResearch Program of Chinese Academy of Sciences, XDA25010401. -0.02 -0.01 0 0.01 0.02 x Semi-analyticnumerical (a) Density ρ -0.02 -0.01 0 0.01 0.02 x v Semi-analyticnumerical (b) Velocity v -0.02 -0.01 0 0.01 0.02 x T Semi-analyticnumerical (c) Temperature T -0.02 -0.01 0 0.01 0.02 x T r Semi-analyticnumerical (d) Radiation temperature T r Figure 8.
Comparison of the results using our AP scheme at time t = 0 .
04 and the semi-analytic solution for Mach number M = 2. For AP scheme, we use ∆ x = 1 /
800 and∆ t = 0 . x . Appendix A. The processes of UGKS for the non-equilibrium regime
Recall the RTE (4.1a) in the slab geometry, and integrating it over [ t s , t s +1 ] and [ x i − , x ] gives I s +1 i − I si ∆ t + 1∆ x ( ζ i + − ζ i − ) = cσ a,i (cid:18) ( T s +1 i ) π − I s +1 i (cid:19) + cσ s,i ε ( J s +1 i − I s +1 i ) + G. The microscopic flux ζ is computed by solving the following initial value problem at the cell boundary x = x i + : I s +1 i − I si ∆ t + 1∆ x ( ζ i + − ζ i − ) = cσ a,i (cid:18) ( T s +1 i ) π − I s +1 i (cid:19) + cσ s,i ε ( J s +1 i − I s +1 i ) + G,I ( t, x, n ) | t = t s = I ( t s , x, n ) . (A.1)In the UGKS scheme, we firstly assume the coefficients σ a and σ s in space and time are piecewise constant,thus the radiative transfer equation in RMHD system (4.1) is equivalent to ddt e µt I (cid:16) t, x + cnε t, n (cid:17) = e µt (cid:16) cσ s ε J (cid:16) t, x + cnε t, n (cid:17) + cσ a π T (cid:16) t, x + cnε t, n (cid:17) + G (cid:17) , (A.2) SPATIAL-TEMPORAL ASYMPTOTIC PRESERVING SCHEME FOR RADIATION MAGNETOHYDRODYNAMICS 23 with µ = ε cσ a + cσ s ε . Solving the equivalent equation (A.2) with the initial value in the system (A.1), whichgives the solution I ( t, x i + , n ) ≈ e − µ i + 12 ( t − t s ) I (cid:16) t s , x i + + cnε ( t − t s ) (cid:17) + 1 − e − µ i + 12 ( t − t s ) µ i + (cid:98) G i + + t (cid:90) t s e − µ i + 12 ( t − t s ) (cid:18) cσ s,i + ε J (cid:16) z, x i + − cnε ( t − z ) (cid:17) + cσ a,i + π T (cid:16) z, x i + − cnε ( t − z ) (cid:17)(cid:19) dz, (A.3)with µ i + , σ a,i + and σ s,i + being the constant value at the corresponding cell boundary for µ , σ a and σ s respectively. In solution (A.3), there remains three terms that need to approximate: the first one is theinitial condition I ( t s ) around x i + , namely the function I (cid:16) t s , x i + + cnε ( t − t s ) (cid:17) ; the second one is the twofunctions J, T localized in the time interval [ t s , t s +1 ] and around the boundary x i + , i.e., the functions J (cid:16) z, x i + − cnε ( t − z ) (cid:17) and T (cid:16) z, x i + − cnε ( t − z ) (cid:17) ; the last one is the term (cid:98) G i + . The first term:
For the first term, a piecewise linear reconstruction function is used to approximate theinitial function I ( t s ) around x i + : I ( t s , x, n ) = (cid:40) I si + δ x I si ( x − x i ) , if x < x i + ,I si +1 + δ x I si +1 ( x − x i +1 ) , if x > x i + . (A.4)Here δ x I si is the limited slope determined by the minmod limiter [26]: δ x I si = min { I si +1 − I si ∆ x , I si +1 − I si − x , I si − I si − ∆ x } , if all are positive , max { I si +1 − I si ∆ x , I si +1 − I si − x , I si − I si − ∆ x } , if all are negative , , otherwise . (A.5) The second term:
For the second term, the approximations of two functions J and T between t s and t s +1 and around x i + are constructed in the same manner, so only the approximation of J is introduced in detail.The function J is treated implicitly in time by using piecewise linear polynomials, so the reconstruction for J reads: J ( t, x ) = J s +1 i + + δ x J s +1 , + i + ( x − x i + ) , if x < x i + ,J s +1 i + + δ x J s +1 , − i + ( x − x i + ) , if x > x i + , (A.6)with the interface value of J defined by J s +1 i + = (cid:68) I s +1 i + ( n ) (cid:69) = 12 (cid:10) I s +1 i + I s +1 i +1 (cid:11) , and with the spatial derivatives given by δ x J s +1 , + i + = J s +1 i + − J s +1 i ∆ x/ , δ x J s +1 , − i + = J s +1 i +1 − J s +1 i + ∆ x/ . The third term:
At last, the approximation for the term (cid:98) G is discussed. Recall the definition of G i : G i =3 εσ a,i nv s +1 x,i (cid:18) ( T s +1 i ) π − J s +1 i (cid:19) + nv s +1 x,i (cid:16) εσ a,i + σ s,i ε (cid:17) (cid:0) J s +1 i + ε ( nR s +1 i + Q si ) (cid:1) − σ s,i v s +1 x,i R s +1 i − ( ε σ a,i − σ s,i )( v s +1 x,i ) c (cid:18) J s +1 i + εK sQ,i (cid:19) , The upwind scheme is used to determine (cid:98) G , i.e., the value of (cid:98) G at the center of the interval [ x i , x i +1 ] isdetermined by the velocity field v at the both sides boundary, If v s +1 x,i ≥ v s +1 x,i +1 : (cid:98) G i + = G i , if ( G i +1 − G i ) / (cid:0) v s +1 x,i +1 − v s +1 x,i (cid:1) > , (cid:98) G i + = G i +1 , if ( G i +1 − G i ) / (cid:0) v s +1 x,i +1 − v s +1 x,i (cid:1) ≤ v s +1 x,i < v s +1 x,i +1 : (cid:98) G i + = G i , if v s +1 x,i > , (cid:98) G i + = G i +1 , if v s +1 x,i +1 < , (cid:98) G i + = 0 , if v s +1 x,i ≤ ≤ v s +1 x,i +1 . (A.7) Appendix B. The full discretization for the equilibrium regime
Here, we will show the full discretization for system (2.11). Applying semi-implicit finite volume schemein system (2.11), the fully discretization of (2.11) is same with (4.3) except the RTE and source terms S re and S rp , which read: J s +1 i − J si ∆ t + ε nR s +1 i + Q s +1 i − nR si − Q si ∆ t + 1∆ x ( (cid:101) ζ i + − (cid:101) ζ i − ) = cσ a,i ε (cid:18) ( T s +1 i ) π − J s +1 i (cid:19) − (cid:16) cσ a,i ε + cεσ s,i (cid:17) ( nR s +1 i + Q s +1 i ) + (cid:101) G i , with the source term: (cid:101) G i =3 σ a,i ε nv s +1 x,i (cid:18) ( T s +1 i ) π − J s +1 i (cid:19) + nv s +1 x,i (cid:16) σ a,i ε + εσ s,i (cid:17) (cid:0) J s +1 i + ε ( nR s +1 i + Q si ) (cid:1) − ε σ s,i v s +1 x,i R s +1 i − ( σ a,i − ε σ s,i )( v s +1 x,i ) c (cid:18) J s +1 i + εK sQ,i (cid:19) , and( (cid:102) S re ) s +1 i = 4 π σ a,i ε (cid:18) ( T s +1 i ) π − J s +1 i (cid:19) + 4 π ( εσ a,i − ε σ s,i ) v s +1 x,i c (cid:20) c R s +1 i − v s +1 x,i (cid:18) J s +1 i + εK sQ,i (cid:19)(cid:21) , ( (cid:103) S rp ) s +1 i = − π ( ε σ s,i + σ a,i ) c (cid:20) c R s +1 i − v s +1 x,i (cid:18) J s +1 i + εK sQ,i (cid:19)(cid:21) + 4 πσ a,i v s +1 x,i c (cid:18) ( T s +1 i ) π − J s +1 i (cid:19) . Moreover, the numerical flux (cid:101) ζ i + is defined by (cid:101) ζ i + = A i + n (cid:16) I s, + i + n> + I s, − i + n< (cid:17) + C i + nJ s +1 i + + C i + π n ( T s +1 i + ) + F i + n ˘ G i + + D i + n (cid:16) δ x J s +1 , + i + n> + δ x J s +1 , − i + n< (cid:17) + D i + π n (cid:16) δ x ( T s +1 , + i + ) n> + δ x ( T s +1 , − i + ) n< (cid:17) , (B.1) SPATIAL-TEMPORAL ASYMPTOTIC PRESERVING SCHEME FOR RADIATION MAGNETOHYDRODYNAMICS 25 where the coefficients in the numerical flux are given by: A = c ∆ tεµ (1 − e − µ ∆ t ) ,C = c σ s ∆ tεµ (cid:18) ∆ t − µ (1 − e − µ ∆ t ) (cid:19) , C = c σ a ∆ tε µ (cid:18) ∆ t − µ (1 − e − µ ∆ t ) (cid:19) ,D = − c σ s ∆ tε µ (cid:18) ∆ t (cid:0) e − µ ∆ t (cid:1) − µ (cid:0) − e − µ ∆ t (cid:1)(cid:19) ,D = − c σ a ∆ tε µ (cid:18) ∆ t (cid:0) e − µ ∆ t (cid:1) − µ (cid:0) − e − µ ∆ t (cid:1)(cid:19) ,F = c ∆ tεµ (cid:18) ∆ t − µ (cid:0) − e − µ ∆ t (cid:1)(cid:19) , (B.2)with µ = cσ a + ε cσ s ε , and ˘ G i + is defined in the same way with (cid:98) G i + . Moreover, the zeroth and first momentsof (B.1) are: 4 π J s +1 i − J si ∆ t + 2 π ∆ x (cid:90) − (cid:101) ζ i + − (cid:101) ζ i − dn = cε ( (cid:102) S re ) s +1 i , (B.3)where (cid:90) − (cid:101) ζ i + dn = (cid:90) − A i + n (cid:16) I s, + i + n> + I s, − i + n< (cid:17) dn + 2 D i + J s +1 i +1 − J s +1 i ∆ x + D i + π ( T s +1 i +1 ) − ( T s +1 i ) ∆ x + F i + σ a,i + ε v s +1 x,i + ( T s +1 i + ) π − J s +1 i + + v s +1 x,i + σ a,i + + ε σ s,i + ε (cid:18) J s +1 i + + 2 εK sQ,i + (cid:19) , πε R s +1 i − R si ∆ t + 2 π ∆ x (cid:90) − n (cid:16)(cid:101) ζ i + − (cid:101) ζ i − (cid:17) dn = cε ( (cid:103) S rp ) s +1 i , (B.4)where (cid:90) − n (cid:101) ζ i + dn = (cid:90) − A i + n (cid:16) I s, + i + n> + I s, − i + n< (cid:17) dn + C i + J s +1 i +1 + J s +1 i ) + C i + π (( T s +1 i +1 ) + ( T s +1 i ) )+ F i + (cid:20) v s +1 x,i + (cid:16) ε σ s,i + + σ a,i + (cid:17) (cid:18) R s +1 i + + (cid:90) − n Q si + dn (cid:19) − ε σ s,i + v s +1 x,i + R s +1 i + − σ a,i + − ε σ s,i + )( v s +1 x,i + ) c (cid:18) J s +1 i + + εK sQ,i + (cid:19) . At last, we will give the diffusion limit of the full discretization of the equilibrium case. Firstly, assuming σ a and σ s are positive, as ε →
0, the leading order of the coefficients in the numerical flux (B.1) are A (∆ t, ε, c, σ a , σ s ) = 0 + O ( ε ) ,εC (∆ t, ε, c, σ a , σ s ) = 0 + O ( ε ) , εC = c + O ( ε ) ,D = 0 + O ( ε ) , D = − cσ a + O ( ε ) , ε F = 1 σ a + O ( ε ) . which means the zero moment and the first moment of the numerical flux, as ε →
0, have the following limit: (cid:90) − (cid:98) ζ i + dn → v s +1 x,i + J s +1 i + − c σ a,i + J s +1 i +1 − J s +1 i ∆ x , (B.5) ε (cid:90) − n (cid:98) ζ i + dn → c J s +1 i +1 + J s +1 i ) . (B.6)In equation (B.3), letting ε → πσ a,i (cid:18) ( T s +1 i ) π − J s +1 i (cid:19) = 0 . Multiplying equation (B.3) by P , adding it to energy equation and substituting limit (B.5) into the equation,then multiplying equation (B.4) by εc P , adding it to momentum equation and substituting limit (B.6) intothe equation, at last letting ε → a ( ρT ) s +1 i − ( ρT ) si ∆ t + ( ρv + B ) s +1 i − ( ρv + B ) si t + P ( T s +1 i ) − ( T si ) ∆ t + 1∆ x ( F s ,i + − F s ,i − )+ 4 P v s +1 x,i + ( T s +1 i + ) − v s +1 x,i − ( T s +1 i − ) ∆ x = P c σ a,i + ∆ x ( T s +1 i +1 ) − ( T s +1 i ) ∆ x − P c σ a,i − ∆ x ( T s +1 i ) − ( T s +1 i − ) ∆ x , (B.7)( ρv x ) s +1 i − ( ρv x ) si ∆ t + 1∆ x ( F s ,i + − F s ,i − ) + P T s +1 i +1 ) − ( T s +1 i − ) ∆ x = 0 , (B.8)which are fully discretization for (2.14c) and (2.14b). References [1] Simon Bolding, Joshua Hansel, Jarrod D. Edwards, Jim E. Morel, and Robert B. Lowrie,
Second-order discretization inspace and time for radiation-hydrodynamics , Journal of Computational Physics (2017), 511 –526.[2] Moysey Brio and Cheng Chin Wu,
An upwind differencing scheme for the equations of ideal magnetohydrodynamics ,Journal of computational physics (1988), no. 2, 400–422.[3] Francis Filbet and Shi Jin, A class of asymptotic-preserving schemes for kinetic equations and related problems with stiffsources , Journal of Computational Physics (2010), no. 20, 7625–7648.[4] Pauline Godillon-Lafitte and Thierry Goudon,
A coupled model for radiative transfer: Doppler effects, equilibrium, andnonequilibrium diffusion asymptotics , Multiscale Modeling & Simulation (2005), no. 4, 1245–1279.[5] Fran¸cois Golse, Shi Jin, and C David Levermore, The convergence of numerical transfer schemes in diffusive regimes i:Discrete-ordinate method , SIAM journal on numerical analysis (1999), no. 5, 1333–1369.[6] Yan-Fei Jiang, James M. Stone, and Shane W. Davis, A godunov method for multidimensional radiation magnetohydrody-namics based on a variable eddington tensor , The Astrophysical Journal Supplement Series (2012), no. 1, 14.[7] Yan-Fei Jiang, James M. Stone, and Shane W. Davis,
An algorithm for radiation magnetohydrodynamics based on solvingthe time-dependent transfer equation , The Astrophysical Journal Supplement Series (2014), no. 1, 7.[8] Shi Jin,
Efficient asymptotic-preserving (ap) schemes for some multiscale kinetic equations , SIAM Journal on ScientificComputing (1999), no. 2, 441–454.[9] Shi Jin, Asymptotic preserving (ap) schemes for multiscale kinetic and hyperbolic equations (2010).[10] Shi Jin and C David Levermore,
Fully-discrete numerical transfer in diffusive regimes , Transport theory and statisticalphysics (1993), no. 6, 739–791.[11] Shi Jin, Lorenzo Pareschi, and Giuseppe Toscani, Diffusive relaxation schemes for multiscale discrete-velocity kineticequations , SIAM Journal on Numerical Analysis (1998), no. 6, 2405–2439.[12] Shi Jin, Min Tang, and Houde Han, A uniformly second order numerical method for the one-dimensional discrete-ordinatetransport equation and its diffusion limit with interface , Networks & Heterogeneous Media (2009), no. 1, 35.[13] Samet Y. Kadioglu, Dana A. Knoll, Robert B. Lowrie, and Rick M. Rauenzahn, A second order self-consistent imex methodfor radiation hydrodynamics , Journal of Computational Physics (2010), no. 22, 8313 –8332.[14] Axel Klar,
An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit , SIAM journal onnumerical analysis (1998), no. 3, 1073–1094.[15] Kerstin K¨upper, Martin Frank, and Shi Jin, An asymptotic preserving two-dimensional staggered grid method for multiscaletransport equations , SIAM Journal on Numerical Analysis (2016), no. 1, 440–461.[16] Edward W Larsen and Jim E Morel, Asymptotic solutions of numerical transport problems in optically thick, diffusiveregimes ii (1989).[17] Edward W Larsen, Jim E Morel, and Warren F Miller Jr,
Asymptotic solutions of numerical transport problems in opticallythick, diffusive regimes , Journal of Computational Physics (1987), no. 2, 283–324. SPATIAL-TEMPORAL ASYMPTOTIC PRESERVING SCHEME FOR RADIATION MAGNETOHYDRODYNAMICS 27 [18] Fengyan Li, Liwei Xu, and Sergey Yakovlev,
Central discontinuous galerkin methods for ideal mhd equations with theexactly divergence-free magnetic field , Journal of Computational Physics (2011), no. 12, 4828–4847.[19] RB Lowrie, JE Morel, and JA Hittinger,
The coupling of radiation and hydrodynamics , The astrophysical journal (1999), no. 1, 432.[20] Robert B Lowrie and Jarrod D Edwards,
Radiative shock solutions with grey nonequilibrium diffusion , Shock Waves (2008), no. 2, 129–143.[21] Ryan G. McClarren, Thomas M. Evans, Robert B. Lowrie, and Jeffery D. Densmore, Semi-implicit time integration forpn thermal radiative transfer , Journal of Computational Physics (2008), no. 16, 7561 –7586.[22] Luc Mieussens,
On the asymptotic preserving property of the unified gas kinetic scheme for the diffusion limit of linearkinetic models , Journal of Computational Physics (2013), 138–156.[23] D Mihalas and BW Mihalas,
Foundations of radiation hydrodynamics oxford university press , New York (1984).[24] J.E. Morel, Todd A. Wareing, and Kenneth Smith,
A linear-discontinuous spatial differencing scheme forsnradiativetransfer calculations , Journal of Computational Physics (1996), no. 2, 445 –462.[25] Gerald C Pomraning,
The equations of radiation hydrodynamics , Courier Corporation, 2005.[26] Philip L Roe,
Characteristic-based schemes for the euler equations , Annual review of fluid mechanics (1986), no. 1,337–365.[27] Michael Sekora and James Stone, A higher-order godunov method for radiation hydrodynamics: Radiation subsystem ,Communications in Applied Mathematics and Computational Science (2009), no. 1, 135–152.[28] Michael D. Sekora and James M. Stone, A hybrid godunov method for radiation hydrodynamics , Journal of ComputationalPhysics (2010), no. 19, 6819 –6852.[29] James M Stone, Thomas A Gardiner, Peter Teuben, John F Hawley, and Jacob B Simon,
Athena: a new code for astro-physical mhd , The Astrophysical Journal Supplement Series (2008), no. 1, 137.[30] Wenjun Sun, Song Jiang, and Kun Xu,
An asymptotic preserving unified gas kinetic scheme for gray radiative transferequations , Journal of Computational Physics (2015), 265–279.[31] Wenjun Sun, Song Jiang, Kun Xu, and Guiyu Cao,
Multiscale simulation for the system of radiation hydrodynamics ,Journal of Scientific Computing (2020), no. 2, 1–24. School of Mathematics, Institute of Natural Sciences and MOE-LSC, Shanghai Jiao Tong University, China
Email address : [email protected] School of Mathematics, Institute of Natural Sciences and MOE-LSC, Shanghai Jiao Tong University, China
Email address : [email protected] School of Mathematics, Institute of Natural Sciences and MOE-LSC, Shanghai Jiao Tong University, China
Email address ::