Merge pull request #202 from ranocha/hr/taal

Taal: add trixi_include and move definitions of abstract types for dispatch

src/Trixi.jl

@@ -24,6 +24,7 @@ using DiffEqBase: ODEProblem, ODESolution, get_du, u_modified!, set_proposed_dt!
 using DiffEqCallbacks: CallbackSet, DiscreteCallback
 using EllipsisNotation # ..
 using HDF5: h5open, attrs
+using LinearMaps: LinearMap
 using StaticArrays: @MVector, @SVector, MVector, MMatrix, MArray, SVector, SMatrix, SArray
 using TimerOutputs: @notimeit, @timeit, @timeit_debug, TimerOutput, print_timer, reset_timer!
 using UnPack: @unpack
@@ -40,12 +41,12 @@ const globals = Dict{Symbol, Any}()
 export globals
 
 
-# TODO: Taal, decide where to define the entry points of our type hierarchy,
-# e.g. AbstractEquations, AbstractSemidiscretization etc.
-# Placing them here will allow us to make use of them for dispatch even for
+# Define the entry points of our type hierarchy, e.g.
+#     AbstractEquations, AbstractSemidiscretization etc.
+# Placing them here allows us to make use of them for dispatch even for
 # other stuff defined very early in our include pipeline, e.g.
-#   IndicatorLöhner(semi::AbstractSemidiscretization)
-
+#     IndicatorLöhner(semi::AbstractSemidiscretization)
+include("abstract_types.jl")
 
 # Include all top-level source files
 include("auxiliary/auxiliary.jl")
@@ -59,7 +60,7 @@ include("amr/amr.jl")
 include("callbacks/callbacks.jl")
 include("semidiscretization_euler_gravity.jl")
 
-# TODO: Taal remove run
+# TODO: Taal refactor, get rid of old run methods, rename the file
 # Include top-level run method
 include("run.jl")
 
@@ -100,7 +101,7 @@ export density, pressure, density_pressure
 
 export entropy, energy_total, energy_kinetic, energy_internal
 
-export examples_dir, get_examples, default_example
+export trixi_include, examples_dir, get_examples, default_example
 
 
 end

src/abstract_types.jl

@@ -0,0 +1,12 @@
+
+# abstract supertype of specific semidiscretizations such as
+# - SemidiscretizationHyperbolic for hyperbolic conservation laws
+# - SemidiscretizationEulerGravity for Euler with self-gravity
+abstract type AbstractSemidiscretization end
+
+
+# abstract supertype of specific equations such as the compressible Euler equations
+abstract type AbstractEquations{NDIMS, NVARS} end
+
+
+# TODO: Taal decide, which abstract types shall be defined here?

src/auxiliary/auxiliary.jl

@@ -127,7 +127,7 @@ end
 
 Return the path to an example parameter file that can be used to quickly see Trixi in action.
 """
-default_example() = joinpath(examples_dir(), "2d", "parameters.toml") # TODO: Taal, use parameters.jl
+default_example() = joinpath(examples_dir(), "2d", "parameters.jl")
 
 
 """

src/callbacks/alive.jl

@@ -4,7 +4,7 @@
 
 Inexpensive callback showing that a simulation is still running by printing
 some information such as the current time to the screen every `alive_interval`
-time steps. If `analysis_interval ≂̸ 0`, sthe output is omitted every
+time steps. If `analysis_interval ≂̸ 0`, the output is omitted every
 `analysis_interval` time steps.
 """
 mutable struct AliveCallback

src/equations/equations.jl

@@ -1,7 +1,4 @@
 
-# Basic abstract types creating the hierarchy
-abstract type AbstractEquations{NDIMS, NVARS} end
-
 # Retrieve number of variables from equation type
 @inline nvariables(::Type{AbstractEquations{NDIMS, NVARS}}) where {NDIMS, NVARS} = NVARS
 

src/run.jl

@@ -1,4 +1,51 @@
-using LinearMaps: LinearMap
+
+# Apply the function `f` to `expr` and all sub-expressions recursively.
+walkexpr(f, expr::Expr) = f(Expr(expr.head, (walkexpr(f, arg) for arg in expr.args)...))
+walkexpr(f, x) = f(x)
+
+# Replace assignments to `key` in `expr` by `key = val` for all `(key,val)` in `kwargs`.
+function replace_assignments(expr; kwargs...)
+  walkexpr(expr) do x
+    if x isa Expr
+      for (key,val) in kwargs
+        if x.head === Symbol("=") && x.args[1] === Symbol(key)
+          x.args[2] = :( $val )
+          # dump(x)
+        end
+      end
+    end
+    return x
+  end
+end
+
+# Note: Wa can't call the method below `Trixi.include` since that is created automatically
+# inside `module Trixi` to `include` source files and evaluate them within the global scope
+# of `Trixi`. However, users will want to evaluate in the global scope of `Main` or something
+# similar to manage dependencies on their own.
+"""
+    trixi_include([mod::Module=Main,] parameters_file; kwargs...)
+
+`include` the file `parameters_file` and evaluate its content in the global scope of module `mod`.
+You can override specific assignments in `parameters_file` by supplying keyword arguments.
+It's basic purpose is to make it easier to modify some parameters while running Trixi from the
+REPL. Additionally, this is used in tests to reduce the computational burden for CI while still
+providing examples with sensible default values for users.
+
+# Examples
+
+```jldoctest
+julia> trixi_include(default_example(), tspan=(0.0, 0.1))
+...
+julia> sol.t[end]
+0.1
+```
+"""
+function trixi_include(mod::Module, parameters_file::AbstractString; kwargs...)
+  Base.include(ex -> replace_assignments(ex; kwargs...), mod, parameters_file)
+end
+
+trixi_include(parameters_file::AbstractString; kwargs...) = trixi_include(Main, parameters_file; kwargs...)
+
 
 
 """

src/semidiscretization.jl

@@ -1,7 +1,4 @@
 
-abstract type AbstractSemidiscretization end
-
-
 """
     ndofs(semi::AbstractSemidiscretization)
 

src/solvers/dg/dg.jl

@@ -321,7 +321,8 @@ struct IndicatorHennemannGassner{RealT<:Real, Variable, Cache} <: AbstractIndica
   cache::Cache
 end
 
-function IndicatorHennemannGassner(equations, basis;
+# this method is used when the indicator is constructed as for shock-capturing volume integrals
+function IndicatorHennemannGassner(equations::AbstractEquations, basis;
                                    alpha_max=0.5,
                                    alpha_min=0.001,
                                    alpha_smooth=true,
@@ -332,7 +333,8 @@ function IndicatorHennemannGassner(equations, basis;
     alpha_max, alpha_min, alpha_smooth, variable, cache)
 end
 
-function IndicatorHennemannGassner(semi;
+# this method is used when the indicator is constructed as for AMR
+function IndicatorHennemannGassner(semi::AbstractSemidiscretization;
                                    alpha_max=0.5,
                                    alpha_min=0.001,
                                    alpha_smooth=true,
@@ -384,12 +386,16 @@ struct IndicatorLöhner{RealT<:Real, Variable, Cache} <: AbstractIndicator
   cache::Cache
 end
 
-function IndicatorLöhner(basis, equations; f_wave=0.2, variable=first)
+# this method is used when the indicator is constructed as for shock-capturing volume integrals
+function IndicatorLöhner(equations::AbstractEquations, basis;
+                         f_wave=0.2, variable=first)
   cache = create_cache(IndicatorLöhner, equations, basis)
   IndicatorLöhner{typeof(f_wave), typeof(variable), typeof(cache)}(f_wave, variable, cache)
 end
 
-function IndicatorLöhner(semi; f_wave=0.2, variable=first)
+# this method is used when the indicator is constructed as for AMR
+function IndicatorLöhner(semi::AbstractSemidiscretization;
+                         f_wave=0.2, variable=first)
   cache = create_cache(IndicatorLöhner, semi)
   IndicatorLöhner{typeof(f_wave), typeof(variable), typeof(cache)}(f_wave, variable, cache)
 end
@@ -415,12 +421,16 @@ struct IndicatorMax{Variable, Cache<:NamedTuple} <: AbstractIndicator
   cache::Cache
 end
 
-function IndicatorMax(basis, equations::AbstractEquations; variable=first)
+# this method is used when the indicator is constructed as for shock-capturing volume integrals
+function IndicatorMax(equations::AbstractEquations, basis;
+                      variable=first)
   cache = create_cache(IndicatorMax, equations, basis)
   IndicatorMax{typeof(variable), typeof(cache)}(variable, cache)
 end
 
-function IndicatorMax(semi; variable=first)
+# this method is used when the indicator is constructed as for AMR
+function IndicatorMax(semi::AbstractSemidiscretization;
+                      variable=first)
   cache = create_cache(IndicatorMax, semi)
   return IndicatorMax{typeof(variable), typeof(cache)}(variable, cache)
 end

src/solvers/dg/dg_2d.jl

@@ -16,7 +16,6 @@ function create_cache(mesh::TreeMesh{2}, equations::AbstractEquations{2},
   interfaces = init_interfaces(leaf_cell_ids, mesh, elements,
                                RealT, nvariables(equations), polydeg(dg))
 
-  # TODO: Taal implement BCs
   boundaries, _ = init_boundaries(leaf_cell_ids, mesh, elements,
                                   RealT, nvariables(equations), polydeg(dg))
 
@@ -266,7 +265,8 @@ function calc_volume_integral!(du::AbstractArray{<:Any,4}, u,
   end
 end
 
-@inline function split_form_kernel!(du, u, nonconservative_terms::Val{false}, equations,
+@inline function split_form_kernel!(du::AbstractArray{<:Any,4}, u,
+                                    nonconservative_terms::Val{false}, equations,
                                     volume_flux, dg::DGSEM, cache,
                                     element, alpha=true)
   # true * [some floating point value] == [exactly the same floating point value]
@@ -309,7 +309,8 @@ end
   end
 end
 
-@inline function split_form_kernel!(du, u, nonconservative_terms::Val{true}, equations,
+@inline function split_form_kernel!(du::AbstractArray{<:Any,4}, u,
+                                    nonconservative_terms::Val{true}, equations,
                                     volume_flux, dg::DGSEM, cache,
                                     element, alpha=true)
   @unpack derivative_split_transpose = dg.basis
@@ -339,6 +340,7 @@ end
 end
 
 
+# TODO: Taal dimension agnostic
 function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, nonconservative_terms, equations,
                                volume_integral::VolumeIntegralShockCapturingHG,
                                dg::DGSEM, cache)

test/test_manual.jl

@@ -183,7 +183,7 @@ isdir(outdir) && rm(outdir, recursive=true)
   @testset "example parameters" begin
     @test basename(examples_dir()) == "examples"
     @test !isempty(get_examples())
-    @test endswith(default_example(), "parameters.toml")
+    @test endswith(default_example(), "parameters.jl")
   end
 
   @testset "DG L2 mortar container debug output" begin

test/test_trixi.jl

@@ -29,44 +29,39 @@ function test_trixi_run(parameters_file; l2=nothing, linf=nothing, atol=200*eps(
 end
 
 
+"""
+    test_trixi_include(parameters_file; l2=nothing, linf=nothing, atol=10*eps(), rtol=0.001, parameters...)
 
-walkexpr(f, expr::Expr) = f(Expr(expr.head, (walkexpr(f, arg) for arg in expr.args)...))
-walkexpr(f, x) = f(x)
-
-function mapexpr(expr; kwargs...)
-  walkexpr(expr) do x
-    if x isa Expr
-      for (key,val) in kwargs
-        if x.head === Symbol("=") && x.args[1] === Symbol(key)
-          x.args[2] = :( $val )
-          # dump(x)
-        end
-      end
-    end
-    return x
-  end
-end
-
+Test Trixi by calling `trixi_include(parameters_file; parameters...)`.
+By default, only the absence of error output is checked.
+If `l2` or `linf` are specified, in addition the resulting L2/Linf errors
+are compared approximately against these reference values, using `atol, rtol`
+as absolute/relative tolerance.
+"""
 function test_trixi_include(parameters_file; l2=nothing, linf=nothing,
                                              atol=200*eps(), rtol=0.001,
                                              kwargs...)
 
   println("#"^80)
   println(parameters_file)
 
-  @test_nowarn include(ex -> mapexpr(ex; kwargs...), parameters_file)
-  l2_measured, linf_measured = analysis_callback(sol)
+  # evaluate examples in the scope of the module they're called from
+  @test_nowarn trixi_include(@__MODULE__, parameters_file; kwargs...)
 
   # If present, compare L2 and Linf errors against reference values
-  if !isnothing(l2)
-    for (l2_expected, l2_actual) in zip(l2, l2_measured)
-      @test isapprox(l2_expected, l2_actual, atol=atol, rtol=rtol)
+  if !isnothing(l2) || !isnothing(linf)
+    l2_measured, linf_measured = analysis_callback(sol)
+
+    if !isnothing(l2)
+      for (l2_expected, l2_actual) in zip(l2, l2_measured)
+        @test isapprox(l2_expected, l2_actual, atol=atol, rtol=rtol)
+      end
     end
-  end
 
-  if !isnothing(linf)
-    for (linf_expected, linf_actual) in zip(linf, linf_measured)
-      @test isapprox(linf_expected, linf_actual, atol=atol, rtol=rtol)
+    if !isnothing(linf)
+      for (linf_expected, linf_actual) in zip(linf, linf_measured)
+        @test isapprox(linf_expected, linf_actual, atol=atol, rtol=rtol)
+      end
     end
   end