Maxwell 1D (Wave Equation) (#1949)

* Maxwell 1D

* test

* another maxwell example

* non-trivial example

* fmt

* remove comment

* test for coverage

* Update src/equations/maxwell_1d.jl

Co-authored-by: Simon Candelaresi <[email protected]>

* Update src/equations/maxwell_1d.jl

* Apply suggestions from code review

* Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <[email protected]>

* c is scalar

* try return

* Abstract Maxwell

* similar constructor as lin euler 1d

* try another constructor

* Update src/equations/maxwell_1d.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

---------

Co-authored-by: Simon Candelaresi <[email protected]>
Co-authored-by: Hendrik Ranocha <[email protected]>

examples/tree_1d_dgsem/elixir_maxwell_E_excitation.jl

@@ -0,0 +1,59 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+equations = MaxwellEquations1D()
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = 0.0
+coordinates_max = 1.0
+
+# Create a uniformly refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 30_000) # set maximum capacity of tree data structure
+
+# Excite the electric field which causes a standing wave.
+# The solution is an undamped exchange between electric and magnetic energy.
+function initial_condition_E_excitation(x, t, equations::MaxwellEquations1D)
+    c = equations.speed_of_light
+    E = -c * sin(2 * pi * x[1])
+    B = 0.0
+
+    return SVector(E, B)
+end
+
+initial_condition = initial_condition_E_excitation
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# As the wave speed is equal to the speed of light which is on the order of 3 * 10^8
+# we consider only a small time horizon.
+ode = semidiscretize(semi, (0.0, 1e-7));
+
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval = 100)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl = 1.6)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()

examples/tree_1d_dgsem/elixir_maxwell_convergence.jl

@@ -0,0 +1,49 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+equations = MaxwellEquations1D()
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_hll)
+
+coordinates_min = 0.0
+coordinates_max = 1.0
+
+# Create a uniformly refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 30_000) # set maximum capacity of tree data structure
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,
+                                    solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# As the wave speed is equal to the speed of light which is on the order of 3 * 10^8
+# we consider only a small time horizon.
+ode = semidiscretize(semi, (0.0, 1e-8));
+
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval = 100)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl = 1.6)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()

src/Trixi.jl

@@ -162,7 +162,8 @@ export AcousticPerturbationEquations2D,
        ShallowWaterEquationsQuasi1D,
        LinearizedEulerEquations1D, LinearizedEulerEquations2D, LinearizedEulerEquations3D,
        PolytropicEulerEquations2D,
-       TrafficFlowLWREquations1D
+       TrafficFlowLWREquations1D,
+       MaxwellEquations1D
 
 export LaplaceDiffusion1D, LaplaceDiffusion2D, LaplaceDiffusion3D,
        CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D,

src/equations/equations.jl

@@ -512,4 +512,8 @@ abstract type AbstractEquationsParabolic{NDIMS, NVARS, GradientVariables} <:
 abstract type AbstractTrafficFlowLWREquations{NDIMS, NVARS} <:
               AbstractEquations{NDIMS, NVARS} end
 include("traffic_flow_lwr_1d.jl")
+
+abstract type AbstractMaxwellEquations{NDIMS, NVARS} <:
+              AbstractEquations{NDIMS, NVARS} end
+include("maxwell_1d.jl")
 end # @muladd

src/equations/maxwell_1d.jl

@@ -0,0 +1,104 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+@doc raw"""
+    MaxwellEquations1D(c = 299_792_458.0)
+
+The Maxwell equations of electro dynamics
+```math
+\frac{\partial}{\partial t}
+\begin{pmatrix}
+E \\ B
+\end{pmatrix}
++ 
+\frac{\partial}{\partial x}
+\begin{pmatrix}
+c^2 B \\ E
+\end{pmatrix}
+=
+\begin{pmatrix}
+0 \\ 0 
+\end{pmatrix}
+```
+in one dimension with speed of light `c = 299792458 m/s` (in vacuum).
+In one dimension the Maxwell equations reduce to a wave equation.
+The orthogonal magnetic (e.g.`B_y`) and electric field (`E_z`) propagate as waves 
+through the domain in `x`-direction.
+For reference, see 
+- e.g. p.15 of Numerical Methods for Conservation Laws: From Analysis to Algorithms 
+  https://doi.org/10.1137/1.9781611975109
+
+- or equation (1) in https://inria.hal.science/hal-01720293/document
+"""
+struct MaxwellEquations1D{RealT <: Real} <: AbstractMaxwellEquations{1, 2}
+    speed_of_light::RealT # c
+
+    function MaxwellEquations1D(c::Real = 299_792_458.0)
+        new{typeof(c)}(c)
+    end
+end
+
+function varnames(::typeof(cons2cons), ::MaxwellEquations1D)
+    ("E", "B")
+end
+function varnames(::typeof(cons2prim), ::MaxwellEquations1D)
+    ("E", "B")
+end
+
+"""
+    initial_condition_convergence_test(x, t, equations::MaxwellEquations1D)
+
+A smooth initial condition used for convergence tests.
+"""
+function initial_condition_convergence_test(x, t, equations::MaxwellEquations1D)
+    c = equations.speed_of_light
+    char_pos = c * t + x[1]
+
+    sin_char_pos = sin(2 * pi * char_pos)
+
+    E = -c * sin_char_pos
+    B = sin_char_pos
+
+    return SVector(E, B)
+end
+
+# Calculate 1D flux for a single point
+@inline function flux(u, orientation::Integer,
+                      equations::MaxwellEquations1D)
+    E, B = u
+    return SVector(equations.speed_of_light^2 * B, E)
+end
+
+# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
+@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Int,
+                                     equations::MaxwellEquations1D)
+    λ_max = equations.speed_of_light
+end
+
+@inline have_constant_speed(::MaxwellEquations1D) = True()
+
+@inline function max_abs_speeds(equations::MaxwellEquations1D)
+    return equations.speed_of_light
+end
+
+@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer,
+                                     equations::MaxwellEquations1D)
+    min_max_speed_davis(u_ll, u_rr, orientation, equations)
+end
+
+@inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer,
+                                     equations::MaxwellEquations1D)
+    λ_min = -equations.speed_of_light
+    λ_max = equations.speed_of_light
+
+    return λ_min, λ_max
+end
+
+# Convert conservative variables to primitive
+@inline cons2prim(u, ::MaxwellEquations1D) = u
+@inline cons2entropy(u, ::MaxwellEquations1D) = u
+end # @muladd

test/test_tree_1d.jl

@@ -51,6 +51,9 @@ isdir(outdir) && rm(outdir, recursive = true)
 
     # Linearized Euler
     include("test_tree_1d_linearizedeuler.jl")
+
+    # Maxwell
+    include("test_tree_1d_maxwell.jl")
 end
 
 # Coverage test for all initial conditions

test/test_tree_1d_maxwell.jl

@@ -0,0 +1,44 @@
+module TestExamples1DMaxwell
+
+using Test
+using Trixi
+
+include("test_trixi.jl")
+
+EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
+
+@testset "Maxwell" begin
+#! format: noindent
+
+@trixi_testset "elixir_maxwell_convergence.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_maxwell_convergence.jl"),
+                        l2=[8933.196486422636, 2.979793603210305e-5],
+                        linf=[21136.527033627033, 7.050386515528029e-5])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
+
+@trixi_testset "elixir_maxwell_E_excitation.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_maxwell_E_excitation.jl"),
+                        l2=[1.8181768208894413e6, 0.09221738723979069],
+                        linf=[2.5804473693440557e6, 0.1304024464192847])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
+end
+
+end # module

test/test_unit.jl

@@ -1598,6 +1598,25 @@ end
                                     min_max_speed_davis(u_ll, u_rr, normal_z,
                                                         equations)))
     end
+
+    @timed_testset "Maxwell 1D" begin
+        equations = MaxwellEquations1D()
+
+        u_values_left = [SVector(1.0, 0.0),
+            SVector(0.0, 1.0),
+            SVector(0.5, -0.5),
+            SVector(-1.2, 0.3)]
+
+        u_values_right = [SVector(1.0, 0.0),
+            SVector(0.0, 1.0),
+            SVector(0.5, -0.5),
+            SVector(-1.2, 0.3)]
+        for u_ll in u_values_left, u_rr in u_values_right
+            @test all(isapprox(x, y)
+                      for (x, y) in zip(min_max_speed_naive(u_ll, u_rr, 1, equations),
+                                        min_max_speed_davis(u_ll, u_rr, 1, equations)))
+        end
+    end
 end
 
 @testset "SimpleKronecker" begin