From 43900f4c0557a0374174594938b53d8283dfcd97 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 5 Jul 2024 17:00:19 +0200 Subject: [PATCH 1/4] Update docstring `AnalysisSurfaceIntegral` (#2003) * Update docstring `AnalysisSurfaceIntegral` * Update src/callbacks_step/analysis_surface_integral_2d.jl --- src/callbacks_step/analysis_surface_integral_2d.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl index d92ac5a2076..e22d8b14e94 100644 --- a/src/callbacks_step/analysis_surface_integral_2d.jl +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -9,9 +9,9 @@ # surface forces """ - AnalysisSurfaceIntegral{Semidiscretization, Variable}(semi, - boundary_symbol_or_boundary_symbols, - variable) + AnalysisSurfaceIntegral{Variable, NBoundaries}(semi, + boundary_symbols::NTuple{NBoundaries, Symbol}, + variable) This struct is used to compute the surface integral of a quantity of interest `variable` alongside the boundary/boundaries associated with particular name(s) given in `boundary_symbol` From fc2a7618c613e11266ca0f5824f161d03d7ee3ba Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 5 Jul 2024 17:02:42 +0200 Subject: [PATCH 2/4] 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 <10759273+SimonCan@users.noreply.github.com> * Update src/equations/maxwell_1d.jl * Apply suggestions from code review * Apply suggestions from code review Co-authored-by: Hendrik Ranocha * 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 --------- Co-authored-by: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Co-authored-by: Hendrik Ranocha --- .../elixir_maxwell_E_excitation.jl | 59 ++++++++++ .../elixir_maxwell_convergence.jl | 49 +++++++++ src/Trixi.jl | 3 +- src/equations/equations.jl | 4 + src/equations/maxwell_1d.jl | 104 ++++++++++++++++++ test/test_tree_1d.jl | 3 + test/test_tree_1d_maxwell.jl | 44 ++++++++ test/test_unit.jl | 19 ++++ 8 files changed, 284 insertions(+), 1 deletion(-) create mode 100644 examples/tree_1d_dgsem/elixir_maxwell_E_excitation.jl create mode 100644 examples/tree_1d_dgsem/elixir_maxwell_convergence.jl create mode 100644 src/equations/maxwell_1d.jl create mode 100644 test/test_tree_1d_maxwell.jl diff --git a/examples/tree_1d_dgsem/elixir_maxwell_E_excitation.jl b/examples/tree_1d_dgsem/elixir_maxwell_E_excitation.jl new file mode 100644 index 00000000000..b5478331a7d --- /dev/null +++ b/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() diff --git a/examples/tree_1d_dgsem/elixir_maxwell_convergence.jl b/examples/tree_1d_dgsem/elixir_maxwell_convergence.jl new file mode 100644 index 00000000000..e62f7215ace --- /dev/null +++ b/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() diff --git a/src/Trixi.jl b/src/Trixi.jl index b8364eef445..1a509ed92d1 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -162,7 +162,8 @@ export AcousticPerturbationEquations2D, ShallowWaterEquationsQuasi1D, LinearizedEulerEquations1D, LinearizedEulerEquations2D, LinearizedEulerEquations3D, PolytropicEulerEquations2D, - TrafficFlowLWREquations1D + TrafficFlowLWREquations1D, + MaxwellEquations1D export LaplaceDiffusion1D, LaplaceDiffusion2D, LaplaceDiffusion3D, CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D, diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 69a54f7faf8..7875954c5f0 100644 --- a/src/equations/equations.jl +++ b/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 diff --git a/src/equations/maxwell_1d.jl b/src/equations/maxwell_1d.jl new file mode 100644 index 00000000000..838e4e7b812 --- /dev/null +++ b/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 diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index 76129f15e07..98d6ad11c7f 100644 --- a/test/test_tree_1d.jl +++ b/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 diff --git a/test/test_tree_1d_maxwell.jl b/test/test_tree_1d_maxwell.jl new file mode 100644 index 00000000000..0d936f703fe --- /dev/null +++ b/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 diff --git a/test/test_unit.jl b/test/test_unit.jl index d7ec2084361..09cf506576c 100644 --- a/test/test_unit.jl +++ b/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 From f4920c64b50f5e1179aaf979e46684a643c413d3 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 5 Jul 2024 17:03:24 +0200 Subject: [PATCH 3/4] set version to v0.8.2 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index f942de89309..a787ffa8d57 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.8.2-DEV" +version = "0.8.2" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From a160481c4094ef8ff9151f34076a4c93daae3c4e Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 5 Jul 2024 17:03:40 +0200 Subject: [PATCH 4/4] set development version to v0.8.3-DEV --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index a787ffa8d57..c59e01b9e63 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.8.2" +version = "0.8.3-DEV" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"