Taal aka Trixi as a library

[ranocha]

Taal: Trixi as a library

This is a draft of a possible approach to restructure Trixi.jl to

• make it more Julian in general (General plan and layout #42)
• make it easier to use without having to modify any internals
• make it more flexible and composable with other Julia packages (Separate implementation of time integration schemes #83)
• finally get rid of TOML parameter files (TOML? Why? #150)

Design principles

Some design principles I had in mind while writing this draft:

• Verbs aka functions are the most important design decisions. Their meaning should be clear
  and they should be re-usable for different types via multiple dispatch.
• Everything evolves around some basic structs which describe the discretization and related functions.
• Structs shouldn't be jack of all trades devices suitable for every possible purpose.
  Instead, they should serve a well-defined purpose wherever orthogonal design decisions are possible.
  Structs can be composed to combine them.
• Only the parts that are needed for some computations are passed to the
  associated functions. Then, it's comparably easy to use multiple dispatch
  and methods do not need to be recompiled that often because something changed.
  For example, by separating the initial_conditions from the basic DG stuff,
  the function computing the volume integral does not need to be recompiled
  when different initial_conditions are chosen. This could help reducing the
  latency and time-to-first-plot.

Draft of the new layout of Trixi

Here is a list of suggested changes.

• Move everything related to Lobatto-Legendre nodes to a separate struct, maybe

  struct LobattoLegendreBasis{RealT<:Real, POLYDEG, DerivativeMatrix, MortarMatrix} <: AbstractSBPOperator{RealT}
    nodes          ::SVector{NNODES, RealT}
    weights        ::SVector{NNODES, RealT}
    inverse_weights::SVector{NNODES, RealT}
  
    inverse_vandermonde_legendre::Matrix{RealT}
    boundary_interpolation      ::Matrix{RealT} # lhat
  
    derivative_matrix   ::DerivativeMatrix # dsplit
    derivative_adjoint  ::DerivativeMatrix # dhat, adjoint wrt the SBP dot product
    derivative_transpose::DerivativeMatrix # dsplit_transposed
  
    mortar_forward_upper::MortarMatrix
    mortar_forward_lower::MortarMatrix
    mortar_reverse_upper::MortarMatrix
    mortar_reverse_lower::MortarMatrix
  end

• We could create new structs for different volume integral types.

  struct VolumeIntegralWeakForm <: AbstractVolumeIntegral end
  
  struct VolumeIntegralSplitForm{VolumeFlux} <: AbstractVolumeIntegral end
    volume_flux_function::VolumeFlux
  end
  
  struct VolumeIntegralShockCapturing{VolumeFluxDG, VolumeFluxFV, ShockIndicatorVariable, RealT<:Real} <: AbstractVolumeIntegral end
    volume_flux_function_dg::VolumeFluxDG # symmetric, e.g. split-form or entropy-conservative
    volume_flux_function_fv::VolumeFluxFV # non-symmetric in general, e.g.entropy-dissipative
    shock_indicator_variable::ShockIndicatorVariable
    shock_alpha_max::RealT
    shock_alpha_min::RealT
    shock_alpha_smooth::Bool
  end

• A DG struct will wrap all of these components.

  struct Dg{RealT, Operator<:AbstractSBPOperator{RealT}, SurfaceFlux, VolumeIntegral} <: AbstractSolver
    operator::Operator
    surface_flux_function::SurfaceFlux
    volume_integral_operator::VolumeIntegral
  
    # TODO: caches...???
  end

• Everything necessary for computing the rate of change of the conserved variables will be wrapped in a semidiscretization struct.

  struct Semidiscretization{Mesh, Equations, InitialConditions, SourceTerms, Solver}
    mesh::Mesh
  
    equations::Equations
    # this guy is a bit messy since we abuse it as some kind of "exact solution"
    # although this doesn't really exist...
    initial_conditions::InitialConditions
    source_terms::SourceTerms
  
    solver::Solver
  
    # TODO: caches...???
  end

  This semidisc::Semidiscretization is basically the ODE right-hand side. It will spit out an ODEProblem which can be passed to solve from the DiffEq suite.
  For example, we could have something like

  const SemidiscretizationDG = Semidiscretization{Mesh, Equations, InitialConditions, SourceTerms, Solver} where {Mesh, Equations, InitialConditions, SourceTerms, Solver<:AbstractDG}
  
  # functor style
  function (semidisc::SemidiscretizationDG)(du, u, params, t)
    # Reset u_t
    # Calculate volume integral
    # Prolong solution to interfaces
    # Calculate interface fluxes
    # Prolong solution to boundaries
    # Calculate boundary fluxes
    # Prolong solution to mortars
    # Calculate mortar fluxes
    # Calculate surface integrals
    # Apply Jacobian from mapping to reference element
    # Calculate source terms
  end
  
  # function with parameters style
  function rhs!(du, u, semidisc::SemidiscretizationDG, t) # do we like rhs! ?
    # Reset u_t
    # Calculate volume integral
    # Prolong solution to interfaces
    # Calculate interface fluxes
    # Prolong solution to boundaries
    # Calculate boundary fluxes
    # Prolong solution to mortars
    # Calculate mortar fluxes
    # Calculate surface integrals
    # Apply Jacobian from mapping to reference element
    # Calculate source terms
  end

• We need to figure out a nice way to store caches, temporary information, intermediate results etc. There are different kinds of such buffers.

  • Something like the thread_cache in DgXD: Only needed to improve the performance, does never change, and will not be of interest for output etc.
  • Something like the interfaces, boundaries, mortars: Really needed to improve the performance, contain some important logic and connectivity information, need to be modified, e.g. by AMR. Not needed for visualization.
  • Something like the blending_factors or amr_indicators: Needed to improve the performance, need to be modified by AMR, interesting for visualization.

  Maybe we can have something like cache = create_cache(eqs, mesh, solver), which will return an appropriate cache by iterating over all components of the solver etc. and asking them for stuff to put into the cache? This cache could be a NamedTuple, initialized as cache = NamedTuple() and updated as cache = (cache..., new_part=new_part).

• All analysis stuff will be moved to new structs and wrapped in callbacks, e.g. a PresetTimeCallback or a PeriodicCallback.

  struct LobattoLegendreAnalysis{RealT<:Real, POLYDEG}
    analysis_nodes  ::SVector{NNODES, RealT}
    analysis_weights::SVector{NNODES, RealT}
    analysis_vandermonde::Matrix{RealT}
  end
  
  struct AnalysisCallback{RealT, AnalysisHelper}
    analysis_total_volume::RealT
    save_analysis::Bool
    analysis_filename::String
    analysis_quantities::Vector{Symbol} # something else?
    analysis_helper::AnalysisHelper
  end
  
  function AnalysisCallback(basis::LobattoLegendreBasis{RealT}; analysis_polydeg=2*polydeg(basis)) where {RealT}
    analysis_helper = LobattoLegendreAnalysis{RealT, analysis_polydeg}()
    # etc.
  end

• All AMR stuff will be moved to new structs wrapped in callbacks, e.g. a PresetTimeCallback or a PeriodicCallback. Resizing the caches in DiffEq can be done similar to https://diffeq.sciml.ai/latest/features/callback_functions/#Example-3:-Growing-Cell-Population.

• CFL-based time stepping is obtained via callbacks. But we should improve the interface (Compute CFL factors a la Fluxo #139).

• Steady state stuff: https://diffeq.sciml.ai/latest/features/callback_library/#TerminateSteadyState

• Get rid of globals like globals, parameters, main_timer.
  Parameters etc. should be provided when constructing the specific objects
  and the timer should be created/passed (if/when desired).

• We should create some convenience functions to construct everything.
  For example, I could imagine to end up with something along the lines of

  using OrdinaryDiffEq
  using Trixi
  using StaticArrays
  
  advection_velocity = SVector(1.0,1 1.0)
  eqs = LinearScalarAdvectionEquation2D(advection_velocity)
  initial_conditions = initial_conditions_convergence_test # from Trixi
  
  coordinates_min = (-1.0, -1.0)
  coordinates_min = ( 1.0,  1.0)
  initial_refinement_level = 4
  mesh = TreeMesh2D(coordinates_min=coordinates_min,
                    coordinates_max=coordinates_max,
                    initial_refinement_level=initial_refinement_level)
  
  solver = DG2D(polydeg=3, surface_flux=flux_lax_friedrichs)
  
  semidisc = Semidiscretization(mesh, eqs, initial_conditions, solver) # source_terms default to nothing
  
  cfl_callback = StepsizeLimiter(calc_max_dt(eqs, mesh, solver, cfl=0.8))
  analysis_callback = AnalysisCallback(solver, extra_analysis_quantities=["entropy", "energy_total"])
  # primitive_variables is something like Trixi.cons2prim
  output_callback = OutputCallback(solution_interval=100, solution_variables=primitive_variables, restart_interval=10)
  callback = CallbackSet(cfl_callback, analysis_callback, output_callback)
  
  ode = init!(semidisc) # some better name...
  
  tspan = (0.0, 1.0)
  sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
              dt=calc_max_dt(ode), save_everystep=false,
              callback=callback)

  from examples/2d/parameters.toml.

• Figure out a good way to test everything. Modifying some parameter settings which are otherwise taken from a file isn't possible anymore. Shall we just write the tests inside every examples/**.jl file?

• Port previous elixirs: compute_linear_structure, compute_jacobian_dg, fluctuation stuff

• Port previous syntactic sugar: convtest

Disclaimer

This is just a rough first draft and needs to be refined. We definitely have
to figure out a lot of details and invest a lot of time to make this work.

Thoughts, ideas, suggestions? I'm open to discussions.

[ranocha]

By the way,

julia> round(Int, 42 * 4.2, RoundUp)
177

Seems to fit to #42 😉

[gregorgassner]

I am curious about the very last file: where does a user get the interface information, when looking at this linear scalar advection problem.

Say, I want to do a new initial condition, but don't mess around with the internals. In principles, I could implement it directly in the elixir...however, I have to know how the interface of initial condition looks like, right? Hence, I have to dig through the Trixi internals anyways...?

[sloede]

First of all: Great work! 👍 I think this is a very good start for getting the ball rolling on this most complex undertaking we have tried so far with Trixi...

Two initial thoughts/questions:

1. Where would elements live? I didn't see it up there. Or would this also end up in the cache?
2. How does multi-physics work? For instance, how would your example file look for, e.g., examples/paper-self-gravitating-gas-dynamics/parameters_sedov_self_gravity.toml?

[ranocha]

Say, I want to do a new initial condition, but don't mess around with the internals. In principles, I could implement it directly in the elixir...however, I have to know how the interface of initial condition looks like, right? Hence, I have to dig through the Trixi internals anyways...?

The interface will be documented. Hence, people can either browse the source code of Trixi or look at the docs - maybe some of the tutorials?

[ranocha]

1. Where would elements live? I didn't see it up there. Or would this also end up in the cache?

Right now, elements is an ElementContainer, which looks like this.

struct ElementContainer2D{NVARS, POLYDEG} <: AbstractContainer
  u::Array{Float64, 4}                   # [variables, i, j, elements]
  u_t::Array{Float64, 4}                 # [variables, i, j, elements]
  u_tmp2::Array{Float64, 4}              # [variables, i, j, elements]
  u_tmp3::Array{Float64, 4}              # [variables, i, j, elements]
  inverse_jacobian::Vector{Float64}      # [elements]
  node_coordinates::Array{Float64, 4}    # [orientation, i, j, elements]
  surface_ids::Matrix{Int}               # [direction, elements]
  surface_flux_values::Array{Float64, 4} # [variables, i, direction, elements]
  cell_ids::Vector{Int}                  # [elements]
end

• u is the current state, u_t it's time derivative. Both will be handled by the time integrator
• u_tmp2, u_tmp3 are caches of the time integrator and will be handled there internally
• inverse_jacobian, node_coordinates, surface_ids, and cell_ids contain structural information. They belong to the second category of cache stuff. I would probably store them in the Semidiscretization.
• surface_flux_values is a temporary buffer that definitely belongs to the cache.

[ranocha]

How does multi-physics work? For instance, how would your example file look for, e.g., examples/paper-self-gravitating-gas-dynamics/parameters_sedov_self_gravity.toml?

That depends on how much we want to hide from the user. For example, we could create a type EulerGravitySemidiscretization, that will contain a parameterized ODE problem for the hyperbolic diffusion part. This ODE problem is re-initialized with the given data and solved to steady state for every RHS evaluation.

[sloede]

Could you maybe give an example of how the Taal script for examples/paper-self-gravitating-gas-dynamics/parameters_sedov_self_gravity.toml could look like, just to get an idea?

[ranocha]

Maybe something like

using OrdinaryDiffEq
using Trixi
using StaticArrays

eqs_euler = CompressibleEulerEquations2D(gamma=1.4)
eqs_gravity = HyperbolicDiffusionEquations2D()
initial_conditions = initial_conditions_sedov_self_gravity # from Trixi

mesh = TreeMesh2D(coordinates_min=(-4.0, -4.0),
                  coordinates_max=( 4.0,  4.0),
                  initial_refinement_level=2,
                  periodicity=false)

volume_integral = VolumeIntegralShockCapturing(volume_flux_function_dg=flux_chandrashekar,
                                               volume_flux_function_fv=flux_hll,
                                               shock_indicator_variable=density_pressure)
solver = DG2D(polydeg=3, surface_flux=flux_hll,
              volume_integral=volume_integral)

# maybe we can also pass a time integration method as parameter, depending on how we would like
# to structure this part
semidisc = EulerGravitySemidiscretization(mesh, eqs_euler, eqs_gravity, initial_conditions, solver,
                                          source_terms=source_terms_harmonic,
                                          resid_tol=1.0e-4, cfl_gravity=1.2,
                                          G=6.674e-8, rho0=0.0)

cfl_callback = StepsizeLimiter(calc_max_dt(eqs, mesh, solver, cfl=0.5))
analysis_callback = AnalysisCallback(solver, extra_analysis_quantities=["entropy", "energy_total","energy_internal"])
amr_callback = AdaptiveMeshRefinementCallback(indicator=amr_indicator_sedov_self_gravity,
                                              interval=1, alpha_max=1.0, alpha_min=0.0)
# primitive_variables is something like Trixi.cons2prim
output_callback = OutputCallback(solution_interval=100, solution_variables=primitive_variables, restart_interval=10)
callback = CallbackSet(cfl_callback, analysis_callback, output_callback)

ode = init!(semidisc) # some better name...

tspan = (0.0, 1.0)
sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=calc_max_dt(ode), save_everystep=false,
            callback=callback)