Created
December 13, 2021 13:35
-
-
Save jwscook/bcdc55b5d17e808caea5face17c3bf5e to your computer and use it in GitHub Desktop.
Trixi Discontinuous Galerkin Maxwell
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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() | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment