Discontinuous Galerkin Methods for anisotropic diffusion Praveen Chandrashekar praveen@tifrbng.res.in Center for Applicable Mathematics Tata Institute of Fundamental Research Bangalore-560065, India http://cpraveen.github.io Oberseminar, Dept. of Mathematics, Univ. of W¨ urzburg 21 July, 2016 Supported by Airbus Foundation Chair at TIFR-CAM, Bangalore http://math.tifrbng.res.in/airbus-chair 1 / 57
2.1. Equations of MHD with Anisotropic Heat Conduction The equations of MHD with the addition of the heat flux, Q , and a vertical gravitational acceleration, g = − g 0 ˆ z are ∂ρ ∂t + ∇ · ( ρ v ) = 0 , (1) � � � � p + B 2 ∂ ( ρ v ) I − BB + ∇ · ρ vv + + ρ g = 0 , (2) ∂t 8 π 4 π ∂ B ∂t + ∇ × ( v × B ) = 0 , (3) � � � � E + p + B 2 ∂E − B ( B · v ) ∂t + ∇ · v + ∇ · Q + ρ g · v = 0 , (4) 8 π 4 π where the symbols have their usual meaning. The total energy E is given as E = � + ρ v · v + B · B 8 π , (5) 2 with the internal energy, � = p/ ( γ − 1). For this paper, we assume γ = 5 / 3 throughout. Parrish & Stone, http://arxiv.org/abs/astro-ph/0507212 2 / 57
The heat flux contains contributions both from electron motions (which are constrained to move primarily along field lines) and from isotropic transport which may arise due to photons or particle collisions which drive cross field di ff usion. Thus, Q = Q C + Q R , where Q C = − χ C ˆ bˆ b · ∇ T, (6) Q R = − χ R ∇ T, (7) 1962), ˆ where χ C is the Spitzer Coulombic conductivity (Spitzer b is a unit vector in the direction of the magnetic field, and χ R is the coe ffi cient of isotropic conductivity, ostensibly due to radiation. We consider both χ C and χ R as free parameters in the problem and will vary both independently. Parrish & Stone, http://arxiv.org/abs/astro-ph/0507212 3 / 57
2 THERMAL CONDUCTION ionization state and the composition of the gas. Since we We start with the energy conservation equation which can are dealing mainly with a high temperature, low metallic- be written as ity plasma of primordial composition we assume that the ρ@u gas is fully ionised and consequently set the mean molecular @t + r · j = 0 , (4) weight to the constant value µ ∼ 0 . 6 m p . The final equation where u is the gas internal energy per unit mass, j is the di- for anisotropic thermal conduction then becomes rectional heat flux and ρ is the gas density. If the conduction @u 1 is isotropic, then the heat flux j is opposite to the direction @t = c v ρ r · [ κ b ( b · r ) u ] , (13) of the temperature gradient and below we focus on studying numerical solutions schemes j = − κ sp ( T ) r T. (5) for this partial differential equation. In the presence of the magnetic fields, electrons prefer- entially move along the magnetic field lines. Therefore the heat flux is modified as j = − κ [ b ( b · r T )] , (10) where b = B / | B | . The coefficient of conduction perpendic- ular to the magnetic fields is set to zero (see Arth et al. 2014 for a discussion on non vanishing perpendicular conduction coefficients). We can also rewrite Equation (10) as j = − κ c v [ b ( b · r u )] , (11) Kannan et al., Accurately simulating by replacing temperature for an ideal gas with the internal anisotropic thermal conduction on a energy per unit mass k B T moving mesh (2015) u = ( γ − 1) µ = c v T, (12) http://arxiv.org/abs/1512.03053 where γ is the adiabatic index, and µ is the mean molec- ular weight. The value of µ depends on the temperature, 4 / 57
Isotropic diffusion model ∂t = χ∂ 2 θ ∂θ ∂x 2 Usual finite difference j + χ ∆ t θ n +1 = θ n h 2 ( θ n j − 1 − 2 θ n j + θ n j +1 ) j Stable in L 2 and L ∞ norms if σ = χ ∆ t ≤ 1 h 2 2 5 / 57
Anisotropic diffusion model ∂θ ∂t + ∇ · q = 0 , in Ω (1) where the heat flux vector is given by [7], [6] q ( θ ) = − ( χ � − χ ⊥ ) b ( b · ∇ θ ) − χ ⊥ ∇ θ The vector field b has unit magnitude and is assumed to be given, and χ ⊥ ≥ 0 , χ � − χ ⊥ > 0 are taken to be constant In many applications χ � ≫ χ ⊥ and the diffusion is mainly along the vector field b which is refered to as anisotropic diffusion, while the term containing χ ⊥ leads to isotropic diffusion. For simplicity of notation, let us define κ = χ � − χ ⊥ > 0 , η = χ ⊥ ≥ 0 6 / 57
Anisotropic diffusion model The heat flux vector can be written as D = ηI + κ bb ⊤ q = − D ∇ θ where The diffusion matrix is positive definite and satisfies η | p | 2 ≤ ( D p , p ) ≤ ( η + κ ) | p | 2 ∀ p ∈ R 2 If η = 0 , then the diffusion matrix is only positive semi-definite. Boundary conditions θ = Θ on Γ d (2) q · n = g on Γ n (3) Total energy � � � � � d 1 2 θ 2 d x + κ | b ·∇ θ | 2 d x + η |∇ θ | 2 d x + ( q ( θ ) · n )Θ d s + gθ d s = 0 d t Ω Ω Ω Γ d Γ n If Θ = 0 , g = 0 , we observe that the energy decreases with time. 7 / 57
Discontinuous solutions: χ � = 1 , χ ⊥ = 0 Consider b = (1 , 0) and initial condition � θ u y < 0 θ ( x, y, 0) = θ 0 ( x, y ) = θ l y > 0 On the surface y = 0 , n = (0 , 1) , heat flux q · n = − ( b · n )( b · ∇ θ ) = 0 since b · n = 0 . Hence θ ( x, y, t ) = θ 0 ( x, y ) 8 / 57
Ring test: Diffusion on a ring ( χ � = 1 , χ ⊥ = 0 ) The initial temperature is given by � 11 12 π < φ < 13 12 , if 0 . 5 < r < 0 . 7 and 12 π θ ( r, φ, 0) = 10 , otherwise The “magnetic field” is b = ( − sin φ, cos φ ) . 9 / 57
Ring test: Galerkin method 7046 cells, 3624 vertices CG(1) using BDF2 and σ = 1 at t = 200 10 / 57
Finite difference/volume methods ( η = 0 ) Heat flux q = ( q x , q y ) q x = − κB x ( B x ∂ x θ + B y ∂ y θ ) , q y = − κB y ( B x ∂ x θ + B y ∂ y θ ) Finite volume method θ n +1 ( q x ) i + 1 2 ,j − ( q x ) i − 1 ( q y ) i,j + 1 2 − ( q y ) i,j − 1 − θ n 2 ,j i,j i,j + + 2 = 0 ∆ t ∆ x ∆ y Approximation of fluxes: centered asymmetric scheme [6] � ∂θ � θ i +1 ,j − θ i,j ( q x ) i + 1 2 ,j = − κB x B x + B y ∆ x ∂y i + 1 2 ,j � �� � ? � ∂θ � 1 = 4∆ y ( θ i,j +1 + θ i +1 ,j +1 − θ i,j − 1 − θ i +1 ,j − 1 ) ∂y i + 1 2 ,j 11 / 57
T =0.1 T =10 i,j+1 i+1,j+1 T =0.1 T =0.1 i,j i+1,j 3 2.5 temperature at (i,j) 2 1.5 1 0.5 0 − 0.5 0 2 4 6 8 10 time Fig. 4. Test problem to show that the asymmetric method can result in negative temperature. Magnetic field lines are along the diagonal p ffiffiffi with b x ¼ � b y ¼ 1 = 2 . With the asymmetric method heat flows out of the third quadrant which is already a temperature minimum, resulting in a negative temperature T i , j . However due to numerical perpendicular di ff usion, at late times the temperature becomes positive again. The temperature at ( i , j ) is shown for di ff erent methods: asymmetric (solid line), symmetric (dotted line), asymmetric and symmetric with slope limiters (dashed line; both give the same result), and symmetric with entropy limiting (dot dashed line). Sharma & Hammett [6] 12 / 57
Mesh and solution space T h = triangulation of the domain by disjoint triangles, Ω = ∪ K ∈T h K . E i = interior edges, E d = Dirichlet edges, E n = Neumann edges h e = length of edge e h K = diameter of K ∈ T h = length of largest side of K Diameter of the mesh h = max K ∈T h h K ρ K = diameter of the largest circle that can be inscribed in triangle K The mesh is shape regular , i.e., there is a constant α > 0 such that h K ≤ α ρ K h min = length of the smallest edge in the triangulation 13 / 57
Mesh and solution space Space of broken polynomials V k h = { ϕ ∈ L 2 (Ω) : ϕ | K ∈ P k ( K ) , ∀ K ∈ T h } and their vector version W k h = V k h × V k h For each interior edge 1 2( q − + q + ) { { q } } = ϕ + − ϕ − � ϕ � = n − ϕ − n − + ϕ + n + � ϕ � n = ϕ − ϕ + On a boundary edge n + { { q } } = q � ϕ � = ϕ � ϕ � n = ϕ n 14 / 57
Inverse inequality 1 h e � ϕ � e ≤ C t � ϕ � K , 2 ∀ ϕ ∈ P k ( K ) � where C t = ( k + 1)( k + 2) . Inverse Poincare-type inequality h K �∇ ϕ � K ≤ C p � ϕ � K , ∀ ϕ ∈ P k ( K ) √ with C p = 2 α , where α is the shape regularity parameter. 15 / 57
Interior Penalty DG (IPDG) n − Multiply the heat equation (1) by a test function ϕ h ∈ V k h and integrate ϕ − on one cell K K � � � h n − = 0 ˆ } · ϕ − ϕ h ∂ t θ h d x + D ∇ θ h · ∇ ϕ h d x − D { {∇ θ h } K K ∂K where ˆ } ⊤ D = ηI + κ { { b } }{ { b } 16 / 57
Interior Penalty DG (IPDG) Sum over all cells K ∈ T h � � � � ˆ ϕ h ∂ t θ h d x + D ∇ θ h · ∇ ϕ h d x − D { {∇ θ h } } · � ϕ h � n d s Ω Ω e e ∈E i � � � � − ( D ∇ θ h · n ) ϕ h d s + gϕ h d s = 0 e e e ∈E d e ∈E n However such a scheme is unstable !!! IPDG scheme : find θ h ∈ V k h such that for all ϕ h ∈ V k h ( ∂ t θ h , ϕ h ) + A ( κ,η ) ( θ h , ϕ h ) = ℓ Θ ( ϕ h ) + ℓ g ( ϕ h ) (4) h 17 / 57
Recommend
More recommend
