diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 7aadc3e320fd68d7e119b8602468a7565a50e6d3..04b8e2a1b89d80c455bb74da4809d2d7257a1ca7 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -22,15 +22,16 @@ unittests:
integrationtests:
script:
- julia --project=test/IntegrationTests -e "import Pkg; Pkg.develop(path=\".\"); Pkg.instantiate()"
- - julia --project=test/IntegrationTests -e "using IntegrationTests; runtests()"
- - mv test/IntegrationTests/tests/tests.log public
+ # using !runtest() because Int(false) = 0 and 0 exit code indicates success in shells
+ - julia --project=test/IntegrationTests -e "using IntegrationTests; exit(!runtests())"
only:
- main
- develop
- merge_requests
artifacts:
+ when: always
paths:
- - public
+ - test/IntegrationTests/tests/tests.log
pages:
script:
diff --git a/Project.toml b/Project.toml
index 04327461d2a0a2fec0911e52cae37de7e78714db..d61af21ff70fecab3ecd7d00e5db073677433a8c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,9 +5,7 @@ version = "1.0.0"
[deps]
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
-Configurations = "5218b696-f38b-4ac9-8b61-a12ec717816d"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
@@ -29,6 +27,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+ToggleableAsserts = "07ecc579-1b30-493c-b961-3180daf6e3ae"
[compat]
julia = "1.7"
diff --git a/experiments/grhd_tov/grhd_tov_1.toml b/experiments/grhd_tov/grhd_tov_1.toml
index 13ff95d455cc293e52f225eb5adefbe63562feb4..9e291be252f3b7bda865296e584ba1008184726c 100644
--- a/experiments/grhd_tov/grhd_tov_1.toml
+++ b/experiments/grhd_tov/grhd_tov_1.toml
@@ -1,19 +1,31 @@
[EquationOfState]
-type = "HybridEoS"
+eos = "hybrid"
-[BoundaryConditions]
-type = "from initial data"
+[GRHD]
+bc = "from_initial_data"
+id = "tov"
+id_filename = "TOV.h5"
+
+[TCI]
+method = "mda"
+
+[HRSC]
+# # type = "BernsteinReconstruction"
+method = "av"
+av_method = "mda"
+av_smoother_order = 1
+
+# Non-project specific parameters
[Mesh]
xrange = [ -20.0, 20.0 ]
-CFL = 0.1
-N = 2
-K = 500
+k = 50
+cfl = 0.1
[Plot]
every_dt = 0
-every_iteration = 0
-pause = false
+every_iteration = 1
+pause = true
[Output]
every_iteration = 1
@@ -21,19 +33,8 @@ every_dt = 0
variables = [ "sD", "sSr", "stau", "D", "Sr", "tau", "rho", "eps", "vr", "p" ]
[Grid]
-type = "LGL_lumped"
+basis = "lgl_lumped"
+n = 8
[Evolution]
-tend = 1.2
-
-[Initialdata]
-type = "tov"
-
-[TCI]
-type = "ModalDecayAverage"
-
-# [HRSC]
-# # type = "BernsteinReconstruction"
-# type = "MDAViscosity"
-# smoother = true
-# order = 1
+tend = 10.2
diff --git a/scripts/movie.jl b/scripts/movie.jl
index dc11671e0d4f2efcfbc5fde9a8e6a5aa3854cacf..d1fb28fde07f1c0a34b6be946f617caefaa0d9e8 100644
--- a/scripts/movie.jl
+++ b/scripts/movie.jl
@@ -50,9 +50,26 @@ function make_movie(dir;
dir = abspath(joinpath(".",dir))
tomlfile = get_parameter_file(dir)
- parameters = FlatDict(TOML.parsefile(joinpath(dir, tomlfile)))
+ pars = TOML.parsefile(joinpath(dir, tomlfile))
+ mesh = if haskey(pars, "Equation")
+ parameters = FlatDict(TOML.parsefile(joinpath(dir, tomlfile)))
+ display(parameters)
+ mesh = dg1d.Mesh(N=parameters["Mesh.N"],
+ K=parameters["Mesh.K"],
+ CFL=parameters["Mesh.CFL"],
+ xrange=parameters["Mesh.xrange"],
+ basis=lowercase(parameters["Grid.type"]))
+ else
+ project, parameters = dg1d.load_parameters(joinpath(dir, tomlfile))
+ mesh_prms = parameters["Mesh"]
+ mesh = dg1d.Mesh(N=mesh_prms["n"],
+ K=mesh_prms["k"],
+ CFL=mesh_prms["cfl"],
+ xrange=mesh_prms["xrange"],
+ basis=lowercase(mesh_prms["basis"]))
+ end
- mesh = make_Mesh(trim_headers(parameters, "Mesh"))
+ # mesh = make_Mesh(trim_headers(parameters, "Mesh"))
# read data
output = Output(joinpath(dir, "output.h5"))
@@ -115,7 +132,7 @@ function make_movie(dir;
prevbt = headergrid[2,1] = Button(fig, label=labels[:prev], font=font, fontsize=35)
idxbox = headergrid[2,2] = Textbox(fig, stored_string="1", validator=Int64,
displayed_string="1", width=Fixed(50),
- halign=:center, reset_on_defocuse=true,
+ halign=:center, #reset_on_defocuse=true,
defocus_on_submit=false)
nextbt = headergrid[2,3] = Button(fig, label=labels[:next], font=font, fontsize=35)
playbwrdbt = headergrid[2,4] = Button(fig, label=labels[:playbwrd], font=font, fontsize=35)
diff --git a/src/bdryconditions.jl b/src/bdryconditions.jl
index 41b5d293f05939bfdc95b96afc574fc0adc552dd..78c1dfca92f52f43efc856d12f6cef5b15b1ceda 100644
--- a/src/bdryconditions.jl
+++ b/src/bdryconditions.jl
@@ -1,5 +1,4 @@
-export AbstractBdryConditions,
- DirichletBC, OutflowBC, InflowBC, PeriodicBC,
+export DirichletBC, OutflowBC, InflowBC, PeriodicBC,
BoundaryConditions, isperiodic,
broadcast_boundaryconditions!
@@ -72,7 +71,7 @@ julia> rhs_bc = OutflowBC()
julia> bdryconds = BoundaryConditions(lhs_bc, rhs_bc, rsolver)
```
"""
-struct BoundaryConditions{T_BC_LHS<:AbstractBC,
+struct BoundaryConditions{T_BC_LHS<:AbstractBC,
T_BC_RHS<:AbstractBC,
T_RSolver<:AbstractRiemannSolver}
lhs_bc::T_BC_LHS
diff --git a/src/cache.jl b/src/cache.jl
index 1726f250740f8e716570cabfe667cee232754cab..56c424f376feb99469540c9fb17bf20a2291ad3d 100644
--- a/src/cache.jl
+++ b/src/cache.jl
@@ -150,7 +150,7 @@ end
function register_variables!(C::Cache;
- dynamic_variablenames = [],
+ dynamic_variablenames = [],
rhs_variablenames = [],
static_variablenames = [],
cell_variablenames = [],
diff --git a/src/callbacks.jl b/src/callbacks.jl
index 120e3ee8cc0cc6a0eb3ae6687302c28933ffeca1..bc4a3838131e1bfc2145a6df48919f00859392e7 100644
--- a/src/callbacks.jl
+++ b/src/callbacks.jl
@@ -147,6 +147,8 @@ end
Base.push!(cbs::CallbackSet, cb::AbstractCallback) = push!(cbs.callbacks, cb)
+Base.push!(cbs::CallbackSet, ::Nothing) = nothing
+Base.append!(cbs::CallbackSet, cbs2::CallbackSet) = append!(cbs.callbacks, cbs2.callbacks)
function Base.show(io::IO, cbs::CallbackSet)
@@ -223,6 +225,9 @@ mutable struct SaveCallback <: AbstractCallback
end
# normalize and sanitize filename
+ if length(filename) == 0
+ error("SaveCallback: Invalid filename = '$filename' ")
+ end
filename = Filesystem.abspath(filename)
if filename[end-2:end] != ".h5"
filename = "$(filename).h5"
diff --git a/src/dg1d.jl b/src/dg1d.jl
index 8aead18e4ce59dc52f3249b24cf11888b2bc1750..f2d8535638546eadabb9606f9b6406f7fb03c11c 100644
--- a/src/dg1d.jl
+++ b/src/dg1d.jl
@@ -2,9 +2,7 @@ module dg1d
using Reexport
import Base: Filesystem, Threads
- using Configurations
using DataFrames
- using DocStringExtensions
using Glob
using InteractiveUtils
using Jacobi
@@ -22,15 +20,15 @@ module dg1d
using StaticArrays
using StructArrays
@reexport using TOML
+ @reexport using ToggleableAsserts
- export Maybe
+ export Maybe, main, Environment
# utils
include("utils.jl")
include("rootfinders.jl")
include("with_signature.jl")
include("units.jl")
- include("parameters.jl")
include("flatdict.jl")
include("output.jl")
@@ -45,9 +43,6 @@ module dg1d
# custom evolution
include("evolve.jl")
- # evolution based on OrdinaryDiffEq pkg
- # using OrdinaryDiffEq
- # include("evolve_ordinarydiffeqs.jl")
# interfaces
include("equations.jl")
@@ -63,6 +58,10 @@ module dg1d
include("segment.jl")
+ include("types.jl")
+ include("projects/projects.jl")
+ include("parameters.jl")
+
## projects
include("projects/TCI/TCI.jl")
@reexport using .TCI
@@ -83,4 +82,6 @@ module dg1d
include("projects/GRHD/GRHD.jl")
+ include("main.jl")
+
end # module
diff --git a/src/evolve.jl b/src/evolve.jl
index 4402c9e118538e80190dcb162f542913ed6d2e55..04da398b05bad875d45d79d21706d4d5e8100614 100644
--- a/src/evolve.jl
+++ b/src/evolve.jl
@@ -3,11 +3,7 @@
#######################################################################
-"""
-$(FUNCTIONNAME)
-$(FIELDS)
-"""
-@option struct Parameters_Evolution <: dg1d.AbstractParameters
+struct Parameters_Evolution
"Time stepping algorithm, choose one out of LSERK, RK4, SSPRK3, ..."
algorithm::String = "LSERK"
"End time of evolution"
@@ -15,12 +11,6 @@ $(FIELDS)
end
-function sanitize(prms::Parameters_Evolution)
- @assert prms.algorithm in [ "LSERK", "RK4", "SSPRK3" ]
- @assert tstart < tend
-end
-
-
#######################################################################
# Butcher-Tablaeus #
#######################################################################
@@ -424,19 +414,19 @@ Generate an Explicit-Runge-Kutta stepper for iterate `u0` that uses scheme `sche
"""
function make_ERK_stepper(u0, scheme::Symbol, callback_substep)
- coeffs, method = if scheme == :MIDPT
+ coeffs, method = if scheme == :midpt
BT_MIDPT(), baseERK_step!
- elseif scheme == :RK4
+ elseif scheme == :rk4
BT_RK4(), baseERK_step!
- elseif scheme == :RK4_twothirds
+ elseif scheme == :rk4_twothirds
BT_RK4_twothirds(), baseERK_step!
- elseif scheme == :RK4_Ralston
+ elseif scheme == :rk4_ralston
BT_RK4_Ralston(), baseERK_step!
- elseif scheme == :SSPRK3
+ elseif scheme == :ssprk3
BT_SSPRK3(), baseERK_step!
- elseif scheme == :LSERK4
+ elseif scheme == :lserk4
COEFFS_LSERK4(), LSERK4_step!
- elseif scheme == :RKF10
+ elseif scheme == :rkf10
BT_RKF10(), baseERK_step!
else
error("Requested ERK scheme not implemented!")
diff --git a/src/evolve_ordinarydiffeqs.jl b/src/evolve_ordinarydiffeqs.jl
deleted file mode 100644
index f2db726a1674a46228769c565841baf6f029096d..0000000000000000000000000000000000000000
--- a/src/evolve_ordinarydiffeqs.jl
+++ /dev/null
@@ -1,25 +0,0 @@
-function evolve!(prms, rhs!, u, p, callbacks::CallbackSet)
- timestep_prms = get_or_error(prms, "Evolution", "Missing section 'Evolution'")
- # TODO write getter that allows to select between parameters, e.g.
- # get(prms, "timespan|tend", missing), where | means error. This should
- # throw an error if both are set. Tbf, the above example is dumb becase we
- # never need tspan
- # TODO Make get function also take types such that it can throw if type does not
- # match
- # TODO Make get function also take types such that it can throw if type does not
- # match.
- # tspan = get(prms, "timespan", missing)
- tend = get_or_error(timestep_prms, "tend", "Missing parameter 'Evolution.tend'")
- algorithm = get(prms, "algorithm", "RK4") # also need the | here, e.g. RK4|SSPRK3|...
- tspan = (0, tend)
- problem = ODEProblem(rhs!, u, tspan, p)
- stepper = if algorithm == "RK4"
- RK4()
- else
- error("Algorithm '$algorihtm' not available yet!")
- end
- solve(problem, stepper, callback=callbacks,
- adapative=false, # turn of build-in adaptive time stepping, use callback instead
- )
- return
-end
diff --git a/src/main.jl b/src/main.jl
new file mode 100644
index 0000000000000000000000000000000000000000..ee047f5251deb740cb78863ced8228c1a57a4928
--- /dev/null
+++ b/src/main.jl
@@ -0,0 +1,219 @@
+"""
+ main(parfile)
+
+Run program specified by a `parfile` in TOML format.
+"""
+function main(parfile)
+
+ # normalize parfile name
+ parfile = abspath(parfile)
+
+ project_name, prms = load_parameters(parfile)
+
+ # TODO Use PrettyTables.jl for summary
+ # TODO Dump parameters also to file
+ # TODO Dump versioninfo() to some file
+ # TODO Dump Pkg.status() to some file
+ # TODO Dump `git status` to some file
+ # TODO Dump gethostname() and other system info
+ # summary = summarize(prms)
+ # println(summary)
+
+ # backup any previous output dir and setup a fresh one
+ # TODO Rename this to make_outputdir
+ outputdir = mk_outputdir(parfile)
+ cp(parfile, joinpath(outputdir, basename(parfile)))
+
+ # Setup the environment
+ mesh = Mesh(; N=prms["Mesh"]["n"], K=prms["Mesh"]["k"], xrange=prms["Mesh"]["xrange"],
+ basis=prms["Mesh"]["basis"], CFL=prms["Mesh"]["cfl"])
+ cache = Cache(mesh)
+ callbacks = CallbackSet()
+ callbacks_substeps = CallbackSet()
+ env = Environment(mesh, cache, callbacks, callbacks_substeps)
+
+ # Setup the project. This involves
+ # - registering variables in cache
+ # - setting up initial data
+ # - setting up boundary conditions
+ # - setting up callbacks
+ # Need to use Val(project_name) for dispatch here, because project_name is a runtime variable.
+ Project = project_type(Val(project_name))
+ # TODO Till I sorted out the issue with where to place bdry condition types, I have to
+ # return them here and forward them below
+ project, bdryconds = Project(env, prms)
+
+ # wrap project interface methods for useage with time stepper
+ function wrapped_rhs!(state_du, state_u, t)
+ wrap_dynamic_variables!(env.cache, state_u)
+ wrap_rhs_variables!(env.cache, state_du)
+ # provided by project interface
+ rhs!(env, project, bdryconds)
+ end
+ function wrapped_timestep!(state_u, t, n)
+ wrap_dynamic_variables!(env.cache, state_u)
+ # provided by project interface
+ timestep!(env, project)
+ end
+
+ # Setup project independent callbacks
+ # For the time being, we let the Projects freely push callbacks into env.callbacks.
+ # This should make no difference as long as we 1) don't run in parallel and 2) don't allow
+ # multiple projects. Might want to move that logic later from Project to here.
+ # TODO Add an AliveCallback that manages a 'status' file in the output folder
+ # TODO Add priority ordering to CallbackSet and callbacks.
+ # TODO Drop unused callbacks from callbackset upon creation
+
+ # Add some 'global' callbacks
+
+ savecb = SaveCallback(Symbol.(prms["Output"]["variables"]), cache, mesh,
+ joinpath(outputdir, "output.h5"),
+ CallbackTiming(every_dt=prms["Output"]["every_dt"],
+ every_iteration=prms["Output"]["every_iteration"]))
+
+ aligned_ts = prms["Output"]["aligned_ts"]
+ timealigned_savecb, wrapped_timestep! = if length(aligned_ts) > 0
+ atscb = TimeAlignedSaveCallback(aligned_ts, wrapped_timestep!, savecb)
+ atscb, atscb.aligned_timestep_fn
+ else
+ nothing, wrapped_timestep!
+ end
+
+ plotfn(u,t) = plot(env, project)
+ plotcb = PlotCallback(plotfn, CallbackTiming(every_dt=prms["Plot"]["every_dt"],
+ every_iteration=prms["Plot"]["every_iteration"]), pause=prms["Plot"]["pause"])
+
+ push!(callbacks, savecb)
+ push!(callbacks, timealigned_savecb)
+ push!(callbacks, plotcb)
+
+ # TODO Setup a logger, maybe want to setup even earlier
+ # TODO summarize current setup once and log to screen
+ # println(summarize(env.mesh))
+ # println(summarize(env.cache))
+ # println(summarize(env.callbacks))
+ # println(summarize(env.callbacks_substep))
+ # println(summarize(project))
+
+ # plot parameter list
+ # TODO Dumb this to a file
+ for (k,v) in prms
+ table = vcat([k ""], hcat(collect(keys(v)), collect(values(v))))
+ pretty_table(
+ table;
+ body_hlines = [1],
+ body_hlines_format = Tuple('─' for _ = 1:4),
+ alignment = [ :l, :l ],
+ noheader = true,
+ tf = tf_borderless,
+ columns_width = 30,
+ linebreaks = true,
+ autowrap = true,
+ crop = :none
+ )
+ println()
+ end
+
+ # Let's go!!!
+ evolve!(env, wrapped_rhs!, wrapped_timestep!, prms)
+
+ return
+end
+
+
+function load_parameters(parfile)
+
+ # read parameters
+ if !isfile(parfile)
+ error("Cannot locate parameter file '$parfile'")
+ end
+
+ user_prms = TOML.parsefile(parfile)
+ # switch to a flat dict where keys of nested dicts are concatenated with dots
+ user_prms = Dict(string(k) => flatten(v) for (k,v) in pairs(user_prms))
+
+ ### sanitize parameter sections
+ # supported top-level headers that aren't user projects
+ default_projects = [ :Mesh, :Evolution, :Plot, :Output, :TCI, :HRSC, :EquationOfState ]
+
+ for (k,v) in user_prms
+ if !(v isa AbstractDict)
+ error("Invalid parameter file: No top-level parameters allowed, found '$k = $v'")
+ end
+ end
+
+ # find a parameter section for which somebody registered a project with that name
+ all_projects = Symbol.(keys(user_prms))
+ maybe_user_projects = filter(all_projects) do k
+ return !(k in default_projects)
+ end
+ registered_projects = filter(isregistered ∘ Symbol, maybe_user_projects)
+ not_registered_projects = setdiff(maybe_user_projects, registered_projects)
+ if length(not_registered_projects) > 0
+ error("Invalid parameter file: Found unknown top-level sections: $(not_registered_projects)")
+ end
+ if length(registered_projects) == 0
+ error("Invalid parameter file: Missing a parameter section for a project, don't know which project to run ...")
+ end
+ if length(registered_projects) > 1
+ error("Invalid parameter file: Only one project parameter section allowed, found multiple: $(registered_projects)")
+ end
+
+ # we found a user project
+ project_name = Symbol(first(registered_projects))
+
+ # load default parameters
+ prms = Dict(string(s) => parameters(Val(s)) for s in default_projects)
+ prms[string(project_name)] = parameters(project_name)
+ # TODO parameters(Val) does not guarantee that dicts returned from parameter are flat.
+ prms = Dict(k => flatten(v) for (k,v) in pairs(prms))
+
+ # merge defaults with with user parameters
+ for (sec, u_prms) in pairs(user_prms)
+ def_prms = prms[sec]
+ for (k,v) in pairs(u_prms)
+ if !haskey(def_prms, k)
+ error("Invalid parameter file: Unsupported parameter found in section '$sec': '$k = $v'. Keep your parameter files clean!")
+ end
+ def_prms[k] = v
+ end
+ end
+
+ # let's keep things organized
+ prms = OrderedDict( k => OrderedDict(v) for (k,v) in pairs(prms) )
+ # sort sections
+ sort!(prms)
+ # sort each section
+ foreach(p -> sort!(p), values(prms))
+
+ return project_name, prms
+end
+
+
+function evolve!(env, rhs_fn, timestep_fn, prms)
+
+ @unpack tend, method = prms["Evolution"]
+ @unpack cfl = prms["Mesh"]
+ @unpack cache, callbacks, callbacks_substep = env
+
+ # sanity checks
+ if tend <= 0
+ error("TimeStepper: tend must be positive, found 'tend = $tend'")
+ end
+ if cfl < 0
+ error("TimeStepper: cfl must be positive, found 'cfl = $cfl'")
+ end
+ methods = [ "midpt", "rk4", "rk4_twothirds", "rk4_ralston", "ssprk3", "lserk4", "rkf10" ]
+ if !(method in methods)
+ error("TimeStepper: Unknown time stepping method '$method' requested, available methods: '$(methods...)'")
+ end
+
+ # TODO This should be hidden behind an interface
+ all_vars = values(cache.dynamic_variables.datadict)
+ u0 = vcat(all_vars...)
+
+ # Let's go!!!
+ evolve_ERK(rhs_fn, u0, timestep_fn, (0.0, tend), Symbol(method),
+ callback_fullstep = callbacks)#, callback_substep = callbacks_substep)
+
+end
diff --git a/src/mesh.jl b/src/mesh.jl
index 4533fc16fb5ab268a52e9648550468a789be587d..d63ddf8153e24f1b6793ddd18b9e6f057bb414ef 100644
--- a/src/mesh.jl
+++ b/src/mesh.jl
@@ -1,5 +1,4 @@
export Mesh, MeshInterpolator,
- make_Mesh, parse_parameters,
layout, n_points, n_cells,
grid_average, cellwise_average, cellwise_inner_product, broken_inner_product,
local_maxabs, locate_point, find_cell_index
@@ -17,7 +16,7 @@ export Mesh, MeshInterpolator,
#######################################################################
-@doc """
+"""
Mesh
Object that bundles grid and operators needed for modeling a Discontinuous Galerkin
@@ -34,73 +33,51 @@ struct Mesh
x::Matrix{Float64} # grid points
invjac::Matrix{Float64} # inverse jacobian of coordinate transformation
# from computational [-1,1] onto physical grid [xmin,xmax]
-
- function Mesh(; N, K, xrange, CFL, quadrature_method)
- element = SpectralElement(N, Symbol(quadrature_method))
- @unpack z = element
- xmin, xmax = xrange
- L = xmax - xmin
- dl = L/K
- a = [ xmin + dl * idx for idx = 0:K-1 ]
- b = a .+ dl
- Npts = length(z)
- x = Matrix{Float64}(undef, Npts, K)
- invjac = similar(x)
- for i=1:Npts, j=1:K
- aj = a[j]
- bj = b[j]
- zi = z[i]
- x[i,j] = 0.5 * ((bj - aj) * zi + (aj + bj))
- invjac[i,j] = 2.0 / (bj - aj)
- end
- # sanitize cell end points to have same x coordinate at common interfaces
- for j=1:K-1
- xl = x[end,j]
- xr = x[1,j+1]
- xmean = 0.5 * (xl + xr)
- x[end,j] = xmean
- x[1,j+1] = xmean
- end
- return new(element, K, xmin, xmax, CFL, L, dl, x, invjac)
- end
end
-#######################################################################
-# Factories #
-#######################################################################
-
-
-make_Mesh(prms) = Mesh(; parse_parameters(Mesh, prms)...)
-make_Mesh(; args...) = make_Mesh(FlatDict(args))
-
-
-######################################################################
-# Parameters #
-######################################################################
-
+function Mesh(; N=5, K=10, xrange=[-1.0,1.0], CFL=0.5, basis="lgl_lumped")
+
+ @toggled_assert N > 0
+ @toggled_assert K > 0
+ @toggled_assert CFL > 0
+ @toggled_assert length(xrange) == 2
+ @toggled_assert xrange[1] < xrange[2]
+ @toggled_assert basis in ["lgl", "glgl", "lgl_lumped"]
+
+ # TODO Remove this after updating parameter names in SpectralElement
+ quadrature_method = Dict("lgl" => :LGL, "glgl" => :GLGL, "lgl_lumped" => :LGL_lumped)[basis]
+
+ element = SpectralElement(N, Symbol(quadrature_method))
+ @unpack z = element
+ xmin, xmax = xrange
+ L = xmax - xmin
+ dl = L/K
+ a = [ xmin + dl * idx for idx = 0:K-1 ]
+ b = a .+ dl
+ Npts = length(z)
+ x = Matrix{Float64}(undef, Npts, K)
+ invjac = similar(x)
+ for i=1:Npts, j=1:K
+ aj = a[j]
+ bj = b[j]
+ zi = z[i]
+ x[i,j] = 0.5 * ((bj - aj) * zi + (aj + bj))
+ invjac[i,j] = 2.0 / (bj - aj)
+ end
+ # sanitize cell end points to have same x coordinate at common interfaces
+ for j=1:K-1
+ xl = x[end,j]
+ xr = x[1,j+1]
+ xmean = 0.5 * (xl + xr)
+ x[end,j] = xmean
+ x[1,j+1] = xmean
+ end
-function parse_parameters(::Type{Mesh}, prms)
- N = query(Int64, prms, "N", default=5)
- K = query(Int64, prms, "K", default=20)
- CFL = query(Float64, prms, "CFL", default=0.5)
- quadrature_method = query(String, prms, "quadrature_method", default="LGL",
- options=["LGL", "GLGL", "LGL_lumped"])
- xrange = query(Vector{Float64}, prms, "xrange", default = [-1.0,1.0])
- # sanitize
- @assert N > 0 "'N > 0' required"
- @assert K > 0 "'K > 0' required"
- @assert CFL > 0 "'CFL > 0' required"
- @assert length(xrange) == 2 "'xrange' must be a tuple"
- @assert xrange[1] < xrange[2] "'xrange' limits must be ordered"
- @assert quadrature_method in ["LGL", "GLGL", "LGL_lumped"] "'quadrature_method' invalid"
- return (; N, K, CFL, quadrature_method, xrange)
+ return Mesh(element, K, xmin, xmax, CFL, L, dl, x, invjac)
end
-@deprecate Mesh(prms) make_Mesh(prms)
-
-
#######################################################################
# Methods #
#######################################################################
@@ -261,7 +238,7 @@ Operator for interpolating data from one mesh to another.
# Examples
```julia
-julia> sample_mesh, reference_mesh = make_Mesh(; N=3, K=4); make_Mesh(; N=10, K=20);
+julia> sample_mesh, reference_mesh = Mesh(; N=3, K=4); Mesh(; N=10, K=20);
julia> reference_data = sin.(reference_mesh.x);
diff --git a/src/mwe.jl b/src/mwe.jl
deleted file mode 100644
index fcfe9a25f5ba68ec4df28ecfb84b751bc9d438c8..0000000000000000000000000000000000000000
--- a/src/mwe.jl
+++ /dev/null
@@ -1,9 +0,0 @@
-using StructArrays
-
-
-x, y, z = [ randn(5) for _ = 1:3 ]
-
-Vf = Vector{Float64}
-soa = StructArray{Tuple{Tuple{Vf,Vf}, Vf}}(((x,y),z))
-
-
diff --git a/src/parameters.jl b/src/parameters.jl
index 317b05b2e09339402d915b4825d4f70a961b73c9..655de8a73d7aa530a1157b84b50d5f47d3fbeb8d 100644
--- a/src/parameters.jl
+++ b/src/parameters.jl
@@ -1,48 +1,560 @@
-export AbstractParameters,
- sanitize
+export parameters, @parameters
-#######################################################################
-# Types #
-#######################################################################
+"""
+ @parameters name ex
+Macro to define parameters and their default values for a symbol `name`.
+If invoked with `MyProject` this generates the following methods for the `dg1d` project interface:
+- `dg1d.parameters(::Val{:MyProject})`
+- `dg1d.check_parameters(::Val{:MyProject}, prms)`
-@doc """
- abstract AbstractParameters
+---
-Supertype for participating in the Parameters interface.
+# Example
-Subtyped structs can implement the following methods to verify that
-provided parameters are sane:
-```
-sanitize(prms::T) where T<:AbstractParameters
+```julia
+@parameters MyProject begin
+
+ "some variable"
+ x = 1.0
+ @check x > 0.0
+
+ "another variable"
+ y = "String"
+
+end
```
"""
-abstract type AbstractParameters end
+macro parameters(name, expr)
+ return esc(_generate_parameters(name, expr))
+end
-#######################################################################
-# Interface methods #
-#######################################################################
+macro check(ex)
+ return ex
+end
+
+
+function _generate_parameters(name, expr)
+
+ if expr.head !== :block
+ error("@parameters: Second argument must be a begin ... end block containing docstrings and parameter defintions!")
+ end
+
+ pars = expr.args
+ pars = filter(p -> !(p isa LineNumberNode), pars)
+
+ docstrings = String[]
+ kvexprs = Expr[]
+ checks = Vector{Expr}[]
+ i_definition = 0
+ for p in pars
+
+ # if p isa Expr
+ # filter!(pp -> !(pp isa LineNumberNode), p.args)
+ # end
+
+ is_definition = (p isa Expr) && p.head === :macrocall &&
+ (p.args[1] === Symbol("@doc") || p.args[1] isa GlobalRef) && length(p.args) == 4
+ is_check = (p isa Expr) && p.head === :macrocall && p.args[1] === Symbol("@check") &&
+ length(p.args) == 3
+
+ is_valid_block = is_definition || is_check
+ if !is_valid_block
+ error("""
+ @parameters: Invalid parameter definition.
+ Expected either a block with a docstring, parameter name and default value
+ or an argument check, e.g.
+ \"this is a docstring for variable x\"
+ x = 1.0
+ @check x > 1.0
+ Instead found
+ $(string(p))
+ """)
+ end
+
+ if is_definition
+ # parse a parameter definition block
+
+ i_definition += 1
+ # extract docstring, mandatory
+ d = p.args[3]
+ if d isa Expr && d.head === :macrocall && d.args[1] === Symbol("@raw_str")
+ d = d.args[3]
+ end
-function sanitize(prms::AbstractParameters) end
+ # extract variable name and default value, mandatory
+ kv = p.args[4]
+ if !(kv isa Expr)
+ error("""
+ @parameters: Invalid parameter definition. Missing default value for
+ $(string(kv))
+ """)
+
+ end
+ if kv.head !== :(=)
+ error("""
+ @parameters: Invalid parameter definition. Expected pair of parameter name and default value, e.g.
+ x = 1.0
+ Instead found
+ $(string(kv))
+ """)
+ end
+
+ k, v = kv.args[1], kv.args[2]
+ if !(k isa Symbol)
+ error("""
+ @parameters: Invalid parameter name: '$k'"
+ """)
+ end
+
+ push!(docstrings, d)
+ push!(kvexprs, kv)
+ push!(checks, Expr[])
+ end
+
+ # parse a check, they are optional
+ if is_check
+
+ # we already verified above that p is a @check block
+ # we also stripped any LineNumberNodes there
+ ck = p.args[3]
+
+ if i_definition == 0
+ error("@parameters: Invalid definition. Found @check before any parameter definition: $ck")
+ end
+
+ # macro expansion of @check is delayed till check_parameters definition below
+ push!(last(checks), ck)
+ end
+
+ end
+
+ keys = [ kv.args[1] for kv in kvexprs ]
+ for i = 1:length(keys)-1
+ refk = keys[i]
+ foreach(keys[i+1:end]) do k
+ if refk === k
+ error("""
+ @parameters: Duplicated parameters found: '$k'
+ """)
+ end
+ end
+ end
+
+ # generate the docs
+ io = IOBuffer()
+ println(io, name)
+ for (d, kv) in zip(docstrings, kvexprs)
+ println(io, "- `", kv, "`\n ", replace(d, '\n' => "\n "))
+ end
+ docs = String(take!(io))
+
+ # compare key = value to "key" => value for Dict constructor
+ pairs = [ Expr(:call, :(=>), string(kv.args[1]), kv.args[2]) for kv in kvexprs ]
+
+ checkscode = Expr[]
+ for cks in checks
+ evalchecks = [ quote
+ result = try
+ $ck
+ catch
+ error("check_parameters for $($(string(name))): Failed to evaluate check '$($(string(ck)))'")
+ end
+ if !(result isa Bool)
+ error("check_parameters for $($(string(name))): Invalid check found '$($(string(ck)))'. Expected to return a Bool, but found $(typeof(result))")
+ end
+ if result == false
+ error("check_parameters for $($(string(name))): Check failed: '$($(string(ck)))'")
+ end
+ end
+ for ck in cks ]
+ code = quote $(evalchecks...) end
+ push!(checkscode, code)
+ end
+ unpack_prms = [ quote $k = prms["$($("$k"))"] end for k in keys ]
+
+ code = quote
+ @doc $docs dg1d.parameters(::Val{$(QuoteNode(name))}) = Dict($(pairs...))
+ function dg1d.check_parameters(::Val{$(QuoteNode(name))}, prms)
+ $(unpack_prms...)
+ $(checkscode...)
+ end
+ end
+
+ return code
+end
#######################################################################
-# Methods #
+# Parameter definitions #
#######################################################################
-function Base.show(io::IO, prms::AbstractParameters)
- if Configurations.is_option(prms)
- print(strip(to_toml(prms)))
- else
- print(prms)
- end
+@parameters Mesh begin
+
+ "Number of cells"
+ k = 10
+ @check k > 0
+
+ "Physical domain extensions in x direction"
+ xrange = [ -1.0, 1.0 ]
+ @check length(xrange) == 2
+ @check xrange[2] - xrange[1] > 0.0
+
+ # TODO Move to Evolution
+ "Courant-Friedrichs-Lewy factor"
+ cfl = 0.5
+ @check cfl > 0.0
+
+ "Polynomial order in each grid"
+ n = 5
+ @check n > 0
+
+ """
+ Nodal basis used for the polynomial approximation. Available options:
+ - `lgl_lumped` ... mass-lumped Legendre-Gauss-Lobatto grid
+ - `lgl` ... Legendre-Gauss-Lobatto grid
+ - `glgl` ... Generalized Legendre-Gauss-Lobatto grid
+ """
+ basis = "lgl_lumped"
+ @check basis in [ "lgl", "glgl", "lgl_lumped" ]
+
end
-function Base.show(io::IO, ::MIME"text/plain", prms::AbstractParameters)
- show(io, prms)
+@parameters Evolution begin
+
+ "Final time"
+ tend = 10.0
+ @check tend > 0.0
+
+ """
+ Time stepping method. Available options
+ - `midpt` ... 2nd order explicit Runge Kutta method -- midpoint formula
+ - `rk4` ... 4th order explicit Runge Kutta method
+ - `rk4_twothirds` ... 4th order explicit Runge Kutta method -- two thirds version
+ - `rk4_ralston` ... 4th order explicit Runge Kutta method -- Ralston's method,
+ - designed for minimal truncation error
+ - `ssprk3` ... 3rd order strong stability preserving Runge Kutta method
+ - `lserk4` ... 4th order low storage explicit Runge Kutta method
+ - `rkf10` ... 10th order Runge Kutta Feagin method
+ """
+ method = "lserk4"
+ @check method in [ "midpt", "lserk4", "rk4", "rk4_twothirds", "rk_ralston",
+ "ssprk3", "rkf10" ]
+
end
+
+
+@parameters Plot begin
+
+ "Call plot callback after every_iteration full time steps"
+ every_iteration = 0
+ @check every_iteration >= 0
+
+ "Call plot callback after every_dt simulation time"
+ every_dt = 0.0
+ @check every_dt >= 0
+
+ """
+ If true, pauses the simulation to display the plot.
+ If set to a positive number, the simulations is paused only for that
+ number of seconds (walltime).
+ """
+ pause = false
+ @check pause isa Bool || pause > 0.0
+
+ "List of variables that should be plotted"
+ variables = String[]
+
+end
+
+
+@parameters Output begin
+
+ "Output data after every_iteration full time steps"
+ every_iteration = 0
+ @check every_iteration >= 0
+
+ "Output data after every_dt simulation time"
+ every_dt = 0.0
+ @check every_dt >= 0
+
+ "List of variables that should be outputted"
+ variables = String[]
+
+ "Filename for .h5 output"
+ filename = ""
+ @check length(filename) > 0
+
+ "List of times at which data should be outputted"
+ aligned_ts = Float64[]
+ @check all(t -> t > 0, aligned_ts)
+ @check issorted(aligned_ts)
+
+end
+
+
+@parameters EquationOfState begin
+
+ """
+ Available options
+ - `polytrope`
+ - `hybrid`
+ - `idealgas`
+ """
+ eos = "polytrope"
+ @check eos in [ "polytrope", "hybrid", "idealgas" ]
+
+ "Polytropic constant"
+ polytrope_k = 1.0
+ @check polytrope_k > 0.0
+
+ "exponent constant"
+ polytrope_gamma = 1.0
+ @check polytrope_gamma > 0.0
+
+ "Thermalization constant"
+ hybrid_gamma = 1.0
+ @check hybrid_gamma > 0.0
+
+ """
+ Available options
+ - `polytrope` ... adjust with `polytrope_k, polytrope_gamma`
+ """
+ hybrid_cold_eos = "polytrope"
+ @check hybrid_cold_eos in [ "polytrope" ]
+
+ "Polytropic constant"
+ idealgas_gamma = 1.0
+ @check idealgas_gamma > 0.0
+
+end
+
+
+@parameters TCI begin
+
+ """
+ Available options:
+ - `none` ... disabled
+ - `diffusion` ... based on gradient of velocity field
+ - `mda` ... based on average of modal decay of solution
+ - `mdh` ... based on modal decay of highest order coefficient of solution
+ - `minmod` ... based on minmod limiter to detect undesired oscillations
+ - `threshold` ... flag if threshold is exceeded
+ - `entropy_production` ... flag if entropy production is exceeded
+ """
+ method = "none"
+ @check method in [ "none", "diffusion", "mda", "mdh", "minmod", "threshold",
+ "entropy_production" ]
+
+ "List of variable names to monitor. If more than one is given `combine_method`
+ controls how to add the TCI results up."
+ variables = String[]
+ @check length(variables) > 0
+
+ """
+ This controls on how to combine TCIs from different variables. Available options:
+ - `max`
+ - `min`
+ - `average`
+ - `sum`
+ """
+ combine_method = "max"
+ @check combine_method in [ "max", "min", "average", "sum" ]
+
+ "Flag if diffusion exceeds this value"
+ diffusion_max = 1.0
+ @check diffusion_max > 0
+
+ "Flag if modal decay rate falls below this value"
+ mda_threshold_min = 1.0
+ @check mda_threshold_min > 0
+
+ "Don't flag if modal decay rate exceeds this value"
+ mda_threshold_max = 3.0
+ @check mda_threshold_max > mda_threshold_min
+
+ "Flag if modal decay rate falls below this value"
+ mdh_threshold_min = 1.0
+ @check mdh_threshold_min > 0
+
+ "Don't flag if modal decay rate exceeds this value"
+ mdh_threshold_max = 3.0
+ @check mdh_threshold_max > mdh_threshold_min
+
+ "TODO"
+ minmod_threshold = 1e-5
+ @check minmod_threshold > 0
+
+ "Flag if monitor exceeds this value"
+ threshold = 1e-5
+ @check threshold > 0
+
+ "Flag if entropy production exceeds this value"
+ entropy_production_max = 1.0
+ @check entropy_production_max > 0
+
+end
+
+
+@parameters HRSC begin
+
+ """
+ Available options:
+ - `none` ... disabled
+ - `weno` ... Weighted Non-Oscillatory approximation
+ - `av` ... Artificial viscosity
+ - `filter`
+ """
+ method = "none"
+ @check method in [ "weno", "av", "filter" ]
+
+ """
+ Available options:
+ - `mda` ... modal decay average
+ - `mdh` ... modal decay based on highest modal coefficient
+ - `entropy` ... entropy production
+ - `db` ... diffusion based
+ - `tci` ... viscosity based on TCI monitor
+ """
+ av_method = "mda"
+ @check method in [ "mda", "mdh", "entropy", "db", "tci" ]
+
+ """
+ TODO
+ """
+ av_smoother_order = 1
+ @check av_smoother_order in [ 1, 2 ]
+
+ "TODO"
+ mda_cmax = 1.0
+ @check mda_cmax > 0
+
+ "TODO"
+ mdh_cmax = 1.0
+ @check mdh_cmax > 0
+
+ "TODO"
+ mdh_ca = 0.5
+ @check mdh_ca > 0
+
+ "TODO"
+ mdh_ck = 0.1
+ @check mdh_ck > 0
+
+ "TODO"
+ entropy_cmax = 1.0
+ @check entropy_cmax > 0
+
+ "TODO"
+ entropy_ce = 0.0
+ @check entropy_ce > 0
+
+ "TODO"
+ db_cmax = 1.0
+ @check db_cmax > 0
+
+ "TODO"
+ db_cb = 0.5
+ @check db_cb > 0
+
+ "TODO"
+ tci_threshold = 1e-5
+ @check tci_threshold > 0
+
+ "TODO"
+ tci_max = 1e-2
+ @check tci_max > 0
+
+ "Order of reconstruction polynomial"
+ weno_order = 3
+ @check weno_order > 0
+
+ """
+ Small number to prevent division by zero in WENO weights computation.
+ """
+ weno_beta = 1e-6
+ @check weno_beta > 0
+
+ """
+ Available options:
+ - `bernstein`
+ """
+ filter_method = "bernstein"
+ @check filter_method in [ "bernstein" ]
+
+end
+
+
+# #######################################################################
+# # Project parameter definitions #
+# #######################################################################
+#
+#
+# @parameters GRHD begin
+#
+# "Set to true to enable an artificial atmosphere in regions of low density"
+# use_atmosphere = false
+#
+# "Cutoff value below which an artifical atmosphere is imposed"
+# atmosphere_cutoff_density = 1e-10
+#
+# "Restmass density of artifical atmosphere"
+# atmosphere_density = 1e-10
+#
+# end
+#
+#
+# @parameters SRHD begin
+#
+# "Set to true to enable an artificial atmosphere in regions of low density"
+# use_atmosphere = false
+#
+# "Cutoff value below which an artifical atmosphere is imposed"
+# atmosphere_cutoff_density = 1e-10
+#
+# "Restmass density of artifical atmosphere"
+# atmosphere_density = 1e-10
+#
+# end
+#
+#
+# @parameters ScalarEquation begin
+#
+# "Available options: advection, burgers"
+# equation = "advection"
+#
+# "Wave speed for advection equation"
+# advection_speed = 1.0
+#
+# # TODO Do we support reflecitve BCs?
+# "Available options: periodic, timedependent, reflective"
+# boundaryconditions = "periodic"
+#
+# # TODO Do we support timedependent BCs?
+# "TODO Options for time dependent BCs"
+# boundaryconditions_timedependent = ""
+#
+# "Available options: sine, rectangle"
+# initialdata = "sine"
+#
+# ""
+# initialdata_sine_period = 2.0 * pi
+#
+# ""
+# initialdata_sine_amplitude = 1.0
+#
+# ""
+# initialdata_rectangle_width = 0.5
+#
+# ""
+# initialdata_rectangle_start = -0.5
+#
+# ""
+# initialdata_rectangle_amplitude = 1.0
+#
+# end
diff --git a/src/projects/EquationOfState/EquationOfState.jl b/src/projects/EquationOfState/EquationOfState.jl
index a3134334fb60c222168bc510a86daa2d19085447..f796ddbf4f01f421b751b0cfaa37db5639720554 100644
--- a/src/projects/EquationOfState/EquationOfState.jl
+++ b/src/projects/EquationOfState/EquationOfState.jl
@@ -200,15 +200,30 @@ include("hybrideos.jl")
function make_EquationOfState(prms)
- type = Symbol(query(prms, "type", default=missing))
- type === :missing && error("No equation of state provided!")
-
- concretes = [ Base.typename(ct).name for ct in concretetypes(AbstractEquationOfState) ]
- @assert type in concretes "'EquationOfState.type = $type' unknown. Available options: $concretes"
- eos_type = eval(type)
- eos_prms = parse_parameters(eos_type, prms)
- eos = eos_type(; eos_prms...)
- return eos
+
+ @unpack eos = prms
+ # TODO Rename hybrideos to hybrid and fall back to concretetypes to figure
+ # out available EoS implementations
+
+ if eos == "hybrid"
+ @unpack hybrid_cold_eos, hybrid_gamma = prms
+ cold_eos = if hybrid_cold_eos == "polytrope"
+ @unpack polytrope_k, polytrope_gamma = prms
+ Polytrope(polytrope_k, polytrope_gamma)
+ else
+ error("EquationOfState: Unknown cold equation of state requested: $hybrid_eos")
+ end
+ return HybridEoS(cold_eos, hybrid_gamma)
+ elseif eos == "idealgas"
+ @unpack idealgas_gamma = prms
+ return IdealGas(idealgas_gamma)
+ elseif eos == "polytrope"
+ @unpack polytrope_k, polytrope_gamma = prms
+ return Polytrope(polytrope_k, polytrope_gamma)
+ else
+ error("EquationOfState: Unknown equation of state requested: $eos")
+ end
+
end
diff --git a/src/projects/EquationOfState/polytrope.jl b/src/projects/EquationOfState/polytrope.jl
index f255a0629a50534d43558e688bcaeec0457465f0..37f1c3892a908fe4e1f5300196240efdc8509419 100644
--- a/src/projects/EquationOfState/polytrope.jl
+++ b/src/projects/EquationOfState/polytrope.jl
@@ -22,8 +22,8 @@ end
function dg1d.parse_parameters(::Type{Polytrope}, prms=Dict())
- K = query(Float64, prms, "K", default=100.0)
- Gamma = query(Float64, prms, "Gamma", default=2.0)
+ K = query(Float64, prms, "polytrope_k", default=100.0)
+ Gamma = query(Float64, prms, "polytrope_gamma", default=2.0)
@assert K > 0 "Polytropic constant must be positive!"
@assert Gamma > 0 "Polytropic index must be positiv!"
return (; K, Gamma)
diff --git a/src/projects/EulerEq/EulerEq.jl b/src/projects/EulerEq/EulerEq.jl
index c91fd5655a5025211d1a1474880b43757e377417..b0ef0845fdcb7dc4c2bfbb1f0335e942a54e2753 100644
--- a/src/projects/EulerEq/EulerEq.jl
+++ b/src/projects/EulerEq/EulerEq.jl
@@ -4,21 +4,52 @@ module EulerEq
using dg1d
using .EquationOfState
-using Match
using Plots
-using InteractiveUtils
-using PlotThemes
-import Base: Filesystem
-import Parameters: @with_kw, @unpack
+include("types.jl")
include("equation.jl")
include("rhs.jl")
include("initialdata.jl")
include("plots.jl")
include("callbacks.jl")
include("boundaryconditions.jl")
-include("main.jl")
+include("setup.jl")
+#######################################################################
+# dg1d project interface definitions #
+#######################################################################
+
+
+dg1d.project_type(::Val{:EulerEq}) = RHS
+
+# TODO Remove bdryconds
+dg1d.rhs!(env::Environment, p::RHS, bdryconds) = rhs!(env, p, p.hrsc, bdryconds...)
+
+dg1d.timestep!(env::Environment, p::RHS) = timestep(env, p, p.hrsc)
+
+dg1d.@parameters EulerEq begin
+
+ """
+ Initial data configuration. Available options:
+ - `smooth_isentropic_flow`
+ - `sod_shock_tube`
+ - `shu_osher_shock_tube_v1`
+ - `shu_osher_shock_tube_v2`
+ """
+ id = "sine"
+ @check id in [ "smooth_isentropic_flow", "sod_shock_tube",
+ "shu_osher_shock_tube_v1", "shu_osher_shock_tube_v2" ]
+
+ """
+ Boundary conditions. Available options:
+ - `periodic`
+ - `from_id`
+ """
+ bc = "periodic"
+ @check bc in [ "periodic", "from_id" ]
+
+end
+
end # module
diff --git a/src/projects/EulerEq/boundaryconditions.jl b/src/projects/EulerEq/boundaryconditions.jl
index 3c4035a1687008ca1127761c34448b246bbb9669..b9951173b2c27b2f469101293c6914148651d463 100644
--- a/src/projects/EulerEq/boundaryconditions.jl
+++ b/src/projects/EulerEq/boundaryconditions.jl
@@ -1,28 +1,35 @@
-function make_BoundaryConditions(rhs::RHS, prms)
- type = query(String, prms, "type")
-
- lhs_bc, rhs_bc = @match type begin
- "periodic" => (PeriodicBC(), PeriodicBC())
- "from initial data" => begin
- @unpack rho, q, E = get_dynamic_variables(rhs.cache)
- lhs_state = (rho[1], q[1], E[1])
- rhs_state = (rho[end], q[end], E[end])
- (DirichletBC(lhs_state), DirichletBC(rhs_state))
- end
- _ => error("Boundary condition type '$type' not implemented!")
+function make_BoundaryConditions(env, eq, rsolver, prms)
+
+ @unpack bc = prms
+
+ lhs_bc, rhs_bc = if bc == "periodic"
+ PeriodicBC(), PeriodicBC()
+ elseif bc == "from_id"
+ @unpack rho, q, E = get_dynamic_variables(env.cache)
+ lhs_state = (rho[1], q[1], E[1])
+ rhs_state = (rho[end], q[end], E[end])
+ (DirichletBC(lhs_state), DirichletBC(rhs_state))
+ else
+ error("Unknown boundary condition requested: '$bc'")
end
- return BoundaryConditions(lhs_bc, rhs_bc, rhs.rsolver)
+ return BoundaryConditions(lhs_bc, rhs_bc, rsolver)
end
-function make_BoundaryConditions_LDG(ldg_rsolver, prms)
+function make_BoundaryConditions_LDG(env, eq, ldg_rsolver, prms)
+
+ @unpack bc = prms
+
isnothing(ldg_rsolver) && return nothing
- type = query(String, prms, "type")
+
lhs_bc, rhs_bc = if type == "periodic"
PeriodicBC(), PeriodicBC()
- else
+ elseif bc == "from_id"
OutflowBC(), OutflowBC()
+ else
+ error("Unknown boundary condition requested: '$bc'")
end
+
return BoundaryConditions(lhs_bc, rhs_bc, ldg_rsolver)
end
diff --git a/src/projects/EulerEq/callbacks.jl b/src/projects/EulerEq/callbacks.jl
index f986eb1389bdb30fb37dccacc6344221acac49c5..90a9065bf09169a0e2a0b3a54e70c69fafc9ccc9 100644
--- a/src/projects/EulerEq/callbacks.jl
+++ b/src/projects/EulerEq/callbacks.jl
@@ -1,11 +1,11 @@
-function make_callback(rhs::RHS, bdryconds::BoundaryConditions)
- cbfn_equation(u, t) = callback_equation(u, t, rhs, isperiodic(bdryconds), rhs.equation)
+function make_callback(env, rhs::RHS, bdryconds::BoundaryConditions)
+ cbfn_equation(u, t) = callback_equation(u, t, env, rhs, isperiodic(bdryconds), rhs.equation)
cb_equation = FunctionCallback(cbfn_equation,
CallbackTiming(every_iteration=1,every_dt=0))
- cbfn_tci(u, t) = callback_tci(u, t, rhs, isperiodic(bdryconds), rhs.tci)
+ cbfn_tci(u, t) = callback_tci(u, t, env, rhs, isperiodic(bdryconds), rhs.tci)
cb_tci = FunctionCallback(cbfn_tci,
CallbackTiming(every_iteration=1, every_dt=0))
- cbfn_hrsc(u, t) = callback_hrsc(u, t, rhs, isperiodic(bdryconds), rhs.hrsc)
+ cbfn_hrsc(u, t) = callback_hrsc(u, t, env, rhs, isperiodic(bdryconds), rhs.hrsc)
cb_hrsc = FunctionCallback(cbfn_hrsc,
CallbackTiming(every_iteration=1,every_dt=0))
callbackset = CallbackSet(cb_equation, cb_tci, cb_hrsc)
@@ -18,8 +18,9 @@ end
#######################################################################
-function callback_equation(state_u, state_t, rhs, isperiodic, eq::EulerEquation)
- @unpack cache, equation = rhs
+function callback_equation(state_u, state_t, env, rhs, isperiodic, eq::EulerEquation)
+ @unpack cache = env
+ @unpack equation = rhs
wrap_dynamic_variables!(cache, state_u)
@unpack Sm1, S, flx_S, flx_Sm1 = get_static_variables(cache)
@unpack t, tm1 = get_global_variables(cache)
@@ -39,20 +40,20 @@ end
#######################################################################
-function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::Nothing) end
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::Nothing) end
-callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::AbstractReconstruction) = nothing
+callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::AbstractReconstruction) = nothing
-function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::Maybe{HRSC.AbstractHRSC})
- error("To be implemented!")
- error("not yet implemented for type $(typeof(hrsc))")
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::Maybe{HRSC.AbstractHRSC})
+ TODO()
end
-function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.AbstractArtificialViscosity)
- @unpack cache, mesh, equation = rhs
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::HRSC.AbstractArtificialViscosity)
+ @unpack cache, mesh = env
+ @unpack equation = rhs
@unpack mu, max_v = get_cell_variables(cache)
@unpack v = get_static_variables(cache)
@unpack rho = get_dynamic_variables(cache)
@@ -69,24 +70,24 @@ function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.AbstractArt
end
-function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.TCIViscosity)
- @unpack cache, mesh = rhs
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::HRSC.TCIViscosity)
+ @unpack cache, mesh = env
@unpack mu, flag = get_cell_variables(cache)
HRSC.compute_viscosity!(mu, flag, hrsc)
end
-function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.SmoothedArtificialViscosity)
- callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc.av)
- @unpack cache = rhs
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::HRSC.SmoothedArtificialViscosity)
+ callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc.av)
+ @unpack cache = env
@unpack mu = get_cell_variables(cache)
@unpack smoothed_mu = get_static_variables(cache)
smoothed_mu .= vec(hrsc.smoother(mu, isperiodic))
end
-function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.EntropyViscosity)
- @unpack cache, mesh = rhs
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::HRSC.EntropyViscosity)
+ @unpack cache, mesh = env
@unpack rho = get_dynamic_variables(cache)
@unpack S, Sm1, flx_S, flx_Sm1 = get_static_variables(cache)
@unpack mu, max_v = get_cell_variables(cache)
@@ -103,17 +104,17 @@ end
#######################################################################
-function callback_tci(state_u, state_t, rhs, isperiodic, tci::Nothing) end
+function callback_tci(state_u, state_t, env, rhs, isperiodic, tci::Nothing) end
-function callback_tci(state_u, state_t, rhs, isperiodic, tci::TCI.Threshold)
+function callback_tci(state_u, state_t, env, rhs, isperiodic, tci::TCI.Threshold)
error("Threshold TCI must be combined with an indicator that is supposed to be zero for " *
"for smooth solutions, but not for shocks. TODO Implement an indicator variable.")
end
-function callback_tci(state_u, state_t, rhs, isperiodic, tci::TCI.EntropyProduction)
- @unpack cache, mesh = rhs
+function callback_tci(state_u, state_t, env, rhs, isperiodic, tci::TCI.EntropyProduction)
+ @unpack cache, mesh = env
@unpack S, Sm1, flx_S, flx_Sm1, EP = get_static_variables(cache)
@unpack flag = get_cell_variables(cache)
@unpack EPjump_lhs, EPjump_rhs = get_bdry_variables(cache)
@@ -124,11 +125,11 @@ function callback_tci(state_u, state_t, rhs, isperiodic, tci::TCI.EntropyProduct
end
-function callback_tci(state_u, state_t, rhs, isperiodic,
+function callback_tci(state_u, state_t, env, rhs, isperiodic,
tci::Union{TCI.Diffusion, TCI.Minmod, TCI.Threshold,
TCI.ModalDecayAverage, TCI.ModalDecayHighest})
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack rho = get_dynamic_variables(cache)
@unpack flag = get_cell_variables(cache)
TCI.compute_indicator!(flag, rho, tci)
diff --git a/src/projects/EulerEq/equation.jl b/src/projects/EulerEq/equation.jl
index d0a2c28ad5501dcd8899f7363440b85b5818da49..8930e4519356984bf70e95fb6b184d5df3d486cf 100644
--- a/src/projects/EulerEq/equation.jl
+++ b/src/projects/EulerEq/equation.jl
@@ -1,6 +1,7 @@
-@with_kw struct EulerEquation{T_EoS<:EquationOfState.AbstractEquationOfState} <: AbstractEquation
- eos::T_EoS
-end
+#######################################################################
+# Flux methods #
+#######################################################################
+
@with_signature function flux(equation::EulerEquation)
diff --git a/src/projects/EulerEq/initialdata.jl b/src/projects/EulerEq/initialdata.jl
index 0d06af079eb7b1ccd0493ed695c5ddedbf49a039..ca96c2c0d702fa132cbb38e6ec950c4bced34ab8 100644
--- a/src/projects/EulerEq/initialdata.jl
+++ b/src/projects/EulerEq/initialdata.jl
@@ -1,18 +1,16 @@
-function initialdata!(rhs::RHS, prms)
- initialdata_equation_of_state!(rhs, rhs.equation.eos, prms)
- initialdata_entropy_variables!(rhs, prms)
- initialdata_hrsc!(rhs, rhs.hrsc, prms)
- initialdata_tci!(rhs, rhs.tci, prms)
+function initialdata!(env, rhs::RHS, prms)
+ initialdata_equation_of_state!(env, rhs, rhs.equation.eos, prms)
+ initialdata_entropy_variables!(env, rhs, prms)
+ initialdata_hrsc!(env, rhs, rhs.hrsc, prms)
+ initialdata_tci!(env, rhs, rhs.tci, prms)
end
-function initialdata_equation_of_state!(rhs::RHS, eos::EquationOfState.IdealGas, parameters)
+function initialdata_equation_of_state!(env, rhs::RHS, eos::EquationOfState.IdealGas, prms)
- type = query(parameters, "type",
- options=["smooth_isentropic_flow", "sod_shock_tube",
- "shu_osher_shock_tube_v1", "shu_osher_shock_tube_v2"])
+ @unpack id = prms
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
xmin, xmax = extrema(mesh)
γ = eos.gamma
@@ -20,40 +18,40 @@ function initialdata_equation_of_state!(rhs::RHS, eos::EquationOfState.IdealGas,
E(p, rho, v) = p / (γ - 1.0) + 0.5 * rho * v^2
- tuple_data = if type == "smooth_isentropic_flow"
+ tuple_data = if id == "smooth_isentropic_flow"
# Glaubitz, PhD thesis 2020, p. 219
- @assert isapprox(xmin, -1.0) && isapprox(xmax, 1.0)
- @assert γ == 3.0
+ @toggled_assert isapprox(xmin, -1.0) && isapprox(xmax, 1.0)
+ @toggled_assert γ == 3.0
Ï0 = @. 1.0 + 0.5 * sin(Ï€ * x)
v0 = @. 0.0 * x
p0 = @. Ï0^γ
q0 = @. Ï0 * v0
E0 = @. E(p0, Ï0, v0)
(Ï0, q0, E0, p0)
- elseif type == "sod_shock_tube"
- @assert isapprox(xmin, -1.0) && isapprox(xmax, 1.0)
- @assert γ == 1.4
+ elseif id == "sod_shock_tube"
+ @toggled_assert isapprox(xmin, -1.0) && isapprox(xmax, 1.0)
+ @toggled_assert γ == 1.4
Ï0 = [ (xi < 0.0) ? 1.0 : 0.125 for xi in x ]
v0 = [ (xi < 0.0) ? 0.0 : 0.0 for xi in x ]
p0 = [ (xi < 0.0) ? 1.0 : 0.1 for xi in x ]
q0 = @. Ï0 * v0
E0 = @. E(p0, Ï0, v0)
(Ï0, q0, E0, p0)
- elseif type == "shu_osher_shock_tube_v1"
+ elseif id == "shu_osher_shock_tube_v1"
# taken from Glaubitz 2020, PhD thesis
# results do not agree, I believe that ID is wrong
- @assert isapprox(xmin, 0.0) && isapprox(xmax, 4.0)
- @assert γ == 1.4
+ @toggled_assert isapprox(xmin, 0.0) && isapprox(xmax, 4.0)
+ @toggled_assert γ == 1.4
Ï0 = [ (xi <= 2.5) ? 3.857143 : 1.0 + 0.2 * sin(10*Ï€*xi) for xi in x ]
v0 = [ (xi <= 2.5) ? 2.629369 : 0.0 for xi in x ]
p0 = [ (xi <= 2.5) ? 10.333333 : 1.0 for xi in x ]
q0 = @. Ï0 * v0
E0 = @. E(p0, Ï0, v0)
(Ï0, q0, E0, p0)
- elseif type == "shu_osher_shock_tube_v2"
+ elseif id == "shu_osher_shock_tube_v2"
# taken from Yu, Hesthaven, 2019
- @assert isapprox(xmin, -5.0) && isapprox(xmax, 5.0)
- @assert γ == 1.4
+ @toggled_assert isapprox(xmin, -5.0) && isapprox(xmax, 5.0)
+ @toggled_assert γ == 1.4
Ï0 = [ (xi <= -4.0) ? 3.857143 : 1.0 + 0.2 * sin(5*xi) for xi in x ]
v0 = [ (xi <= -4.0) ? 2.629369 : 0.0 for xi in x ]
p0 = [ (xi <= -4.0) ? 10.333333 : 1.0 for xi in x ]
@@ -61,7 +59,7 @@ function initialdata_equation_of_state!(rhs::RHS, eos::EquationOfState.IdealGas,
E0 = @. E(p0, Ï0, v0)
(Ï0, q0, E0, p0)
else
- error("Initial data of type $type not implemented!")
+ error("Initial data of type $id not implemented!")
end
# pack data
@@ -76,8 +74,9 @@ function initialdata_equation_of_state!(rhs::RHS, eos::EquationOfState.IdealGas,
end
-function initialdata_entropy_variables!(rhs::RHS, prms)
- @unpack cache, equation = rhs
+function initialdata_entropy_variables!(env, rhs::RHS, prms)
+ @unpack cache = env
+ @unpack equation = rhs
@unpack t, tm1 = get_global_variables(cache)
@unpack S, Sm1, flx_S, flx_Sm1 = get_static_variables(cache)
tm1 .= 0.0
@@ -89,20 +88,20 @@ function initialdata_entropy_variables!(rhs::RHS, prms)
end
-function initialdata_hrsc!(rhs::RHS, hrsc::Maybe{HRSC.AbstractHRSC}, prms) end
+function initialdata_hrsc!(env, rhs::RHS, hrsc::Maybe{HRSC.AbstractHRSC}, prms) end
-function initialdata_hrsc!(rhs::RHS, hrsc::HRSC.AbstractArtificialViscosity, prms)
- @unpack mu = get_cell_variables(rhs.cache)
+function initialdata_hrsc!(env, rhs::RHS, hrsc::HRSC.AbstractArtificialViscosity, prms)
+ @unpack mu = get_cell_variables(env.cache)
@. mu = 0.0
end
-function initialdata_tci!(rhs::RHS, tci::Nothing, prms) end
+function initialdata_tci!(env, rhs::RHS, tci::Nothing, prms) end
-function initialdata_tci!(rhs::RHS, tci::TCI.EntropyProduction, prms)
- @unpack cache = rhs
+function initialdata_tci!(env, rhs::RHS, tci::TCI.EntropyProduction, prms)
+ @unpack cache = env
@unpack EP = get_static_variables(cache)
@unpack EPjump_lhs, EPjump_rhs = get_bdry_variables(cache)
@. EP = 0.0
@@ -111,7 +110,7 @@ function initialdata_tci!(rhs::RHS, tci::TCI.EntropyProduction, prms)
end
-function initialdata_tci!(rhs::RHS, tci::TCI.AbstractTCI, prms)
- @unpack flag = get_cell_variables(rhs.cache)
+function initialdata_tci!(env, rhs::RHS, tci::TCI.AbstractTCI, prms)
+ @unpack flag = get_cell_variables(env.cache)
@. flag = 0.0
end
diff --git a/src/projects/EulerEq/main.jl b/src/projects/EulerEq/main.jl
deleted file mode 100644
index 862feee8d919f81ea4574c31a387f3c434dcbec2..0000000000000000000000000000000000000000
--- a/src/projects/EulerEq/main.jl
+++ /dev/null
@@ -1,137 +0,0 @@
-using OrderedCollections
-"""
- main(parameters_filename)
-
-Run program specified by TOML file `parameters_filename`.
-
-The following sections in the parameter file are recognized:
-- [`Mesh`](@ref)
-- [`EquationOfState`](@ref)
-- [`Initialdata`](@ref)
-- [`Timestep`](@ref)
-- [`HRSC`](@ref)
-- [`Output`](@ref)
-"""
-function main(parameters_filename=""; override::AbstractDict=Dict())
-
- override_prms = FlatDict(override)
-
- parameters_filename = Filesystem.abspath(parameters_filename)
- outputdir = mk_outputdir(parameters_filename)
- Filesystem.cp(parameters_filename, joinpath(outputdir, Filesystem.basename(parameters_filename)))
-
- # read parameters
- if !isfile(parameters_filename)
- error("Cannot locate parameter file '$parameters_filename'")
- end
- prms = FlatDict(TOML.parsefile(parameters_filename))
- prms = merge(prms, override_prms)
- ordered_prms = OrderedDict(prms)
- display(sort(ordered_prms))
-
- default_eos_prms = Dict("type"=>"IdealGas", "gamma"=>3.0)
- eos_prms = trim_headers(prms, "EquationOfState")
- amend!(eos_prms, default_eos_prms)
- diff = diffkeys(eos_prms, default_eos_prms)
- @assert length(diff) == 0 "Section 'Plot': Unknown options '$diff' encountered!"
- eos = EquationOfState.make_EquationOfState(eos_prms)
-
- # setup RHS
- mesh = make_Mesh(trim_headers(prms, "Mesh"))
- cache = Cache(mesh)
- equation = EulerEquation(eos)
- tci = TCI.make_TCI(trim_headers(prms, "TCI"), mesh)
- hrsc = HRSC.make_HRSC(trim_headers(prms, "HRSC"), mesh)
- rsolver = ApproxRiemannSolver(flux, speed, equation)
- ldg_rsolver, av_rsolver = if hrsc isa AbstractArtificialViscosity
- ApproxRiemannSolver(ldg_flux, ldg_speed, equation),
- ApproxRiemannSolver(av_flux, av_flux, equation)
- else
- nothing, nothing
- end
- rhs = RHS(cache=cache, mesh=mesh, equation=equation, tci=tci, hrsc=hrsc,
- rsolver=rsolver, ldg_rsolver=ldg_rsolver, av_rsolver=av_rsolver)
-
- _register_variables!(rhs)
- register_variables!(cache, rsolver)
- !isnothing(ldg_rsolver) && register_variables!(cache, ldg_rsolver)
- !isnothing(av_rsolver) && register_variables!(cache, av_rsolver, dont_register=true)
- display(cache)
-
- initialdata!(rhs, trim_headers(prms, "Initialdata"))
-
- default_boundaryconditions_prms = Dict("type"=>"periodic")
- boundaryconditions_prms = trim_headers(prms, "BoundaryConditions")
- amend!(boundaryconditions_prms, default_boundaryconditions_prms)
- diff = diffkeys(boundaryconditions_prms, default_boundaryconditions_prms)
- @assert length(diff) == 0 "Section 'Plot': Unknown options '$diff' encountered!"
- bdryconds = make_BoundaryConditions(rhs, boundaryconditions_prms)
- ldg_bdryconds, av_bdryconds = if hrsc isa AbstractArtificialViscosity
- make_BoundaryConditions_LDG(ldg_rsolver, boundaryconditions_prms),
- make_BoundaryConditions_LDG(av_rsolver, boundaryconditions_prms)
- else
- nothing, nothing
- end
-
- rhs_fn!(du, u, t) = rhs!(du, u, t, rhs, bdryconds, ldg_bdryconds, av_bdryconds)
- timestep_fn(u, t, n) = timestep(u, t, rhs)
-
- ### Callbacks
- rhscb = make_callback(rhs, bdryconds)
-
- # Plotting
- default_plot_prms = Dict("every_iteration"=>0, "every_dt"=>0.0, "pause"=>false)
- plot_prms = trim_headers(prms, "Plot")
- dg1d.amend!(plot_prms, default_plot_prms)
- diff = diffkeys(plot_prms, default_plot_prms)
- @assert length(diff) == 0 "Section 'Plot': Unknown options '$diff' encountered!"
- theme(:dark, markerstrokewidth=0.1, marker=:circle)
- plotfn(u,t) = plot(rhs)
- plotcb = PlotCallback(plotfn,
- CallbackTiming(every_dt=plot_prms["every_dt"],
- every_iteration=plot_prms["every_iteration"]),
- pause=plot_prms["pause"])
- default_output_prms = Dict("every_iteration"=>0,
- "every_dt"=>0.5,
- "variables"=>["rho", "q", "E"],
- "filename"=>joinpath(outputdir, "output"),
- "aligned_times"=>[])
- output_prms = trim_headers(prms, "Output")
- amend!(output_prms, default_output_prms)
- diff = diffkeys(output_prms, default_output_prms)
- @assert length(diff) == 0 "Section 'Output': Unknown options '$diff' encountered!"
- savecb = SaveCallback(Symbol.(output_prms["variables"]), cache, mesh,
- output_prms["filename"],
- CallbackTiming(every_dt=output_prms["every_dt"],
- every_iteration=output_prms["every_iteration"])
- )
-
- # Combine callbacks
- callbacks = CallbackSet(rhscb.callbacks..., savecb, plotcb)
-
- aligned_times = get(output_prms, "aligned_times", [])
- timealigned_savecb, timestep_fn = if length(aligned_times) > 0
- atscb = TimeAlignedSaveCallback(aligned_times, timestep_fn, savecb)
- atscb, atscb.aligned_timestep_fn
- else
- nothing, timestep_fn
- end
-
- !isnothing(timealigned_savecb) && push!(callbacks, timealigned_savecb)
-
- default_evolve_prms = Dict("algorithm"=>"LSERK4", "tend"=>1.0)
- evolve_prms = trim_headers(prms, "Evolution")
- amend!(evolve_prms, default_evolve_prms)
- diff = diffkeys(evolve_prms, default_evolve_prms)
- @assert length(diff) == 0 "Section 'Evolution': Unknown options '$diff' encountered!"
-
- @unpack rho, q, E = get_dynamic_variables(cache)
- state_u = vcat(rho, q, E)
-
- state_u = dg1d.evolve_ERK(rhs_fn!, state_u, timestep_fn,
- (0.0, Float64(evolve_prms["tend"])),
- Symbol(evolve_prms["algorithm"]),
- callback_fullstep=callbacks)
-
- return
-end
diff --git a/src/projects/EulerEq/plots.jl b/src/projects/EulerEq/plots.jl
index 93a6b75ab2e21a6dfc769eda6c52a6c36ba465ab..a57b2d70c0a805d83b199213e8f0e7e66de0a252 100644
--- a/src/projects/EulerEq/plots.jl
+++ b/src/projects/EulerEq/plots.jl
@@ -1,8 +1,8 @@
-function plot(rhs::RHS)
- plt_eq = plot_equation(rhs)
+function plot(env, rhs::RHS)
+ plt_eq = plot_equation(env, rhs)
plt_rhs = nothing # plot_rhs(rhs)
- plt_hrsc = plot_hrsc(rhs, rhs.hrsc)
- plt_tci = plot_tci(rhs, rhs.tci)
+ plt_hrsc = plot_hrsc(env, rhs, rhs.hrsc)
+ plt_tci = plot_tci(env, rhs, rhs.tci)
# drop missing plots and figure out layout
list_sub_plts = [ plt_eq, plt_rhs, plt_hrsc, plt_tci ]
@@ -11,7 +11,7 @@ function plot(rhs::RHS)
sub_plt = Plots.plot(list_sub_plts..., layout=layout)
# combine plots
- @unpack cache = rhs
+ @unpack cache = env
@unpack t = get_global_variables(cache)
plt = Plots.plot(sub_plt, layout=layout, size=(1400,1000), plot_title="t=$(t[1])",
left_margin=10Plots.mm)
@@ -24,9 +24,9 @@ end
#######################################################################
-function plot_equation(rhs)
+function plot_equation(env, rhs)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack rho, q, E = get_dynamic_variables(cache)
@unpack p, eps, S = get_static_variables(cache)
@unpack t = get_global_variables(cache)
@@ -56,7 +56,7 @@ function plot_equation(rhs)
end
plt = Plots.plot(plts..., layout=(3,2))
-
+
return plt
end
@@ -66,12 +66,12 @@ end
#######################################################################
-function plot_hrsc(rhs, hrsc::Nothing) end
+function plot_hrsc(env, rhs, hrsc::Nothing) end
-function plot_hrsc(rhs, hrsc::HRSC.AbstractArtificialViscosity)
+function plot_hrsc(env, rhs, hrsc::HRSC.AbstractArtificialViscosity)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack mu = get_cell_variables(cache)
@unpack x = mesh
@@ -93,11 +93,11 @@ function plot_hrsc(rhs, hrsc::HRSC.AbstractArtificialViscosity)
end
-function plot_hrsc(rhs, hrsc::HRSC.SmoothedArtificialViscosity)
+function plot_hrsc(env, rhs, hrsc::HRSC.SmoothedArtificialViscosity)
- plt_mu = plot_hrsc(rhs, hrsc.av)
+ plt_mu = plot_hrsc(env, rhs, hrsc.av)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack smoothed_mu = get_static_variables(cache)
@unpack x = mesh
@@ -122,7 +122,7 @@ function plot_hrsc(rhs, hrsc::HRSC.SmoothedArtificialViscosity)
end
-plot_hrsc(rhs, hrsc::AbstractReconstruction) = nothing
+plot_hrsc(env, rhs, hrsc::AbstractReconstruction) = nothing
#######################################################################
@@ -130,12 +130,12 @@ plot_hrsc(rhs, hrsc::AbstractReconstruction) = nothing
#######################################################################
-function plot_tci(rhs, tci::Nothing) end
+function plot_tci(env, rhs, tci::Nothing) end
-function plot_tci(rhs, tci::TCI.AbstractTCI)
+function plot_tci(env, rhs, tci::TCI.AbstractTCI)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack flag = get_cell_variables(cache)
@unpack x = mesh
@@ -165,9 +165,9 @@ function plot_tci(rhs, tci::TCI.AbstractTCI)
end
-function plot_tci(rhs, tci::TCI.EntropyProduction)
+function plot_tci(env, rhs, tci::TCI.EntropyProduction)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack EP = get_static_variables(cache)
@unpack flag = get_cell_variables(cache)
@unpack x = mesh
@@ -203,9 +203,9 @@ end
#######################################################################
-function plot_rhs(rhs)
+function plot_rhs(env, rhs)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack rhs_rho, rhs_q, rhs_E = get_rhs_variables(cache)
@unpack t = get_global_variables(cache)
@unpack x = mesh
diff --git a/src/projects/EulerEq/rhs.jl b/src/projects/EulerEq/rhs.jl
index 0895a2ca1d28c096ee0be4c8741c25c8d5f58d92..c5f8061acbcf1e868e2abf0b30f3105118957ad6 100644
--- a/src/projects/EulerEq/rhs.jl
+++ b/src/projects/EulerEq/rhs.jl
@@ -1,97 +1,11 @@
-Base.@kwdef mutable struct RHS{T_Equation <:EulerEquation,
- T_RSolver <:AbstractRiemannSolver,
- T_HRSC <:Maybe{HRSC.AbstractHRSC},
- T_TCI <:Maybe{TCI.AbstractTCI},
- T_LDG_RSolver<:Maybe{AbstractRiemannSolver},
- T_AV_RSolver<:Maybe{AbstractRiemannSolver}}
-
- cache::Cache
- mesh::Mesh
- equation::T_Equation
- rsolver::T_RSolver
- hrsc::T_HRSC = nothing
- tci::T_TCI = nothing
- ldg_rsolver::T_LDG_RSolver = nothing
- av_rsolver::T_AV_RSolver = nothing
-end
-
-
-#######################################################################
-# Register variables #
-#######################################################################
-
-
-# Julia limitation/bug: Overloading imported functions without their Module prefix shadows their
-# definition in this module without a warning or error. This does not happen for methods
-# defined in Base.
-# Because of that we prefix with _ below.
-# Cf. https://github.com/JuliaLang/julia/issues/15510
-function _register_variables!(rhs::RHS)
- register_variables!(rhs.cache,
- dynamic_variablenames = (:rho, :q, :E,), # density, momentum density, energy
- rhs_variablenames = (:rhs_rho, :rhs_q, :rhs_E),
- static_variablenames = (:flx_rho, :flx_q, :flx_E,
- :p, :eps, :v, # pressure, internal energy, speed
- :S, :Sm1, # entropy
- :flx_S, :flx_Sm1), # entropy flux
- cell_variablenames = (:max_v,),
- global_variablenames = (:t, :tm1) # recent time stamps
- )
- _register_variables!(rhs, rhs.hrsc)
- _register_variables!(rhs, rhs.tci)
-end
-
-
-_register_variables!(rhs, tci_or_hrsc::Nothing) = nothing
-
-_register_variables!(rhs, hrsc::AbstractReconstruction) = nothing
-
-
-function _register_variables!(rhs, hrsc::HRSC.AbstractArtificialViscosity)
- register_variables!(rhs.cache,
- static_variablenames = (:ldg_rho, :ldg_q, :ldg_E, # local DG variable
- :flx_ldg_rho, :flx_ldg_q, :flx_ldg_E), # local DG flux
- cell_variablenames = (:mu,) # 'one viscosity to rule them all'
- )
-end
-
-
-function _register_variables!(rhs, hrsc::HRSC.SmoothedArtificialViscosity)
- _register_variables!(rhs, hrsc.av)
- register_variables!(rhs.cache,
- static_variablenames = (:smoothed_mu,)
- )
-end
-
-
-function _register_variables!(rhs, tci::TCI.AbstractTCI)
- register_variables!(rhs.cache, cell_variablenames = (:flag,))
-end
-
-
-function _register_variables!(rhs, tci::TCI.EntropyProduction)
- register_variables!(rhs.cache,
- cell_variablenames = (:flag,),
- static_variablenames = (:EP,), # entropy production
- bdry_variablenames = (:EPjump_lhs, :EPjump_rhs)
- # entropy production due to jumps at interfaces
- )
-end
-
-
#######################################################################
# Timestep #
#######################################################################
-function timestep(state_u, t, rhs::RHS)
- wrap_dynamic_variables!(rhs.cache, state_u)
- return timestep(rhs, rhs.hrsc)
-end
-
-
-function timestep(rhs::RHS, hrsc::Maybe{HRSC.AbstractHRSC})
- @unpack cache, mesh, equation = rhs
+function timestep(env, rhs::RHS, hrsc::Maybe{HRSC.AbstractHRSC})
+ @unpack cache, mesh = env
+ @unpack equation = rhs
broadcast_volume!(speed, equation, cache)
v = get_variable(cache, :v)
@@ -110,8 +24,9 @@ function timestep(rhs::RHS, hrsc::Maybe{HRSC.AbstractHRSC})
end
-function timestep(rhs::RHS, hrsc::HRSC.AbstractArtificialViscosity)
- @unpack cache, mesh, equation = rhs
+function timestep(env, rhs::RHS, hrsc::HRSC.AbstractArtificialViscosity)
+ @unpack cache, mesh = env
+ @unpack equation = rhs
broadcast_volume!(speed, equation, cache)
v = get_variable(cache, :v)
@@ -142,20 +57,10 @@ end
#######################################################################
-function rhs!(state_du, state_u, t, rhs::RHS,
- bdryconds::BoundaryConditions, ldg_bdryconds::Maybe{BoundaryConditions},
- av_bdryconds::Maybe{BoundaryConditions})
-
- @unpack cache = rhs
- wrap_dynamic_variables!(cache, state_u)
- wrap_rhs_variables!(cache, state_du)
- rhs!(rhs, t, rhs.hrsc, bdryconds, ldg_bdryconds, av_bdryconds)
-end
-
-
-function rhs!(rhs::RHS, t, hrsc::Nothing, bdryconds, ldg_bdryconds, av_bdryconds)
+function rhs!(env, rhs::RHS, hrsc::Nothing, bdryconds, ldg_bdryconds, av_bdryconds)
- @unpack cache, mesh, equation, rsolver = rhs
+ @unpack cache, mesh = env
+ @unpack equation, rsolver = rhs
@unpack rho, q, E = get_dynamic_variables(cache)
@unpack flx_rho, flx_q, flx_E, p, eps = get_static_variables(cache)
@unpack rhs_rho, rhs_q, rhs_E = get_rhs_variables(cache)
@@ -165,7 +70,7 @@ function rhs!(rhs::RHS, t, hrsc::Nothing, bdryconds, ldg_bdryconds, av_bdryconds
broadcast_volume!(flux, equation, cache)
broadcast_faces!(lax_friedrich_flux, rsolver, cache, mesh)
- broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, t)
+ broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, 0.0)
compute_rhs_weak_form!(rhs_rho, flx_rho, lhs_numflx_rho, rhs_numflx_rho, mesh)
compute_rhs_weak_form!(rhs_q, flx_q, lhs_numflx_q, rhs_numflx_q, mesh)
@@ -175,9 +80,10 @@ function rhs!(rhs::RHS, t, hrsc::Nothing, bdryconds, ldg_bdryconds, av_bdryconds
end
-function rhs!(rhs::RHS, t, hrsc::AbstractReconstruction, bdryconds, ldg_bdryconds, av_bdryconds)
+function rhs!(env, rhs::RHS, hrsc::AbstractReconstruction, bdryconds, ldg_bdryconds, av_bdryconds)
- @unpack cache, mesh, equation, rsolver = rhs
+ @unpack cache, mesh = env
+ @unpack equation, rsolver = rhs
@unpack rho, q, E = get_dynamic_variables(cache)
@unpack flx_rho, flx_q, flx_E, p, eps = get_static_variables(cache)
@unpack flag = get_cell_variables(cache)
@@ -192,7 +98,7 @@ function rhs!(rhs::RHS, t, hrsc::AbstractReconstruction, bdryconds, ldg_bdrycond
broadcast_volume!(flux, equation, cache)
broadcast_faces!(lax_friedrich_flux, rsolver, cache, mesh)
- broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, t)
+ broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, 0.0)
compute_rhs_weak_form!(rhs_rho, flx_rho, lhs_numflx_rho, rhs_numflx_rho, mesh)
compute_rhs_weak_form!(rhs_q, flx_q, lhs_numflx_q, rhs_numflx_q, mesh)
@@ -202,14 +108,15 @@ function rhs!(rhs::RHS, t, hrsc::AbstractReconstruction, bdryconds, ldg_bdrycond
end
-function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractArtificialViscosity,
+function rhs!(env, rhs::RHS, hrsc::HRSC.AbstractArtificialViscosity,
bdryconds, ldg_bdryconds, av_bdryconds)
- @unpack cache, mesh, equation, rsolver, ldg_rsolver, av_rsolver = rhs
+ @unpack cache, mesh = env
+ @unpack equation, rsolver, ldg_rsolver, av_rsolver = rhs
@unpack rho, q, E = get_dynamic_variables(cache)
@unpack flx_rho, flx_q, flx_E, p, eps,
- ldg_rho, ldg_q, ldg_E,
+ ldg_rho, ldg_q, ldg_E,
flx_ldg_rho, flx_ldg_q, flx_ldg_E = get_static_variables(cache)
@unpack rhs_rho, rhs_q, rhs_E = get_rhs_variables(cache)
@unpack lhs_numflx_rho, rhs_numflx_rho,
@@ -222,7 +129,7 @@ function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractArtificialViscosity,
## solve auxiliary equation: q + ∂x u = 0
broadcast_volume!(ldg_flux, equation, cache)
broadcast_faces!(central_flux, ldg_rsolver, cache, mesh)
- broadcast_boundaryconditions!(central_flux, ldg_bdryconds, cache, mesh, t)
+ broadcast_boundaryconditions!(central_flux, ldg_bdryconds, cache, mesh, 0.0)
compute_rhs_weak_form!(ldg_rho, flx_ldg_rho, lhs_numflx_ldg_rho, rhs_numflx_ldg_rho, mesh)
compute_rhs_weak_form!(ldg_q, flx_ldg_q, lhs_numflx_ldg_q, rhs_numflx_ldg_q, mesh)
compute_rhs_weak_form!(ldg_E, flx_ldg_E, lhs_numflx_ldg_E, rhs_numflx_ldg_E, mesh)
@@ -230,10 +137,10 @@ function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractArtificialViscosity,
## compute rhs of regularized equation: ∂t u + ∂x f + ∂x μ q = 0
broadcast_volume!(av_flux, equation, cache)
broadcast_faces!(mean_flux, av_rsolver, cache, mesh)
- broadcast_boundaryconditions!(mean_flux, av_bdryconds, cache, mesh, t)
+ broadcast_boundaryconditions!(mean_flux, av_bdryconds, cache, mesh, 0.0)
broadcast_volume!(flux, equation, cache)
broadcast_faces!(lax_friedrich_flux, rsolver, cache, mesh)
- broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, t)
+ broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, 0.0)
@. flx_rho += flx_ldg_rho
@. flx_q += flx_ldg_q
diff --git a/src/projects/EulerEq/setup.jl b/src/projects/EulerEq/setup.jl
new file mode 100644
index 0000000000000000000000000000000000000000..32ba0bcf8dd16f0ce9ef8787d1f20d0ef1c59a89
--- /dev/null
+++ b/src/projects/EulerEq/setup.jl
@@ -0,0 +1,115 @@
+function RHS(env::Environment, prms)
+
+ # construct RHS
+ eos = EquationOfState.make_EquationOfState(prms["EquationOfState"])
+ equation = EulerEquation(eos)
+ tci = TCI.make_TCI(env.mesh, prms["TCI"])
+ hrsc = HRSC.make_HRSC(env.mesh, prms["HRSC"])
+ rsolver = ApproxRiemannSolver(flux, speed, equation)
+ ldg_rsolver, av_rsolver = if hrsc isa AbstractArtificialViscosity
+ ApproxRiemannSolver(ldg_flux, ldg_speed, equation),
+ ApproxRiemannSolver(av_flux, av_flux, equation)
+ else
+ nothing, nothing
+ end
+
+ rhs = RHS(equation, rsolver, hrsc, tci, ldg_rsolver, av_rsolver)
+
+ # register variables
+ # TODO add a ::Nothing overload for register_variables!
+ # TODO Somehow replace _register_variables! with register_variables!
+ @unpack cache = env
+ _register_variables!(cache, rhs)
+ register_variables!(cache, rsolver)
+ !isnothing(ldg_rsolver) && register_variables!(cache, ldg_rsolver)
+ !isnothing(av_rsolver) && register_variables!(cache, av_rsolver, dont_register=true)
+ display(cache)
+
+ # setup boundary conditions
+ bdryconds = make_BoundaryConditions(env, equation, rsolver, prms["EulerEq"])
+ ldg_bdryconds, av_bdryconds = if hrsc isa AbstractArtificialViscosity
+ make_BoundaryConditions_LDG(env, equation, ldg_rsolver, prms["EulerEq"]),
+ make_BoundaryConditions_LDG(env, equation, av_rsolver, prms["EulerEq"])
+ else
+ nothing, nothing
+ end
+
+ # TODO Bdry conditions should be moved to project and should only be just a 'lightweight'
+ # descriptor. Furthermore, they should also register variables in cache themselves.
+ # Note: Atm, only the 'from_initial_data' option would require a place to store the data
+ # somewhere. For periodic BCs we don't need extra storage at all.
+ # But what other types of boundary conditions will be needed in the future?
+
+ # setup initial data
+ initialdata!(env, rhs, prms["EulerEq"])
+
+ # setup callbacks
+ rhscb = make_callback(env, rhs, bdryconds)
+
+ append!(env.callbacks, CallbackSet(rhscb.callbacks...))
+
+ return rhs, (bdryconds, ldg_bdryconds, av_bdryconds)
+end
+
+
+#######################################################################
+# Register variables #
+#######################################################################
+
+
+# Julia limitation/bug: Overloading imported functions without their Module prefix shadows their
+# definition in this module without a warning or error. This does not happen for methods
+# defined in Base.
+# Because of that we prefix with _ below.
+# Cf. https://github.com/JuliaLang/julia/issues/15510
+function _register_variables!(cache, rhs::RHS)
+ register_variables!(cache,
+ dynamic_variablenames = (:rho, :q, :E,), # density, momentum density, energy
+ rhs_variablenames = (:rhs_rho, :rhs_q, :rhs_E),
+ static_variablenames = (:flx_rho, :flx_q, :flx_E,
+ :p, :eps, :v, # pressure, internal energy, speed
+ :S, :Sm1, # entropy
+ :flx_S, :flx_Sm1), # entropy flux
+ cell_variablenames = (:max_v,),
+ global_variablenames = (:t, :tm1) # recent time stamps
+ )
+ _register_variables!(cache, rhs.hrsc)
+ _register_variables!(cache, rhs.tci)
+end
+
+
+_register_variables!(cache, tci_or_hrsc::Nothing) = nothing
+
+_register_variables!(cache, hrsc::AbstractReconstruction) = nothing
+
+
+function _register_variables!(cache, hrsc::HRSC.AbstractArtificialViscosity)
+ register_variables!(rhs.cache,
+ static_variablenames = (:ldg_rho, :ldg_q, :ldg_E, # local DG variable
+ :flx_ldg_rho, :flx_ldg_q, :flx_ldg_E), # local DG flux
+ cell_variablenames = (:mu,) # 'one viscosity to rule them all'
+ )
+end
+
+
+function _register_variables!(cache, hrsc::HRSC.SmoothedArtificialViscosity)
+ _register_variables!(cache, hrsc.av)
+ register_variables!(cache,
+ static_variablenames = (:smoothed_mu,)
+ )
+end
+
+
+function _register_variables!(cache, tci::TCI.AbstractTCI)
+ register_variables!(cache, cell_variablenames = (:flag,))
+end
+
+
+function _register_variables!(cache, tci::TCI.EntropyProduction)
+ register_variables!(cache,
+ cell_variablenames = (:flag,),
+ static_variablenames = (:EP,), # entropy production
+ bdry_variablenames = (:EPjump_lhs, :EPjump_rhs)
+ # entropy production due to jumps at interfaces
+ )
+end
diff --git a/src/projects/EulerEq/types.jl b/src/projects/EulerEq/types.jl
new file mode 100644
index 0000000000000000000000000000000000000000..171e9d91e72303b2a8ed23bccc5ec5cb6fb3353c
--- /dev/null
+++ b/src/projects/EulerEq/types.jl
@@ -0,0 +1,19 @@
+struct EulerEquation{T_EoS<:EquationOfState.AbstractEquationOfState} <: AbstractEquation
+ eos::T_EoS
+end
+
+
+mutable struct RHS{T_Equation <:EulerEquation,
+ T_RSolver <:AbstractRiemannSolver,
+ T_HRSC <:Maybe{HRSC.AbstractHRSC},
+ T_TCI <:Maybe{TCI.AbstractTCI},
+ T_LDG_RSolver<:Maybe{AbstractRiemannSolver},
+ T_AV_RSolver<:Maybe{AbstractRiemannSolver}}
+
+ equation::T_Equation
+ rsolver::T_RSolver
+ hrsc::T_HRSC
+ tci::T_TCI
+ ldg_rsolver::T_LDG_RSolver
+ av_rsolver::T_AV_RSolver
+end
diff --git a/src/projects/GRHD/GRHD.jl b/src/projects/GRHD/GRHD.jl
index 765fa8c3c1ec5654f8b34bdb4428dbce4a4e0b86..903c7005f8b2bd64bc1c1e59b9f40749310c862f 100644
--- a/src/projects/GRHD/GRHD.jl
+++ b/src/projects/GRHD/GRHD.jl
@@ -6,18 +6,74 @@ using dg1d
using dg1d.EquationOfState
using Interpolations
using LinearAlgebra
-using Match
using Plots
using OrderedCollections
+include("types.jl")
include("equation.jl")
include("rhs.jl")
include("initialdata.jl")
include("callbacks.jl")
include("plots.jl")
include("boundaryconditions.jl")
-include("main.jl")
+include("setup.jl")
+
+
+#######################################################################
+# dg1d project interface definitions #
+#######################################################################
+
+
+dg1d.project_type(::Val{:GRHD}) = RHS
+
+# TODO Remove bdryconds
+dg1d.rhs!(env::Environment, p::RHS, bdryconds) = rhs!(env, p, p.hrsc, bdryconds...)
+
+dg1d.timestep!(env::Environment, p::RHS) = timestep(env, p, p.hrsc)
+
+dg1d.@parameters GRHD begin
+
+ """
+ Multiplicative factor to determine restmass density for the artificial atmosphere
+ relative to the maximum of restmass density over the grid.
+ """
+ atm_factor = 1e-10
+
+ """
+ Multiplicative factor to determine the cutoff restmass density below which
+ we impose artificial atmosphere.
+ """
+ atm_threshold_factor = 100.0
+
+ """
+ Available options:
+ - `periodic`
+ - `from_initia_data`
+ """
+ bc = "periodic"
+
+ """
+ Available options:
+ - `tov`
+ - `bondi_accretion`
+ - `atmosphere`
+ """
+ id = "atmosphere"
+
+ """
+ Path to an .h5 file containing initial data. Supported initial data types:
+ - `tov`
+ - `bondi_accretion`
+ """
+ id_filename = ""
+
+ """
+ Restmass density for `atmosphere` initial data.
+ """
+ id_atmosphere_density = 1e-11
+
+end
end # module
diff --git a/src/projects/GRHD/boundaryconditions.jl b/src/projects/GRHD/boundaryconditions.jl
index 516f30955d00d0c675db1d486bbf2436f91ed525..1d3872258f5f88fdb199f662ea7acb7bef77aa0f 100644
--- a/src/projects/GRHD/boundaryconditions.jl
+++ b/src/projects/GRHD/boundaryconditions.jl
@@ -1,32 +1,40 @@
-function make_BoundaryConditions(rhs::RHS, prms)
- type = query(String, prms, "type")
-
- xrange = extrema(rhs.mesh.x)
- broadcast_volume!(scaledcons2cons, rhs.equation, rhs.cache)
-
- lhs_bc, rhs_bc = @match type begin
- "periodic" => (PeriodicBC(), PeriodicBC())
- "from initial data" => begin
- @unpack sD, sSr, stau = get_dynamic_variables(rhs.cache)
- @unpack D, Sr, tau, A, B, r, rho, vr, eps, p = get_static_variables(rhs.cache)
- lhs_state = (sD[1], sSr[1], stau[1], p[1], rho[1], vr[1], eps[1], A[1], B[1], r[1])
- rhs_state = (sD[end], sSr[end], stau[end], p[end], rho[end], vr[end], eps[end], A[end], B[end], r[end])
- (DirichletBC(lhs_state), DirichletBC(rhs_state))
- end
- _ => error("Boundary condition type '$type' not implemented!")
+function make_BoundaryConditions(env, equation, rsolver, prms)
+ @unpack bc = prms["GRHD"]
+ @unpack cache = env
+
+ broadcast_volume!(scaledcons2cons, equation, cache)
+
+ lhs_bc, rhs_bc = if bc == "periodic"
+ PeriodicBC(), PeriodicBC()
+ elseif bc == "from_initial_data"
+ @unpack sD, sSr, stau = get_dynamic_variables(cache)
+ @unpack D, Sr, tau, A, B, r, rho, vr, eps, p = get_static_variables(cache)
+ lhs_state = (sD[1], sSr[1], stau[1], p[1], rho[1], vr[1], eps[1],
+ A[1], B[1], r[1])
+ rhs_state = (sD[end], sSr[end], stau[end], p[end], rho[end], vr[end], eps[end],
+ A[end], B[end], r[end])
+ DirichletBC(lhs_state), DirichletBC(rhs_state)
+ else
+ error("GRHD: Boundary condition '$type' not implemented!")
end
- return BoundaryConditions(lhs_bc, rhs_bc, rhs.rsolver)
+ return BoundaryConditions(lhs_bc, rhs_bc, rsolver)
end
-function make_BoundaryConditions_LDG(ldg_rsolver, prms)
+function make_BoundaryConditions_LDG(env, equation, ldg_rsolver, prms)
+
+ @unpack bc = prms["GRHD"]
+
isnothing(ldg_rsolver) && return nothing
- type = query(String, prms, "type")
- lhs_bc, rhs_bc = if type == "periodic"
+
+ lhs_bc, rhs_bc = if bc == "periodic"
PeriodicBC(), PeriodicBC()
- else
+ elseif bc == "from_initial_data"
OutflowBC(), OutflowBC()
+ else
+ error("GRHD: Boundary condition '$type' not implemented!")
end
+
return BoundaryConditions(lhs_bc, rhs_bc, ldg_rsolver)
end
diff --git a/src/projects/GRHD/callbacks.jl b/src/projects/GRHD/callbacks.jl
index 23af3017509151f1a8859751f145adaf15f5c781..5ae5fc65ad804921f4c8e13fd6ac9bb0c1281169 100644
--- a/src/projects/GRHD/callbacks.jl
+++ b/src/projects/GRHD/callbacks.jl
@@ -1,11 +1,11 @@
-function make_callback(rhs::RHS, bdryconds::BoundaryConditions)
- cbfn_equation(u, t) = callback_equation(u, t, rhs, isperiodic(bdryconds), rhs.equation)
+function make_callback(env, rhs::RHS, bdryconds::BoundaryConditions)
+ cbfn_equation(u, t) = callback_equation(u, t, env, rhs, isperiodic(bdryconds), rhs.equation)
cb_equation = FunctionCallback(cbfn_equation,
CallbackTiming(every_iteration=1,every_dt=0))
- cbfn_tci(u, t) = callback_tci(u, t, rhs, isperiodic(bdryconds), rhs.tci)
+ cbfn_tci(u, t) = callback_tci(u, t, env, rhs, isperiodic(bdryconds), rhs.tci)
cb_tci = FunctionCallback(cbfn_tci,
CallbackTiming(every_iteration=1, every_dt=0))
- cbfn_hrsc(u, t) = callback_hrsc(u, t, rhs, isperiodic(bdryconds), rhs.hrsc)
+ cbfn_hrsc(u, t) = callback_hrsc(u, t, env, rhs, isperiodic(bdryconds), rhs.hrsc)
cb_hrsc = FunctionCallback(cbfn_hrsc,
CallbackTiming(every_iteration=1,every_dt=0))
callbackset = CallbackSet(cb_equation, cb_tci, cb_hrsc)
@@ -18,8 +18,8 @@ end
#######################################################################
-function callback_equation(state_u, state_t, rhs, isperiodic, eq)
- @unpack cache = rhs
+function callback_equation(state_u, state_t, env, rhs, isperiodic, eq)
+ @unpack cache = env
wrap_dynamic_variables!(cache, state_u)
@unpack t, tm1 = get_global_variables(cache)
tm1 .= t
@@ -33,17 +33,18 @@ end
#######################################################################
-callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::Nothing) = nothing
-callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::AbstractReconstruction) = nothing
+callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::Nothing) = nothing
+callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::AbstractReconstruction) = nothing
-function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::Maybe{HRSC.AbstractHRSC})
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::Maybe{HRSC.AbstractHRSC})
TODO(callback_hrsc)
end
-function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.AbstractArtificialViscosity)
- @unpack cache, mesh, equation = rhs
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::HRSC.AbstractArtificialViscosity)
+ @unpack cache, mesh = env
+ @unpack equation = rhs
@unpack max_v = get_static_variables(cache)
@unpack mu, cellmax_v = get_cell_variables(cache)
@unpack sD = get_dynamic_variables(cache)
@@ -60,20 +61,20 @@ function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.AbstractArt
end
-function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.TCIViscosity)
- @unpack cache, mesh = rhs
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::HRSC.TCIViscosity)
+ @unpack cache, mesh = env
@unpack mu, flag = get_cell_variables(cache)
HRSC.compute_viscosity!(mu, flag, hrsc)
end
-function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.SmoothedArtificialViscosity)
- callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc.av)
- @unpack cache = rhs
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::HRSC.SmoothedArtificialViscosity)
+ callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc.av)
+ @unpack cache, mesh = env
@unpack mu = get_cell_variables(cache)
indices = findall(muk -> muk > 0.0, mu)
tmp = zeros(size(mu))
- _, K = layout(rhs.mesh)
+ _, K = layout(mesh)
# artificially broaden viscosity
for k in indices
k_lhs, k_rhs = periodic_index(k-1, K), periodic_index(k+1, K)
@@ -88,9 +89,9 @@ function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.SmoothedArt
end
-function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.EntropyViscosity)
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::HRSC.EntropyViscosity)
TODO(callback_hrsc)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack rho = get_dynamic_variables(cache)
@unpack S, Sm1, flx_S, flx_Sm1 = get_static_variables(cache)
@unpack mu, cellmax_v = get_cell_variables(cache)
@@ -107,19 +108,19 @@ end
#######################################################################
-function callback_tci(state_u, state_t, rhs, isperiodic, tci::Nothing) end
+function callback_tci(state_u, state_t, env, rhs, isperiodic, tci::Nothing) end
-function callback_tci(state_u, state_t, rhs, isperiodic, tci::TCI.Threshold)
+function callback_tci(state_u, state_t, env, rhs, isperiodic, tci::TCI.Threshold)
error("Threshold TCI must be combined with an indicator that is supposed to be zero for " *
"for smooth solutions, but not for shocks. TODO Implement an indicator variable.")
end
-function callback_tci(state_u, state_t, rhs, isperiodic, tci::TCI.EntropyProduction)
+function callback_tci(state_u, state_t, env, rhs, isperiodic, tci::TCI.EntropyProduction)
@warn "Figure out entropy fluxes for special relativistic euler equation!"
TODO(callback_tci)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack E, Em1, flx_E, flx_Em1, EP = get_static_variables(cache)
@unpack flag = get_cell_variables(cache)
@unpack EPjump_lhs, EPjump_rhs = get_bdry_variables(cache)
@@ -130,11 +131,11 @@ function callback_tci(state_u, state_t, rhs, isperiodic, tci::TCI.EntropyProduct
end
-function callback_tci(state_u, state_t, rhs, isperiodic,
- tci::Union{TCI.Diffusion, TCI.Minmod, TCI.Threshold,
+function callback_tci(state_u, state_t, env, rhs, isperiodic,
+ tci::Union{TCI.Diffusion, TCI.Minmod, TCI.Threshold,
TCI.ModalDecayAverage, TCI.ModalDecayHighest})
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack sD, sSr, stau = get_dynamic_variables(cache)
@unpack flag, sD_flag, sSr_flag, stau_flag = get_cell_variables(cache)
TCI.compute_indicator!(sD_flag, sD, tci)
diff --git a/src/projects/GRHD/equation.jl b/src/projects/GRHD/equation.jl
index 77b64dc1e180b6f04d8666b1e9ea97c3b9d1b135..8fff0cd6cf9b5bc65df145218a69f59530b94097 100644
--- a/src/projects/GRHD/equation.jl
+++ b/src/projects/GRHD/equation.jl
@@ -1,46 +1,4 @@
-#######################################################################
-# Types #
-#######################################################################
-
-
-Base.@kwdef struct Atmosphere
- factor::Float64
- threshold::Float64
-end
-
-
-function parse_parameters(::Type{Atmosphere}, prms)
- factor = query(Float64, prms, "factor", default=1e-7)
- threshold = query(Float64, prms, "threshold", default=100.0)
-
- @assert factor > 0
- @assert threshold > 0
-
- return (; factor, threshold)
-end
-
-
-Atmosphere(prms::AbstractDict) = Atmosphere(; parse_parameters(Atmosphere, prms)...)
-
-
-"""
- struct Equation <: AbstractEquationOfState
-
-Represents the equations of General Relativistic Hydrodynamics (GRHD).
-"""
-Base.@kwdef struct Equation{T_EoS<:AbstractEquationOfState} <: AbstractEquation
- eos::T_EoS
- atmosphere::Atmosphere
- cons2prim_newtonsolver::NewtonSolver
- cons2prim_bisection::Bisection
-end
-
-
-#######################################################################
-# Implementation #
-#######################################################################
-
-
+### TODO Move include to GRHD.jl or merge with this file?
include("cons2prim_implementations.jl")
diff --git a/src/projects/GRHD/initialdata.jl b/src/projects/GRHD/initialdata.jl
index 8cc2a6ca088b8d8dc72a2064b51ebd5452c648b4..248c5219bf0248010010caa8f72d56a09c0698e6 100644
--- a/src/projects/GRHD/initialdata.jl
+++ b/src/projects/GRHD/initialdata.jl
@@ -1,16 +1,14 @@
-struct InitialData{T} end
-
-
-function initialdata!(rhs::RHS, parameters)
- type = query(parameters, "type")
- type = Symbol(lowercase(replace(type, " "=>"_")))
- id_type = InitialData{type}()
-
- initialdata_equation!(id_type, rhs, parameters)
+function initialdata!(env, rhs::RHS, prms)
+ @unpack id = prms
+ let
+ @unpack r = get_static_variables(env.cache)
+ @. r = env.mesh.x[:]
+ end
+ initialdata_equation!(Val(Symbol(id)), env, rhs, prms)
return
- # initialdata_entropy_variables!(id_type, rhs, parameters)
- initialdata_hrsc!(id_type, rhs, rhs.hrsc, parameters)
- initialdata_tci!(id_type, rhs, rhs.tci, parameters)
+ # initialdata_entropy_variables!(id, rhs, prms)
+ initialdata_hrsc!(Val(Symbol(id)), env, rhs, rhs.hrsc, prms)
+ initialdata_tci!(Val(Symbol(id)), env, rhs, rhs.tci, prms)
end
@@ -19,24 +17,29 @@ end
#######################################################################
-function initialdata_equation!(id::InitialData, rhs, parameters)
- error_not_implemented("initialdata_equation!", typeof(id))
+function initialdata_equation!(::Val{T}, env, rhs, parameters) where T
+ error_not_implemented("initialdata_equation!", string(T))
end
-function initialdata_equation!(::InitialData{:bondi_accretion}, rhs, parameters)
+function initialdata_equation!(::Val{:bondi_accretion}, env, rhs, prms)
TODO(initialdata_equation!)
end
-function initialdata_equation!(::InitialData{:tov}, rhs, parameters)
+function initialdata_equation!(::Val{:tov}, env, rhs, prms)
- filename = query(parameters, "filename")
+ filename = prms["id_filename"]
+ if length(filename) == 0
+ error("GRHD: Missing filename for TOV output, include with GRHD.id_filename = 'tov'")
+ end
+ filename = joinpath(get_output_dir(), filename)
# read TOV result from file
- # only contains solution in star interior
- output_dir = dg1d.get_output_dir()
- data = TOV.read_data(joinpath(output_dir, "TOV.h5"))
+ # should only contain solution in star interior
+ # output_dir = dg1d.get_output_dir()
+ # data = TOV.read_data(joinpath(output_dir, "TOV.h5"))
+ data = TOV.read_data(filename)
Ï•int = data.Ï•
Rint = data.R
rint = data.r
@@ -49,8 +52,9 @@ function initialdata_equation!(::InitialData{:tov}, rhs, parameters)
# interpolate data onto our grid
# we don't assume a special grid layout for the loaded data
- @unpack mesh = rhs
+ @unpack cache, mesh = env
@unpack Npts = mesh.element
+
intrp_R = vec(mesh.x) # DG coordinates
NN = length(intrp_R)
@@ -142,7 +146,7 @@ function initialdata_equation!(::InitialData{:tov}, rhs, parameters)
@assert all(map(X->all(isfinite.(X)), [intrp_A, intrp_B, intrp_A_r, intrp_B_r])) == true ""
"Initial data interpolation failed!"
- @unpack A, B, A_r, B_r = get_static_variables(rhs.cache)
+ @unpack A, B, A_r, B_r = get_static_variables(cache)
A .= intrp_A
B .= intrp_B
A_r .= intrp_A_r
@@ -150,17 +154,17 @@ function initialdata_equation!(::InitialData{:tov}, rhs, parameters)
### setup hydrodynamic variables
@unpack equation = rhs
- @unpack atmosphere = equation
+ # @unpack atmosphere = equation
eos = Polytrope(K=K, Gamma=Γ)
# apply atmosphere treatment
pmax = maximum(pint)
Ïmax = eos(Density, Pressure, pmax)
- Ïatm = atmosphere.factor * Ïmax
+ Ïatm = equation.atm_factor * Ïmax
patm = eos(Pressure, Density, Ïatm)
pint[@.(pint < 0)] .= 0.0 # just to be sure, they will be replace with patm below anyways
Ïint = @. eos(Density, Pressure, pint)
- idxatm_treatment = @. Ïint < Ïatm * atmosphere.threshold
+ idxatm_treatment = @. Ïint < Ïatm * equation.atm_threshold
pint[idxatm_treatment] .= patm
pint[idxatm_treatment] .= Ïatm
vint = zeros(size(Rint))
@@ -190,17 +194,17 @@ function initialdata_equation!(::InitialData{:tov}, rhs, parameters)
intrp_ϵ = mirror(intrp_ϵ, exclude_origin)
end
- @unpack p, vr, rho, eps = get_static_variables(rhs.cache)
- @unpack sD, sSr, stau = get_dynamic_variables(rhs.cache)
- @unpack D, Sr, tau, = get_static_variables(rhs.cache)
+ @unpack p, vr, rho, eps = get_static_variables(cache)
+ @unpack sD, sSr, stau = get_dynamic_variables(cache)
+ @unpack D, Sr, tau, = get_static_variables(cache)
p .= intrp_p
rho .= intrp_Ï
vr .= intrp_vr
eps .= intrp_ϵ
- broadcast_volume!(prim2cons, equation, rhs.cache)
- broadcast_volume!(cons2scaledcons, equation, rhs.cache)
+ broadcast_volume!(prim2cons, equation, cache)
+ broadcast_volume!(cons2scaledcons, equation, cache)
@assert all(map(X->all(isfinite.(X)),
[p, vr, rho, eps, D, Sr, tau, sD, sSr, stau])) == true "Initial data interpolation failed!"
@@ -243,9 +247,10 @@ function reflect(A, exclude_bdry::Bool)
end
-function initialdata_equation!(::InitialData{:atmosphere}, rhs, parameters)
+function initialdata_equation!(::Val{:atmosphere}, env, rhs, prms)
- @unpack mesh, equation = rhs
+ @unpack cache, mesh = env
+ @unpack equation = rhs
@unpack Npts = mesh.element
@unpack eos = equation
@unpack A, B, A_r, B_r = get_static_variables(rhs.cache)
@@ -258,7 +263,8 @@ function initialdata_equation!(::InitialData{:atmosphere}, rhs, parameters)
B_r .= zeros(size(R))
# low density atmosphere
- rho_atm = query(parameters, "density", default=1e-11)
+ # rho_atm = query(prms, "density", default=1e-11)
+ rho_atm = prms["id_atmosphere_density"]
@assert rho_atm > 0.0
@assert eos isa Polytrope "Require model for cold atmosphere"
@@ -285,7 +291,7 @@ end
#######################################################################
-function initialdata_entropy_variables!(::InitialData, rhs, parameters)
+function initialdata_entropy_variables!(::Val, rhs, prms)
TODO(initialdata_entropy_variables!)
end
@@ -293,7 +299,7 @@ end
# InitialData HRSC #
#######################################################################
-function initialdata_hrsc!(::InitialData, rhs, hrsc, parameters)
+function initialdata_hrsc!(::Val, rhs, hrsc, prms)
TODO(initialdata_hrsc!)
end
@@ -301,6 +307,6 @@ end
# InitialData TCI #
#######################################################################
-function initialdata_tci!(::InitialData, rhs, tci, prms)
+function initialdata_tci!(::Val, rhs, tci, prms)
TODO(initialdata_tci!)
end
diff --git a/src/projects/GRHD/main.jl b/src/projects/GRHD/main.jl
deleted file mode 100644
index 48490526626735c4015890ca5c8d28e2afb0f97f..0000000000000000000000000000000000000000
--- a/src/projects/GRHD/main.jl
+++ /dev/null
@@ -1,156 +0,0 @@
-"""
- main(parameters_filename)
-
-Run program specified by TOML file `parameters_filename`.
-
-The following sections in the parameter file are recognized:
-- [`Mesh`](@ref)
-- [`EquationOfState`](@ref)
-- [`Initialdata`](@ref)
-- [`Timestep`](@ref)
-- [`HRSC`](@ref)
-- [`Output`](@ref)
-"""
-function main(parameters_filename=""; override::AbstractDict=Dict())
-
- override_prms = FlatDict(override)
-
- parameters_filename = abspath(parameters_filename)
- outputdir = mk_outputdir(parameters_filename)
- cp(parameters_filename, joinpath(outputdir, basename(parameters_filename)))
-
- # read parameters
- if !isfile(parameters_filename)
- error("Cannot locate parameter file '$parameters_filename'")
- end
- prms = FlatDict(TOML.parsefile(parameters_filename))
- prms = merge(prms, override_prms)
- ordered_prms = OrderedDict(prms)
- display(sort(ordered_prms))
-
- default_eos_prms = Dict("type"=>"IdealGas", "gamma"=>5/3)
- eos_prms = trim_headers(prms, "EquationOfState")
- amend!(eos_prms, default_eos_prms)
- diff = diffkeys(eos_prms, default_eos_prms)
- @assert length(diff) == 0 "Section 'Plot': Unknown options '$diff' encountered!"
- eos = EquationOfState.make_EquationOfState(eos_prms)
- atm = Atmosphere(; parse_parameters(Atmosphere, trim_headers(prms, "Atmosphere"))...)
- cons2prim_newtonsolver = NewtonSolver(approx_derivative=true)
- cons2prim_bisection = Bisection()
-
- # setup RHS
- mesh = make_Mesh(trim_headers(prms, "Mesh"))
- cache = Cache(mesh)
- equation = Equation(eos, atm, cons2prim_newtonsolver, cons2prim_bisection)
- tci = TCI.make_TCI(trim_headers(prms, "TCI"), mesh)
- hrsc = HRSC.make_HRSC(trim_headers(prms, "HRSC"), mesh)
- rsolver = ApproxRiemannSolver(flux, speed, equation)
- ldg_rsolver, av_rsolver = if hrsc isa AbstractArtificialViscosity
- ApproxRiemannSolver(ldg_flux, ldg_speed, equation),
- ApproxRiemannSolver(av_flux, av_flux, equation)
- else
- nothing, nothing
- end
- rhs = RHS(cache=cache, mesh=mesh, equation=equation, tci=tci, hrsc=hrsc,
- rsolver=rsolver, ldg_rsolver=ldg_rsolver, av_rsolver=av_rsolver)
-
- _register_variables!(rhs)
- register_variables!(cache, rsolver)
- !isnothing(ldg_rsolver) && register_variables!(cache, ldg_rsolver)
- !isnothing(av_rsolver) && register_variables!(cache, av_rsolver, dont_register=true)
- display(cache)
-
- initialdata!(rhs, trim_headers(prms, "Initialdata"))
-
- default_boundaryconditions_prms = Dict("type"=>"periodic")
- boundaryconditions_prms = trim_headers(prms, "BoundaryConditions")
- amend!(boundaryconditions_prms, default_boundaryconditions_prms)
- diff = diffkeys(boundaryconditions_prms, default_boundaryconditions_prms)
- @assert length(diff) == 0 "Section 'Plot': Unknown options '$diff' encountered!"
- bdryconds = make_BoundaryConditions(rhs, boundaryconditions_prms)
- ldg_bdryconds, av_bdryconds = if hrsc isa AbstractArtificialViscosity
- make_BoundaryConditions_LDG(ldg_rsolver, boundaryconditions_prms),
- make_BoundaryConditions_LDG(av_rsolver, boundaryconditions_prms)
- else
- nothing, nothing
- end
-
- rhs_fn!(du, u, t) = rhs!(du, u, t, rhs, bdryconds, ldg_bdryconds, av_bdryconds)
- timestep_fn(u, t, n) = timestep(u, t, rhs)
-
- ### Callbacks
- rhscb = make_callback(rhs, bdryconds)
-
- # Plotting
- default_plot_prms = Dict("every_iteration"=>0, "every_dt"=>0.0, "pause"=>false)
- plot_prms = trim_headers(prms, "Plot")
- dg1d.amend!(plot_prms, default_plot_prms)
- diff = diffkeys(plot_prms, default_plot_prms)
- @assert length(diff) == 0 "Section 'Plot': Unknown options '$diff' encountered!"
- theme(:dark, markerstrokewidth=0.1, marker=:circle)
- plotfn(u,t) = plot(rhs)
- plotcb = PlotCallback(plotfn,
- CallbackTiming(every_dt=plot_prms["every_dt"],
- every_iteration=plot_prms["every_iteration"]),
- pause=plot_prms["pause"])
- default_output_prms = Dict("every_iteration"=>0,
- "every_dt"=>0.5,
- "variables"=>["rho", "q", "E"],
- "filename"=>joinpath(outputdir, "output"))
- output_prms = trim_headers(prms, "Output")
- amend!(output_prms, default_output_prms)
- diff = diffkeys(output_prms, default_output_prms)
- @assert length(diff) == 0 "Section 'Output': Unknown options '$diff' encountered!"
- savecb = SaveCallback(Symbol.(output_prms["variables"]), cache, mesh,
- output_prms["filename"],
- CallbackTiming(every_dt=output_prms["every_dt"],
- every_iteration=output_prms["every_iteration"])
- )
-
- # artificially smooth initial data
- smoother_called = false
- function id_smoother_fn(u, t)
- if !smoother_called
- println("Smoothing initial data...")
- @unpack sD, sSr, stau = get_dynamic_variables(cache)
- @unpack flag = get_cell_variables(cache)
- id_hrsc = HRSC.BernsteinReconstruction(mesh)
- HRSC.reconstruct!(sD, flag, id_hrsc, isperiodic=isperiodic(bdryconds))
- HRSC.reconstruct!(sSr, flag, id_hrsc, isperiodic=isperiodic(bdryconds))
- HRSC.reconstruct!(stau, flag, id_hrsc, isperiodic=isperiodic(bdryconds))
- smoother_called = true
- broadcast_volume!(cons2prim, equation, cache)
- end
- end
- id_smoother_cb = FunctionCallback(id_smoother_fn,
- CallbackTiming(every_iteration=1))
-
- # Combine callbacks
- callbacks = CallbackSet(rhscb.callbacks..., id_smoother_cb, savecb, plotcb)
-
- aligned_times = get(output_prms, "aligned_times", [])
- timealigned_savecb, timestep_fn = if length(aligned_times) > 0
- atscb = TimeAlignedSaveCallback(aligned_times, timestep_fn, savecb)
- atscb, atscb.aligned_timestep_fn
- else
- nothing, timestep_fn
- end
-
- !isnothing(timealigned_savecb) && push!(callbacks, timealigned_savecb)
-
- default_evolve_prms = Dict("algorithm"=>"LSERK4", "tend"=>1.0)
- evolve_prms = trim_headers(prms, "Evolution")
- amend!(evolve_prms, default_evolve_prms)
- diff = diffkeys(evolve_prms, default_evolve_prms)
- @assert length(diff) == 0 "Section 'Evolution': Unknown options '$diff' encountered!"
-
- @unpack sD, sSr, stau = get_dynamic_variables(cache)
- state_u = vcat(sD, sSr, stau)
-
- state_u = dg1d.evolve_ERK(rhs_fn!, state_u, timestep_fn,
- (0.0, Float64(evolve_prms["tend"])),
- Symbol(evolve_prms["algorithm"]);
- callback_fullstep=callbacks)
-
- return
-end
diff --git a/src/projects/GRHD/parameters.jl b/src/projects/GRHD/parameters.jl
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/projects/GRHD/plots.jl b/src/projects/GRHD/plots.jl
index c6977a732712df4c35ca4d66040960ee6be7b7e3..aaee1f000ed67bfc2516abc6a73d2f138548dde7 100644
--- a/src/projects/GRHD/plots.jl
+++ b/src/projects/GRHD/plots.jl
@@ -1,8 +1,8 @@
-function plot(rhs::RHS)
- plt_eq = plot_equation(rhs)
- plt_rhs = plot_rhs(rhs)
- plt_hrsc = plot_hrsc(rhs, rhs.hrsc)
- plt_tci = plot_tci(rhs, rhs.tci)
+function plot(env, rhs::RHS)
+ plt_eq = plot_equation(env, rhs)
+ plt_rhs = plot_rhs(env, rhs)
+ plt_hrsc = plot_hrsc(env, rhs, rhs.hrsc)
+ plt_tci = plot_tci(env, rhs, rhs.tci)
list_sub_plts = [ plt_eq, #=plt_rhs,=# plt_hrsc, plt_tci ]
filter!(!isnothing, list_sub_plts)
@@ -13,7 +13,7 @@ function plot(rhs::RHS)
sub_plt = Plots.plot(list_sub_plts..., layout=grid(layout..., heights=plt_heights))
# combine plots
- @unpack cache = rhs
+ @unpack cache = env
@unpack t = get_global_variables(cache)
plt = Plots.plot(sub_plt, layout=layout, size=(1400,1000), plot_title="t=$(t[1])",
left_margin=10Plots.mm)
@@ -26,9 +26,9 @@ end
#######################################################################
-function plot_equation(rhs)
+function plot_equation(env, rhs)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack sD, sSr, stau = get_dynamic_variables(cache)
@unpack D, Sr, tau,
rho, p, eps, vr, A, B = get_static_variables(cache)
@@ -75,13 +75,13 @@ end
#######################################################################
-function plot_hrsc(rhs, hrsc::Nothing) end
-function plot_hrsc(rhs, hrsc::AbstractReconstruction) end
+function plot_hrsc(env, rhs, hrsc::Nothing) end
+function plot_hrsc(env, rhs, hrsc::AbstractReconstruction) end
-function plot_hrsc(rhs, hrsc::HRSC.AbstractArtificialViscosity)
+function plot_hrsc(env, rhs, hrsc::HRSC.AbstractArtificialViscosity)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack mu = get_cell_variables(cache)
@unpack x = mesh
@@ -103,11 +103,11 @@ function plot_hrsc(rhs, hrsc::HRSC.AbstractArtificialViscosity)
end
-function plot_hrsc(rhs, hrsc::HRSC.SmoothedArtificialViscosity)
+function plot_hrsc(env, rhs, hrsc::HRSC.SmoothedArtificialViscosity)
- plt_mu = plot_hrsc(rhs, hrsc.av)
+ plt_mu = plot_hrsc(env, rhs, hrsc.av)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack smoothed_mu = get_static_variables(cache)
@unpack x = mesh
@@ -137,12 +137,12 @@ end
#######################################################################
-function plot_tci(rhs, tci::Nothing) end
+function plot_tci(env, rhs, tci::Nothing) end
-function plot_tci(rhs, tci::TCI.AbstractTCI)
+function plot_tci(env, rhs, tci::TCI.AbstractTCI)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack flag = get_cell_variables(cache)
@unpack x = mesh
@@ -172,9 +172,9 @@ function plot_tci(rhs, tci::TCI.AbstractTCI)
end
-function plot_tci(rhs, tci::TCI.EntropyProduction)
+function plot_tci(env, rhs, tci::TCI.EntropyProduction)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack EP = get_static_variables(cache)
@unpack flag = get_cell_variables(cache)
@unpack x = mesh
@@ -211,9 +211,9 @@ end
#######################################################################
-function plot_rhs(rhs)
+function plot_rhs(env, rhs)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack rhs_sD, rhs_sSr, rhs_stau = get_rhs_variables(cache)
@unpack t = get_global_variables(cache)
@unpack x = mesh
diff --git a/src/projects/GRHD/rhs.jl b/src/projects/GRHD/rhs.jl
index 457208bb50bb96fc9ce7f3df5a867b9fbfb01620..662095151c8df2ba265c400e83c5e8d4288411c2 100644
--- a/src/projects/GRHD/rhs.jl
+++ b/src/projects/GRHD/rhs.jl
@@ -1,106 +1,11 @@
-Base.@kwdef struct RHS{T_RSolver <:AbstractRiemannSolver,
- T_HRSC <:Maybe{HRSC.AbstractHRSC},
- T_TCI <:Maybe{TCI.AbstractTCI},
- T_LDG_RSolver<:Maybe{AbstractRiemannSolver},
- T_AV_RSolver<:Maybe{AbstractRiemannSolver}}
-
- cache::Cache
- mesh::Mesh
- equation::Equation
- rsolver::T_RSolver
- hrsc::T_HRSC = nothing
- tci::T_TCI = nothing
- ldg_rsolver::T_LDG_RSolver = nothing
- av_rsolver::T_AV_RSolver = nothing
-end
-
-
-#######################################################################
-# Register variables #
-#######################################################################
-
-
-# Julia limitation/bug: Overloading imported functions without their Module prefix shadows their
-# definition in this module without a warning or error. This does not happen for methods
-# defined in Base.
-# Because of that we prefix with _ below.
-# Cf. https://github.com/JuliaLang/julia/issues/15510
-function _register_variables!(rhs::RHS)
- register_variables!(rhs.cache,
- dynamic_variablenames = (:sD, :sSr, :stau,), # 'scaled' conservatives
- rhs_variablenames = (:rhs_sD, :rhs_sSr, :rhs_stau), # rhs of conservatives
- static_variablenames = (:flx_sD, :flx_sSr, :flx_stau, # conservatives' fluxes
- :src_sD, :src_sSr, :src_stau, # sources
- :D, :Sr, :tau, # true hydrodynamic conservatives
- :rho, :vr, :p, :eps, # primitives
- :max_v, :r, # maximum charactersitic speed, grid
- :E, :Em1, # entropy
- :flx_E, :flx_Em1, # entropy flux
- :A, :B, :A_r, :B_r), # metric fields and derivatives
- cell_variablenames = (:cellmax_v,),
- global_variablenames = (:t, :tm1) # time steps
- )
- _register_variables!(rhs, rhs.hrsc)
- _register_variables!(rhs, rhs.tci)
- @unpack r = get_static_variables(rhs.cache)
- @. r = rhs.mesh.x[:]
-end
-
-
-function _register_variables!(rhs, tci_or_hrsc::Nothing) end
-
-
-_register_variables!(rhs, hrsc::AbstractHRSC) =
- error("_register_variables: Not implemented for hrsc of type '$(typeof(hrsc))'")
-
-_register_variables!(rhs, hrsc::HRSC.AbstractReconstruction) = nothing
-
-
-function _register_variables!(rhs, hrsc::HRSC.AbstractArtificialViscosity)
- register_variables!(rhs.cache,
- static_variablenames = (:ldg_sD, :ldg_sSr, :ldg_stau, # local DG variable
- :flx_ldg_sD, :flx_ldg_sSr, :flx_ldg_stau), # local DG flux
- cell_variablenames = (:mu,), # 'one viscosity to rule them all'
- )
-end
-
-
-function _register_variables!(rhs, hrsc::HRSC.SmoothedArtificialViscosity)
- _register_variables!(rhs, hrsc.av)
- register_variables!(rhs.cache,
- static_variablenames = (:smoothed_mu,)
- )
-end
-
-
-function _register_variables!(rhs, tci::TCI.AbstractTCI)
- register_variables!(rhs.cache, cell_variablenames = (:flag,:sD_flag,:sSr_flag,:stau_flag))
-end
-
-
-function _register_variables!(rhs, tci::TCI.EntropyProduction)
- register_variables!(rhs.cache,
- cell_variablenames = (:flag,),
- static_variablenames = (:EP,), # entropy production
- bdry_variablenames = (:EPjump_lhs, :EPjump_rhs)
- # entropy production due to jumps at interfaces
- )
-end
-
-
#######################################################################
# Timestep #
#######################################################################
-function timestep(state_u, t, rhs::RHS)
- wrap_dynamic_variables!(rhs.cache, state_u)
- return timestep(rhs, rhs.hrsc)
-end
-
-
-function timestep(rhs::RHS, hrsc::Maybe{HRSC.AbstractHRSC})
- @unpack cache, mesh, equation = rhs
+function timestep(env, rhs::RHS, hrsc::Maybe{HRSC.AbstractHRSC})
+ @unpack cache, mesh = env
+ @unpack equation = rhs
@unpack max_v = get_static_variables(cache)
broadcast_volume!(cons2prim, equation, cache)
@@ -120,8 +25,9 @@ function timestep(rhs::RHS, hrsc::Maybe{HRSC.AbstractHRSC})
end
-function timestep(rhs::RHS, hrsc::HRSC.AbstractArtificialViscosity)
- @unpack cache, mesh, equation = rhs
+function timestep(env, rhs::RHS, hrsc::HRSC.AbstractArtificialViscosity)
+ @unpack cache, mesh = env
+ @unpack equation = rhs
broadcast_volume!(cons2prim, equation, cache)
broadcast_volume!(speed, equation, cache)
@@ -153,19 +59,11 @@ end
#######################################################################
-function rhs!(state_du, state_u, t, rhs::RHS, bdryconds::BoundaryConditions,
- ldg_bdryconds::Maybe{BoundaryConditions}, av_bdryconds::Maybe{BoundaryConditions})
-
- @unpack cache = rhs
- wrap_dynamic_variables!(cache, state_u)
- wrap_rhs_variables!(cache, state_du)
- rhs!(rhs, t, rhs.hrsc, bdryconds, ldg_bdryconds, av_bdryconds)
-end
-
+function rhs!(env, rhs::RHS, hrsc::Nothing, bdryconds, ldg_bdryconds, av_bdryconds)
-function rhs!(rhs::RHS, t, hrsc::Nothing, bdryconds, ldg_bdryconds, av_bdryconds)
+ @unpack cache, mesh = env
+ @unpack equation, rsolver = rhs
- @unpack cache, mesh, equation, rsolver = rhs
@unpack sD, sSr, stau = get_dynamic_variables(cache)
@unpack flx_sD, flx_sSr, flx_stau = get_static_variables(cache)
@unpack rhs_sD, rhs_sSr, rhs_stau = get_rhs_variables(cache)
@@ -179,7 +77,7 @@ function rhs!(rhs::RHS, t, hrsc::Nothing, bdryconds, ldg_bdryconds, av_bdryconds
broadcast_volume!(flux, equation, cache)
broadcast_volume!(sources, equation, cache)
broadcast_faces!(lax_friedrich_flux, rsolver, cache, mesh)
- broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, t)
+ broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, 0.0)
@unpack p, vr, rho, eps = get_static_variables(cache)
@@ -191,10 +89,12 @@ function rhs!(rhs::RHS, t, hrsc::Nothing, bdryconds, ldg_bdryconds, av_bdryconds
end
-function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractReconstruction,
+function rhs!(env, rhs::RHS, hrsc::HRSC.AbstractReconstruction,
bdryconds, ldg_bdryconds, av_bdryconds)
- @unpack cache, mesh, equation, rsolver = rhs
+ @unpack cache, mesh = env
+ @unpack equation, rsolver = rhs
+
@unpack sD, sSr, stau = get_dynamic_variables(cache)
@unpack flx_sD, flx_sSr, flx_stau = get_static_variables(cache)
@unpack rhs_sD, rhs_sSr, rhs_stau = get_rhs_variables(cache)
@@ -222,11 +122,12 @@ function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractReconstruction,
end
-function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractArtificialViscosity,
+function rhs!(env, rhs::RHS, hrsc::HRSC.AbstractArtificialViscosity,
bdryconds, ldg_bdryconds, av_bdryconds)
- @unpack cache, mesh, equation, rsolver,
- ldg_rsolver, av_rsolver = rhs
+ @unpack cache, mesh = env
+ @unpack equation, rsolver, ldg_rsolver, av_rsolver = rhs
+
@unpack sD, sSr, stau = get_dynamic_variables(cache)
@unpack flx_sD, flx_sSr, flx_stau,
ldg_sD, ldg_sSr, ldg_stau,
@@ -244,7 +145,7 @@ function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractArtificialViscosity,
## solve auxiliary equation: q + ∂x u = 0
broadcast_volume!(ldg_flux, equation, cache)
broadcast_faces!(mean_flux, ldg_rsolver, cache, mesh)
- broadcast_boundaryconditions!(mean_flux, ldg_bdryconds, cache, mesh, t)
+ broadcast_boundaryconditions!(mean_flux, ldg_bdryconds, cache, mesh, 0.0)
compute_rhs_weak_form!(ldg_sD, flx_ldg_sD, lhs_numflx_ldg_sD, rhs_numflx_ldg_sD, mesh)
compute_rhs_weak_form!(ldg_sSr, flx_ldg_sSr, lhs_numflx_ldg_sSr, rhs_numflx_ldg_sSr, mesh)
compute_rhs_weak_form!(ldg_stau, flx_ldg_stau, lhs_numflx_ldg_stau, rhs_numflx_ldg_stau, mesh)
@@ -252,14 +153,14 @@ function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractArtificialViscosity,
## compute rhs of regularized equation: ∂t u + ∂x f + ∂x μ q = 0
broadcast_volume!(av_flux, equation, cache)
broadcast_faces!(mean_flux, av_rsolver, cache, mesh)
- broadcast_boundaryconditions!(mean_flux, av_bdryconds, cache, mesh, t)
+ broadcast_boundaryconditions!(mean_flux, av_bdryconds, cache, mesh, 0.0)
broadcast_volume!(scaledcons2cons, equation, cache)
postprocess_scaledcons2cons!(cache, equation, mesh)
broadcast_volume!(flux, equation, cache)
broadcast_volume!(sources, equation, cache)
broadcast_faces!(lax_friedrich_flux, rsolver, cache, mesh)
- broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, t)
+ broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, 0.0)
@. flx_sD += flx_ldg_sD
@. flx_sSr += flx_ldg_sSr
diff --git a/src/projects/GRHD/setup.jl b/src/projects/GRHD/setup.jl
new file mode 100644
index 0000000000000000000000000000000000000000..f384e068382b9e7f2494e7cbe535374d677233db
--- /dev/null
+++ b/src/projects/GRHD/setup.jl
@@ -0,0 +1,161 @@
+function RHS(env::Environment, prms)
+
+ # construct RHS
+ eos = EquationOfState.make_EquationOfState(prms["EquationOfState"])
+ equation = Equation(eos, prms["GRHD"])
+ # TODO Update TCI, HRSC to new interface
+ tci = TCI.make_TCI(env.mesh, prms["TCI"])
+ hrsc = HRSC.make_HRSC(env.mesh, prms["HRSC"])
+ # tci = TCI.make_TCI(prms["TCI"], env.mesh)
+ # hrsc = HRSC.make_HRSC(prms["HRSC"], env.mesh)
+ rsolver = ApproxRiemannSolver(flux, speed, equation)
+ ldg_rsolver, av_rsolver = if hrsc isa AbstractArtificialViscosity
+ ApproxRiemannSolver(ldg_flux, ldg_speed, equation),
+ ApproxRiemannSolver(av_flux, av_flux, equation)
+ else
+ nothing, nothing
+ end
+
+ rhs = RHS(equation, hrsc, tci, rsolver, ldg_rsolver, av_rsolver)
+
+ # register variables
+ # TODO add a ::Nothing overload for register_variables!
+ # TODO Somehow replace _register_variables! with register_variables!
+ @unpack cache = env
+ _register_variables!(cache, rhs)
+ register_variables!(cache, rsolver)
+ !isnothing(ldg_rsolver) && register_variables!(cache, ldg_rsolver)
+ !isnothing(av_rsolver) && register_variables!(cache, av_rsolver, dont_register=true)
+ display(cache)
+
+ # setup boundary conditions
+ bdryconds = make_BoundaryConditions(env, equation, rsolver, prms)
+ ldg_bdryconds, av_bdryconds = if hrsc isa AbstractArtificialViscosity
+ make_BoundaryConditions_LDG(env, equation, ldg_rsolver, prms),
+ make_BoundaryConditions_LDG(env, equation, av_rsolver, prms)
+ else
+ nothing, nothing
+ end
+
+ # TODO Bdry conditions should be moved to project and should only be just a 'lightweight'
+ # descriptor. Furthermore, they should also register variables in cache themselves.
+ # Note: Atm, only the 'from_initial_data' option would require a place to store the data
+ # somewhere. For periodic BCs we don't need extra storage at all.
+ # But what other types of boundary conditions will be needed in the future?
+
+ # setup initial data
+ initialdata!(env, rhs, prms["GRHD"])
+
+ # setup callbacks
+ rhscb = make_callback(env, rhs, bdryconds)
+
+ # artificially smooth initial data
+ # TODO Move this into make_callback
+ smoother_called = false
+ function id_smoother_fn(u, t)
+ if !smoother_called
+ println("Smoothing initial data...")
+ @unpack sD, sSr, stau = get_dynamic_variables(cache)
+ @unpack flag = get_cell_variables(cache)
+ id_hrsc = HRSC.BernsteinReconstruction(env.mesh)
+ HRSC.reconstruct!(sD, flag, id_hrsc, isperiodic=isperiodic(bdryconds))
+ HRSC.reconstruct!(sSr, flag, id_hrsc, isperiodic=isperiodic(bdryconds))
+ HRSC.reconstruct!(stau, flag, id_hrsc, isperiodic=isperiodic(bdryconds))
+ smoother_called = true
+ broadcast_volume!(cons2prim, equation, cache)
+ end
+ end
+ id_smoother_cb = FunctionCallback(id_smoother_fn,
+ CallbackTiming(every_iteration=1))
+
+ append!(env.callbacks, CallbackSet(rhscb.callbacks..., id_smoother_cb))
+
+ return rhs, (bdryconds, ldg_bdryconds, av_bdryconds)
+end
+
+
+function Equation(eos, prms)
+ @unpack atm_factor, atm_threshold_factor = prms
+ if atm_factor < 0
+ error("GRHD: atm_factor must be positive, found 'atm_factor = $atm_factor'")
+ end
+ if atm_threshold_factor < 0
+ error("GRHD: atm_threshold_factor must be positive, found 'atm_threshold_factor = $atm_threshold_factor'")
+ end
+ cons2prim_newtonsolver = NewtonSolver(approx_derivative=true)
+ cons2prim_bisection = Bisection()
+ return Equation(eos, atm_factor, atm_threshold_factor,
+ cons2prim_newtonsolver, cons2prim_bisection)
+end
+
+
+#######################################################################
+# Register variables #
+#######################################################################
+
+
+# Julia limitation/bug: Overloading imported functions without their Module prefix shadows their
+# definition in this module without a warning or error. This does not happen for methods
+# defined in Base.
+# Because of that we prefix with _ below.
+# Cf. https://github.com/JuliaLang/julia/issues/15510
+function _register_variables!(cache, rhs::RHS)
+ register_variables!(cache,
+ dynamic_variablenames = (:sD, :sSr, :stau,), # 'scaled' conservatives
+ rhs_variablenames = (:rhs_sD, :rhs_sSr, :rhs_stau), # rhs of conservatives
+ static_variablenames = (:flx_sD, :flx_sSr, :flx_stau, # conservatives' fluxes
+ :src_sD, :src_sSr, :src_stau, # sources
+ :D, :Sr, :tau, # true hydrodynamic conservatives
+ :rho, :vr, :p, :eps, # primitives
+ :max_v, :r, # maximum charactersitic speed, grid
+ :E, :Em1, # entropy
+ :flx_E, :flx_Em1, # entropy flux
+ :A, :B, :A_r, :B_r), # metric fields and derivatives
+ cell_variablenames = (:cellmax_v,),
+ global_variablenames = (:t, :tm1) # time steps
+ )
+ _register_variables!(cache, rhs.hrsc)
+ _register_variables!(cache, rhs.tci)
+ # @unpack r = get_static_variables(cache)
+end
+
+
+function _register_variables!(rhs, tci_or_hrsc::Nothing) end
+
+
+_register_variables!(rhs, hrsc::AbstractHRSC) =
+ error("_register_variables: Not implemented for hrsc of type '$(typeof(hrsc))'")
+
+_register_variables!(rhs, hrsc::HRSC.AbstractReconstruction) = nothing
+
+
+function _register_variables!(cache, hrsc::HRSC.AbstractArtificialViscosity)
+ register_variables!(cache,
+ static_variablenames = (:ldg_sD, :ldg_sSr, :ldg_stau, # local DG variable
+ :flx_ldg_sD, :flx_ldg_sSr, :flx_ldg_stau), # local DG flux
+ cell_variablenames = (:mu,), # 'one viscosity to rule them all'
+ )
+end
+
+
+function _register_variables!(cache, hrsc::HRSC.SmoothedArtificialViscosity)
+ _register_variables!(cache, hrsc.av)
+ register_variables!(cache,
+ static_variablenames = (:smoothed_mu,)
+ )
+end
+
+
+function _register_variables!(cache, tci::TCI.AbstractTCI)
+ register_variables!(cache, cell_variablenames = (:flag,:sD_flag,:sSr_flag,:stau_flag))
+end
+
+
+function _register_variables!(cache, tci::TCI.EntropyProduction)
+ register_variables!(cache,
+ cell_variablenames = (:flag,),
+ static_variablenames = (:EP,), # entropy production
+ bdry_variablenames = (:EPjump_lhs, :EPjump_rhs)
+ # entropy production due to jumps at interfaces
+ )
+end
diff --git a/src/projects/GRHD/types.jl b/src/projects/GRHD/types.jl
new file mode 100644
index 0000000000000000000000000000000000000000..9d3d95fb2b4a11a8445c377b64a9bfec852cc0d7
--- /dev/null
+++ b/src/projects/GRHD/types.jl
@@ -0,0 +1,22 @@
+struct Equation{T_EoS<:AbstractEquationOfState} <: AbstractEquation
+ eos::T_EoS
+ atm_factor::Float64
+ atm_threshold::Float64
+ cons2prim_newtonsolver::dg1d.NewtonSolver
+ cons2prim_bisection::dg1d.Bisection
+end
+
+
+struct RHS{T_HRSC <:Maybe{HRSC.AbstractHRSC},
+ T_TCI <:Maybe{TCI.AbstractTCI},
+ T_RSolver <:AbstractRiemannSolver,
+ T_LDG_RSolver <:Maybe{AbstractRiemannSolver},
+ T_AV_RSolver <:Maybe{AbstractRiemannSolver}} <: dg1d.AbstractProject
+
+ equation::Equation
+ hrsc::T_HRSC
+ tci::T_TCI
+ rsolver::T_RSolver
+ ldg_rsolver::T_LDG_RSolver
+ av_rsolver::T_AV_RSolver
+end
diff --git a/src/projects/HRSC/interface.jl b/src/projects/HRSC/interface.jl
index 412de66f4632cca3115b4e71e86f01bb777ebe1b..2e0021b8e7f71d96ed3d430fd0c6af191af2e298 100644
--- a/src/projects/HRSC/interface.jl
+++ b/src/projects/HRSC/interface.jl
@@ -1,20 +1,46 @@
abstract type AbstractHRSC end
-function make_HRSC(parameters, mesh)
- type = Symbol(query(parameters, "type", default=missing))
- concretes = [ Base.typename(ct).name for ct in concretetypes(AbstractHRSC) ]
- if type === :missing
+function make_HRSC(mesh, prms)
+ @unpack method = prms
+
+ if method == "none"
return nothing
+ elseif method == "wenoz"
+ TODO()
+ elseif method == "av"
+ @unpack av_method, av_smoother_order = prms
+ hrsc = if av_method == "mda"
+ MDAViscosity(mesh, prms["mda_cmax"])
+ elseif av_method == "mdh"
+ MDAViscosity(mesh, prms["mdh_cmax"], prms["mdh_ca"], prms["mdh_ck"])
+ elseif av_method == "entropy"
+ MDAViscosity(mesh, prms["entropy_cmax"], prms["entropy_ce"])
+ else
+ error("HRSC: Unknown method requested: 'av_method = $av_method'")
+ end
+
+ if av_smoother_order < 0
+ error("HRSC: 'av_smoother_order' must be positive, found $av_smoother_order")
+ elseif av_smoother_order > 0
+ return SmoothedArtificialViscosity(mesh, hrsc, smoother_order=av_smoother_order)
+ else
+ return hrsc
+ end
+
+ elseif method == "filter"
+ @unpack filter_method = prms
+ if filter_method == "bernstein"
+ return BernsteinReconstruction(mesh)
+ else
+ error("HRSC: Unknown method requested: 'filter_method = $filter_method'")
+ end
+ elseif method == "fv-hybrid"
+ TODO()
+ elseif method == "fd-hybrid"
+ TODO()
+ else
+ error("HRSC: Unknown method requested: '$method'")
end
- @assert type in concretes "HRSC parameter 'type = $type' unknown. Available options: $concretes"
- hrsc_type = eval(type)
- hrsc_parameters = parse_parameters(hrsc_type, parameters)
- hrsc = hrsc_type(; mesh, hrsc_parameters...)
- smoothav = query(parameters, "smoother", default=missing)
- if hrsc isa AbstractArtificialViscosity && smoothav !== missing
- order = query(Int64, parameters, "smoother_order", default=1)
- return SmoothedArtificialViscosity(mesh, hrsc, smoother_order=order)
- end
- return hrsc
+
end
diff --git a/src/projects/SREulerEq/SREulerEq.jl b/src/projects/SREulerEq/SREulerEq.jl
index 0983e6d08ef85bed654c5220d3c0d7634c03520d..3457b9236d1eb0bfea7ed215580e91facb0d4ecb 100644
--- a/src/projects/SREulerEq/SREulerEq.jl
+++ b/src/projects/SREulerEq/SREulerEq.jl
@@ -4,17 +4,72 @@ module SREulerEq
using dg1d
using Plots
-using Match
-import Parameters: @unpack
+include("types.jl")
include("equation.jl")
include("rhs.jl")
include("initialdata.jl")
include("callbacks.jl")
include("plots.jl")
include("boundaryconditions.jl")
-include("main.jl")
+include("setup.jl")
+#######################################################################
+# dg1d project interface definitions #
+#######################################################################
+
+
+dg1d.project_type(::Val{:SREulerEq}) = RHS
+
+# TODO Remove bdryconds
+dg1d.rhs!(env::Environment, p::RHS, bdryconds) = rhs!(env, p, p.hrsc, bdryconds...)
+
+dg1d.timestep!(env::Environment, p::RHS) = timestep(env, p, p.hrsc)
+
+dg1d.@parameters SREulerEq begin
+
+ """
+ Density of artificial atmosphere.
+ """
+ atm_density = 1e-10
+ @check atm_density > 0.0
+
+ """
+ Multiplicative factor to determine the cutoff restmass density below which
+ we impose artificial atmosphere.
+ """
+ atm_threshold_factor = 100.0
+ @check atm_threshold_factor > 0.0
+
+ """
+ Initial data configuration. Available options:
+ - `id_smooth`
+ - `id_periodic_relativistic_shocktube`
+ - `id_periodic_relativistic_shocktube_2`
+ - `id_periodic_shocktube_pressure_jump`
+ - `id_periodic_shocktube_less_extreme`
+ - `id_periodic_sod_shocktube`
+ - `id_sod_shocktube`
+ """
+ id = "sine"
+ @check id in [ "id_smooth",
+ "id_periodic_relativistic_shocktube",
+ "id_periodic_relativistic_shocktube_2",
+ "id_periodic_shocktube_pressure_jump",
+ "id_periodic_shocktube_less_extreme",
+ "id_periodic_sod_shocktube",
+ "id_sod_shocktube" ]
+
+ """
+ Boundary conditions. Available options:
+ - `periodic`
+ - `from_id`
+ """
+ bc = "periodic"
+ @check bc in [ "periodic", "from_id" ]
+
+end
+
end # module
diff --git a/src/projects/SREulerEq/boundaryconditions.jl b/src/projects/SREulerEq/boundaryconditions.jl
index 4ab76543c310fead4a3ad434c5a8ca95fa798f54..d3b547e53960241e75ac528c44116a34ece52e26 100644
--- a/src/projects/SREulerEq/boundaryconditions.jl
+++ b/src/projects/SREulerEq/boundaryconditions.jl
@@ -1,37 +1,46 @@
-function make_BoundaryConditions(rhs::RHS, prms)
- type = query(String, prms, "type")
-
- xrange = extrema(rhs.mesh.x)
-
- lhs_bc, rhs_bc = @match type begin
- "periodic" => (PeriodicBC(), PeriodicBC())
- "from initial data" => begin
- @unpack D, S, tau = get_dynamic_variables(rhs.cache)
- p0 = 0.0
- DStau_lhs = (D[1], S[1], tau[1])
- (p_lhs,) = find_pressure_insane((DStau_lhs..., p0, xrange[1]), rhs.equation)
- (rho_lhs, v_lhs, eps_lhs) = primitives((DStau_lhs..., p_lhs), rhs.equation)
- lhs_state = (DStau_lhs..., p_lhs, rho_lhs, v_lhs, eps_lhs)
- DStau_rhs = (D[end], S[end], tau[end])
- (p_rhs,) = find_pressure_insane((DStau_rhs..., p0, xrange[end]), rhs.equation)
- (rho_rhs, v_rhs, eps_rhs) = primitives((DStau_rhs..., p_rhs), rhs.equation)
- rhs_state = (DStau_rhs..., p_rhs, rho_rhs, v_rhs, eps_rhs)
- (DirichletBC(lhs_state), DirichletBC(rhs_state))
- end
- _ => error("Boundary condition type '$type' not implemented!")
+function make_BoundaryConditions(env, eq, rsolver, prms)
+
+ @unpack mesh = env
+ @unpack bc = prms
+ @unpack xmin, xmax = mesh
+
+ lhs_bc, rhs_bc = if bc == "periodic"
+ PeriodicBC(), PeriodicBC()
+ elseif bc == "from_id"
+ @unpack D, S, tau = get_dynamic_variables(rhs.cache)
+ p0 = 0.0
+ DStau_lhs = (D[1], S[1], tau[1])
+ (p_lhs,) = find_pressure_insane((DStau_lhs..., p0, xmin), rhs.equation)
+ (rho_lhs, v_lhs, eps_lhs) = primitives((DStau_lhs..., p_lhs), rhs.equation)
+ lhs_state = (DStau_lhs..., p_lhs, rho_lhs, v_lhs, eps_lhs)
+ DStau_rhs = (D[end], S[end], tau[end])
+ (p_rhs,) = find_pressure_insane((DStau_rhs..., p0, xmax), rhs.equation)
+ (rho_rhs, v_rhs, eps_rhs) = primitives((DStau_rhs..., p_rhs), rhs.equation)
+ rhs_state = (DStau_rhs..., p_rhs, rho_rhs, v_rhs, eps_rhs)
+ (DirichletBC(lhs_state), DirichletBC(rhs_state))
+ else
+ error("Unknown boundary condition requested: '$bc'")
end
- return BoundaryConditions(lhs_bc, rhs_bc, rhs.rsolver)
+ return BoundaryConditions(lhs_bc, rhs_bc, rsolver)
end
-function make_BoundaryConditions_LDG(ldg_rsolver, prms)
+function make_BoundaryConditions_LDG(env, eq, ldg_rsolver, prms)
+
+ @unpack bc = prms
+
isnothing(ldg_rsolver) && return nothing
- type = query(String, prms, "type")
- lhs_bc, rhs_bc = if type == "periodic"
+
+ lhs_bc, rhs_bc = if bc == "periodic"
PeriodicBC(), PeriodicBC()
- else
+ elseif bc == "no_inflow"
+ OutflowBC(), OutflowBC()
+ elseif bc == "from_id"
OutflowBC(), OutflowBC()
+ else
+ error("Unknown boundary condition requested: '$bc'")
end
+
return BoundaryConditions(lhs_bc, rhs_bc, ldg_rsolver)
end
diff --git a/src/projects/SREulerEq/callbacks.jl b/src/projects/SREulerEq/callbacks.jl
index add288d1047a0b54c36c4ab81a9cc9e1684771b1..4b21c5553c93b16eb01e7d6b4f0634ab55119a97 100644
--- a/src/projects/SREulerEq/callbacks.jl
+++ b/src/projects/SREulerEq/callbacks.jl
@@ -1,11 +1,11 @@
-function make_callback(rhs::RHS, bdryconds::BoundaryConditions)
- cbfn_equation(u, t) = callback_equation(u, t, rhs, isperiodic(bdryconds), rhs.equation)
+function make_callback(env, rhs::RHS, bdryconds::BoundaryConditions)
+ cbfn_equation(u, t) = callback_equation(u, t, env, rhs, isperiodic(bdryconds), rhs.equation)
cb_equation = FunctionCallback(cbfn_equation,
CallbackTiming(every_iteration=1,every_dt=0))
- cbfn_tci(u, t) = callback_tci(u, t, rhs, isperiodic(bdryconds), rhs.tci)
+ cbfn_tci(u, t) = callback_tci(u, t, env, rhs, isperiodic(bdryconds), rhs.tci)
cb_tci = FunctionCallback(cbfn_tci,
CallbackTiming(every_iteration=1, every_dt=0))
- cbfn_hrsc(u, t) = callback_hrsc(u, t, rhs, isperiodic(bdryconds), rhs.hrsc)
+ cbfn_hrsc(u, t) = callback_hrsc(u, t, env, rhs, isperiodic(bdryconds), rhs.hrsc)
cb_hrsc = FunctionCallback(cbfn_hrsc,
CallbackTiming(every_iteration=1,every_dt=0))
callbackset = CallbackSet(cb_equation, cb_tci, cb_hrsc)
@@ -18,8 +18,8 @@ end
#######################################################################
-function callback_equation(state_u, state_t, rhs, isperiodic, eq::SREulerEquation)
- @unpack cache = rhs
+function callback_equation(state_u, state_t, env, rhs, isperiodic, eq::SREulerEquation)
+ @unpack cache = env
wrap_dynamic_variables!(cache, state_u)
@unpack t, tm1 = get_global_variables(cache)
tm1 .= t
@@ -33,17 +33,18 @@ end
#######################################################################
-callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::Nothing) = nothing
-callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::AbstractReconstruction) = nothing
+callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::Nothing) = nothing
+callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::AbstractReconstruction) = nothing
-function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::Maybe{HRSC.AbstractHRSC})
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::Maybe{HRSC.AbstractHRSC})
TODO(callback_hrsc)
end
-function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.AbstractArtificialViscosity)
- @unpack cache, mesh, equation = rhs
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::HRSC.AbstractArtificialViscosity)
+ @unpack cache, mesh = env
+ @unpack equation = rhs
@unpack max_v = get_static_variables(cache)
@unpack mu, cellmax_v = get_cell_variables(cache)
@unpack D = get_dynamic_variables(cache)
@@ -60,16 +61,16 @@ function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.AbstractArt
end
-function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.TCIViscosity)
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::HRSC.TCIViscosity)
@unpack cache, mesh = rhs
@unpack mu, flag = get_cell_variables(cache)
HRSC.compute_viscosity!(mu, flag, hrsc)
end
-function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.SmoothedArtificialViscosity)
- callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc.av)
- @unpack cache = rhs
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::HRSC.SmoothedArtificialViscosity)
+ callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc.av)
+ @unpack cache = env
@unpack mu = get_cell_variables(cache)
indices = findall(muk -> muk > 0.0, mu)
tmp = zeros(size(mu))
@@ -88,15 +89,15 @@ function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.SmoothedArt
end
-function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.EntropyViscosity)
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::HRSC.EntropyViscosity)
TODO()
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack rho = get_dynamic_variables(cache)
@unpack S, Sm1, flx_S, flx_Sm1 = get_static_variables(cache)
@unpack mu, cellmax_v = get_cell_variables(cache)
@unpack t, tm1 = get_global_variables(cache)
HRSC.compute_viscosity!(
- mu,
+ mu,
rho, cellmax_v, S, Sm1, flx_S, flx_Sm1, t[1], tm1[1],
hrsc)
end
@@ -107,19 +108,19 @@ end
#######################################################################
-function callback_tci(state_u, state_t, rhs, isperiodic, tci::Nothing) end
+function callback_tci(state_u, state_t, env, rhs, isperiodic, tci::Nothing) end
-function callback_tci(state_u, state_t, rhs, isperiodic, tci::TCI.Threshold)
+function callback_tci(state_u, state_t, env, rhs, isperiodic, tci::TCI.Threshold)
error("Threshold TCI must be combined with an indicator that is supposed to be zero for " *
"for smooth solutions, but not for shocks. TODO Implement an indicator variable.")
end
-function callback_tci(state_u, state_t, rhs, isperiodic, tci::TCI.EntropyProduction)
+function callback_tci(state_u, state_t, env, rhs, isperiodic, tci::TCI.EntropyProduction)
@warn "Figure out entropy fluxes for special relativistic euler equation!"
TODO()
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack E, Em1, flx_E, flx_Em1, EP = get_static_variables(cache)
@unpack flag = get_cell_variables(cache)
@unpack EPjump_lhs, EPjump_rhs = get_bdry_variables(cache)
@@ -130,11 +131,11 @@ function callback_tci(state_u, state_t, rhs, isperiodic, tci::TCI.EntropyProduct
end
-function callback_tci(state_u, state_t, rhs, isperiodic,
- tci::Union{TCI.Diffusion, TCI.Minmod, TCI.Threshold,
+function callback_tci(state_u, state_t, env, rhs, isperiodic,
+ tci::Union{TCI.Diffusion, TCI.Minmod, TCI.Threshold,
TCI.ModalDecayAverage, TCI.ModalDecayHighest})
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack D, S, tau = get_dynamic_variables(cache)
@unpack flag, D_flag, S_flag, tau_flag = get_cell_variables(cache)
TCI.compute_indicator!(D_flag, D, tci)
diff --git a/src/projects/SREulerEq/equation.jl b/src/projects/SREulerEq/equation.jl
index 805f717e34b8bb75091db1afa840e87ccdd7cd17..519eb3ae5b58d083faa992da115bcbe42ba31ab9 100644
--- a/src/projects/SREulerEq/equation.jl
+++ b/src/projects/SREulerEq/equation.jl
@@ -1,38 +1,9 @@
-#######################################################################
-# Types #
-#######################################################################
-
-
-Base.@kwdef struct Atmosphere
- rho::Float64
- threshold::Float64
-end
-
-
-function parse_parameters(::Type{Atmosphere}, prms)
- rho = query(Float64, prms, "rho", default=1e-10)
- threshold = query(Float64, prms, "threshold", default=100.0)
-
- @assert rho > 0
- @assert threshold > 0
-
- return (; rho, threshold)
-end
-
-
-Atmosphere(prms::AbstractDict) = Atmosphere(; parse_parameters(Atmosphere, prms)...)
-
-
-"""
- SREulerEquation
-
-Special relativistic Euler equations in cartesian coordinates and 1D.
-"""
-Base.@kwdef struct SREulerEquation{T_EoS<:EquationOfState.AbstractEquationOfState} <: AbstractEquation
- eos::T_EoS
- atmosphere::Atmosphere
- cons2prim_newtonsolver::NewtonSolver
- cons2prim_bisection::Bisection
+function make_SREulerEquation(eos, prms)
+ @unpack atm_density, atm_threshold_factor = prms
+ atm = Atmosphere(atm_density, atm_threshold_factor)
+ cons2prim_newtonsolver = NewtonSolver(approx_derivative=true)
+ cons2prim_bisection = Bisection()
+ return SREulerEquation(eos, atm, cons2prim_newtonsolver, cons2prim_bisection)
end
@@ -201,6 +172,7 @@ end
@returns p
end
+
#######################################################################
# LDG #
#######################################################################
diff --git a/src/projects/SREulerEq/initialdata.jl b/src/projects/SREulerEq/initialdata.jl
index 7b087af4e8c90031004765758ffbb06486ee7a90..49ec858f7a283216fd229b27c1c166336db05e69 100644
--- a/src/projects/SREulerEq/initialdata.jl
+++ b/src/projects/SREulerEq/initialdata.jl
@@ -1,31 +1,24 @@
-function initialdata!(rhs::RHS, prms)
- initialdata_equation_of_state!(rhs, rhs.equation.eos, prms)
- # initialdata_entropy_variables!(rhs, prms)
- initialdata_hrsc!(rhs, rhs.hrsc, prms)
- initialdata_tci!(rhs, rhs.tci, prms)
+function initialdata!(env, rhs::RHS, prms)
+ initialdata_equation_of_state!(env, rhs, rhs.equation.eos, prms)
+ # initialdata_entropy_variables!(env, rhs, prms)
+ initialdata_hrsc!(env, rhs, rhs.hrsc, prms)
+ initialdata_tci!(env, rhs, rhs.tci, prms)
return
end
-struct InitialData{T} end
-
-
#######################################################################
# Initial data equation of state #
#######################################################################
-function initialdata_equation_of_state!(rhs::RHS, eos::EquationOfState.IdealGas, parameters)
+function initialdata_equation_of_state!(env, rhs::RHS, eos::EquationOfState.IdealGas, prms)
- # TODO Dropped options from query call in favor of using dispatch for setting up initial data
- # this makes it more attractive for implementing, but less so when trying to figure out
- # which initial data are available
- type = query(parameters, "type")
- type = Symbol(lowercase(replace(type, " "=>"_")))
- @unpack cache, mesh = rhs
+ @unpack id = prms
+ @unpack cache, mesh = env
x = vec(mesh.x)
- Ï0, v0, p0 = initialdata_equation_of_state(InitialData{type}(), eos, x)
+ Ï0, v0, p0 = initialdata_equation_of_state(Val(Symbol(id)), eos, x)
ϵ0 = @. eos(InternalEnergy, Density, Ï0, Pressure, p0)
@@ -36,76 +29,77 @@ function initialdata_equation_of_state!(rhs::RHS, eos::EquationOfState.IdealGas,
@. S = Formulae.S(Ï0, v0, ϵ0, p0)
@. tau = Formulae.Ï„(Ï0, v0, ϵ0, p0)
- broadcast_volume!(cons2prim, rhs.equation, rhs.cache)
+ broadcast_volume!(cons2prim, rhs.equation, cache)
return
end
-function initialdata_equation_of_state(
- ::InitialData{type}, eos, x) where type
+function initialdata_equation_of_state(::Val{type}, eos, x) where type
error("initialdata_equation_of_state: Unknown initial data type '$type'")
end
-function initialdata_equation_of_state(
- ::InitialData{:smooth}, eos::EquationOfState.IdealGas, x)
+function initialdata_equation_of_state(::Val{:smooth}, eos::EquationOfState.IdealGas, x)
xmin, xmax = extrema(x)
@unpack gamma = eos
# from Bugner, PhD thesis 2018
- @assert isapprox(xmin, -1.0) && isapprox(xmax, 1.0)
- @assert gamma == 5/3
+ @toggled_assert isapprox(xmin, -1.0) && isapprox(xmax, 1.0)
+ @toggled_assert gamma == 5/3
Ï0 = @. 1.0 + 0.2 * sin(2.0 * Ï€ * x)
v0 = @. 0.2 * x^0
p0 = @. 1.0 * x^0
+
(Ï0, v0, p0)
end
function initialdata_equation_of_state(
- ::InitialData{:periodic_relativistic_shocktube}, eos::EquationOfState.IdealGas, x)
+ ::Val{:periodic_relativistic_shocktube}, eos::EquationOfState.IdealGas, x)
xmin, xmax = extrema(x)
# inspired by Bugner, PhD thesis 2018; Marti, Mueller, JCP 1996
- @assert isapprox(xmin, -2.0) && isapprox(xmax, 2.0)
- @assert eos.gamma == 5/3
+ @toggled_assert isapprox(xmin, -2.0) && isapprox(xmax, 2.0)
+ @toggled_assert eos.gamma == 5/3
x0 = 1.0
Ï0 = [ (xi < -x0 || xi > x0) ? 10.0 : 1.0 for xi in x ]
v0 = [ (xi < -x0 || xi > x0) ? 0.0 : 0.0 for xi in x ]
p0 = [ (xi < -x0 || xi > x0) ? 13.33 : 1e-7 for xi in x ]
+
(Ï0, v0, p0)
end
function initialdata_equation_of_state(
- ::InitialData{:periodic_relativistic_shocktube_2}, eos::EquationOfState.IdealGas, x)
+ ::Val{:periodic_relativistic_shocktube_2}, eos::EquationOfState.IdealGas, x)
xmin, xmax = extrema(x)
# inspired by Bugner, PhD thesis 2018; Marti, Mueller, JCP 1996
- @assert isapprox(xmin, -2.0) && isapprox(xmax, 2.0)
- @assert eos.gamma == 5/3
+ @toggled_assert isapprox(xmin, -2.0) && isapprox(xmax, 2.0)
+ @toggled_assert eos.gamma == 5/3
x0 = 1.0
Ï0 = [ (xi < -x0 || xi > x0) ? 1.0 : 1.0 for xi in x ]
v0 = [ (xi < -x0 || xi > x0) ? 0.0 : 0.0 for xi in x ]
p0 = [ (xi < -x0 || xi > x0) ? 1000.0 : 0.01 for xi in x ]
+
(Ï0, v0, p0)
end
function initialdata_equation_of_state(
- ::InitialData{:periodic_shocktube_pressure_jump}, eos::EquationOfState.IdealGas, x)
+ ::Val{:periodic_shocktube_pressure_jump}, eos::EquationOfState.IdealGas, x)
xmin, xmax = extrema(x)
γ = eos.gamma
# selfmade, less extreme case than periodic_relativistic_shocktube with only a jump
# in pressure
- @assert isapprox(xmin, -1.0) && isapprox(xmax, 1.0)
- @assert γ == 5/3
+ @toggled_assert isapprox(xmin, -1.0) && isapprox(xmax, 1.0)
+ @toggled_assert γ == 5/3
x0 = 0.5
Ï0 = [ (xi < -x0 || xi > x0) ? 10.0 : 10.0 for xi in x ]
v0 = [ (xi < -x0 || xi > x0) ? 0.0 : 0.0 for xi in x ]
@@ -115,47 +109,50 @@ end
function initialdata_equation_of_state(
- ::InitialData{:periodic_shocktube_less_extreme}, eos::EquationOfState.IdealGas, x)
+ ::Val{:periodic_shocktube_less_extreme}, eos::EquationOfState.IdealGas, x)
xmin, xmax = extrema(x)
# selfmade, less extreme case than periodic_relativistic_shocktube
- @assert isapprox(xmin, -1.0) && isapprox(xmax, 1.0)
- @assert gamma == 5/3
+ @toggled_assert isapprox(xmin, -1.0) && isapprox(xmax, 1.0)
+ @toggled_assert gamma == 5/3
x0 = 0.5
Ï0 = [ (xi < -x0 || xi > x0) ? 10.0 : 0.5 for xi in x ]
v0 = [ (xi < -x0 || xi > x0) ? 0.0 : 0.0 for xi in x ]
p0 = [ (xi < -x0 || xi > x0) ? 8.0 : 1.5 for xi in x ]
+
(Ï0, v0, p0)
end
function initialdata_equation_of_state(
- ::InitialData{:periodic_sod_shocktube}, eos::EquationOfState.IdealGas, x)
+ ::Val{:periodic_sod_shocktube}, eos::EquationOfState.IdealGas, x)
xmin, xmax = extrema(x)
- @assert isapprox(xmin, -2.0) && isapprox(xmax, 2.0)
- @assert eos.gamma == 5/3
+ @toggled_assert isapprox(xmin, -2.0) && isapprox(xmax, 2.0)
+ @toggled_assert eos.gamma == 5/3
x0 = 1.0
Ï0 = [ (xi < -x0 || xi > x0) ? 1.0 : 0.125 for xi in x ]
v0 = [ (xi < -x0 || xi > x0) ? 0.0 : 0.0 for xi in x ]
p0 = [ (xi < -x0 || xi > x0) ? 1.0 : 0.1 for xi in x ]
+
(Ï0, v0, p0)
end
function initialdata_equation_of_state(
- ::InitialData{:sod_shocktube}, eos::EquationOfState.IdealGas, x)
+ ::Val{:sod_shocktube}, eos::EquationOfState.IdealGas, x)
xmin, xmax = extrema(x)
- @assert isapprox(xmin, -1.0) && isapprox(xmax, 1.0)
- @assert eos.gamma == 5/3
+ @toggled_assert isapprox(xmin, -1.0) && isapprox(xmax, 1.0)
+ @toggled_assert eos.gamma == 5/3
x0 = 0.0
Ï0 = [ (xi < x0) ? 1.0 : 0.125 for xi in x ]
v0 = [ (xi < x0) ? 0.0 : 0.0 for xi in x ]
p0 = [ (xi < x0) ? 1.0 : 0.1 for xi in x ]
+
(Ï0, v0, p0)
end
@@ -165,9 +162,9 @@ end
#######################################################################
-function initialdata_entropy_variables!(rhs::RHS, prms)
+function initialdata_entropy_variables!(env, rhs::RHS, prms)
@TODO
- @unpack cache = rhs
+ @unpack cache = env
@unpack t, tm1 = get_global_variables(cache)
@unpack E, Em1, flx_E, flx_Em1 = get_static_variables(cache)
tm1 .= 0.0
@@ -183,18 +180,18 @@ end
#######################################################################
-function initialdata_hrsc!(rhs::RHS, hrsc::Maybe{HRSC.AbstractHRSC}, prms) end
+function initialdata_hrsc!(env, rhs::RHS, hrsc::Maybe{HRSC.AbstractHRSC}, prms) end
-function initialdata_hrsc!(rhs::RHS, hrsc::HRSC.AbstractArtificialViscosity, prms)
- @unpack mu = get_cell_variables(rhs.cache)
+function initialdata_hrsc!(env, rhs::RHS, hrsc::HRSC.AbstractArtificialViscosity, prms)
+ @unpack mu = get_cell_variables(env.cache)
@. mu = 0.0
end
-function initialdata_hrsc!(rhs::RHS, hrsc::HRSC.SmoothedArtificialViscosity, prms)
- initialdata_hrsc!(rhs, hrsc.av, prms)
- @unpack smoothed_mu = get_static_variables(rhs.cache)
+function initialdata_hrsc!(env, rhs::RHS, hrsc::HRSC.SmoothedArtificialViscosity, prms)
+ initialdata_hrsc!(env, rhs, hrsc.av, prms)
+ @unpack smoothed_mu = get_static_variables(env.cache)
@. smoothed_mu = 0.0
end
@@ -204,11 +201,11 @@ end
#######################################################################
-function initialdata_tci!(rhs::RHS, tci::Nothing, prms) end
+function initialdata_tci!(env, rhs::RHS, tci::Nothing, prms) end
-function initialdata_tci!(rhs::RHS, tci::TCI.EntropyProduction, prms)
- @unpack cache = rhs
+function initialdata_tci!(env, rhs::RHS, tci::TCI.EntropyProduction, prms)
+ @unpack cache = env
@unpack EP = get_static_variables(cache)
@unpack EPjump_lhs, EPjump_rhs = get_bdry_variables(cache)
@. EP = 0.0
@@ -217,7 +214,7 @@ function initialdata_tci!(rhs::RHS, tci::TCI.EntropyProduction, prms)
end
-function initialdata_tci!(rhs::RHS, tci::TCI.AbstractTCI, prms)
- @unpack flag = get_cell_variables(rhs.cache)
+function initialdata_tci!(env, rhs::RHS, tci::TCI.AbstractTCI, prms)
+ @unpack flag = get_cell_variables(env.cache)
@. flag = 0.0
end
diff --git a/src/projects/SREulerEq/main.jl b/src/projects/SREulerEq/main.jl
deleted file mode 100644
index 1e674a5102afd7b993e00ed09ee48f851e952136..0000000000000000000000000000000000000000
--- a/src/projects/SREulerEq/main.jl
+++ /dev/null
@@ -1,157 +0,0 @@
-using OrderedCollections
-"""
- main(parameters_filename)
-
-Run program specified by TOML file `parameters_filename`.
-
-The following sections in the parameter file are recognized:
-- [`Mesh`](@ref)
-- [`EquationOfState`](@ref)
-- [`Initialdata`](@ref)
-- [`Timestep`](@ref)
-- [`HRSC`](@ref)
-- [`Output`](@ref)
-"""
-function main(parameters_filename=""; override::AbstractDict=Dict())
-
- override_prms = FlatDict(override)
-
- parameters_filename = abspath(parameters_filename)
- outputdir = mk_outputdir(parameters_filename)
- cp(parameters_filename, joinpath(outputdir, basename(parameters_filename)))
-
- # read parameters
- if !isfile(parameters_filename)
- error("Cannot locate parameter file '$parameters_filename'")
- end
- prms = FlatDict(TOML.parsefile(parameters_filename))
- prms = merge(prms, override_prms)
- ordered_prms = OrderedDict(prms)
- display(sort(ordered_prms))
-
- default_eos_prms = Dict("type"=>"IdealGas", "gamma"=>3.0)
- eos_prms = trim_headers(prms, "EquationOfState")
- amend!(eos_prms, default_eos_prms)
- diff = diffkeys(eos_prms, default_eos_prms)
- @assert length(diff) == 0 "Section 'Plot': Unknown options '$diff' encountered!"
- eos = EquationOfState.make_EquationOfState(eos_prms)
- atm = Atmosphere(; parse_parameters(Atmosphere, trim_headers(prms, "Atmosphere"))...)
- cons2prim_newtonsolver = NewtonSolver(approx_derivative=true)
- cons2prim_bisection = Bisection()
-
- # setup RHS
- mesh = make_Mesh(trim_headers(prms, "Mesh"))
- cache = Cache(mesh)
- equation = SREulerEquation(eos, atm, cons2prim_newtonsolver, cons2prim_bisection)
- tci = TCI.make_TCI(trim_headers(prms, "TCI"), mesh)
- hrsc = HRSC.make_HRSC(trim_headers(prms, "HRSC"), mesh)
- rsolver = ApproxRiemannSolver(flux, speed, equation)
- ldg_rsolver, av_rsolver = if hrsc isa AbstractArtificialViscosity
- ApproxRiemannSolver(ldg_flux, ldg_speed, equation),
- ApproxRiemannSolver(av_flux, av_flux, equation)
- else
- nothing, nothing
- end
- rhs = RHS(cache=cache, mesh=mesh, equation=equation, tci=tci, hrsc=hrsc,
- rsolver=rsolver, ldg_rsolver=ldg_rsolver, av_rsolver=av_rsolver)
-
- _register_variables!(rhs)
- register_variables!(cache, rsolver)
- !isnothing(ldg_rsolver) && register_variables!(cache, ldg_rsolver)
- !isnothing(av_rsolver) && register_variables!(cache, av_rsolver, dont_register=true)
- display(cache)
-
- initialdata!(rhs, trim_headers(prms, "Initialdata"))
-
- default_boundaryconditions_prms = Dict("type"=>"periodic")
- boundaryconditions_prms = trim_headers(prms, "BoundaryConditions")
- amend!(boundaryconditions_prms, default_boundaryconditions_prms)
- diff = diffkeys(boundaryconditions_prms, default_boundaryconditions_prms)
- @assert length(diff) == 0 "Section 'Plot': Unknown options '$diff' encountered!"
- bdryconds = make_BoundaryConditions(rhs, boundaryconditions_prms)
- ldg_bdryconds, av_bdryconds = if hrsc isa AbstractArtificialViscosity
- make_BoundaryConditions_LDG(ldg_rsolver, boundaryconditions_prms),
- make_BoundaryConditions_LDG(av_rsolver, boundaryconditions_prms)
- else
- nothing, nothing
- end
-
- rhs_fn!(du, u, t) = rhs!(du, u, t, rhs, bdryconds, ldg_bdryconds, av_bdryconds)
- timestep_fn(u, t, n) = timestep(u, t, rhs)
-
- ### Callbacks
- rhscb = make_callback(rhs, bdryconds)
-
- # Plotting
- default_plot_prms = Dict("every_iteration"=>0, "every_dt"=>0.0, "pause"=>false)
- plot_prms = trim_headers(prms, "Plot")
- dg1d.amend!(plot_prms, default_plot_prms)
- diff = diffkeys(plot_prms, default_plot_prms)
- @assert length(diff) == 0 "Section 'Plot': Unknown options '$diff' encountered!"
- theme(:dark, markerstrokewidth=0.1, marker=:circle)
- plotfn(u,t) = plot(rhs)
- plotcb = PlotCallback(plotfn,
- CallbackTiming(every_dt=plot_prms["every_dt"],
- every_iteration=plot_prms["every_iteration"]),
- pause=plot_prms["pause"])
- default_output_prms = Dict("every_iteration"=>0,
- "every_dt"=>0.5,
- "variables"=>["rho", "q", "E"],
- "filename"=>joinpath(outputdir, "output"))
- output_prms = trim_headers(prms, "Output")
- amend!(output_prms, default_output_prms)
- diff = diffkeys(output_prms, default_output_prms)
- @assert length(diff) == 0 "Section 'Output': Unknown options '$diff' encountered!"
- savecb = SaveCallback(Symbol.(output_prms["variables"]), cache, mesh,
- output_prms["filename"],
- CallbackTiming(every_dt=output_prms["every_dt"],
- every_iteration=output_prms["every_iteration"])
- )
-
- # artificially smooth initial data
- smoother_called = false
- function id_smoother_fn(u, t)
- if !smoother_called
- println("Smoothing initial data...")
- @unpack D, S, tau = get_dynamic_variables(cache)
- @unpack flag = get_cell_variables(cache)
- id_hrsc = HRSC.BernsteinReconstruction(mesh)
- HRSC.reconstruct!(D, flag, id_hrsc, isperiodic=isperiodic(bdryconds), m=0.0)
- HRSC.reconstruct!(S, flag, id_hrsc, isperiodic=isperiodic(bdryconds))
- HRSC.reconstruct!(tau, flag, id_hrsc, isperiodic=isperiodic(bdryconds), m=0.0)
- smoother_called = true
- broadcast_volume!(cons2prim, equation, cache)
- end
- end
- id_smoother_cb = FunctionCallback(id_smoother_fn,
- CallbackTiming(every_iteration=1))
-
- # Combine callbacks
- callbacks = CallbackSet(rhscb.callbacks..., id_smoother_cb, savecb, plotcb)
-
- aligned_times = get(output_prms, "aligned_times", [])
- timealigned_savecb, timestep_fn = if length(aligned_times) > 0
- atscb = TimeAlignedSaveCallback(aligned_times, timestep_fn, savecb)
- atscb, atscb.aligned_timestep_fn
- else
- nothing, timestep_fn
- end
-
- !isnothing(timealigned_savecb) && push!(callbacks, timealigned_savecb)
-
- default_evolve_prms = Dict("algorithm"=>"LSERK4", "tend"=>1.0)
- evolve_prms = trim_headers(prms, "Evolution")
- amend!(evolve_prms, default_evolve_prms)
- diff = diffkeys(evolve_prms, default_evolve_prms)
- @assert length(diff) == 0 "Section 'Evolution': Unknown options '$diff' encountered!"
-
- @unpack D, S, tau = get_dynamic_variables(cache)
- state_u = vcat(D, S, tau)
-
- state_u = dg1d.evolve_ERK(rhs_fn!, state_u, timestep_fn,
- (0.0, Float64(evolve_prms["tend"])),
- Symbol(evolve_prms["algorithm"]),
- callback_fullstep=callbacks)
-
- return
-end
diff --git a/src/projects/SREulerEq/plots.jl b/src/projects/SREulerEq/plots.jl
index 01079f07fa5b7d157cf980eed900e1691d9a17f1..57f53915d55673a027b8a10f1f9496269bcc26b4 100644
--- a/src/projects/SREulerEq/plots.jl
+++ b/src/projects/SREulerEq/plots.jl
@@ -1,8 +1,8 @@
-function plot(rhs::RHS)
- plt_eq = plot_equation(rhs)
- plt_rhs = nothing # plot_rhs(rhs)
- plt_hrsc = plot_hrsc(rhs, rhs.hrsc)
- plt_tci = plot_tci(rhs, rhs.tci)
+function plot(env, rhs::RHS)
+ plt_eq = plot_equation(env, rhs)
+ plt_rhs = nothing # plot_rhs(env, rhs)
+ plt_hrsc = plot_hrsc(env, rhs, rhs.hrsc)
+ plt_tci = plot_tci(env, rhs, rhs.tci)
list_sub_plts = [ plt_eq, plt_rhs, plt_hrsc, plt_tci ]
filter!(!isnothing, list_sub_plts)
@@ -11,7 +11,7 @@ function plot(rhs::RHS)
sub_plt = Plots.plot(list_sub_plts..., layout=grid(layout..., heights=plt_heights))
# combine plots
- @unpack cache = rhs
+ @unpack cache = env
@unpack t = get_global_variables(cache)
plt = Plots.plot(sub_plt, layout=layout, size=(1400,1000), plot_title="t=$(t[1])",
left_margin=10Plots.mm)
@@ -24,9 +24,9 @@ end
#######################################################################
-function plot_equation(rhs)
+function plot_equation(env, rhs)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack D, S, tau = get_dynamic_variables(cache)
@unpack rho, p, eps, v = get_static_variables(cache)
@unpack t = get_global_variables(cache)
@@ -66,13 +66,13 @@ end
#######################################################################
-function plot_hrsc(rhs, hrsc::Nothing) end
-function plot_hrsc(rhs, hrsc::AbstractReconstruction) end
+function plot_hrsc(env, rhs, hrsc::Nothing) end
+function plot_hrsc(env, rhs, hrsc::AbstractReconstruction) end
-function plot_hrsc(rhs, hrsc::HRSC.AbstractArtificialViscosity)
+function plot_hrsc(env, rhs, hrsc::HRSC.AbstractArtificialViscosity)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack mu = get_cell_variables(cache)
@unpack x = mesh
@@ -94,11 +94,11 @@ function plot_hrsc(rhs, hrsc::HRSC.AbstractArtificialViscosity)
end
-function plot_hrsc(rhs, hrsc::HRSC.SmoothedArtificialViscosity)
+function plot_hrsc(env, rhs, hrsc::HRSC.SmoothedArtificialViscosity)
- plt_mu = plot_hrsc(rhs, hrsc.av)
+ plt_mu = plot_hrsc(env, rhs, hrsc.av)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack smoothed_mu = get_static_variables(cache)
@unpack x = mesh
@@ -128,12 +128,12 @@ end
#######################################################################
-function plot_tci(rhs, tci::Nothing) end
+function plot_tci(env, rhs, tci::Nothing) end
-function plot_tci(rhs, tci::TCI.AbstractTCI)
+function plot_tci(env, rhs, tci::TCI.AbstractTCI)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack flag = get_cell_variables(cache)
@unpack x = mesh
@@ -163,9 +163,9 @@ function plot_tci(rhs, tci::TCI.AbstractTCI)
end
-function plot_tci(rhs, tci::TCI.EntropyProduction)
+function plot_tci(env, rhs, tci::TCI.EntropyProduction)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack EP = get_static_variables(cache)
@unpack flag = get_cell_variables(cache)
@unpack x = mesh
@@ -202,9 +202,9 @@ end
#######################################################################
-function plot_rhs(rhs)
+function plot_rhs(env, rhs)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack rhs_rho, rhs_q, rhs_E = get_rhs_variables(cache)
@unpack t = get_global_variables(cache)
@unpack x = mesh
diff --git a/src/projects/SREulerEq/rhs.jl b/src/projects/SREulerEq/rhs.jl
index 9fcf39384e38343f81ffb6e95c1008711d305cf9..aca7ed492caa6f503a06f46e5e38178c4f761a53 100644
--- a/src/projects/SREulerEq/rhs.jl
+++ b/src/projects/SREulerEq/rhs.jl
@@ -1,104 +1,11 @@
-Base.@kwdef struct RHS{T_Equation <:SREulerEquation,
- T_RSolver <:AbstractRiemannSolver,
- T_HRSC <:Maybe{HRSC.AbstractHRSC},
- T_TCI <:Maybe{TCI.AbstractTCI},
- T_LDG_RSolver<:Maybe{AbstractRiemannSolver},
- T_AV_RSolver<:Maybe{AbstractRiemannSolver}}
-
- cache::Cache
- mesh::Mesh
- equation::T_Equation
- rsolver::T_RSolver
- hrsc::T_HRSC = nothing
- tci::T_TCI = nothing
- ldg_rsolver::T_LDG_RSolver = nothing
- av_rsolver::T_AV_RSolver = nothing
-end
-
-
-#######################################################################
-# Register variables #
-#######################################################################
-
-
-# Julia limitation/bug: Overloading imported functions without their Module prefix shadows their
-# definition in this module without a warning or error. This does not happen for methods
-# defined in Base.
-# Because of that we prefix with _ below.
-# Cf. https://github.com/JuliaLang/julia/issues/15510
-function _register_variables!(rhs::RHS)
- register_variables!(rhs.cache,
- dynamic_variablenames = (:D, :S, :tau,), # conservatives
- rhs_variablenames = (:rhs_D, :rhs_S, :rhs_tau), # rhs of conservatives
- static_variablenames = (:flx_D, :flx_S, :flx_tau, # conservatives' fluxes
- :rho, :v, :p, :eps, # primitives
- :max_v, :x, # maximum charactersitic speed, grid
- :E, :Em1, # entropy
- :flx_E, :flx_Em1), # entropy flux
- cell_variablenames = (:cellmax_v,),
- global_variablenames = (:t, :tm1) # time steps
- )
- _register_variables!(rhs, rhs.hrsc)
- _register_variables!(rhs, rhs.tci)
- @unpack x = get_static_variables(rhs.cache)
- @. x = rhs.mesh.x[:]
-end
-
-
-function _register_variables!(rhs, tci_or_hrsc::Nothing) end
-
-
-_register_variables!(rhs, hrsc::AbstractHRSC) =
- error("_register_variables: Not implemented for hrsc of type '$(typeof(hrsc))'")
-
-_register_variables!(rhs, hrsc::HRSC.AbstractReconstruction) = nothing
-
-
-function _register_variables!(rhs, hrsc::HRSC.AbstractArtificialViscosity)
- register_variables!(rhs.cache,
- static_variablenames = (:ldg_D, :ldg_S, :ldg_tau, # local DG variable
- :flx_ldg_D, :flx_ldg_S, :flx_ldg_tau), # local DG flux
- cell_variablenames = (:mu,), # 'one viscosity to rule them all'
- )
-end
-
-
-function _register_variables!(rhs, hrsc::HRSC.SmoothedArtificialViscosity)
- _register_variables!(rhs, hrsc.av)
- register_variables!(rhs.cache,
- static_variablenames = (:smoothed_mu,)
- )
-end
-
-
-function _register_variables!(rhs, tci::TCI.AbstractTCI)
- register_variables!(rhs.cache, cell_variablenames = (:flag,:D_flag,:S_flag,:tau_flag))
-end
-
-
-function _register_variables!(rhs, tci::TCI.EntropyProduction)
- register_variables!(rhs.cache,
- cell_variablenames = (:flag,),
- static_variablenames = (:EP,), # entropy production
- bdry_variablenames = (:EPjump_lhs, :EPjump_rhs)
- # entropy production due to jumps at interfaces
- )
-end
-
-
#######################################################################
# Timestep #
#######################################################################
-function timestep(state_u, t, rhs::RHS)
- wrap_dynamic_variables!(rhs.cache, state_u)
- return timestep(rhs, rhs.hrsc)
-end
-
-
-function timestep(rhs::RHS, hrsc::Maybe{HRSC.AbstractHRSC})
- @unpack cache, mesh, equation = rhs
+function timestep(env, rhs::RHS, hrsc::Maybe{HRSC.AbstractHRSC})
+ @unpack cache, mesh = env
+ @unpack equation = rhs
@unpack max_v = get_static_variables(cache)
broadcast_volume!(cons2prim, equation, cache)
@@ -119,8 +26,9 @@ function timestep(rhs::RHS, hrsc::Maybe{HRSC.AbstractHRSC})
end
-function timestep(rhs::RHS, hrsc::HRSC.AbstractArtificialViscosity)
- @unpack cache, mesh, equation = rhs
+function timestep(env, rhs::RHS, hrsc::HRSC.AbstractArtificialViscosity)
+ @unpack cache, mesh = env
+ @unpack equation = rhs
broadcast_volume!(cons2prim, equation, cache)
broadcast_volume!(speed, equation, cache)
@@ -152,19 +60,10 @@ end
#######################################################################
-function rhs!(state_du, state_u, t, rhs::RHS, bdryconds::BoundaryConditions,
- ldg_bdryconds::Maybe{BoundaryConditions}, av_bdryconds::Maybe{BoundaryConditions})
-
- @unpack cache = rhs
- wrap_dynamic_variables!(cache, state_u)
- wrap_rhs_variables!(cache, state_du)
- rhs!(rhs, t, rhs.hrsc, bdryconds, ldg_bdryconds, av_bdryconds)
-end
-
-
-function rhs!(rhs::RHS, t, hrsc::Nothing, bdryconds, ldg_bdryconds, av_bdryconds)
+function rhs!(env, rhs::RHS, hrsc::Nothing, bdryconds, ldg_bdryconds, av_bdryconds)
- @unpack cache, mesh, equation, rsolver = rhs
+ @unpack cache, mesh = env
+ @unpack equation, rsolver = rhs
@unpack D, S, tau = get_dynamic_variables(cache)
@unpack flx_D, flx_S, flx_tau = get_static_variables(cache)
@unpack rhs_D, rhs_S, rhs_tau = get_rhs_variables(cache)
@@ -175,7 +74,7 @@ function rhs!(rhs::RHS, t, hrsc::Nothing, bdryconds, ldg_bdryconds, av_bdryconds
broadcast_volume!(cons2prim, equation, cache)
broadcast_volume!(flux, equation, cache)
broadcast_faces!(lax_friedrich_flux, rsolver, cache, mesh)
- broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, t)
+ broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, 0.0)
compute_rhs_weak_form!(rhs_D, flx_D, lhs_numflx_D, rhs_numflx_D, mesh)
compute_rhs_weak_form!(rhs_S, flx_S, lhs_numflx_S, rhs_numflx_S, mesh)
@@ -185,10 +84,11 @@ function rhs!(rhs::RHS, t, hrsc::Nothing, bdryconds, ldg_bdryconds, av_bdryconds
end
-function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractReconstruction,
+function rhs!(env, rhs::RHS, hrsc::HRSC.AbstractReconstruction,
bdryconds, ldg_bdryconds, av_bdryconds)
- @unpack cache, mesh, equation, rsolver = rhs
+ @unpack cache, mesh = env
+ @unpack equation, rsolver = rhs
@unpack D, S, tau = get_dynamic_variables(cache)
@unpack flx_D, flx_S, flx_tau = get_static_variables(cache)
@unpack rhs_D, rhs_S, rhs_tau = get_rhs_variables(cache)
@@ -213,11 +113,12 @@ function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractReconstruction,
end
-function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractArtificialViscosity,
+function rhs!(env, rhs::RHS, hrsc::HRSC.AbstractArtificialViscosity,
bdryconds, ldg_bdryconds, av_bdryconds)
- @unpack cache, mesh, equation, rsolver,
- ldg_rsolver, av_rsolver = rhs
+ @unpack cache, mesh = env
+ @unpack equation, rsolver, ldg_rsolver,
+ av_rsolver = rhs
@unpack D, S, tau = get_dynamic_variables(cache)
@unpack flx_D, flx_S, flx_tau,
ldg_D, ldg_S, ldg_tau,
@@ -233,7 +134,7 @@ function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractArtificialViscosity,
## solve auxiliary equation: q + ∂x u = 0
broadcast_volume!(ldg_flux, equation, cache)
broadcast_faces!(mean_flux, ldg_rsolver, cache, mesh)
- broadcast_boundaryconditions!(mean_flux, ldg_bdryconds, cache, mesh, t)
+ broadcast_boundaryconditions!(mean_flux, ldg_bdryconds, cache, mesh, 0.0)
compute_rhs_weak_form!(ldg_D, flx_ldg_D, lhs_numflx_ldg_D, rhs_numflx_ldg_D, mesh)
compute_rhs_weak_form!(ldg_S, flx_ldg_S, lhs_numflx_ldg_S, rhs_numflx_ldg_S, mesh)
compute_rhs_weak_form!(ldg_tau, flx_ldg_tau, lhs_numflx_ldg_tau, rhs_numflx_ldg_tau, mesh)
@@ -241,12 +142,12 @@ function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractArtificialViscosity,
## compute rhs of regularized equation: ∂t u + ∂x f + ∂x μ q = 0
broadcast_volume!(av_flux, equation, cache)
broadcast_faces!(mean_flux, av_rsolver, cache, mesh)
- broadcast_boundaryconditions!(mean_flux, av_bdryconds, cache, mesh, t)
+ broadcast_boundaryconditions!(mean_flux, av_bdryconds, cache, mesh, 0.0)
broadcast_volume!(cons2prim, equation, cache)
broadcast_volume!(flux, equation, cache)
broadcast_faces!(lax_friedrich_flux, rsolver, cache, mesh)
- broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, t)
+ broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, 0.0)
@. flx_D += flx_ldg_D
@. flx_S += flx_ldg_S
diff --git a/src/projects/SREulerEq/setup.jl b/src/projects/SREulerEq/setup.jl
new file mode 100644
index 0000000000000000000000000000000000000000..117ea492f101218e03f6360b24de4551e33a1439
--- /dev/null
+++ b/src/projects/SREulerEq/setup.jl
@@ -0,0 +1,122 @@
+function RHS(env::Environment, prms)
+
+ # construct RHS
+ eos = EquationOfState.make_EquationOfState(prms["EquationOfState"])
+ equation = make_SREulerEquation(eos, prms["SREulerEq"])
+ tci = TCI.make_TCI(env.mesh, prms["TCI"])
+ hrsc = HRSC.make_HRSC(env.mesh, prms["HRSC"])
+ rsolver = ApproxRiemannSolver(flux, speed, equation)
+ ldg_rsolver, av_rsolver = if hrsc isa AbstractArtificialViscosity
+ ApproxRiemannSolver(ldg_flux, ldg_speed, equation),
+ ApproxRiemannSolver(av_flux, av_flux, equation)
+ else
+ nothing, nothing
+ end
+
+ rhs = RHS(equation, rsolver, hrsc, tci, ldg_rsolver, av_rsolver)
+
+ # register variables
+ # TODO add a ::Nothing overload for register_variables!
+ # TODO Somehow replace _register_variables! with register_variables!
+ @unpack cache = env
+ _register_variables!(cache, rhs)
+ register_variables!(cache, rsolver)
+ !isnothing(ldg_rsolver) && register_variables!(cache, ldg_rsolver)
+ !isnothing(av_rsolver) && register_variables!(cache, av_rsolver, dont_register=true)
+ display(cache)
+
+ # setup boundary conditions
+ bdryconds = make_BoundaryConditions(env, equation, rsolver, prms["SREulerEq"])
+ ldg_bdryconds, av_bdryconds = if hrsc isa AbstractArtificialViscosity
+ make_BoundaryConditions_LDG(env, equation, ldg_rsolver, prms["SREulerEq"]),
+ make_BoundaryConditions_LDG(env, equation, av_rsolver, prms["SREulerEq"])
+ else
+ nothing, nothing
+ end
+
+ # TODO Bdry conditions should be moved to project and should only be just a 'lightweight'
+ # descriptor. Furthermore, they should also register variables in cache themselves.
+ # Note: Atm, only the 'from_initial_data' option would require a place to store the data
+ # somewhere. For periodic BCs we don't need extra storage at all.
+ # But what other types of boundary conditions will be needed in the future?
+
+ # setup initial data
+ initialdata!(env, rhs, prms["SREulerEq"])
+
+ # setup callbacks
+ rhscb = make_callback(env, rhs, bdryconds)
+
+ append!(env.callbacks, CallbackSet(rhscb.callbacks...))
+
+ return rhs, (bdryconds, ldg_bdryconds, av_bdryconds)
+end
+
+
+#######################################################################
+# Register variables #
+#######################################################################
+
+
+# Julia limitation/bug: Overloading imported functions without their Module prefix shadows their
+# definition in this module without a warning or error. This does not happen for methods
+# defined in Base.
+# Because of that we prefix with _ below.
+# Cf. https://github.com/JuliaLang/julia/issues/15510
+function _register_variables!(cache, rhs::RHS)
+ register_variables!(cache,
+ dynamic_variablenames = (:D, :S, :tau,), # conservatives
+ rhs_variablenames = (:rhs_D, :rhs_S, :rhs_tau), # rhs of conservatives
+ static_variablenames = (:flx_D, :flx_S, :flx_tau, # conservatives' fluxes
+ :rho, :v, :p, :eps, # primitives
+ :max_v, :x, # maximum charactersitic speed, grid
+ :E, :Em1, # entropy
+ :flx_E, :flx_Em1), # entropy flux
+ cell_variablenames = (:cellmax_v,),
+ global_variablenames = (:t, :tm1) # time steps
+ )
+ _register_variables!(cache, rhs.hrsc)
+ _register_variables!(cache, rhs.tci)
+ @unpack x = get_static_variables(cache)
+ # @. x = env.mesh.x[:]
+end
+
+
+function _register_variables!(cache, tci_or_hrsc::Nothing) end
+
+
+_register_variables!(cache, hrsc::AbstractHRSC) =
+ error("_register_variables: Not implemented for hrsc of type '$(typeof(hrsc))'")
+
+_register_variables!(cache, hrsc::HRSC.AbstractReconstruction) = nothing
+
+
+function _register_variables!(cache, hrsc::HRSC.AbstractArtificialViscosity)
+ register_variables!(cache,
+ static_variablenames = (:ldg_D, :ldg_S, :ldg_tau, # local DG variable
+ :flx_ldg_D, :flx_ldg_S, :flx_ldg_tau), # local DG flux
+ cell_variablenames = (:mu,), # 'one viscosity to rule them all'
+ )
+end
+
+
+function _register_variables!(cache, hrsc::HRSC.SmoothedArtificialViscosity)
+ _register_variables!(cache, hrsc.av)
+ register_variables!(cache,
+ static_variablenames = (:smoothed_mu,)
+ )
+end
+
+
+function _register_variables!(cache, tci::TCI.AbstractTCI)
+ register_variables!(cache, cell_variablenames = (:flag,:D_flag,:S_flag,:tau_flag))
+end
+
+
+function _register_variables!(cache, tci::TCI.EntropyProduction)
+ register_variables!(cache,
+ cell_variablenames = (:flag,),
+ static_variablenames = (:EP,), # entropy production
+ bdry_variablenames = (:EPjump_lhs, :EPjump_rhs)
+ # entropy production due to jumps at interfaces
+ )
+end
diff --git a/src/projects/SREulerEq/types.jl b/src/projects/SREulerEq/types.jl
new file mode 100644
index 0000000000000000000000000000000000000000..8a4c2e7c17c857cc3eebdb797c31e384748d18d6
--- /dev/null
+++ b/src/projects/SREulerEq/types.jl
@@ -0,0 +1,28 @@
+Base.@kwdef struct Atmosphere
+ rho::Float64
+ threshold::Float64
+end
+
+
+struct SREulerEquation{T_EoS<:EquationOfState.AbstractEquationOfState} <: AbstractEquation
+ eos::T_EoS
+ atmosphere::Atmosphere
+ cons2prim_newtonsolver::NewtonSolver
+ cons2prim_bisection::Bisection
+end
+
+
+struct RHS{T_Equation <:SREulerEquation,
+ T_RSolver <:AbstractRiemannSolver,
+ T_HRSC <:Maybe{HRSC.AbstractHRSC},
+ T_TCI <:Maybe{TCI.AbstractTCI},
+ T_LDG_RSolver<:Maybe{AbstractRiemannSolver},
+ T_AV_RSolver<:Maybe{AbstractRiemannSolver}}
+
+ equation::T_Equation
+ rsolver::T_RSolver
+ hrsc::T_HRSC
+ tci::T_TCI
+ ldg_rsolver::T_LDG_RSolver
+ av_rsolver::T_AV_RSolver
+end
diff --git a/src/projects/ScalarEquation/ScalarEquation.jl b/src/projects/ScalarEquation/ScalarEquation.jl
index 337d04d6ed5183325bb5ac701570c0c348dcc5de..07f4d4d039b42047ef4b4fc4789c27964fba231d 100644
--- a/src/projects/ScalarEquation/ScalarEquation.jl
+++ b/src/projects/ScalarEquation/ScalarEquation.jl
@@ -3,18 +3,75 @@ module ScalarEquation
using dg1d
-using Match
-import Parameters: @unpack, @with_kw
using Plots
+include("types.jl")
include("equation.jl")
include("rhs.jl")
include("initialdata.jl")
include("plots.jl")
include("callbacks.jl")
include("boundaryconditions.jl")
-include("main.jl")
+include("setup.jl")
+
+
+#######################################################################
+# dg1d project interface definitions #
+#######################################################################
+
+
+dg1d.project_type(::Val{:ScalarEquation}) = RHS
+
+# TODO Remove bdryconds
+dg1d.rhs!(env::Environment, p::RHS, bdryconds) = rhs!(env, p, p.hrsc, bdryconds...)
+
+dg1d.timestep!(env::Environment, p::RHS) = timestep(env, p, p.hrsc)
+
+dg1d.@parameters ScalarEquation begin
+
+ """
+ Available options:
+ - `advection`
+ - `burgers`
+ """
+ equation = "advection"
+ @check equation in [ "advection", "burgers" ]
+
+ "Advection speed"
+ advection_v = 1.0
+
+ """
+ Initial data configuration. Available options:
+ - `sine`
+ - `bump`
+ - `gaussian`
+ """
+ id = "sine"
+ @check id in [ "sine", "bump", "gaussian" ]
+
+ "Period of sine initial data"
+ id_sine_period = 1
+
+ "Phase of sine initial data"
+ id_sine_phase = 0.0
+
+ "RMS width of gaussian"
+ id_gaussian_width = 1.0
+
+ "Amplitude of gaussian"
+ id_gaussian_amplitude = 1.0
+
+ """
+ Boundary conditions. Available options:
+ - `periodic`
+ - `no_inflow`
+ - `from_id`
+ """
+ bc = "periodic"
+ @check bc in [ "periodic", "no_inflow", "from_id" ]
+
+end
end # ScalarEquation
diff --git a/src/projects/ScalarEquation/boundaryconditions.jl b/src/projects/ScalarEquation/boundaryconditions.jl
index c0d06a4ca0c371f6707f11800d2693e831dcae99..9ee123823495542f3295f3de0278b6d5549a199c 100644
--- a/src/projects/ScalarEquation/boundaryconditions.jl
+++ b/src/projects/ScalarEquation/boundaryconditions.jl
@@ -1,30 +1,39 @@
-function make_BoundaryConditions_equation(rhs::RHS, prms)
- type = query(String, prms, "type")
-
- lhs_bc, rhs_bc = @match type begin
- "periodic" => (PeriodicBC(), PeriodicBC())
- "no inflow" => (DirichletBC((0.0,)), OutflowBC())
- "initial data" => begin
- @unpack u = get_dynamic_variables(rhs.cache)
- u_lhs, u_rhs = deepcopy(u[1]), deepcopy(u[end])
- lhs_bdry = t -> (u_lhs,)
- rhs_bdry = t -> (u_rhs,)
- DirichletBC(lhs_bdry), DirichletBC(rhs_bdry)
- end
- _ => error("Boundary condition type '$type' not implemented!")
+function make_BoundaryConditions(env, equation, rsolver, prms)
+ @unpack bc = prms
+
+ lhs_bc, rhs_bc = if bc == "periodic"
+ PeriodicBC(), PeriodicBC()
+ elseif bc == "no_inflow"
+ DirichletBC((0.0,)), OutflowBC()
+ elseif bc == "from_id"
+ @unpack u = get_dynamic_variables
+ u_lhs, u_rhs = deepcopy(u[1]), deepcopy(u[end])
+ lhs_bdry = t -> (u_lhs,)
+ rhs_bdry = t -> (u_rhs,)
+ DirichletBC(lhs_bdry), DirichletBC(rhs_bdry)
+ else
+ error("Unknown boundary condition requested: '$bc'")
end
-
- return BoundaryConditions(lhs_bc, rhs_bc, rhs.rsolver)
+
+ return BoundaryConditions(lhs_bc, rhs_bc, rsolver)
end
-function make_BoundaryConditions_LDG(ldg_rsolver, prms)
+function make_BoundaryConditions_LDG(env, equation, ldg_rsolver, prms)
+
+ @unpack bc = prms
+
isnothing(ldg_rsolver) && return nothing
- type = query(String, prms, "type")
- lhs_bc, rhs_bc = if type == "periodic"
+
+ lhs_bc, rhs_bc = if bc == "periodic"
PeriodicBC(), PeriodicBC()
- else
+ elseif bc == "no_inflow"
+ OutflowBC(), OutflowBC()
+ elseif bc == "from_id"
OutflowBC(), OutflowBC()
+ else
+ error("Unknown boundary condition requested: '$bc'")
end
+
return BoundaryConditions(lhs_bc, rhs_bc, ldg_rsolver)
end
diff --git a/src/projects/ScalarEquation/callbacks.jl b/src/projects/ScalarEquation/callbacks.jl
index 9d11c0e368d0a4be77d5175975a8b7fee0002a78..5bde11a0faab87b12aa5025648c91e1e0b0586f9 100644
--- a/src/projects/ScalarEquation/callbacks.jl
+++ b/src/projects/ScalarEquation/callbacks.jl
@@ -1,11 +1,11 @@
-function make_callback(rhs::RHS, bdryconds::BoundaryConditions)
- cbfn_equation(u, t) = callback_equation(u, t, rhs, isperiodic(bdryconds), rhs.equation)
+function make_callback(env, rhs::RHS, bdryconds::BoundaryConditions)
+ cbfn_equation(u, t) = callback_equation(u, t, env, rhs, isperiodic(bdryconds), rhs.equation)
cb_equation = FunctionCallback(cbfn_equation,
CallbackTiming(every_iteration=1,every_dt=0))
- cbfn_tci(u, t) = callback_tci(u, t, rhs, isperiodic(bdryconds), rhs.tci)
+ cbfn_tci(u, t) = callback_tci(u, t, env, rhs, isperiodic(bdryconds), rhs.tci)
cb_tci = FunctionCallback(cbfn_tci,
CallbackTiming(every_iteration=1, every_dt=0))
- cbfn_hrsc(u, t) = callback_hrsc(u, t, rhs, isperiodic(bdryconds), rhs.hrsc)
+ cbfn_hrsc(u, t) = callback_hrsc(u, t, env, rhs, isperiodic(bdryconds), rhs.hrsc)
cb_hrsc = FunctionCallback(cbfn_hrsc,
CallbackTiming(every_iteration=1,every_dt=0))
callbackset = CallbackSet(cb_equation, cb_tci, cb_hrsc)
@@ -18,8 +18,9 @@ end
#######################################################################
-function callback_equation(state_u, state_t, rhs, isperiodic, eq::AbstractScalarEquation)
- @unpack cache, mesh, equation = rhs
+function callback_equation(state_u, state_t, env, rhs, isperiodic, eq::AbstractScalarEquation)
+ @unpack cache, mesh = env
+ @unpack equation = rhs
wrap_dynamic_variables!(cache, state_u)
@unpack u = get_dynamic_variables(cache)
@unpack E, Em1, F, Fm1 = get_static_variables(cache)
@@ -41,19 +42,20 @@ end
#######################################################################
-callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::Nothing) = nothing
+callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::Nothing) = nothing
-callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.AbstractHRSC) =
+callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::HRSC.AbstractHRSC) =
error("not yet implemented for type $(typeof(hrsc))")
-function callback_hrsc(state_u, state_t, rhs, isperiodic,
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic,
hrsc::Union{HRSC.BernsteinReconstruction, HRSC.WENOReconstruction}) end
-function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.AbstractArtificialViscosity)
- @unpack cache, mesh, equation = rhs
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::HRSC.AbstractArtificialViscosity)
+ @unpack cache, mesh = env
+ @unpack equation = rhs
@unpack v = get_static_variables(cache)
@unpack mu, max_v = get_cell_variables(cache)
@unpack u = get_dynamic_variables(cache)
@@ -70,24 +72,26 @@ function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.AbstractArt
end
-function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.TCIViscosity)
- @unpack cache, mesh = rhs
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::HRSC.TCIViscosity)
+ @unpack cache, mesh = env
@unpack mu, flag = get_cell_variables(cache)
HRSC.compute_viscosity!(mu, flag, hrsc)
end
-function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.SmoothedArtificialViscosity)
- callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc.av)
- @unpack cache, rsolver = rhs
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::HRSC.SmoothedArtificialViscosity)
+ callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc.av)
+ @unpack cache = env
+ @unpack rsolver = rhs
@unpack mu = get_cell_variables(cache)
@unpack smoothed_mu = get_static_variables(cache)
smoothed_mu .= vec(hrsc.smoother(mu,isperiodic))
end
-function callback_hrsc(state_u, state_t, rhs, isperiodic, hrsc::HRSC.EntropyViscosity)
- @unpack cache, mesh, equation = rhs
+function callback_hrsc(state_u, state_t, env, rhs, isperiodic, hrsc::HRSC.EntropyViscosity)
+ @unpack cache, mesh = env
+ @unpack equation = rhs
@unpack u = get_dynamic_variables(cache)
@unpack E, F, Em1, Fm1, v = get_static_variables(cache)
@unpack mu, max_v = get_cell_variables(cache)
@@ -110,15 +114,15 @@ end
#######################################################################
-callback_tci(state_u, state_t, rhs, isperiodic, tci::Nothing) = nothing
+callback_tci(state_u, state_t, env, rhs, isperiodic, tci::Nothing) = nothing
-callback_tci(state_u, state_t, rhs, isperiodic, tci::TCI.AbstractTCI) =
+callback_tci(state_u, state_t, env, rhs, isperiodic, tci::TCI.AbstractTCI) =
error("not yet implemented for type $(typeof(tci))")
-function callback_tci(state_u, state_t, rhs, isperiodic, tci::TCI.EntropyProduction)
- @unpack cache, mesh = rhs
+function callback_tci(state_u, state_t, env, rhs, isperiodic, tci::TCI.EntropyProduction)
+ @unpack cache, mesh = env
@unpack u = get_dynamic_variables(cache)
@unpack E, Em1, F, Fm1, EP = get_static_variables(cache)
@unpack flag = get_cell_variables(cache)
@@ -130,10 +134,10 @@ function callback_tci(state_u, state_t, rhs, isperiodic, tci::TCI.EntropyProduct
end
-function callback_tci(state_u, state_t, rhs, isperiodic,
+function callback_tci(state_u, state_t, env, rhs, isperiodic,
tci::Union{TCI.Diffusion, TCI.Minmod, TCI.Threshold,
TCI.ModalDecayAverage, TCI.ModalDecayHighest})
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack u = get_dynamic_variables(cache)
@unpack flag = get_cell_variables(cache)
TCI.compute_indicator!(flag, u, tci)
diff --git a/src/projects/ScalarEquation/equation.jl b/src/projects/ScalarEquation/equation.jl
index c710a7c1a80b335182c66a88304415b567ce0209..81b9571b92bbc068f6e8d25d7c6dcc4bbbf48ebc 100644
--- a/src/projects/ScalarEquation/equation.jl
+++ b/src/projects/ScalarEquation/equation.jl
@@ -1,40 +1,23 @@
-"""
- abstract type AbstractScalarEquation <: dg1d.AbstractEquation end
-
-Supertype for ScalarEquation problems.
-New problems should subtype it and also implement the following methods:
-- `make_Equation(prms, ::Type{NewProblem})`: factory for an instance of type `NewProblem`.
-The following methods have to be defined with the `@with_signature` interface:
-- `speed(args, ::NewProblem)`
-- `flux(args, ::NewProblem)`
-- `entropy(args, ::NewProblem)`
-- `entropy_flux(args, ::NewProblem)`
-"""
-abstract type AbstractScalarEquation <: dg1d.AbstractEquation end
-
-
#######################################################################
-# Equation Factory #
+# Factory #
#######################################################################
function make_Equation(prms)
- available_types = concretetypes(AbstractScalarEquation)
- all_options = last.(split.(string.(available_types), "."))
- type = query(prms, "type", options=all_options)
- expr = quote $(Symbol(type)) end
- return make_Equation(prms, eval(expr))
+ @unpack equation = prms
+ if equation == "advection"
+ return Advection(prms["advection_v"])
+ elseif equation == "burgers"
+ return Burgers()
+ else
+ error("Unknown equation requested: '$equation'")
+ end
end
-
#######################################################################
# Advection #
#######################################################################
-struct Advection <: AbstractScalarEquation
- v::Float64
-end
-
function make_Equation(prms, ::Type{Advection})
v = query(Float64, prms, "v", default=1.0)
return Advection(v)
@@ -68,8 +51,6 @@ end
# Burgers #
#######################################################################
-struct Burgers <: AbstractScalarEquation end
-
make_Equation(prms, ::Type{Burgers}) = Burgers()
make_Burgers(prms) = Burgers()
diff --git a/src/projects/ScalarEquation/initialdata.jl b/src/projects/ScalarEquation/initialdata.jl
index 1807e8869cc6ed73a002386d9edf4c0bb105c0e5..089a8e1d55733c6782c7d399c188fbebeea8e056 100644
--- a/src/projects/ScalarEquation/initialdata.jl
+++ b/src/projects/ScalarEquation/initialdata.jl
@@ -1,57 +1,34 @@
-function initialdata!(rhs::RHS, prms)
- initialdata_equation!(rhs, rhs.equation, prms)
- initialdata_hrsc!(rhs, rhs.hrsc, prms)
- initialdata_tci!(rhs, rhs.tci, prms)
+function initialdata!(env, rhs::RHS, prms)
+ initialdata_equation!(env, rhs, rhs.equation, prms)
+ initialdata_hrsc!(env, rhs, rhs.hrsc, prms)
+ initialdata_tci!(env, rhs, rhs.tci, prms)
end
-function initialdata_equation!(rhs::RHS, equation::AbstractScalarEquation, parameters)
- @unpack cache, mesh = rhs
- @unpack u = get_dynamic_variables(cache)
- @unpack L, x = mesh
- x = vec(x)
- tstart = 0.0 # always
- id_fn = make_function_initialdata_equation(parameters)
- @. u = id_fn(x)
- return
-end
+function initialdata_equation!(env, rhs::RHS, equation::AbstractScalarEquation, prms)
-
-"""
- make_function_initialdata(parameters)
-
-Generate a function `id_fn(x)` to compute initial data on grid `x`.
-`parameters` must contain sections `Initialdata, Mesh`.
-"""
-function make_function_initialdata_equation(parameters)
-
- type = query(parameters, "type", options=["sine", "bump", "bugner gaussian", "gaussian"])
- xrange = query(Vector{Float64}, parameters, "xrange", default = [-1.0,1.0])
- L = first(diff(xrange))
-
- id_fn = if type == "sine"
- period = query(Int64, parameters, "period", default=1) * 2.0 * π / L
- phase = query(Float64, parameters, "phase", default=0.0)
- (x) -> sin(period * x + phase)
- elseif type == "bump"
- xmin, xmax = xrange
- height = query(Float64, parameters, "height", default=1)
- width = query(Float64, parameters, "width", default=(xmax-xmin)/2)
- @assert width < L "Initialdata 'bump': Bump width exceeds x domain '$xrange'"
- xl, xr = xmin+(L-width)/2, xmin+(L+width)/2
- (x) -> (x < xl || x > xr) ? 0.5 : height
- elseif type == "gaussian"
- amplitude = query(Float64, parameters, "height", default=1)
- width = query(Float64, parameters, "width", default=1)
- x0 = (xrange[2] + xrange[1]) / 2
- display(x0)
- (x) -> amplitude * exp(-(x-x0)^2 / width^2)
- elseif type == "bugner gaussian"
- xmin, xmax = xrange
- @assert xmin == -3.0 && xmax == 3.0
- A, sigma = 1.0, 0.4
- (x) -> A * exp(-x^2 / sigma^2)
- elseif type == "smooth bump"
+ @unpack cache, mesh = env
+ @unpack u = get_dynamic_variables(cache)
+ @unpack L, x, xmin, xmax = mesh
+
+ @unpack id = prms
+ x = x[:]
+
+ id_fn = if id == "sine"
+ @unpack id_sine_period, id_sine_phase = prms
+ u .= @. sin.(2 * pi * id_sine_period * x/L + id_sine_phase)
+ elseif id == "bump"
+ @unpack id_bump_height, id_bump_width = prms
+ @toggled_assert id_bump_width < L "Initialdata 'bump': Bump width exceeds x domain '$xrange'"
+ xl, xr = xmin+(L-id_bump_width)/2, xmin+(L+id_bump_width)/2
+ u .= @. (x < xl || x > xr) ? 0.5 : id_bump_height
+ elseif id == "gaussian"
+ A = prms["id_gaussian_amplitude"]
+ w = prms["id_gaussian_width"]
+ @toggled_assert w > 0.0
+ x0 = (xmax + xmin) / 2
+ u .= @. A * exp(-(x-x0)^2 / w^2)
+ elseif id == "smooth bump"
error("TODO Implement 'smooth bump' initial data")
# elseif type === :smooth_bump
# xl = xmin + L * 0.10
@@ -70,25 +47,25 @@ function make_function_initialdata_equation(parameters)
# end
# end
else
- error("Requested initialdata type $id_type not implemented!")
+ error("Unknown initialdata requested: '$id'")
end
- return id_fn
+ return
end
-function initialdata_hrsc!(rhs::RHS, hrsc::Maybe{HRSC.AbstractHRSC}, prms) end
+function initialdata_hrsc!(env, rhs::RHS, hrsc::Maybe{HRSC.AbstractHRSC}, prms) end
-function initialdata_hrsc!(rhs::RHS, hrsc::HRSC.AbstractArtificialViscosity, prms)
- @unpack mu = get_cell_variables(rhs.cache)
+function initialdata_hrsc!(env, rhs::RHS, hrsc::HRSC.AbstractArtificialViscosity, prms)
+ @unpack mu = get_cell_variables(env.cache)
@. mu = 0.0
end
-function initialdata_hrsc!(rhs::RHS, hrsc::HRSC.EntropyViscosity, prms)
- @unpack E, Em1, F, Fm1 = get_static_variables(rhs.cache)
- @unpack mu = get_cell_variables(rhs.cache)
+function initialdata_hrsc!(env, rhs::RHS, hrsc::HRSC.EntropyViscosity, prms)
+ @unpack E, Em1, F, Fm1 = get_static_variables(env.cache)
+ @unpack mu = get_cell_variables(env.cache)
@. E = 0.0
@. F = 0.0
@. Em1 = 0.0
@@ -97,10 +74,10 @@ function initialdata_hrsc!(rhs::RHS, hrsc::HRSC.EntropyViscosity, prms)
end
-function initialdata_tci!(rhs::RHS, tci::Nothing, prms) end
+function initialdata_tci!(env, rhs::RHS, tci::Nothing, prms) end
-function initialdata_tci!(rhs::RHS, tci::TCI.AbstractTCI, prms)
- @unpack flag = get_cell_variables(rhs.cache)
+function initialdata_tci!(env, rhs::RHS, tci::TCI.AbstractTCI, prms)
+ @unpack flag = get_cell_variables(env.cache)
@. flag = 0.0
end
diff --git a/src/projects/ScalarEquation/main.jl b/src/projects/ScalarEquation/main.jl
deleted file mode 100644
index 6e7c794dd03fd8b332e367f632006d1a2ad30781..0000000000000000000000000000000000000000
--- a/src/projects/ScalarEquation/main.jl
+++ /dev/null
@@ -1,199 +0,0 @@
-using InteractiveUtils
-using PlotThemes
-
-
-function section(prms, name)
- subprms = get(prms, name, missing)
- ismissing(subprms) && error("Missing parameter section '$name'")
- return subprms
-end
-
-
-@doc """
- main(parameters_filename)
-
-Run program specified by TOML file `parameters_filename`.
-
-The following sections in the parameter file are recognized:
-- [`Mesh`](@ref)
-- [`Equation`](@ref)
-- [`Initialdata`](@ref)
-- [`Timestep`](@ref)
-- [`Analysis`](@ref)
-- [`HRSC`](@ref)
-- [`Output`](@ref)
-
-TODO The above should generate links to the parameter lists for the
-corresponding sections.
-
-TODO Add backup::Bool flag to main.
-"""
-function main(parameters_filename="")
-
- parameters_filename = abspath(parameters_filename)
- outputdir = mk_outputdir(parameters_filename)
- cp(parameters_filename, joinpath(outputdir, basename(parameters_filename)))
-
- # read parameters
- if !isfile(parameters_filename)
- error("Cannot locate parameter file '$parameters_filename'")
- end
- prms = FlatDict(TOML.parsefile(parameters_filename))
-
- # setup RHS
- mesh = make_Mesh(trim_headers(prms, "Mesh"))
- cache = Cache(mesh)
- equation = make_Equation(trim_headers(prms, "Equation"))
- tci = TCI.make_TCI(trim_headers(prms, "TCI"), mesh)
- hrsc = HRSC.make_HRSC(trim_headers(prms, "HRSC"), mesh)
- rsolver = ApproxRiemannSolver(flux, speed, equation)
- ldg_rsolver, av_rsolver = if hrsc isa AbstractArtificialViscosity
- ApproxRiemannSolver(ldg_flux, ldg_speed, equation),
- ApproxRiemannSolver(av_flux, av_speed, equation)
- else
- nothing, nothing
- end
- rhs = RHS(cache=cache, mesh=mesh, equation=equation, tci=tci, hrsc=hrsc,
- rsolver=rsolver, ldg_rsolver=ldg_rsolver, av_rsolver=av_rsolver)
-
- _register_variables!(rhs)
- register_variables!(cache, rsolver)
- !isnothing(ldg_rsolver) && register_variables!(cache, ldg_rsolver)
- !isnothing(av_rsolver) && register_variables!(cache, av_rsolver, dont_register=true)
- display(cache)
-
- initialdata!(rhs, trim_headers(prms, "Initialdata", "Mesh"))
-
- default_boundaryconditions_prms = Dict("type"=>"periodic")
- boundaryconditions_prms = trim_headers(prms, "BoundaryConditions")
- amend!(boundaryconditions_prms, default_boundaryconditions_prms)
- diff = diffkeys(boundaryconditions_prms, default_boundaryconditions_prms)
- @assert length(diff) == 0 "Section 'BoundaryConditions': Unknown options '$diff' encountered!"
- bdryconds = make_BoundaryConditions_equation(rhs, boundaryconditions_prms)
- ldg_bdryconds = make_BoundaryConditions_LDG(ldg_rsolver, boundaryconditions_prms)
- av_bdryconds = make_BoundaryConditions_LDG(av_rsolver, boundaryconditions_prms)
-
- rhs_fn!(du, u, t) = rhs!(du, u, t, rhs, bdryconds, ldg_bdryconds, av_bdryconds)
- timestep_fn(u, t, n) = timestep(u, t, rhs)
-
- ### Callbacks
- rhscb = make_callback(rhs, bdryconds)
-
- # Plotting
- default_plot_prms = Dict("every_iteration"=>0,
- "every_dt"=>0.0,
- "pause"=>false)
- plot_prms = trim_headers(prms, "Plot")
- dg1d.amend!(plot_prms, default_plot_prms)
- diff = diffkeys(plot_prms, default_plot_prms)
- theme(:dark, markerstrokewidth=0.1, marker=:circle)
- plotfn(u,t) = plot(rhs)
- @assert length(diff) == 0 "Section 'Plot': Unknown options '$diff' encountered!"
- plotcb = PlotCallback(plotfn,
- CallbackTiming(every_dt=plot_prms["every_dt"],
- every_iteration=plot_prms["every_iteration"]),
- pause=plot_prms["pause"])
- default_output_prms = Dict("every_iteration"=>0,
- "every_dt"=>0.5,
- "variables"=>["u"],
- "filename"=>joinpath(outputdir, "output"),
- "aligned_times"=>[])
- output_prms = trim_headers(prms, "Output")
- amend!(output_prms, default_output_prms)
- diff = diffkeys(output_prms, default_output_prms)
- @assert length(diff) == 0 "Section 'Output': Unknown options '$diff' encountered!"
- savecb = SaveCallback(Symbol.(output_prms["variables"]), cache, mesh,
- output_prms["filename"],
- CallbackTiming(every_dt=output_prms["every_dt"],
- every_iteration=output_prms["every_iteration"])
- )
-
- # Combine callbacks
- callbacks = CallbackSet(rhscb.callbacks..., savecb, plotcb)
-
- aligned_times = get(output_prms, "aligned_times", [])
- timealigned_savecb, timestep_fn = if length(aligned_times) > 0
- atscb = TimeAlignedSaveCallback(aligned_times, timestep_fn, savecb)
- atscb, atscb.aligned_timestep_fn
- else
- nothing, timestep_fn
- end
-
- !isnothing(timealigned_savecb) && push!(callbacks, timealigned_savecb)
- # callback_fn(u, t, i) = callbacks(u, t, i)
-
- default_evolve_prms = Dict("algorithm"=>"LSERK4", "tend"=>1.0)
- evolve_prms = trim_headers(prms, "Evolution")
- amend!(evolve_prms, default_evolve_prms)
- diff = diffkeys(evolve_prms, default_evolve_prms)
- @assert length(diff) == 0 "Section 'Evolution': Unknown options '$diff' encountered!"
-
- @unpack u = get_dynamic_variables(cache)
-
- u = dg1d.evolve_ERK(rhs_fn!, u, timestep_fn,
- (0.0, Float64(evolve_prms["tend"])),
- Symbol(evolve_prms["algorithm"]),
- callback_fullstep=callbacks)
-
- return
-
- # default_output_filename = Filesystem.basename(Filesystem.splitext(parameters_filename)[1])
- # default_output_prms = Dict("every_iteration"=>10,
- # "every_dt"=>0.0,
- # "variables"=>["u"],
- # "filename"=>joinpath(outputdir, "output"))
- # output_prms = pop!(prms, "Output", Dict())
- # amend!(output_prms, default_output_prms)
- # diff = diffkeys(output_prms, default_output_prms)
- # @assert length(diff) == 0 "Section 'Output': Unknown options '$diff' encountered!"
- #
- # default_plot_prms = Dict("every_iteration"=>10,
- # "every_dt"=>0.0,
- # "pause"=>true)
- # plot_prms = pop!(prms, "Plot", Dict())
- # dg1d.amend!(plot_prms, default_plot_prms)
- # diff = diffkeys(plot_prms, default_plot_prms)
- # @assert length(diff) == 0 "Section 'Plot': Unknown options '$diff' encountered!"
- #
- # savecb = SaveCallback(Symbol.(output_prms["variables"]), D, output_prms["filename"],
- # CallbackTiming(every_dt=output_prms["every_dt"],
- # every_iteration=output_prms["every_iteration"]))
- # theme(:dark, markerstrokewidth=0.1)
- # plotfn(u,t) = plot(I, u, t)
- # plotcb = PlotCallback(plotfn,
- # CallbackTiming(every_dt=plot_prms["every_dt"],
- # every_iteration=plot_prms["every_iteration"]),
- # pause=plot_prms["pause"])
- # hrscfn(u,t) = update_hrsc!(I, u, t)
- # hrsccb = FunctionCallback(hrscfn, CallbackTiming(every_iteration=1, every_dt=0))
- # tcicb = TCI.make_callback(tci, I, CallbackTiming(every_iteration=1, every_dt=0))
- # # tcifn(u, t) = TCI.compute_indicator!(D, tci)
- # # tcicb = FunctionCallback(tcifn, CallbackTiming(every_iteration=1, every_dt=0))
- #
- # # callbacks = CallbackSet(plotcb, savecb)
- # callbacks = CallbackSet(tcicb, hrsccb, plotcb, savecb)
- # # callbacks = CallbackSet(hrsccb, plotcb)
- # callback_fn(u, t, i) = callbacks(u, t, i)
- # # callback_fn(u, t, i) = plotcb(u, t, i)
- #
- #
- # ##### Evolution
- #
- # rhs_fn!(du, u, t) = rhs!(du, u, t, I)
- # timestep_fn(u, t) = timestep(u, t, I)
- #
- # default_evolve_prms = Dict("algorithm"=>"LSERK4",
- # "tend"=>1.0)
- # evolve_prms = pop!(prms, "Evolution", Dict())
- # amend!(evolve_prms, default_evolve_prms)
- # diff = diffkeys(evolve_prms, default_evolve_prms)
- # @assert length(diff) == 0 "Section 'Evolution': Unknown options '$diff' encountered!"
- #
- # # dg1d.evolve_ERK(myrhs!, u, timestep, (0.0, 1.0), :RK4, post_hook=post_hook)
- # u = dg1d.evolve_ERK(rhs_fn!, u, timestep_fn,
- # (0.0, Float64(evolve_prms["tend"])),
- # Symbol(evolve_prms["algorithm"]),
- # post_hook=callback_fn)
-
- return
-end
diff --git a/src/projects/ScalarEquation/plots.jl b/src/projects/ScalarEquation/plots.jl
index c2a7ac3aa8432f1381810468cbae0abf6021e393..f5d54cf879067b5d4aa8540e633e569be88b73b7 100644
--- a/src/projects/ScalarEquation/plots.jl
+++ b/src/projects/ScalarEquation/plots.jl
@@ -1,7 +1,7 @@
-function plot(rhs::RHS)
- plt_eq = plot_equation(rhs, rhs.equation)
- plt_hrsc = plot_hrsc(rhs, rhs.hrsc)
- plt_tci = plot_tci(rhs, rhs.tci)
+function plot(env, rhs::RHS)
+ plt_eq = plot_equation(env, rhs, rhs.equation)
+ plt_hrsc = plot_hrsc(env, rhs, rhs.hrsc)
+ plt_tci = plot_tci(env, rhs, rhs.tci)
# drop missing plots and figure out layout
plts = [ plt_eq, plt_hrsc, plt_tci ]
@@ -9,7 +9,7 @@ function plot(rhs::RHS)
layout = (length(plts), 1)
# combine plots
- @unpack cache = rhs
+ @unpack cache = env
@unpack t = get_global_variables(cache)
plt = Plots.plot(plts..., layout=layout, size=(1400,1000), plot_title="t=$(t[1])",
left_margin=10Plots.mm)
@@ -22,9 +22,9 @@ end
#######################################################################
-function plot_equation(rhs, eq::AbstractScalarEquation)
+function plot_equation(env, rhs, eq::AbstractScalarEquation)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack u = get_dynamic_variables(cache)
@unpack x = mesh
@@ -51,12 +51,12 @@ end
#######################################################################
-function plot_hrsc(rhs, hrsc::Nothing) end
+function plot_hrsc(env, rhs, hrsc::Nothing) end
-function plot_hrsc(rhs, hrsc::HRSC.AbstractArtificialViscosity)
+function plot_hrsc(env, rhs, hrsc::HRSC.AbstractArtificialViscosity)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack x = mesh
xlabel = "x"
@@ -83,7 +83,7 @@ function plot_hrsc(rhs, hrsc::HRSC.AbstractArtificialViscosity)
end
-plot_hrsc(rhs, hrsc::Union{HRSC.BernsteinReconstruction,
+plot_hrsc(env, rhs, hrsc::Union{HRSC.BernsteinReconstruction,
HRSC.WENOReconstruction}) = nothing
@@ -92,12 +92,12 @@ plot_hrsc(rhs, hrsc::Union{HRSC.BernsteinReconstruction,
#######################################################################
-function plot_tci(rhs, tci::Nothing) end
+function plot_tci(env, rhs, tci::Nothing) end
-function plot_tci(rhs, tci::TCI.AbstractTCI)
+function plot_tci(env, rhs, tci::TCI.AbstractTCI)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack flag = get_cell_variables(cache)
@unpack x = mesh
@@ -126,9 +126,9 @@ function plot_tci(rhs, tci::TCI.AbstractTCI)
end
-function plot_tci(rhs, tci::TCI.EntropyProduction)
+function plot_tci(env, rhs, tci::TCI.EntropyProduction)
- @unpack cache, mesh = rhs
+ @unpack cache, mesh = env
@unpack EP = get_static_variables(cache)
@unpack flag = get_cell_variables(cache)
@unpack x = mesh
diff --git a/src/projects/ScalarEquation/rhs.jl b/src/projects/ScalarEquation/rhs.jl
index efe61315d12d7d7c75cfd518e17337db972130dc..4884320c5c04900d554fe9981095bde816009405 100644
--- a/src/projects/ScalarEquation/rhs.jl
+++ b/src/projects/ScalarEquation/rhs.jl
@@ -1,102 +1,11 @@
-# Alternative names:
-# - Discretization ... too long
-# - Approximation or Approx ... the second one might fit, but is also unspecific
-# - Implementation or Impl ... too generic
-# - Project ... too generic
-# - State ... is at really a state?
-Base.@kwdef mutable struct RHS{T_Equation <:AbstractScalarEquation,
- T_RSolver <:AbstractRiemannSolver,
- T_HRSC <:Maybe{HRSC.AbstractHRSC},
- T_TCI <:Maybe{TCI.AbstractTCI},
- T_LDG_RSolver<:Maybe{AbstractRiemannSolver},
- T_AV_RSolver<:Maybe{AbstractRiemannSolver}}
-
- cache::Cache
- mesh::Mesh
- equation::T_Equation
- rsolver::T_RSolver
- hrsc::T_HRSC = nothing
- tci::T_TCI = nothing
- ldg_rsolver::T_LDG_RSolver = nothing
- av_rsolver::T_AV_RSolver = nothing
-end
-
-
-#######################################################################
-# Register variables #
-#######################################################################
-
-
-
-# Julia limitation/bug: Overloading imported functions without their Module prefix shadows their
-# definition in this module without a warning or error. This does not happen for methods
-# defined in Base.
-# Because of that we prefix with _ below.
-# Cf. https://github.com/JuliaLang/julia/issues/15510
-# Perhaps this is the import vs. using confusion people talk on the forums.
-function _register_variables!(rhs::RHS)
- _register_variables!(rhs, rhs.equation)
- _register_variables!(rhs, rhs.hrsc)
- _register_variables!(rhs, rhs.tci)
-end
-
-
-function _register_variables!(rhs, eq::AbstractScalarEquation)
- register_variables!(rhs.cache,
- dynamic_variablenames = (:u,),
- rhs_variablenames = (:rhs_u,),
- static_variablenames = (:flx_u, :E, :Em1, :F, :Fm1, :v),
- cell_variablenames = (:max_v,),
- global_variablenames = (:t, :tm1))
-end
-
-
-function _register_variables!(rhs, hrsc::HRSC.AbstractArtificialViscosity)
- register_variables!(rhs.cache,
- static_variablenames = (:q, :flx_q),
- cell_variablenames = (:mu,))
-end
-
-
-function _register_variables!(rhs, hrsc::HRSC.SmoothedArtificialViscosity)
- _register_variables!(rhs, hrsc.av)
- register_variables!(rhs.cache,
- static_variablenames = (:smoothed_mu,))
-end
-
-
-function _register_variables!(rhs, hrsc::Maybe{HRSC.AbstractHRSC}) end
-
-
-function _register_variables!(rhs, tci::Nothing) end
-
-
-function _register_variables!(rhs, tci::TCI.AbstractTCI)
- register_variables!(rhs.cache, cell_variablenames = (:flag,))
-end
-
-
-function _register_variables!(rhs, tci::TCI.EntropyProduction)
- register_variables!(rhs.cache,
- cell_variablenames = (:flag,),
- static_variablenames = (:EP,),
- bdry_variablenames = (:EPjump_lhs, :EPjump_rhs))
-end
-
-
#######################################################################
# Timestep #
#######################################################################
-function timestep(state_u, t, rhs::RHS)
- wrap_dynamic_variables!(rhs.cache, state_u)
- return timestep(t, rhs, rhs.hrsc)
-end
-
-
-function timestep(t, rhs::RHS, hrsc::Maybe{HRSC.AbstractHRSC})
- @unpack cache, mesh, equation = rhs
+function timestep(env, rhs::RHS, hrsc::Maybe{HRSC.AbstractHRSC})
+ @unpack cache, mesh = env
+ @unpack equation = rhs
@unpack v = get_static_variables(cache)
broadcast_volume!(speed, equation, cache)
@@ -120,24 +29,14 @@ end
#######################################################################
-function rhs!(state_du, state_u, t, rhs::RHS,
- bdryconds::BoundaryConditions,
- ldg_bdryconds::Maybe{BoundaryConditions},
- av_bdryconds::Maybe{BoundaryConditions})
-
- wrap_dynamic_variables!(rhs.cache, state_u)
- wrap_rhs_variables!(rhs.cache, state_du)
- return rhs!(rhs, t, rhs.hrsc, bdryconds, ldg_bdryconds, av_bdryconds)
-end
-
-
-function rhs!(rhs::RHS, t, hrsc::Nothing, bdryconds, ldg_bdryconds::Nothing, av_bdryconds::Nothing)
+function rhs!(env, rhs::RHS, hrsc::Nothing, bdryconds, ldg_bdryconds::Nothing, av_bdryconds::Nothing)
- @unpack cache, mesh, equation, rsolver = rhs
+ @unpack cache, mesh = env
+ @unpack equation, rsolver = rhs
broadcast_volume!(flux, equation, cache)
broadcast_faces!(lax_friedrich_flux, rsolver, cache, mesh)
- broadcast_boundaryconditions!(central_flux, bdryconds, cache, mesh, t)
+ broadcast_boundaryconditions!(central_flux, bdryconds, cache, mesh, 0.0)
@unpack flx_u, v = get_static_variables(cache)
@unpack u = get_dynamic_variables(cache)
@@ -149,10 +48,11 @@ function rhs!(rhs::RHS, t, hrsc::Nothing, bdryconds, ldg_bdryconds::Nothing, av_
end
-function rhs!(rhs::RHS, t, hrsc::Union{BernsteinReconstruction, WENOReconstruction},
+function rhs!(env, rhs::RHS, hrsc::Union{BernsteinReconstruction, WENOReconstruction},
bdryconds, ldg_bdryconds::Nothing, av_bdryconds::Nothing)
- @unpack cache, mesh, equation, rsolver = rhs
+ @unpack cache, mesh = env
+ @unpack equation, rsolver = rhs
@unpack flag = get_cell_variables(cache)
@unpack flx_u = get_static_variables(cache)
@unpack u = get_dynamic_variables(cache)
@@ -163,7 +63,7 @@ function rhs!(rhs::RHS, t, hrsc::Union{BernsteinReconstruction, WENOReconstructi
broadcast_volume!(flux, equation, cache)
broadcast_faces!(lax_friedrich_flux, rsolver, cache, mesh)
- broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, t)
+ broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, 0.0)
compute_rhs_weak_form!(rhs_u, flx_u, lhs_numflx_u, rhs_numflx_u, mesh)
@@ -171,10 +71,11 @@ function rhs!(rhs::RHS, t, hrsc::Union{BernsteinReconstruction, WENOReconstructi
end
-function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractArtificialViscosity,
+function rhs!(env, rhs::RHS, hrsc::HRSC.AbstractArtificialViscosity,
bdryconds, ldg_bdryconds, av_bdryconds)
- @unpack cache, mesh, equation, rsolver, ldg_rsolver, av_rsolver = rhs
+ @unpack cache, mesh = env
+ @unpack equation, rsolver, ldg_rsolver, av_rsolver = rhs
@unpack flx_u, flx_q, q, v = get_static_variables(cache)
@unpack u = get_dynamic_variables(cache)
@@ -186,7 +87,7 @@ function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractArtificialViscosity,
## solve auxiliary equation: q + ∂x u = 0
broadcast_volume!(ldg_flux, equation, cache)
broadcast_faces!(central_flux, ldg_rsolver, cache, mesh)
- broadcast_boundaryconditions!(central_flux, ldg_bdryconds, cache, mesh, t)
+ broadcast_boundaryconditions!(central_flux, ldg_bdryconds, cache, mesh, 0.0)
compute_rhs_weak_form!(q, flx_q, lhs_numflx_q, rhs_numflx_q, mesh)
## compute rhs of regularized equation: ∂t u + ∂x f + ∂x μ q = 0
@@ -194,10 +95,10 @@ function rhs!(rhs::RHS, t, hrsc::HRSC.AbstractArtificialViscosity,
# @. flx_q = smoothed_mu * q
broadcast_volume!(av_flux, equation, cache)
broadcast_faces!(mean_flux, av_rsolver, cache, mesh)
- broadcast_boundaryconditions!(mean_flux, av_bdryconds, cache, mesh, t)
+ broadcast_boundaryconditions!(mean_flux, av_bdryconds, cache, mesh, 0.0)
broadcast_volume!(flux, equation, cache)
broadcast_faces!(lax_friedrich_flux, rsolver, cache, mesh)
- broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, t)
+ broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, 0.0)
# # add up fluxs and numerical fluxes
@. flx_u += flx_q
diff --git a/src/projects/ScalarEquation/setup.jl b/src/projects/ScalarEquation/setup.jl
new file mode 100644
index 0000000000000000000000000000000000000000..7d827a637c2daa12a2473aa40597312ceb90b35a
--- /dev/null
+++ b/src/projects/ScalarEquation/setup.jl
@@ -0,0 +1,111 @@
+function RHS(env::Environment, prms)
+
+ # construct RHS
+ equation = make_Equation(prms["ScalarEquation"])
+ tci = TCI.make_TCI(env.mesh, prms["TCI"])
+ hrsc = HRSC.make_HRSC(env.mesh, prms["HRSC"])
+ rsolver = ApproxRiemannSolver(flux, speed, equation)
+ ldg_rsolver, av_rsolver = if hrsc isa AbstractArtificialViscosity
+ ApproxRiemannSolver(ldg_flux, ldg_speed, equation),
+ ApproxRiemannSolver(av_flux, av_flux, equation)
+ else
+ nothing, nothing
+ end
+
+ rhs = RHS(equation, rsolver, hrsc, tci, ldg_rsolver, av_rsolver)
+
+ # register variables
+ # TODO add a ::Nothing overload for register_variables!
+ # TODO Somehow replace _register_variables! with register_variables!
+ @unpack cache = env
+ _register_variables!(cache, rhs)
+ register_variables!(cache, rsolver)
+ !isnothing(ldg_rsolver) && register_variables!(cache, ldg_rsolver)
+ !isnothing(av_rsolver) && register_variables!(cache, av_rsolver, dont_register=true)
+ display(cache)
+
+ # setup boundary conditions
+ bdryconds = make_BoundaryConditions(env, equation, rsolver, prms["ScalarEquation"])
+ ldg_bdryconds, av_bdryconds = if hrsc isa AbstractArtificialViscosity
+ make_BoundaryConditions_LDG(env, equation, ldg_rsolver, prms["ScalarEquation"]),
+ make_BoundaryConditions_LDG(env, equation, av_rsolver, prms["ScalarEquation"])
+ else
+ nothing, nothing
+ end
+
+ # TODO Bdry conditions should be moved to project and should only be just a 'lightweight'
+ # descriptor. Furthermore, they should also register variables in cache themselves.
+ # Note: Atm, only the 'from_initial_data' option would require a place to store the data
+ # somewhere. For periodic BCs we don't need extra storage at all.
+ # But what other types of boundary conditions will be needed in the future?
+
+ # setup initial data
+ initialdata!(env, rhs, prms["ScalarEquation"])
+
+ # setup callbacks
+ rhscb = make_callback(env, rhs, bdryconds)
+
+ append!(env.callbacks, CallbackSet(rhscb.callbacks...))
+
+ return rhs, (bdryconds, ldg_bdryconds, av_bdryconds)
+end
+
+
+#######################################################################
+# Register variables #
+#######################################################################
+
+
+# Julia limitation/bug: Overloading imported functions without their Module prefix shadows their
+# definition in this module without a warning or error. This does not happen for methods
+# defined in Base.
+# Because of that we prefix with _ below.
+# Cf. https://github.com/JuliaLang/julia/issues/15510
+# Perhaps this is the import vs. using confusion people talk on the forums.
+function _register_variables!(cache, rhs::RHS)
+ _register_variables!(cache, rhs.equation)
+ _register_variables!(cache, rhs.hrsc)
+ _register_variables!(cache, rhs.tci)
+end
+
+
+function _register_variables!(cache, eq::AbstractScalarEquation)
+ register_variables!(cache,
+ dynamic_variablenames = (:u,),
+ rhs_variablenames = (:rhs_u,),
+ static_variablenames = (:flx_u, :E, :Em1, :F, :Fm1, :v),
+ cell_variablenames = (:max_v,),
+ global_variablenames = (:t, :tm1))
+end
+
+
+function _register_variables!(cache, hrsc::HRSC.AbstractArtificialViscosity)
+ register_variables!(cache,
+ static_variablenames = (:q, :flx_q),
+ cell_variablenames = (:mu,))
+end
+
+
+function _register_variables!(cache, hrsc::HRSC.SmoothedArtificialViscosity)
+ _register_variables!(cache, hrsc.av)
+ register_variables!(cache, static_variablenames = (:smoothed_mu,))
+end
+
+
+function _register_variables!(cache, hrsc::Maybe{HRSC.AbstractHRSC}) end
+
+
+function _register_variables!(cache, tci::Nothing) end
+
+
+function _register_variables!(cache, tci::TCI.AbstractTCI)
+ register_variables!(cache, cell_variablenames = (:flag,))
+end
+
+
+function _register_variables!(cache, tci::TCI.EntropyProduction)
+ register_variables!(cache,
+ cell_variablenames = (:flag,),
+ static_variablenames = (:EP,),
+ bdry_variablenames = (:EPjump_lhs, :EPjump_rhs))
+end
diff --git a/src/projects/ScalarEquation/types.jl b/src/projects/ScalarEquation/types.jl
new file mode 100644
index 0000000000000000000000000000000000000000..d29cb5e085e5e5cefa4f1f7f31f8fe768b51a548
--- /dev/null
+++ b/src/projects/ScalarEquation/types.jl
@@ -0,0 +1,37 @@
+"""
+ abstract type AbstractScalarEquation <: dg1d.AbstractEquation end
+
+Supertype for ScalarEquation problems.
+New problems should subtype it and also implement the following methods:
+- `make_Equation(prms, ::Type{NewProblem})`: factory for an instance of type `NewProblem`.
+The following methods have to be defined with the `@with_signature` interface:
+- `speed(args, ::NewProblem)`
+- `flux(args, ::NewProblem)`
+- `entropy(args, ::NewProblem)`
+- `entropy_flux(args, ::NewProblem)`
+"""
+abstract type AbstractScalarEquation <: dg1d.AbstractEquation end
+
+
+struct Advection <: AbstractScalarEquation
+ v::Float64
+end
+
+
+struct Burgers <: AbstractScalarEquation end
+
+
+struct RHS{T_Equation <:AbstractScalarEquation,
+ T_RSolver <:AbstractRiemannSolver,
+ T_HRSC <:Maybe{HRSC.AbstractHRSC},
+ T_TCI <:Maybe{TCI.AbstractTCI},
+ T_LDG_RSolver <:Maybe{AbstractRiemannSolver},
+ T_AV_RSolver <:Maybe{AbstractRiemannSolver}} <: dg1d.AbstractProject
+
+ equation::T_Equation
+ rsolver::T_RSolver
+ hrsc::T_HRSC
+ tci::T_TCI
+ ldg_rsolver::T_LDG_RSolver
+ av_rsolver::T_AV_RSolver
+end
diff --git a/src/projects/TCI/interface.jl b/src/projects/TCI/interface.jl
index ede8647e72b08f1e8560503e3e562cac701572ef..c361aa16b77199f0dc21957ca489f83133e02900 100644
--- a/src/projects/TCI/interface.jl
+++ b/src/projects/TCI/interface.jl
@@ -21,14 +21,25 @@ abstract type AbstractTCI end
#######################################################################
-function make_TCI(parameters, mesh::Mesh)
- type = Symbol(get(parameters, "type", missing))
- concretes = [ Base.typename(ct).name for ct in concretetypes(AbstractTCI) ]
- if type === :missing
+function make_TCI(mesh::Mesh, prms)
+ @unpack method = prms
+
+ if method == "none"
return nothing
+ elseif method == "mda"
+ return ModalDecayAverage(mesh, prms["mda_threshold_min"], prms["mda_threshold_max"])
+ elseif method == "mdh"
+ return ModalDecayHighest(mesh, prms["mdh_threshold_min"], prms["mdh_threshold_max"])
+ elseif method == "diffusion"
+ return Diffusion(mesh, prms["diffusion_max"])
+ elseif method == "minmod"
+ return Minmod(mesh, prms["minmod_threshold"])
+ elseif method == "threshold"
+ return Threshold(mesh, prms["threshold"])
+ elseif method == "entropy_production"
+ return EntropyProduction(mesh, prms["entropy_production_max"])
+ else
+ error("TCI: Unknown method requested: '$method'")
end
- @assert type in concretes "TCI parameter 'type = $type' unknown. Available options: $concretes"
- tci_type = eval(type)
- tci_prms = parse_parameters(tci_type, parameters)
- return tci_type(; mesh, tci_prms...)
+
end
diff --git a/src/projects/projects.jl b/src/projects/projects.jl
new file mode 100644
index 0000000000000000000000000000000000000000..d4227c6e53c2b3b1d755cee7206b3b7bede324ac
--- /dev/null
+++ b/src/projects/projects.jl
@@ -0,0 +1,178 @@
+#######################################################################
+# Interface #
+#######################################################################
+
+
+"""
+ abstract type AbstractProject end
+
+The supertype for the Project interface of `dg1d`.
+
+Subtypes must implement the constructor
+```
+function (::Type{<:AbstractProject})(env::Environment, prms)
+ # ...
+end
+```
+
+---
+
+# Project interface
+
+How to add a project called `MyProject`:
+
+0. Setup a submodule in `src/projects/MyProject/MyProject.jl` called `MyProject`.
+1. `include` the project in `src/projects/projects.jl`.
+2. Define a subtype of `AbstractProject`. This will serve as the project's dispatcher.
+3. Overload the following methods:
+ - `dg1d.project_type`
+ - `dg1d.parameters`
+ - `dg1d.check_parameters`
+ - `dg1d.rhs!`
+ - `dg1d.timestep!`
+
+---
+
+# Example
+
+```
+# src/project/MyProject/MyProject.jl
+
+module MyProject
+
+using dg1d
+
+# subtype AbstractProject
+struct Project <: dg1d.AbstractProject
+ # ...
+end
+
+function MyProject.Project(env::Environment, prms)
+ # This constructor should
+ # - register variables in cache(env)
+ # - set up initial data for this project
+ # - set up boundary conditions for this project
+ # - set up callbacks for this project
+end
+
+dg1d.project_type(::Val{:MyModule}) = MyProject.Project
+
+dg1d.@parameters MyModule begin
+ ...
+end
+
+function dg1d.rhs!(env::dg1d.Environment, p::MyProject.Project)
+ # compute the rhs inplace in `env`
+end
+
+function dg1d.timestep!(env::dg1d.Environment, p::MyProject.Project)
+ # compute a timestep based on the state of `env`
+end
+
+end # MyModule
+```
+"""
+abstract type AbstractProject end
+
+function (::Type{T})(env::Environment, prms) where T<:AbstractProject
+ error("Interface: Project type '$T' does not implement a constructor!")
+end
+
+
+"""
+ project_type(::Val{T}) where T -> Type{S} where S<:AbstractProject
+
+Register a module `T` with a dispatch type `S<:AbstractProject` in the dg1d
+interface. `T` will be the also the name for the parameter section in parameter files.
+"""
+function project_type(::Val{T}) where T
+ error("Interface: Project '$T' does not provide a dispatch type")
+end
+
+
+"""
+ parameters(::Val{T}) where T -> Dict{String,Any}
+
+Return a dictionary with the default parameters for a project `T`.
+
+Do not directly overload this method. Instead use `dg1d.@parameters` to combine
+docstring generation with parameter definitions.
+
+Parameters defined here can be modified from parameter files.
+E.g. if project `MyProject` defines a parameter called `name`, you can modify it by
+```
+# parameter_file.toml
+[MyProject]
+name = 1
+```
+
+See: @parameters
+---
+
+# List of all parameters and their defaults available for inclusion in parameter files.
+
+"""
+function parameters(::Val{T}) where T
+ error("Interface: Project '$T' does not provide default parameters")
+end
+parameters(s::Symbol) = parameters(Val(s))
+
+
+"""
+ check_parameters(::Val{T}, prms) where T -> nothing
+
+Check consistency of parameters from `parameters(:MyProject)` and throw and
+error if any check is violated.
+
+Do not directly overload this method. Instead use `dg1d.@parameters` to combine
+docstring generation with parameter definitions and optional checks.
+
+See: @parameters
+"""
+function check_parameters(::Val{T}, prms) where T end
+check_parameters(s::Symbol, prms) = check_parameters(Val(s), prms)
+
+
+"""
+ rhs!(env::Environment, p::AbstractProject)
+
+Compute the rhs for project `p` inplace in `cache(env)`.
+"""
+function rhs!(env::Environment, p::AbstractProject, bdryconds)
+ # TODO Remove bdryconds from this interface!!!
+ error("Interface: Project type '$(typeof(p))' does not implement dg1d.rhs!")
+end
+
+
+# TODO Should a timestep computation be allowed to modify `cache(env)`?
+"""
+ timestep!(env::Environment, p::AbstractProject) -> Float64
+
+Return a time step for project `p`.
+
+The function is marked as inplace (note the `!`), because time step computations are allowed
+to modify `cache(env)`.
+"""
+function timestep!(env::Environment, p::AbstractProject)
+ error("Interface: Project type '$(typeof(p))' does not implement dg1d.timestep!")
+end
+
+
+function plot(env::Environment, p::AbstractProject)
+ nothing
+end
+
+
+#######################################################################
+# Helpers #
+#######################################################################
+
+
+function isregistered(p::Symbol)
+ try
+ project_type(Val(p))
+ return true
+ catch
+ return false
+ end
+end
diff --git a/src/riemannsolver.jl b/src/riemannsolver.jl
index 1d857c01ce021a0aa8e0c3aff0137d52ed0a763b..5721d40653a948193f161f37305c83fdc08c55a4 100644
--- a/src/riemannsolver.jl
+++ b/src/riemannsolver.jl
@@ -9,7 +9,7 @@ export AbstractRiemannSolver,
"""
abstract AbstractRiemannSolver
-
+
Supertype for `RiemannSolver` interface.
Any Subtype `Solver <: RiemannSolver` must implement the following methods
- `register_variables!(cache, rsolver::Solver)`,
@@ -40,7 +40,7 @@ Assumes an `equation::SomeEquation` that implements two `@with_signature` functi
# Example
```julia
-julia> mesh = make_Mesh(); cache = Cache(mesh);
+julia> mesh = Mesh(); cache = Cache(mesh);
julia> rsolver = ApproxRiemannSolver(flux, speed, equation)
diff --git a/src/types.jl b/src/types.jl
new file mode 100644
index 0000000000000000000000000000000000000000..0a9f57e35c978a56bd0142c87eeba5f49d1633d2
--- /dev/null
+++ b/src/types.jl
@@ -0,0 +1,9 @@
+struct Environment
+ mesh::Mesh
+ cache::Cache
+ callbacks::CallbackSet
+ callbacks_substep::CallbackSet
+end
+
+
+const Maybe{T} = Union{Nothing,T}
diff --git a/src/utils.jl b/src/utils.jl
index 3303d905fcbe8533644facfc716cb6deccad74d4..3c1eb202764c2d6db7ae077436a9ef95b4af8367 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -519,6 +519,16 @@ end
to_dict(df::dg1d.DataFrameRow) = to_dict(Dict, df)
+"""
+ to_dict([::Type{T}], pairs::Vector{Pair}) where T<:AbstractDict
+
+Convert a `pairs` into a Dictionary of optional type `T`.
+"""
+function to_dict(::Type{T}, pairs::Vector{<:Pair}) where T<:AbstractDict
+ return T(string(p.first) => p.second for p in pairs)
+end
+to_dict(pairs::Vector{Pair}) = to_dict(Dict, pairs)
+
"""
dirname_path(path)
diff --git a/test/IntegrationTests/.gitignore b/test/IntegrationTests/.gitignore
index c291ddf0d5a1db13c47d75f407d0b0157df2e6b2..e6c08d869a7e35f78ad1b541c846bcd40fb9ba2f 100644
--- a/test/IntegrationTests/.gitignore
+++ b/test/IntegrationTests/.gitignore
@@ -1,4 +1,4 @@
tests
Manifest.toml
-refs/*.toml
-refs/*.previous
+refs/*
+refs/*/*
diff --git a/test/IntegrationTests/refs/advection_bugner_gaussian/advection_bugner_gaussian.toml b/test/IntegrationTests/refs/advection_bugner_gaussian/advection_bugner_gaussian.toml
deleted file mode 100644
index c4349a96c07f89c581d3ab6f9fc537b64845516e..0000000000000000000000000000000000000000
--- a/test/IntegrationTests/refs/advection_bugner_gaussian/advection_bugner_gaussian.toml
+++ /dev/null
@@ -1,27 +0,0 @@
-[Equation]
-type = "Advection"
-v = 1.0
-
-[Mesh]
-xrange = [-3.0, 3.0]
-CFL = 0.5
-N = 3
-K = 128
-
-[Output]
-variables = [ "u" ]
-aligned_times = [ 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
- 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0 ]
-
-[Grid]
-type = "LGL_lumped"
-
-[Evolution]
-tend = 2.0
-
-[Initialdata]
-type = "bugner gaussian"
-
-[TCI]
-
-[HRSC]
diff --git a/test/IntegrationTests/refs/advection_bugner_gaussian/output.h5 b/test/IntegrationTests/refs/advection_bugner_gaussian/output.h5
deleted file mode 100644
index 7f0832defa24f7f5ca40207d8b6008d3565b9315..0000000000000000000000000000000000000000
Binary files a/test/IntegrationTests/refs/advection_bugner_gaussian/output.h5 and /dev/null differ
diff --git a/test/IntegrationTests/refs/advection_sine/advection_sine.toml b/test/IntegrationTests/refs/advection_sine/advection_sine.toml
index 5146418634e54080d98117d00d0c09f4954225bb..6ca8efc3908c78caeb04e1375fbc605d3a6acc46 100644
--- a/test/IntegrationTests/refs/advection_sine/advection_sine.toml
+++ b/test/IntegrationTests/refs/advection_sine/advection_sine.toml
@@ -1,31 +1,17 @@
-[Equation]
-type = "Advection"
-v = 1.0
+[ScalarEquation]
+equation = "advection"
+id = "sine"
[Mesh]
xrange = [ 0.0, 1.0 ]
-CFL = 0.5
-N = 4
-K = 40
-
-[Plot]
-every_dt = 0
-every_iteration = 0
-pause = false
+cfl = 0.5
+n = 4
+k = 40
+basis = "lgl_lumped"
[Output]
-variables = [ "u" ]
-aligned_times = [ 0.05, 0.1, 0.15, 0.2 ]
-
-[Grid]
-type = "LGL_lumped"
+variables = [ "u" ]
+aligned_ts = [ 0.05, 0.1, 0.15, 0.2 ]
[Evolution]
tend = 2.2
-
-[Initialdata]
-type = "sine"
-
-[TCI]
-
-[HRSC]
diff --git a/test/IntegrationTests/refs/advection_sine/output.h5 b/test/IntegrationTests/refs/advection_sine/output.h5
index 2674833a63fad397ff887c9b57ce5d20ac29a63e..bc5587cdc75673cb14118d15cf8a2a77f3920896 100644
Binary files a/test/IntegrationTests/refs/advection_sine/output.h5 and b/test/IntegrationTests/refs/advection_sine/output.h5 differ
diff --git a/test/IntegrationTests/refs/advection_sine_no_inflow/advection_sine_no_inflow.toml b/test/IntegrationTests/refs/advection_sine_no_inflow/advection_sine_no_inflow.toml
index ff85ba000cc3f335a8dad6d224034660a770b285..456f2fd5313b1a1175a61872555943bb27246a35 100644
--- a/test/IntegrationTests/refs/advection_sine_no_inflow/advection_sine_no_inflow.toml
+++ b/test/IntegrationTests/refs/advection_sine_no_inflow/advection_sine_no_inflow.toml
@@ -1,34 +1,23 @@
-[Equation]
-type = "Advection"
-v = 1.0
+[ScalarEquation]
+equation = "advection"
+advection_v = 1.0
+id = "sine"
+bc = "no_inflow"
[Mesh]
xrange = [ 0.0, 1.0 ]
-CFL = 0.5
-N = 3
-K = 80
-
-[Plot]
-every_dt = 0
-every_iteration = 0
-pause = false
+cfl = 0.5
+n = 3
+k = 80
+basis = "lgl"
[Output]
every_iteration = 0
variables = [ "u" ]
-aligned_times = [ 0.05, 0.1, 0.15, 0.2 ]
-
-[Grid]
-type = "LGL_lumped"
+aligned_ts = [ 0.05, 0.1, 0.15, 0.2 ]
[Evolution]
tend = 0.2
-[Initialdata]
-type = "sine"
-
-[TCI]
-
-[HRSC]
-[BoundaryConditions]
-type = "no inflow"
+# [HRSC]
+# method = "av"
diff --git a/test/IntegrationTests/refs/burgers_gaussian/burgers_gaussian.toml b/test/IntegrationTests/refs/burgers_gaussian/burgers_gaussian.toml
index 8a2f1d1009978e55575c9b6d36887769cf1752c8..69e6c63f2de4fdabd864718f164e0377026df186 100644
--- a/test/IntegrationTests/refs/burgers_gaussian/burgers_gaussian.toml
+++ b/test/IntegrationTests/refs/burgers_gaussian/burgers_gaussian.toml
@@ -1,34 +1,24 @@
-[Equation]
-type = "Burgers"
+[ScalarEquation]
+equation = "burgers"
+id = "gaussian"
+id_gaussian_width = 0.4
+id_gaussian_amplitude = 1
[Mesh]
xrange = [ -1.4, 1.4 ]
-CFL = 0.10
-N = 5
-K = 41
-
-[Plot]
-every_dt = 0
-every_iteration = 0
-pause = false
+cfl = 0.10
+n = 5
+k = 41
+basis = "lgl"
[Output]
every_iteration = 0
every_dt = 0
variables = [ "u" ]
-aligned_times = [ 1e-10, 0.02, 0.04, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16, 0.18, 0.20 ]
-
-[Grid]
-type = "LGL_lumped"
+aligned_ts = [ 1e-10, 0.02, 0.04, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16, 0.18, 0.20 ]
[Evolution]
tend = 0.2
-[Initialdata]
-type = "gaussian"
-width = 0.4
-amplitude = 1
-
-[TCI]
-
[HRSC]
+method = "none"
diff --git a/test/IntegrationTests/refs/burgers_sine/burgers_sine.toml b/test/IntegrationTests/refs/burgers_sine/burgers_sine.toml
index a469505eb9891ce68f28dc497b6798273189eecd..16a69fec6ea9aeea6ad3b44add57067e9b6907fd 100644
--- a/test/IntegrationTests/refs/burgers_sine/burgers_sine.toml
+++ b/test/IntegrationTests/refs/burgers_sine/burgers_sine.toml
@@ -1,33 +1,20 @@
-[Equation]
-type = "Burgers"
+[ScalarEquation]
+equation = "burgers"
+id = "sine"
[Mesh]
xrange = [ 0, 1 ]
-CFL = 0.10
-N = 2
-K = 160
-
-[Plot]
-every_dt = 0
-every_iteration = 0
-pause = false
+cfl = 0.10
+n = 2
+k = 160
+basis = "lgl"
[Output]
every_iteration = 0
every_dt = 0
variables = [ "u" ]
-aligned_times = [ 0.02, 0.04, 0.06, 0.08, 0.10, 0.12, 0.14 ]
-
-[Grid]
-type = "LGL"
+aligned_ts = [ 0.02, 0.04, 0.06, 0.08, 0.10, 0.12, 0.14 ]
[Evolution]
# stop before discontinuity appears
tend = 0.14
-
-[Initialdata]
-type = "sine"
-
-[TCI]
-
-[HRSC]
diff --git a/test/IntegrationTests/refs/euler_isentropic_flow/euler_isentropic_flow.toml b/test/IntegrationTests/refs/euler_isentropic_flow/euler_isentropic_flow.toml
index 982431ecf1c740630608d7e244b772e6a5d4e360..f45d16a01160b611e35efb2ee286f2228322b3ea 100644
--- a/test/IntegrationTests/refs/euler_isentropic_flow/euler_isentropic_flow.toml
+++ b/test/IntegrationTests/refs/euler_isentropic_flow/euler_isentropic_flow.toml
@@ -1,32 +1,21 @@
+[EulerEq]
+id = "smooth_isentropic_flow"
+
[EquationOfState]
-type = "IdealGas"
-gamma = 3.0
+eos = "idealgas"
+idealgas_gamma = 3.0
[Mesh]
xrange = [ -1.0, 1.0 ]
-CFL = 0.10
-N = 3
-K = 21
-
-[Plot]
-every_dt = 0
-every_iteration = 0
-pause = false
+cfl = 0.10
+n = 3
+k = 21
+basis = "lgl"
[Output]
every_iteration = 0
every_dt = 0.02
variables = [ "rho", "q", "E" ]
-[Grid]
-type = "LGL_lumped"
-
[Evolution]
tend = 0.25
-
-[Initialdata]
-type = "smooth_isentropic_flow"
-
-[TCI]
-
-[HRSC]
diff --git a/test/IntegrationTests/refs/euler_sod_shock_tube/euler_sod_shock_tube.toml b/test/IntegrationTests/refs/euler_sod_shock_tube/euler_sod_shock_tube.toml
index c327fe44c5e8e901c53f8de86aa08fdce0075851..5693d7c15c7e8634202789a17f626abd683ce04a 100644
--- a/test/IntegrationTests/refs/euler_sod_shock_tube/euler_sod_shock_tube.toml
+++ b/test/IntegrationTests/refs/euler_sod_shock_tube/euler_sod_shock_tube.toml
@@ -1,20 +1,17 @@
[EquationOfState]
-type = "IdealGas"
-gamma = 1.4
+eos = "idealgas"
+idealgas_gamma = 1.4
-[BoundaryConditions]
-type = "from initial data"
+[EulerEq]
+bc = "from_id"
+id = "sod_shock_tube"
[Mesh]
xrange = [ -1.0, 1.0 ]
-CFL = 0.10
-N = 6
-K = 51
-
-[Plot]
-every_dt = 0
-every_iteration = 0
-pause = false
+cfl = 0.10
+n = 6
+k = 51
+basis = "lgl"
[Output]
# atm we are happy when it does not crash
@@ -22,17 +19,12 @@ every_iteration = 0
every_dt = 0
variables = [ "rho", "q", "E" ]
-[Grid]
-type = "LGL_lumped"
-
[Evolution]
tend = 0.6
-[Initialdata]
-type = "sod_shock_tube"
-
[TCI]
-type = "ModalDecayAverage"
+method = "mda"
[HRSC]
-type = "BernsteinReconstruction"
+method = "filter"
+filter_method = "bernstein"
diff --git a/test/IntegrationTests/refs/sreuler_smooth/sreuler_smooth.toml b/test/IntegrationTests/refs/sreuler_smooth/sreuler_smooth.toml
index 9e9cba0fa6e32befee1430bafafffe305769bc7e..9f295fdfec732befe49b81b54b82edf2009a778e 100644
--- a/test/IntegrationTests/refs/sreuler_smooth/sreuler_smooth.toml
+++ b/test/IntegrationTests/refs/sreuler_smooth/sreuler_smooth.toml
@@ -1,36 +1,24 @@
[EquationOfState]
-type = "IdealGas"
-gamma = 1.6666666666666666666667
-
-[Equation]
+eos = "idealgas"
+idealgas_gamma = 1.6666666666666666666667
+[SREulerEq]
+id = "smooth"
[Mesh]
xrange = [ -1.0, 1.0 ]
-CFL = 0.10
-N = 4
-K = 67
-
-[Plot]
-every_dt = 0
-every_iteration = 0
-pause = false
+cfl = 0.10
+n = 4
+k = 67
+basis = "lgl"
[Output]
every_iteration = 0
every_dt = 0.05
variables = [ "D", "S", "tau" ]
-[Grid]
-type = "LGL_lumped"
-
[Evolution]
tend = 0.4
-[Initialdata]
-type = "smooth"
-
[TCI]
-type = "ModalDecayAverage"
-
-[HRSC]
+method = "mda"
diff --git a/test/IntegrationTests/refs/testconfig.toml b/test/IntegrationTests/refs/testconfig.toml
index c0a4ddae70c5ba1bb1f49fb32cdd8628ce8662af..50c7bb7684d583eb62514b630452db5ce5b0b3a4 100644
--- a/test/IntegrationTests/refs/testconfig.toml
+++ b/test/IntegrationTests/refs/testconfig.toml
@@ -12,9 +12,6 @@ variables = [ "u" ]
reltol = 1e-6
abstol = 1e-4
-[advection_bugner_gaussian]
-variables = [ "u" ]
-
[burgers_gaussian]
# stops before discontinuity appears
variables = [ "u" ]
diff --git a/test/IntegrationTests/src/IntegrationTests.jl b/test/IntegrationTests/src/IntegrationTests.jl
index 473f22237296035ea3116a9a22b07e3d0b9663c9..48ff8a669956b97bf471c3faa194a92015ed6167 100644
--- a/test/IntegrationTests/src/IntegrationTests.jl
+++ b/test/IntegrationTests/src/IntegrationTests.jl
@@ -60,32 +60,16 @@ end
function run_pars(pars)
- # TODO This is a dirty hack needed because projects are not yet structured.
- # Will be fixed soon.
- modules = Dict("grhd" => :GRHD,
- "advection" => :ScalarEquation,
- "burgers" => :ScalarEquation,
- "euler" => :EulerEq,
- "sreuler" => :SREulerEq,
- "testy" => :GRHD)
success = Bool[]
dirs = String[]
for p in pars
@info "Running parameter file '$p' ..."
- m = match(r"(.*?)_.*\.toml", basename(p))
- @assert !isnothing(m)
- prefix = m.captures[1]
- mod = modules[prefix]
- @assert isdefined(dg1d, mod)
- cmd = quote
- try
- $mod.main($p)
- true
- catch
- false
- end
+ s = try
+ dg1d.main(p)
+ true
+ catch
+ false
end
- s = @eval $cmd
push!(success, s)
end
@@ -117,8 +101,8 @@ function compare_outputs(io, refoutput, output,
if length(refts) != length(ts)
println(io, """
FAIL: Mismatch in number of time steps for which outputs were recorded
- reference $(length(refts))
- test $(length(ts))
+ reference $(@sprintf("%10d",length(refts))) : $refts
+ test $(@sprintf("%10d", length(ts))) : $ts
""")
return false
end
@@ -157,7 +141,7 @@ function compare_outputs(io, refoutput, output,
if !all(zip(refvi, vi)) do (rv, v)
isapprox(rv, v, rtol=reltol, atol=abstol)
end
- println(io, "FAIL: Mismatch in variable '$var' @ (t,idx) = ($(@sprintf("%10.10e",t)),$idx)")
+ println(io, "FAIL: Mismatch in variable '$var' @ (t,i) = ($(@sprintf("%10.10e",ts[i])),$i), norm(ref-test,1) = $(norm(refvi.-vi,1))")
success = false
end
end
@@ -168,14 +152,11 @@ end
"""
- runtests()
runtests(tests::Vector{String})
Run integration tests for which there is reference data
in `refs/` and a configuration entry in `refs/testconfig.toml`.
"""
-runtests() = runtests(String[])
-
function runtests(tests::Vector{String}=String[])
allrefdirs, allrefpars, allrefoutputs = discover_output_dirs(REFDIR)
@@ -202,7 +183,6 @@ function runtests(tests::Vector{String}=String[])
display(testconfig)
for (registered_test, pars) in testconfig
- display(pars)
if !(pars isa AbstractDict)
error("Invalid testconfig.toml: Found unsupported parameter '$registered_test'")
end
@@ -256,9 +236,9 @@ function runtests(tests::Vector{String}=String[])
println(logs, repeat('+', n_header), "\n\nTest: $test\n")
if !run_success[i]
- println(log, "FAIL: Failed to run test!")
+ println(logs, "FAIL: Failed to run test!")
push!(comparsion_success, missing)
- println(logs, "\nRESULT: Integration test failed\n")
+ println(logs, "\nResult: FAILED\n")
continue
end
@@ -294,9 +274,9 @@ function runtests(tests::Vector{String}=String[])
Summary
-Tests: $(@sprintf("%3d/%-3d = %3.1f%%", n_runs_succeeded, n_runs, pct_runs_succeeded))
-Comparsions: $(@sprintf("%3d/%-3d = %3.1f%%", n_compars_succeeded, n_compars, pct_compars_succeeded))
-Final result: $(final_result ? "PASSED" : "FAILED")
+Tests succeded: $(@sprintf("%3d/%-3d = %3.1f%%", n_runs_succeeded, n_runs, pct_runs_succeeded))
+Comparsions succeded: $(@sprintf("%3d/%-3d = %3.1f%%", n_compars_succeeded, n_compars, pct_compars_succeeded))
+Final result: $(final_result ? "PASSED" : "FAILED")
""")
logfilename = joinpath(TESTDIR, "tests.log")
diff --git a/test/dg1dTests/src/test_EquationOfState.jl b/test/dg1dTests/src/test_EquationOfState.jl
index 2adb5920191769024bb0ede9d7d9a3804ab08a43..9cbec77c27fe09f3379710e14f4a8ac9a1bc2691 100644
--- a/test/dg1dTests/src/test_EquationOfState.jl
+++ b/test/dg1dTests/src/test_EquationOfState.jl
@@ -6,8 +6,8 @@
# Call all methods we expect an eqaution of state to implement at least once here.
for EOS in concretetypes(AbstractEquationOfState)
- parameters = parse_parameters(EOS)
- eos = EOS(parameters...)
+ parameters = dg1d.parameters(:EquationOfState)
+ eos = EquationOfState.make_EquationOfState(parameters)
rho = 1.0
eps = 1.0
diff --git a/test/dg1dTests/src/test_cache.jl b/test/dg1dTests/src/test_cache.jl
index d798c6780842ed10b59082f133e634f4b319771b..cd2fe738f9e25fd0093d7f76dbad5b952bf5d343 100644
--- a/test/dg1dTests/src/test_cache.jl
+++ b/test/dg1dTests/src/test_cache.jl
@@ -4,7 +4,7 @@
@testset "Constructor and getters" begin
# constructor
- mesh = make_Mesh()
+ mesh = Mesh()
C = dg1d.Cache(mesh)
# getters
@@ -20,7 +20,7 @@ end
@testset "Getters" begin
- mesh = make_Mesh()
+ mesh = Mesh()
C = dg1d.Cache(mesh)
dynamic_variablenames = (:u, :v)
rhs_variablenames = (:rhs_u, :rhs_v)
@@ -68,7 +68,7 @@ end
@testset "Registering" begin
- mesh = make_Mesh()
+ mesh = Mesh()
C = dg1d.Cache(mesh)
dynamic_variablenames = (:u, :v)
@@ -99,12 +99,12 @@ end
@testset "Wrapping" begin
- mesh = make_Mesh()
+ mesh = Mesh()
# register some variables
dynamic_variablenames = (:var1, :var2)
rhs_variablenames = (:rhs_var1, :rhs_var2)
- mesh = make_Mesh()
+ mesh = Mesh()
C = dg1d.Cache(mesh)
register_variables!(C; dynamic_variablenames = dynamic_variablenames,
rhs_variablenames = rhs_variablenames)
@@ -129,7 +129,7 @@ end
# hence, the order in which variables are registered matters!
dynamic_variablenames = (:var2, :var1)
rhs_variablenames = (:rhs_var1, :rhs_var2)
- mesh = make_Mesh()
+ mesh = Mesh()
C = dg1d.Cache(mesh)
register_variables!(C; dynamic_variablenames = dynamic_variablenames,
rhs_variablenames = rhs_variablenames)
@@ -154,7 +154,7 @@ end
@testset "Printing" begin
- mesh = make_Mesh()
+ mesh = Mesh()
# register some variables
dynamic_variablenames = (:var1, :var2)
@@ -163,7 +163,7 @@ end
cell_variablenames = (:cell_speed, :cell_mean, :cell_avg)
static_variablenames = (:entropy, :entropy_flux)
global_variablenames = (:t,)
- mesh = make_Mesh()
+ mesh = Mesh()
C = dg1d.Cache(mesh)
register_variables!(C; dynamic_variablenames = dynamic_variablenames,
rhs_variablenames = rhs_variablenames,
@@ -203,7 +203,7 @@ end
end
# setup cache
- mesh = make_Mesh(Npts=5,K=5)
+ mesh = Mesh(N=5,K=5)
cache = Cache(mesh)
dg1d.register_variables!(cache,
dynamic_variablenames = (:u, :v),
diff --git a/test/dg1dTests/src/test_callbacks.jl b/test/dg1dTests/src/test_callbacks.jl
index b78be19e8fc83e10769e9f4914a38fba4440dd33..a3f8c55e59a7f7800cdfae172c99875cb9f8f71f 100644
--- a/test/dg1dTests/src/test_callbacks.jl
+++ b/test/dg1dTests/src/test_callbacks.jl
@@ -109,7 +109,7 @@ end
ispath(mockdir) && rm(mockdir, recursive=true)
mkdir(mockdir)
- mesh = make_Mesh()
+ mesh = Mesh()
cache = Cache(mesh)
dg1d.register_variables!(cache, static_variablenames=[:u, :v], cell_variablenames=[:a,:b])
@@ -167,7 +167,7 @@ end
ispath(mockdir) && rm(mockdir, recursive=true)
mkdir(mockdir)
- mesh = make_Mesh()
+ mesh = Mesh()
cache = Cache(mesh)
dg1d.register_variables!(cache, dynamic_variablenames=[:u])
@@ -190,7 +190,7 @@ end
cbs = CallbackSet(scb, atscb)
# evolve
- dg1d.evolve_ERK(rhs!, u0, aligned_timestep_fn, (0.0, 5.0), :RK4, callback_fullstep=cbs)
+ dg1d.evolve_ERK(rhs!, u0, aligned_timestep_fn, (0.0, 5.0), :rk4, callback_fullstep=cbs)
# test output
@test isfile(filename)
diff --git a/test/dg1dTests/src/test_cons2prim.jl b/test/dg1dTests/src/test_cons2prim.jl
index 94bbf1967552707ca9baf3dfb606617943b0060a..c82c54b6b66c40e6cda53109ea096bac8bc96de3 100644
--- a/test/dg1dTests/src/test_cons2prim.jl
+++ b/test/dg1dTests/src/test_cons2prim.jl
@@ -2,8 +2,8 @@
F = SREulerEq.Formulae
eos = IdealGas(gamma = 5/3)
- atms = SREulerEq.Atmosphere(Dict())
- equation = SREulerEq.SREulerEquation(eos, atms, NewtonSolver(), Bisection())
+ prms = dg1d.parameters(:SREulerEq)
+ equation = SREulerEq.make_SREulerEquation(eos, prms)
# initial state of relativistic shock tube problem from Bugner's thesis
rhol, vl, pl = 10.0, 0.0, 13.33
diff --git a/test/dg1dTests/src/test_evolve.jl b/test/dg1dTests/src/test_evolve.jl
index 8b5b5605605ae953c0cca9041db864434b8af6cd..b3bd081a2f936a88f66aab223d54c373ecc2e4e5 100644
--- a/test/dg1dTests/src/test_evolve.jl
+++ b/test/dg1dTests/src/test_evolve.jl
@@ -15,7 +15,7 @@
u_ana = u(tend)
- schemes = [ :MIDPT, :RK4, :RK4_twothirds, :RK4_Ralston, :SSPRK3, :LSERK4, :RKF10 ]
+ schemes = [ :midpt, :rk4, :rk4_twothirds, :rk4_ralston, :ssprk3, :lserk4, :rkf10 ]
orders = [ 2, 4, 4, 4, 3, 4, 10 ]
tspan = (t0, tend)
@@ -63,7 +63,7 @@ end
u_ana = u(tend)
- schemes = [ :RK4, :RK4_twothirds, :RK4_Ralston, :SSPRK3, :LSERK4 ]
+ schemes = [ :rk4, :rk4_twothirds, :rk4_ralston, :ssprk3, :lserk4 ]
tspan = (t0, tend)
dts = [0.1, 0.05, 0.01]
@@ -119,7 +119,7 @@ end
called = true
end
- scheme = :RK4
+ scheme = :rk4
tspan = (t0, tend)
u_num = dg1d.evolve_ERK(rhs!, u0, dt, tspan, scheme, callback_fullstep=callback_fullstep)
@@ -151,7 +151,7 @@ end
# dynamic time stepping
dt(u, t, n) = t < 0.5 ? 0.1 : 0.05
- schemes = [ :RK4, :RK4_twothirds, :RK4_Ralston, :SSPRK3, :LSERK4 ]
+ schemes = [ :rk4, :rk4_twothirds, :rk4_ralston, :ssprk3, :lserk4 ]
tspan = (t0, tend)
for scheme in schemes
@@ -197,7 +197,7 @@ end
return base_dt / (n_fails+1)
end
- schemes = [ :RK4, :RK4_twothirds, :RK4_Ralston, :SSPRK3, :LSERK4 ]
+ schemes = [ :rk4, :rk4_twothirds, :rk4_ralston, :ssprk3, :lserk4 ]
tspan = (t0, tend)
for scheme in schemes
test_n_failed = 0 # reset counter
diff --git a/test/dg1dTests/src/test_mesh.jl b/test/dg1dTests/src/test_mesh.jl
index 9b4caa6ad08a1ec3debb732d9a4f229e6056ed9d..a3edc778d27ed063eefb9ff7d8e813e8aee92576 100644
--- a/test/dg1dTests/src/test_mesh.jl
+++ b/test/dg1dTests/src/test_mesh.jl
@@ -3,15 +3,11 @@
@testset "Factories" begin
- # use dictionary
- prms = Dict("N"=>4, "K"=>10)
- mesh = make_Mesh(prms)
-
# or keyword arguments
- mesh = make_Mesh(N=4, K=10)
+ mesh = Mesh(N=4, K=10)
# or just the defaults
- mesh = make_Mesh()
+ mesh = Mesh()
end
@@ -19,7 +15,7 @@ end
@testset "Mesh quadrature" begin
# ∫(0≤x≤2π) sin(2 π x) dx = 0
- mesh = make_Mesh(N=10, K=100)
+ mesh = Mesh(N=10, K=100)
u = @. sin(2.0*Ï€*mesh.x)
grid_avg = grid_average(u, mesh)
result = 0.0
@@ -59,9 +55,8 @@ end
@testset "Utilities" begin
- prms = parse_parameters(Mesh, Dict())
- mesh = Mesh(; prms...)
- @unpack N, K = prms
+ N, K = 5, 10
+ mesh = Mesh(N=N, K=K)
#### dimensions, length, layout
npts, ncells = layout(mesh)
@@ -71,7 +66,7 @@ end
@test len == (N+1)*K
#### find cells given nodes on mesh
- mesh = make_Mesh(; N=2, K=10) # =^= equidistant grid with doubled nodes at interfaces
+ mesh = Mesh(N=2, K=10) # =^= equidistant grid with doubled nodes at interfaces
@unpack L, K = mesh
# as a sample we place a node into to center of every cell
dL = L / K
@@ -122,8 +117,8 @@ end
@testset "MeshInterpolator" begin
- smpl_mesh = make_Mesh(; N=4, K=6)
- ref_mesh = make_Mesh(; N=10, K=20)
+ smpl_mesh = Mesh(N=4, K=6)
+ ref_mesh = Mesh(N=10, K=20)
ref_data = sin.(ref_mesh.x)
meshintrp = MeshInterpolator(smpl_mesh, ref_mesh)
observed = meshintrp(ref_data)
@@ -137,8 +132,8 @@ end
#### edge cases
# sample mesh lies (partly) outside of reference mesh
- smpl_mesh_outside = make_Mesh(; N=4, K=6, xrange=[-10.0,0.0])
- ref_mesh = make_Mesh(; N=10, K=20)
+ smpl_mesh_outside = Mesh(N=4, K=6, xrange=[-10.0,0.0])
+ ref_mesh = Mesh(N=10, K=20)
@test_throws ErrorException MeshInterpolator(smpl_mesh_outside, ref_mesh)
# reference data does not match reference mesh
diff --git a/test/dg1dTests/src/test_output.jl b/test/dg1dTests/src/test_output.jl
index eb7053f8562f525b18917e0c2acf14b0d1e7bcfc..49d78a36b5f853bde5989429a339d23cd7d4b8bf 100644
--- a/test/dg1dTests/src/test_output.jl
+++ b/test/dg1dTests/src/test_output.jl
@@ -6,7 +6,7 @@
ispath(mockdir) && rm(mockdir, recursive=true)
mkdir(mockdir)
- mesh = make_Mesh(; N=4,K=6)
+ mesh = Mesh(; N=4,K=6)
cache = Cache(mesh)
dg1d.register_variables!(cache, static_variablenames=[:u], cell_variablenames=[:a])
diff --git a/test/dg1dTests/src/test_parameters.jl b/test/dg1dTests/src/test_parameters.jl
index 9d45892a6d8275b423510dc20bf52fa75e33d36e..28ddd06edc7fdb4597e1761b1b08b9722a90c74a 100644
--- a/test/dg1dTests/src/test_parameters.jl
+++ b/test/dg1dTests/src/test_parameters.jl
@@ -1,58 +1,85 @@
-abstract type AbstractParameters end
-
-
-dg1d.@option "ParametersA" struct DummyParametersA <: dg1d.AbstractParameters
- "Parameter a"
- a::Int64 = 5
- "Parameter b"
- b::Float64 = 2.4
- "Parameter c"
- c::String = "method c"
-end
-
-
-function sanitize(prms::DummyParametersA)
- @assert prms.a > 0 && prms.a < 6
- @assert prms.b > 1.3
- @assert prms.c in [ "method c", "method cc", "method ccc" ]
-end
-
-
-dg1d.@option "ParametersB" struct DummyParametersB <: dg1d.AbstractParameters
- "Parameters A"
- prmsA::DummyParametersA = DummyParametersA(a = 3, b = 3.4)
-end
-
-
-function sanitize(prms::DummyParametersB)
- sanitize(prms.prmsA)
-end
-
-
@testset "Parameters" begin
-
-@testset "Sanitizing" begin
-
- # sane values
- prmsA = DummyParametersA()
- sanitize(prmsA)
-
- # insane values
- @test_throws AssertionError sanitize(DummyParametersA(a = -1))
- @test_throws AssertionError sanitize(DummyParametersA(b = 0.0))
- @test_throws AssertionError sanitize(DummyParametersA(c = "method xyz"))
-
- # sane values nested parameters
- prmsB = DummyParametersB()
- sanitize(prmsB)
-
- # insane values nested parameters
- prmsA = DummyParametersA(a = -1)
- prmsB = DummyParametersB(prmsA=prmsA)
- @test_throws AssertionError sanitize(prmsB)
-
-end
+ # a working example
+ dg1d.@parameters MyType begin
+ "some parameter"
+ x = 1.0
+ @check x > 0.0
+ @check x > 0.5
+
+ "another parameter"
+ y = "sers"
+
+ "and another parameter"
+ z = 123
+ end
+ prms = dg1d.parameters(:MyType)
+ @test prms["x"] == 1.0
+ @test prms["y"] == "sers"
+ @test prms["z"] == 123
+
+ # verify check functionality
+ dg1d.check_parameters(:MyType, prms)
+ # trigger first check
+ prms["x"] = 0.2
+ @test_throws ErrorException dg1d.check_parameters(:MyType, prms)
+ # trigger second check
+ prms["x"] = -1.0
+ @test_throws ErrorException dg1d.check_parameters(:MyType, prms)
+
+ # test wrong usage
+
+ # missing docstring
+ @test_throws LoadError @eval dg1d.@parameters MyType begin
+ x = 1.0
+ end
+
+ # missing default value
+ @test_throws LoadError @eval dg1d.@parameters MyType begin
+ "some parameter"
+ x
+ end
+
+ # duplicated parameter
+ @test_throws LoadError @eval dg1d.@parameters MyType begin
+ "some parameter"
+ x = 1.0
+ "same parameter again"
+ x = 2.0
+ end
+
+ # duplicated parameter
+ @test_throws LoadError @eval dg1d.@parameters MyType begin
+ "some parameter"
+ x = 1.0
+ "same parameter again"
+ x = 2.0
+ end
+
+ # checks can't appear at the start
+ @test_throws LoadError @eval dg1d.@parameters MyType begin
+ @check x > 0.0
+ "some parameter"
+ x = 1.0
+ "same parameter again"
+ x = 2.0
+ end
+
+ # no other macros or expressions allowed
+ @test_throws LoadError @eval dg1d.@parameters MyType begin
+ "some parameter"
+ x = 1.0
+ @test x > 0.0 # should be @check
+ "same parameter again"
+ x = 2.0
+ end
+ @test_throws LoadError @eval dg1d.@parameters MyType begin
+ "some parameter"
+ x = 1.0
+ x > 2.0 # forgot to add @check
+ "same parameter again"
+ x = 2.0
+ end
end
diff --git a/test/dg1dTests/src/test_riemannsolver.jl b/test/dg1dTests/src/test_riemannsolver.jl
index 598f8f47a95b8b3bd02420196fda22be2753f426..07128003f8b90c8bf92aa196aaf9f622466395c0 100644
--- a/test/dg1dTests/src/test_riemannsolver.jl
+++ b/test/dg1dTests/src/test_riemannsolver.jl
@@ -39,7 +39,7 @@ end
@testset "ApproxRiemannSolver" begin
- mesh = make_Mesh()
+ mesh = Mesh()
cache = Cache(mesh)
equation = MyDummyEquation()
register_variables!(cache, equation)
@@ -72,7 +72,7 @@ end
@testset "BoundaryConditions" begin
- mesh = make_Mesh()
+ mesh = Mesh()
cache = Cache(mesh)
equation = MyDummyEquation()
register_variables!(cache, equation)
diff --git a/test/dg1dTests/src/test_segment.jl b/test/dg1dTests/src/test_segment.jl
index e0d6b0b8b944088e9edff4f0707b04eb74091d4f..910cd7078c020c3b497c68bb6ae9e27fae1e3bdb 100644
--- a/test/dg1dTests/src/test_segment.jl
+++ b/test/dg1dTests/src/test_segment.jl
@@ -1,6 +1,6 @@
@testset "Segment" begin
- mesh = make_Mesh(N=3, K=5)
+ mesh = Mesh(N=3, K=5)
NN = length(mesh)
N, K = layout(mesh)
vars = Variables(NN, :a, :b, :c)