Structure-preserving numerical methods in relativity Douglas N. - PowerPoint PPT Presentation

Structure-preserving numerical methods in relativity Douglas N. Arnold, University of Minnesota Advances and Challenges in Computational Relativity September 14, 2020 • ICERM

Structure-preserving discretization Structure-preserving discretizations are numerical methods which preserve, on the discrete level, key geometric, topological, and algebraic structures possessed by the original continuous system. Classical examples are symplectic integrators for Hamiltonian ODEs which preserve the symplectic 2-form associated to the Hamiltonian. Here the Kepler problem is integrated over 4 periods using 200,000 timesteps. Euler’s method symplectic Euler x n + 1 − x n x n + 1 − x n = v n = v n h h v n + 1 − v n v n + 1 − v n = − x n + 1 = − x n | x n + 1 | 3 h | x n | 3 h 1 / 32

Higher order methods 1,000 periods, 1,000,000 time steps Classical 4-stage Runge–Kutta 4-stage Gauss–Legendre (symplectic) 2 / 32

Higher order methods 1,000 periods, 1,000,000 time steps Classical 4-stage Runge–Kutta 4-stage Gauss–Legendre (symplectic) Ruth 1983, Feng Kang 1985, Sanz-Serna 1990, Leimkuhler-Reich 2004 2 / 32

Symplectic integration of the solar system In a famous 2009 paper in Nature , Laskar and Gastineau used a symplectic integrator to simulate the evolution of the solar system for the next 5 Gyr! In fact, they did it 2,500 times, varying the initial position of Mercury by 0.38 mm each time. 1% of the simulations resulted in unstable or collisional orbits. 3 / 32

Structure-preserving discretization for PDEs: Finite Element Exterior Calculus

Hilbert complexes FEEC applies to PDEs which are related to a Hilbert complex . closed, densely-defined domain of d operator Hilbert space · · · → W k − 1 d k − 1 , V k − 1 → W k d k , V k → W k + 1 → · · · − − − − − − − d ◦ d = 0 The sequence must be a complex: Further we require that the Hilbert complex is closed , meaning that the range of each d is closed. (This can be difficult to prove.) Example: grad, H 1 ( Ω ) curl, H ( curl ) · · · → L 2 ( Ω ) → L 2 ( Ω , R 3 ) → L 2 ( Ω , R 3 ) → · · · − − − − − − − − − − − − 4 / 32

Structure of Hilbert complexes · · · → W k − 1 d k − 1 , V k − 1 → W k d k , V k → W k + 1 → · · · − − − − − − − − Each closed Hilbert complex has: cohomology and harmonic forms: N ( d k ) / R ( d k − 1 ) ≈ H k d ∗ k , V ∗ d ∗ k + 1 , V ∗ − W k + 1 ← · · · · · · ← W k − 1 − W k k + 1 k ← − − ← − − − − dual complex: N ( d k ) ⊥ = R ( d ∗ k + 1 ) duality: W k = R ( d ) ⊕ H ⊕ R ( d ∗ ) Hodge decomposition: � u � ≤ c � du � ∀ u ⊥ N ( d ) Poincar´ e inequality: and more . . . 5 / 32

The Hodge Laplacian d d W k − 1 W k W k + 1 d ∗ d ∗ For every Hilbert complex and for each k , there is associated the Hodge Laplacian operator dd ∗ + d ∗ d : W k → W k The Hodge-Laplace equation ( dd ∗ + d ∗ d ) u = f is always well-posed up to harmonic forms = nullspace. It has the obvious weak formulation � d ∗ u , d ∗ v � + � du , dv � = � f , v � ∀ v ∈ V ∩ V ∗ . The mixed weak formulation turns out to be better for discretization: ( σ , u ) ∈ V k − 1 × V k Find s.t. τ ∈ V k − 1 , � σ , τ � − � u , d τ � = 0, v ∈ V k . � d σ , v � + � du , dv � = � f , v � , 6 / 32

The de Rham complex on a domain in 3D grad → L 2 ( Ω , R 3 ) curl → L 2 ( Ω , R 3 ) div 0 → L 2 ( Ω ) → L 2 ( Ω ) → 0 − − − − − This is a special case of the de Rham complex on any manifold in n -D: 0 → Λ 0 ( Ω ) d → Λ 1 ( Ω ) d → · · · d → Λ n ( Ω ) → 0 − − − It is a closed Hilbert complex for any bounded domain Ω ⊂ R n with Lipschitz boundary. 7 / 32

The de Rham complex on a domain in 3D grad → L 2 ( Ω , R 3 ) curl → L 2 ( Ω , R 3 ) div 0 → L 2 ( Ω ) → L 2 ( Ω ) → 0 − − − − − This is a special case of the de Rham complex on any manifold in n -D: 0 → Λ 0 ( Ω ) d → Λ 1 ( Ω ) d → · · · d → Λ n ( Ω ) → 0 − − − It is a closed Hilbert complex for any bounded domain Ω ⊂ R n with Lipschitz boundary. In 3D: The 0-form Hodge Laplacian = − ∆ (standard Laplacian) The 1-form Hodge Laplacian = curl curl − grad div (vector Laplacian) The 2-form Hodge Laplacian = curl curl − grad div with different boundary conditions and different weak formulation. The 3-form Hodge Laplacian = − ∆ with the mixed weak formulation. 7 / 32

The Hodge wave equation U + ( dd ∗ + d ∗ d ) U = 0, ¨ ˙ U ( 0 ) = U 0 , U ( 0 ) = U 1 Then σ : = d ∗ U , ρ : = dU , u : = ˙ U satisfy       − d ∗ σ ˙ 0 0 σ strong  + d ∗  = 0 u ˙ d 0 u     − d ρ ˙ 0 0 ρ ( σ , u , ρ ) : [ 0, T ] → V 0 × V 1 × W 2 Find s.t. τ ∈ V 0 , � ˙ σ , τ � − � u , d τ � = 0, weak v ∈ V 1 , � ˙ u , v � + � d σ , v � + � ρ , dv � = 0, η ∈ W 2 . � ˙ ρ , η � − � du , η � = 0, T HEOREM For any initial data ( σ 0 , u 0 , ρ 0 ) ∈ V 0 × V 1 × W 2 , there exists a unique solution ( σ , u , ρ ) ∈ C 0 ([ 0, T ] , V 0 × V 1 × W 2 ) ∩ C 1 ([ 0, T ] , W 0 × W 1 × W 2 ) . 8 / 32

Example: Maxwell’s equations ˙ ˙ D = curl H B = − curl E div D = 0 div B = 0 D = ǫ E B = µ H W 0 = L 2 ( Ω ) Find ( σ , E , B ) : [ 0, T ] → H 1 × H ( curl ) × L 2 s.t. W 1 = L 2 ( Ω , R 3 , ǫ dx ) � ˙ σ , τ �−� ǫ E , grad τ � = 0 ∀ τ , W 2 = L 2 ( Ω , R 3 , µ − 1 dx ) � ǫ ˙ E , F � + � ǫ grad σ , F � − � µ − 1 B , curl E � = 0 ∀ F , W 0 grad, H 1 � µ − 1 ˙ B , C � + � µ − 1 curl E , C � = 0 ∀ C . → W 1 − curl, H ( curl ) → W 2 − − − − − − − − − − − − T HEOREM For any given initial data, there is a unique solution σ , E, and B to this Hodge wave equation. If the initial data satisfies the constraints ( σ , div D, and div B vanish at t = 0 ), and we set D = ǫ E, H = µ − 1 B, then σ vanishes identically, and E, B, D, H satisfy Maxwell’s equations. 9 / 32

Discretization of the Hodge Laplacian and Hodge wave eq Discretize the mixed weak formulation by Galerkin’s method using finite dimensional space V k h ⊂ V k . The discretization is stable and convergent under two crucial conditions: h ⊂ V k + 1 dV k Subcomplex property: h Bounded cochain projections: d d → V k + 1 → · · · · · · → V k − 1 → V k − − − − − − − −     � π k − 1   � π k � π k + 1 h h h d d · · · → V k − 1 → V k + 1 → V k − − − − − − − − → · · · h h h When these hold, the discrete spaces themselves form a Hilbert complex, and the discretization of the Hodge Laplacian is exactly the Hodge Laplacian associated to this discrete Hilbert complex. This leads to stability and convergence. Structure-preserving discretization! 10 / 32

Finite element Galerkin subspaces A key challenge is to construct finite element spaces to use as the V k h . This means specifying: element shape (simplex, cube, . . . ) space of shape functions on each element (scalar or vector polynomials up to a certain degree, perhaps trimmed or enriched . . . ) P 2 Λ 1 ( R 2 ) degrees of freedom (associated to faces, unisolvent) and ensuring the subcomplex and bounded cochain projection properties. Remarks: 1. The degrees of freedom are crucial for cochain projections. 2. The construction of such structure-preserving finite element spaces is specific to each particular complex. 3. For the de Rham complex the FE spaces to use for each order of differential forms are well understood. 11 / 32

http://z.umn.edu/femtable

Motivating examples for FEEC

Eigenvalues of 1 -form Laplacian ( d ∗ d + dd ∗ ) u = ( curl curl − grad div ) u = λ u , u · n = 0, curl u = 0 on bdry 13 / 32

Eigenvalues of 1 -form Laplacian ( d ∗ d + dd ∗ ) u = ( curl curl − grad div ) u = λ u , u · n = 0, curl u = 0 on bdry Find nonzero u ∈ Λ 1 such that ( du , dv ) + ( d ∗ u , d ∗ v ) = λ ( u , v ) ∀ v ∈ Λ 1 std P 1 FEM λ 1 = 1.94 λ 2 = 2.02 λ 3 = 2.26 13 / 32

Eigenvalues of 1 -form Laplacian ( d ∗ d + dd ∗ ) u = ( curl curl − grad div ) u = λ u , u · n = 0, curl u = 0 on bdry Find nonzero u ∈ Λ 1 such that ( du , dv ) + ( d ∗ u , d ∗ v ) = λ ( u , v ) ∀ v ∈ Λ 1 std P 1 FEM λ 1 = 1.94 λ 2 = 2.02 λ 3 = 2.26 FEEC λ 1 = 0 λ 1 = 0.617 λ 2 = 0.658 13 / 32

Convergence plots for the 1 -form Laplacian, P 1 vs FEEC P 1 finite elements 14 / 32

Convergence plots for the 1 -form Laplacian, P 1 vs FEEC P 1 finite elements true eigenvalues 14 / 32

Convergence plots for the 1 -form Laplacian, P 1 vs FEEC P 1 finite elements FEEC true eigenvalues 14 / 32

Convergence plots for the 1 -form Laplacian, P 1 vs FEEC P 1 finite elements FEEC true eigenvalues The standard P 1 FEM is not convergent for this problem. Nightmare! Fortunately, FEEC saves the day! 14 / 32

The Maxwell eigenvalue problem with std P 1 FE Find nonzero u such that ( du , dv ) = λ ( u , v ) ∀ v (here d = curl) Brezzi–Boffi–Gastaldi ’99 15 / 32