Closed ranocha closed 3 years ago
By the way,
julia> round(Int, 42 * 4.2, RoundUp)
177
Seems to fit to #42 :wink:
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...?
First of all: Great work! :+1: 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:
elements
live? I didn't see it up there. Or would this also end up in the cache
?examples/paper-self-gravitating-gas-dynamics/parameters_sedov_self_gravity.toml
?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?
- Where would
elements
live? I didn't see it up there. Or would this also end up in thecache
?
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 integratoru_tmp2
, u_tmp3
are caches of the time integrator and will be handled there internallyinverse_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.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.
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?
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)
Taal: Trixi as a library
This is a draft of a possible approach to restructure Trixi.jl to
Design principles
Some design principles I had in mind while writing this draft:
initial_conditions
from the basic DG stuff, the function computing the volume integral does not need to be recompiled when differentinitial_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.
[x] Move everything related to Lobatto-Legendre nodes to a separate struct, maybe
[x] We could create new structs for different volume integral types.
[x] A DG struct will wrap all of these components.
[x] Everything necessary for computing the rate of change of the conserved variables will be wrapped in a semidiscretization struct.
This
semidisc::Semidiscretization
is basically the ODE right-hand side. It will spit out anODEProblem
which can be passed tosolve
from the DiffEq suite. For example, we could have something like[x] We need to figure out a nice way to store caches, temporary information, intermediate results etc. There are different kinds of such buffers.
thread_cache
inDgXD
: Only needed to improve the performance, does never change, and will not be of interest for output etc.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.blending_factors
oramr_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 thesolver
etc. and asking them for stuff to put into the cache? Thiscache
could be aNamedTuple
, initialized ascache = NamedTuple()
and updated ascache = (cache..., new_part=new_part)
.[x] All analysis stuff will be moved to new structs and wrapped in callbacks, e.g. a PresetTimeCallback or a PeriodicCallback.
globals
,parameters
,main_timer
. Parameters etc. should be provided when constructing the specific objects and the timer should be created/passed (if/when desired).[x] We should create some convenience functions to construct everything. For example, I could imagine to end up with something along the lines of
from
examples/2d/parameters.toml
.examples/**.jl
file?compute_linear_structure
,compute_jacobian_dg
, fluctuation stuffconvtest
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.