Skip to content

Commit

Permalink
WIP: tutorial on semidiscretizations
Browse files Browse the repository at this point in the history
  • Loading branch information
ranocha committed Sep 12, 2023
1 parent f5f2f14 commit 168c786
Show file tree
Hide file tree
Showing 2 changed files with 116 additions and 2 deletions.
116 changes: 115 additions & 1 deletion docs/literate/src/files/custom_semidiscretization.jl
Original file line number Diff line number Diff line change
Expand Up @@ -24,11 +24,125 @@
# passed to [`semidiscretize`](@ref) as parameter.
# For a [`SemidiscretizationHyperbolic`](@ref), the `ODEProblem` wraps
# `Trixi.rhs!` as ODE RHS.
#
# For a [`SemidiscretizationHyperbolicParabolic`](@ref), Trixi.jl
# uses a `SplitODEProblem` combining `Trixi.rhs_parabolic!` for the
# (potentially) stiff part and `Trixi.rhs!` for the other part.


# ## Standard Trixi.jl setup

# In this tutorial, we will consider the linear advection equation
# with source term
# ```math
# \partial_t u(t,x) + \partial_x u(t,x) = -\exp(-t) sin\bigl(\pi (x - t) \bigr)
# ```
# with periodic boundary conditions in the domain `[-1, 1]` as a
# model problem.
# The initial condition is
# ```math
# u(0,x) = \sin(\pi x).
# ```
# The source term results in some damping and the analytical solution
# ```math
# u(t,x) = \exp(-t) \sin\bigl(\pi (x - t) \bigr).
# ```
# First, we discretize this equation using the standard functionality
# of Trixi.jl.

using Trixi, OrdinaryDiffEq, Plots

# The linear scalar advection equation is already implemented in
# Trixi.jl as [`LinearScalarAdvectionEquation1D`](@ref). We construct
# it with an advection velocity `1.0`.

equations = LinearScalarAdvectionEquation1D(1.0)

# Next, we use a standard [`DGSEM`](@ref) solver.

solver = DGSEM(polydeg = 3)

# We create a simple [`TreeMesh`](@ref) in 1D.

coordinates_min = (-1.0,)
coordinates_max = (+1.0,)
mesh = TreeMesh(coordinates_min, coordinates_max;
initial_refinement_level = 4,
n_cells_max = 10^4,
periodicity = true)

# We wrap everything in in a semidiscretization and pass the source
# terms as a standard Julia function. Please note that Trixi.jl uses
# `SVector`s from
# [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl)
# to store the conserved variables `u`. Thus, the return value of the
# source terms must be wrapped in an `SVector` - even if we consider
# just a scalar problem.

function initial_condition(x, t, equations)
return SVector(exp(-t) * sinpi(x[1] - t))
end

function source_terms_standard(u, x, t, equations)
return -initial_condition(x, t, equations)
end

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition,
solver;
source_terms = source_terms_standard)

# Now, we can create the `ODEProblem`, solve the resulting ODE
# using a time integration method from
# [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl),
# and visualize the numerical solution at the final time using
# [Plots.jl](https://github.com/JuliaPlots/Plots.jl).

tspan = (0.0, 3.0)
ode = semidiscretize(semi, tspan)

sol = solve(ode, RDPK3SpFSAL49(); ode_default_options()...)

plot(sol; label = "numerical sol.")

# We can also plot the analytical solution for comparison.
# Since Trixi.jl uses `SVector`s for the variables, we take their `first`
# (and only) component to get the scalar value for manual plotting.

let
x = range(-1.0, 1.0; length = 200)
plot!(x, first.(initial_condition.(x, sol.t[end], equations)),
label = "analytical sol.", linestyle = :dash, legend = :topleft)
end

# We can also add the initial condition to the plot.

plot!(sol.u[1], semi, label = "u0", linestyle = :dot, legend = :topleft)

# You can of course also use some
# [callbacks](https://trixi-framework.github.io/Trixi.jl/stable/callbacks/)
# provided by Trixi.jl as usual.

summary_callback = SummaryCallback()
analysis_interval = 100
analysis_callback = AnalysisCallback(semi; interval = analysis_interval)
alive_callback = AliveCallback(; analysis_interval)
callbacks = CallbackSet(summary_callback,
analysis_callback,
alive_callback)

sol = solve(ode, RDPK3SpFSAL49();
ode_default_options()..., callback = callbacks)
summary_callback()


# ## Using a custom ODE right-hand side function

# TODO


# ## Setting up a custom semidiscretization

# TODO

# Required
# - `Trixi.rhs!(du_ode, u_ode, semi::SemidiscretizationEulerGravity, t)`
# - `Trixi.mesh_equations_solver_cache(semi::SemidiscretizationEulerGravity)`
Expand Down
2 changes: 1 addition & 1 deletion docs/src/callbacks.md
Original file line number Diff line number Diff line change
Expand Up @@ -106,7 +106,7 @@ will yield the following plot:
the automated performance measurements, including an output of the recorded timers after a simulation.
* The [`VisualizationCallback`](@ref) can be used for in-situ visualization. See
[Visualizing results during a simulation](@ref).
* The [`TrivialCallback`](@ref) does nothing and can be used to to easily disable some callbacks
* The [`TrivialCallback`](@ref) does nothing and can be used to easily disable some callbacks
via [`trixi_include`](@ref).

### Equation-specific callbacks
Expand Down

0 comments on commit 168c786

Please sign in to comment.