Trixi Discontinuous Galerkin Maxwell

elixir_maxwell.jl

using OrdinaryDiffEq
using LinearAlgebra
using StaticArrays
using Trixi
import Trixi: AbstractEquations, varnames, flux

abstract type AbstractMaxwellEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
struct MaxwellEquations2D{RealT<:Real} <: AbstractMaxwellEquations{2, 6}
  μ₀::Float64
  ϵ₀::Float64
  α::RealT
  function MaxwellEquations2D(μ₀=1.0, ϵ₀=1.0, α::T=1.0) where T<:Real
    0.0 <= α <= 1.0 || throw(ArgumentError("Upwindedness parameter, α, must be ∈ [0,1]."))
    return new{T}(μ₀, ϵ₀, α)
  end
end

varnames(_, MaxwellEquations2D) = ("Ex","Ey","Ez","Bx","By","Bz")#,"Jx","Jy","Jz","rho","μ", "ϵ")
electricfield(u, eq::MaxwellEquations2D) = (@view u[1:3])
magneticfield(u, eq::MaxwellEquations2D) = (@view u[4:6])
#currentfield(u, eq::MaxwellEquations2D) = (@view u[7:9])
#chargefield(u, eq::MaxwellEquations2D) = u[10]
#permeability(u, eq::MaxwellEquations2D) = u[11]
#permitivity(u, eq::MaxwellEquations2D) = u[12]

# Calculate 1D flux for a single point
@inline function flux(u, normal::AbstractVector, eq::MaxwellEquations2D)
  # Calculate flux for conservative state variables
  FE = -cross(normal, magneticfield(u, eq)) / (eq.μ₀ * eq.ϵ₀) # * permeability(u, eq)
  FH =  cross(normal, electricfield(u, eq))
  return @SVector [FE[1], FE[2], FE[3], FH[1], FH[2], FH[3]]#, 0, 0, 0, 0, 0, 0]
end

flux(u, i::Int, eq::MaxwellEquations2D) = flux(u, (@SVector [j == i for j ∈ 1:3]), eq)

# Calculate 1D flux for a single point
@inline function flux_upwind(u⁻, u⁺, i::Int, eq::MaxwellEquations2D)
  n̂ = @SVector [j == i for j ∈ 1:3] 
  return flux_upwind(u⁻, u⁺, n̂, eq)
end
@inline function flux_upwind(u⁻::T, u⁺::T, n::AbstractArray, eq::MaxwellEquations2D)::T where {T}
  H⁻ = magneticfield(u⁻, eq) / eq.μ₀ # * permeability(u⁻, eq)
  H⁺ = magneticfield(u⁺, eq) / eq.μ₀ # * permeability(u⁺, eq)
  E⁻ = electricfield(u⁻, eq)
  E⁺ = electricfield(u⁺, eq)
  Z⁻ = 1 #sqrt(permeability(u⁻, eq) / permitivity(u⁻, eq))
  Z⁺ = 1 #sqrt(permeability(u⁺, eq) / permitivity(u⁺, eq))
  Z̄ = Z⁺ + Z⁻
  FE = cross(n,  (H⁺ - H⁻) * Z⁺ - cross(n, E⁺ - E⁻)) / Z̄
  FH = cross(n, -(E⁺ - E⁻) / Z⁺ - cross(n, H⁺ - H⁻)) * Z̄
  return T([FE[1], FE[2], FE[3], FH[1], FH[2], FH[3]])#, 0, 0, 0, 0, 0, 0])
end


function initial_condition(x, t, equations::MaxwellEquations2D)
  kx = ky = pi
  Ex =  sin(kx*x[1]) * cos(ky*x[2]) * kx * 0
  Ey = -cos(kx*x[1]) * sin(ky*x[2]) * ky * 0
  Ez = Bx = By = 0.0
  Bz = sin(kx * x[1]) * sin(ky * x[2])
  return @SVector [Ex, Ey, Ez, Bx, By, Bz]
end

function run()
  # const volume_flux  = (flux_central, flux) # ? how use
  # const surface_flux = (flux_lax_friedrichs, flux) # ? how use
  
  solver = DGSEM(; RealT=Float64, polydeg=3,
    surface_flux=flux_upwind, # ? is this right
    surface_integral=SurfaceIntegralWeakForm(flux_upwind),# ? is this right
    volume_integral=VolumeIntegralWeakForm())# ? is this right
  
  eqns = MaxwellEquations2D()
  ncells1D = 128
  mesh = TreeMesh((-1.0, -1.0), (1.0, 1.0), # min/max coordinates
                  initial_refinement_level=4,
                  n_cells_max=ncells1D^2)
  
  semi = SemidiscretizationHyperbolic(mesh, eqns, initial_condition, solver)
  
  tspan = (0.0, 2.0)
  ode = semidiscretize(semi, tspan)
  summary_callback = SummaryCallback()
  analysis_interval = 100
  
  save_solution = SaveSolutionCallback(
    interval=10000,
    save_initial_solution=true,
    save_final_solution=true,
    solution_variables=(u, eq) -> u)
  
  callbacks = CallbackSet(summary_callback, save_solution)
  
  sol = solve(ode,
              RK4(),#CarpenterKennedy2N54(williamson_condition=false),
              dt=1/ncells1D / 10, #  how small?
              save_everystep=false,
              callback=callbacks);
  summary_callback() # print the timer summary
end
run()