Support multiple dimensions

[ranocha]

We can track ideas/suggestions to extend Trixi to multiple dimensions (1D, 2D, and 3D) here. Here is a short summary of our meeting on 2020-08-21.

• [x] Use multiple dispatch and type parameters (equations, mesh, solver) This could be either along the lines of

  struct LinearScalarAdvectionEquation{DIMS} <: LinearScalarAdvectionEquation{DIMS, 2}
  sources::String
  advectionvelocity::SVector{DIMS, Float64}
  end
  const LinearScalarAdvectionEquation2D = LinearScalarAdvectionEquation{2}

  or

  abstract type LinearScalarAdvectionEquation{DIMS} <: AbstractEquation{DIMS,1} end
  struct LinearScalarAdvectionEquation2D <: LinearScalarAdvectionEquation{2}
  sources::String
  advectionvelocity::SVector{2, Float64}
  end

  Multiple dispatch will work accordingly without changing the function names, e.g.

  function flux_lax_friedrichs(u_ll, u_rr, orientation, equation::CompressibleEulerEquations2D)

• [x] Introduce Base.ndims(::AbstractEquation{DIMS}) where DIMS = DIMS and replace const ndim

• [x] Introduce a folder structure like

  # solvers
  # | dg # general stuff: struct Dg, rhs!, interpolation.jl etc.
  # | | dg1D # specific stuff for 1D
  # | | dg2D # specific stuff for 2D
  
  # equations
  # | linear_scalar_advection.jl # general stuff, includes dimension-specific stuff
  # | compressible_euler.jl
  #
  # | 2d
  # | | compressible_euler.jl
  # | | linear_scalar_advection.jl
  # | 3d
  # | | compressible_euler.jl
  #
  # | compressible_euler_2d.jl
  # | linear_scalar_advection_2d.jl
  # | compressible_euler_3d.jl

• [x] Introduce parameter dims in all parameter files along the lines of

  # Simulation
  dims = 2
  equations = "LinearScalarAdvectionEquation"

  Query this parameter when creating the mesh and equations - the solver will get to know it this way

• [x] linear advection, uniform Cartesian mesh
• [x] compressible Euler, uniform Cartesian mesh
• [x] non-conforming interfaces, L^2 mortars
• [x] AMR in 3D
• [x] visualization
• [x] split forms
• [x] TGV
• [x] shock capturing (with Sedov blast setup like for self-gravitating Sedov blast)
• [x] hyperbolic diffusion
• [x] Euler + gravity (Sedov with self-gravity)
• [x] everything with # FIXME ndims
• [x] MHD
• [x] various numerical fluxes
• [x] clean up...
• [ ] update docs and README.md
• Morning and afternoon sessions every day next week
  • 09:00 - 13:00
  • 14:00 - 17:00

[ranocha]

Some inspiration for the first task:

julia> abstract type AbstractEquation{NDIMS, V} end

julia> struct AbstractCompressibleEuler{NDIMS, V, T} <: AbstractEquation{NDIMS, V}
           gamma::T
       end

julia> function AbstractCompressibleEuler{NDIMS}(gamma::T) where {NDIMS, T}
           AbstractCompressibleEuler{NDIMS, NDIMS+2, T}(gamma)
       end

julia> AbstractCompressibleEuler{2}(1.4)
AbstractCompressibleEuler{2,4,Float64}(1.4)

julia> @code_warntype AbstractCompressibleEuler{2}(1.4)
Variables
  #self#::Type{AbstractCompressibleEuler{2,V,T} where T where V}
  gamma::Float64

Body::AbstractCompressibleEuler{2,4,Float64}
1 ─ %1 = $(Expr(:static_parameter, 1))::Core.Compiler.Const(2, false)
│   %2 = ($(Expr(:static_parameter, 1)) + 2)::Core.Compiler.Const(4, false)
│   %3 = Core.apply_type(Main.AbstractCompressibleEuler, %1, %2, $(Expr(:static_parameter, 2)))::Core.Compiler.Const(AbstractCompressibleEuler{2,4,Float64}, false)
│   %4 = (%3)(gamma)::AbstractCompressibleEuler{2,4,Float64}
└──      return %4

julia> CompressibleEuler1D = AbstractCompressibleEuler{1}
AbstractCompressibleEuler{1,V,T} where T where V

julia> eq = CompressibleEuler1D(1.4)
AbstractCompressibleEuler{1,3,Float64}(1.4)

julia> @code_warntype CompressibleEuler1D(1.4)
Variables
  #self#::Type{AbstractCompressibleEuler{1,V,T} where T where V}
  gamma::Float64

Body::AbstractCompressibleEuler{1,3,Float64}
1 ─ %1 = $(Expr(:static_parameter, 1))::Core.Compiler.Const(1, false)
│   %2 = ($(Expr(:static_parameter, 1)) + 2)::Core.Compiler.Const(3, false)
│   %3 = Core.apply_type(Main.AbstractCompressibleEuler, %1, %2, $(Expr(:static_parameter, 2)))::Core.Compiler.Const(AbstractCompressibleEuler{1,3,Float64}, false)
│   %4 = (%3)(gamma)::AbstractCompressibleEuler{1,3,Float64}
└──      return %4

julia> function foo(eq::CompressibleEuler1D)
           eq.gamma
       end
foo (generic function with 1 method)

julia> foo(eq)
1.4

julia> @code_warntype foo(eq)
Variables
  #self#::Core.Compiler.Const(foo, false)
  eq::AbstractCompressibleEuler{1,3,Float64}

Body::Float64
1 ─ %1 = Base.getproperty(eq, :gamma)::Float64
└──      return %1

[sloede]

I think this looks very good! However, in the spirit of one framework, multiple solvers, I suggest to make AbstractCompressibleEulerEquations really abstract and instead duplicate the gamma part for each dimension. This way, everything for one dimension is really inside a single folder, including the struct.

[ranocha]

Okay, that's also fine with me. Then we would have something along the lines of

julia> abstract type AbstractEquation{NDIMS, V, T<:Real} end

julia> abstract type CompressibleEulerEquations{NDIMS, V, T} end

julia> struct CompressibleEulerEquations1D{T} <: CompressibleEulerEquations{1, 3, T}
           gamma::T
       end

julia> function CompressibleEulerEquations{NDIMS}(gamma) where {NDIMS}
           # this is something that is needed in the current framework if we don't specify the
           # equations and mesh completely but have to read an additional parameter ndims
           if NDIMS == 1
               CompressibleEulerEquations1D(gamma)
           elseif NDIMS == 2
               CompressibleEulerEquations2D(gamma)
           elseif NDIMS == 3
               CompressibleEulerEquations1D(gamma)
           else
               throw(ArgumentError("$NDIMS dimensions not supported."))
           end
       end

julia> @code_warntype CompressibleEulerEquations{1}(1.4)
Variables
  #self#::Type{CompressibleEulerEquations{1,V,T} where T where V}
  gamma::Float64

Body::CompressibleEulerEquations1D{Float64}
1 ─ %1 = ($(Expr(:static_parameter, 1)) == 1)::Core.Compiler.Const(true, false)
│        %1
│   %3 = Main.CompressibleEulerEquations1D(gamma)::CompressibleEulerEquations1D{Float64}
└──      return %3
2 ─      Core.Compiler.Const(:($(Expr(:static_parameter, 1)) == 2), false)
│        Core.Compiler.Const(:(%5), false)
│        Core.Compiler.Const(:(Main.CompressibleEulerEquations2D(gamma)), false)
│        Core.Compiler.Const(:(return %7), false)
│        Core.Compiler.Const(:($(Expr(:static_parameter, 1)) == 3), false)
│        Core.Compiler.Const(:(%9), false)
│        Core.Compiler.Const(:(Main.CompressibleEulerEquations1D(gamma)), false)
│        Core.Compiler.Const(:(return %11), false)
│        Core.Compiler.Const(:(Base.string($(Expr(:static_parameter, 1)), " dimensions not supported.")), false)
│        Core.Compiler.Const(:(Main.ArgumentError(%13)), false)
│        Core.Compiler.Const(:(Main.throw(%14)), false)
└──      Core.Compiler.Const(:(return %15), false)

julia> Base.ndims(::AbstractEquation{NDIMS, V, T}) where {NDIMS, V, T} = NDIMS

julia> nvariables(::AbstractEquation{NDIMS, V, T}) where {NDIMS, V, T} = V # Shall we rename V to NVARS?
nvariables (generic function with 1 method)

julia> Base.eltype(::AbstractEquation{NDIMS, V, T}) where {NDIMS, V, T} = T

[sloede]

This looks good to me! Is it really required to have

julia> function CompressibleEulerEquations{NDIMS}(gamma) where {NDIMS}
...

three times or is this just an artifact from copy-paste?

[ranocha]

It's indeed an artifact - I've fixed that, thanks!

[sloede]

julia> nvariables(::AbstractEquation{NDIMS, V, T}) where {NDIMS, V, T} = V # Shall we rename V to NVARS?

Yes!