1.2: First steps in Trixi.jl: Create first setup

In this part of the introductory guide, we will create a first Trixi.jl setup as an extension of elixir_advection_basic.jl. Since Trixi.jl has a common basic structure for the setups, you can create your own by extending and modifying the following example.

Let's consider the linear advection equation for a state $u = u(x, y, t)$ on the two-dimensional spatial domain $[-1, 1] \times [-1, 1]$ with a source term

\[\frac{\partial}{\partial t}u + \frac{\partial}{\partial x} (0.2 u) - \frac{\partial}{\partial y} (0.7 u) = - 2 e^{-t} \sin\bigl(2 \pi (x - t) \bigr) \sin\bigl(2 \pi (y - t) \bigr),\]

with the initial condition

\[u(x, y, 0) = \sin\bigl(\pi x \bigr) \sin\bigl(\pi y \bigr),\]

and periodic boundary conditions.

The first step is to create and open a file with the .jl extension. You can do this with your favorite text editor (if you do not have one, we recommend VS Code). In this file you will create your setup.

To be able to use functionalities of Trixi.jl, you always need to load Trixi.jl itself and the OrdinaryDiffEq.jl package.

using Trixi
using OrdinaryDiffEq

The next thing to do is to choose an equation that is suitable for your problem. To see all the currently implemented equations, take a look at src/equations. If you are interested in adding a new physics model that has not yet been implemented in Trixi.jl, take a look at the tutorials Adding a new scalar conservation law or Adding a non-conservative equation.

The linear scalar advection equation in two spatial dimensions

\[\frac{\partial}{\partial t}u + \frac{\partial}{\partial x} (a_1 u) + \frac{\partial}{\partial y} (a_2 u) = 0\]

is already implemented in Trixi.jl as LinearScalarAdvectionEquation2D, for which we need to define a two-dimensional parameter advection_velocity describing the parameters $a_1$ and $a_2$. Appropriate for our problem is (0.2, -0.7).

advection_velocity = (0.2, -0.7)
equations = LinearScalarAdvectionEquation2D(advection_velocity)
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ LinearScalarAdvectionEquation2D                                                                  │
│ ═══════════════════════════════                                                                  │
│ #variables: ………………………………………………… 1                                                                │
│ │ variable 1: …………………………………………… scalar                                                           │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

To solve our problem numerically using Trixi.jl, we have to discretize the spatial domain, for which we set up a mesh. One of the most used meshes in Trixi.jl is the TreeMesh. The spatial domain used is $[-1, 1] \times [-1, 1]$. We set an initial number of elements in the mesh using initial_refinement_level, which describes the initial number of hierarchical refinements. In this simple case, the total number of elements is 2^initial_refinement_level throughout the simulation. The variable n_cells_max is used to limit the number of elements in the mesh, which cannot be exceeded when using adaptive mesh refinement.

All minimum and all maximum coordinates must be combined into Tuples.

coordinates_min = (-1.0, -1.0)
coordinates_max = ( 1.0,  1.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level = 4,
                n_cells_max = 30_000)
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ TreeMesh{2, Trixi.SerialTree{2}}                                                                 │
│ ════════════════════════════════                                                                 │
│ center: …………………………………………………………… [0.0, 0.0]                                                       │
│ length: …………………………………………………………… 2.0                                                              │
│ periodicity: ……………………………………………… (true, true)                                                     │
│ current #cells: ……………………………………… 341                                                              │
│ #leaf-cells: ……………………………………………… 256                                                              │
│ maximum #cells: ……………………………………… 30000                                                            │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

To approximate the solution of the defined model, we create a DGSEM solver. The solution in each of the recently defined mesh elements will be approximated by a polynomial of degree polydeg. For more information about discontinuous Galerkin methods, check out the Introduction to DG methods tutorial.

solver = DGSEM(polydeg=3)
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ DG{Float64}                                                                                      │
│ ═══════════                                                                                      │
│ basis: ……………………………………………………………… LobattoLegendreBasis{Float64}(polydeg=3)                         │
│ mortar: …………………………………………………………… LobattoLegendreMortarL2{Float64}(polydeg=3)                      │
│ surface integral: ………………………………… SurfaceIntegralWeakForm                                          │
│ │ surface flux: ……………………………………… flux_central                                                     │
│ volume integral: …………………………………… VolumeIntegralWeakForm                                           │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

Now we need to define an initial condition for our problem. All the already implemented initial conditions for LinearScalarAdvectionEquation2D can be found in src/equations/linear_scalar_advection_2d.jl. If you want to use, for example, a Gaussian pulse, it can be used as follows:

initial_conditions = initial_condition_gauss

But to show you how an arbitrary initial condition can be implemented in a way suitable for Trixi.jl, we define our own initial conditions.

\[u(x, y, 0) = \sin\bigl(\pi x \bigr) \sin\bigl(\pi y \bigr).\]

The initial conditions function must take spatial coordinates, time and equation as arguments and returns an initial condition as a statically sized vector SVector. Following the same structure, you can define your own initial conditions. The time variable t can be unused in the initial condition, but might also be used to describe an analytical solution if known. If you use the initial condition as analytical solution, you can analyze your numerical solution by computing the error, see also the section about analyzing the solution.

function initial_condition_sinpi(x, t, equations::LinearScalarAdvectionEquation2D)
    scalar = sinpi(x[1]) * sinpi(x[2])
    return SVector(scalar)
end
initial_condition = initial_condition_sinpi
initial_condition_sinpi (generic function with 1 method)

The next step is to define a function of the source term corresponding to our problem.

\[f(u, x, y, t) = - 2 e^{-t} \sin\bigl(2 \pi (x - t) \bigr) \sin\bigl(2 \pi (y - t) \bigr)\]

This function must take the state variable, the spatial coordinates, the time and the equation itself as arguments and returns the source term as a static vector SVector.

function source_term_exp_sinpi(u, x, t, equations::LinearScalarAdvectionEquation2D)
    scalar = - 2 * exp(-t) * sinpi(2*(x[1] - t)) * sinpi(2*(x[2] - t))
    return SVector(scalar)
end
source_term_exp_sinpi (generic function with 1 method)

Now we collect all the information that is necessary to define a spatial discretization, which leaves us with an ODE problem in time with a span from 0.0 to 1.0. This approach is commonly referred to as the method of lines.

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
                                    source_terms = source_term_exp_sinpi)
tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan);

At this point, our problem is defined. We will use the solve function defined in OrdinaryDiffEq.jl to get the solution. OrdinaryDiffEq.jl gives us the ability to customize the solver using callbacks without actually modifying it. Trixi.jl already has some implemented Callbacks. The most widely used callbacks in Trixi.jl are step control callbacks that are activated at the end of each time step to perform some actions, e.g. to print statistics. We will show you how to use some of the common callbacks.

To print a summary of the simulation setup at the beginning and to reset timers we use the SummaryCallback. When the returned callback is executed directly, the current timer values are shown.

summary_callback = SummaryCallback()
SummaryCallback

We also want to analyze the current state of the solution in regular intervals. The AnalysisCallback outputs some useful statistical information during the solving process every interval time steps.

analysis_callback = AnalysisCallback(semi, interval = 5)
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ AnalysisCallback                                                                                 │
│ ════════════════                                                                                 │
│ interval: ……………………………………………………… 5                                                                │
│ analyzer: ……………………………………………………… LobattoLegendreAnalyzer{Float64}(polydeg=6)                      │
│ │ error 1: …………………………………………………… l2_error                                                         │
│ │ error 2: …………………………………………………… linf_error                                                       │
│ │ integral 1: …………………………………………… entropy_timederivative                                           │
│ save analysis to file: …………………… no                                                               │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

It is also possible to control the time step size using the StepsizeCallback if the time integration method isn't adaptive itself. To get more details, look at CFL based step size control.

stepsize_callback = StepsizeCallback(cfl = 1.6)
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ StepsizeCallback                                                                                 │
│ ════════════════                                                                                 │
│ CFL number: ………………………………………………… 1.6                                                              │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

To save the current solution in regular intervals we use the SaveSolutionCallback. We would like to save the initial and final solutions as well. The data will be saved as HDF5 files located in the out folder. Afterwards it is possible to visualize a solution from saved files using Trixi2Vtk.jl and ParaView, which is described below in the section Visualize the solution.

save_solution = SaveSolutionCallback(interval = 5,
                                     save_initial_solution = true,
                                     save_final_solution = true)
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SaveSolutionCallback                                                                             │
│ ════════════════════                                                                             │
│ interval: ……………………………………………………… 5                                                                │
│ solution variables: …………………………… cons2prim                                                        │
│ save initial solution: …………………… yes                                                              │
│ save final solution: ………………………… yes                                                              │
│ output directory: ………………………………… /home/runner/work/Trixi.jl/Trixi.jl/docs/out                     │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

Alternatively, we have the option to print solution files at fixed time intervals.

save_solution = SaveSolutionCallback(dt = 0.1,
                                     save_initial_solution = true,
                                     save_final_solution = true)

Another useful callback is the SaveRestartCallback. It saves information for restarting in regular intervals. We are interested in saving a restart file for the final solution as well. To perform a restart, you need to configure the restart setup in a special way, which is described in the section Restart simulation.

save_restart = SaveRestartCallback(interval = 100, save_final_restart = true)
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SaveRestartCallback                                                                              │
│ ═══════════════════                                                                              │
│ interval: ……………………………………………………… 100                                                              │
│ save final solution: ………………………… yes                                                              │
│ output directory: ………………………………… /home/runner/work/Trixi.jl/Trixi.jl/docs/out                     │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

Create a CallbackSet to collect all callbacks so that they can be passed to the solve function.

callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback, save_solution,
                        save_restart)
CallbackSet{Tuple{}, Tuple{DiscreteCallback{typeof(Trixi.summary_callback), typeof(Trixi.summary_callback), Trixi.var"#initialize#1356"{Bool}, typeof(SciMLBase.FINALIZE_DEFAULT)}, DiscreteCallback{Trixi.var"#1362#1369"{Int64}, AnalysisCallback{Trixi.LobattoLegendreAnalyzer{Float64, 7, SVector{7, Float64}, Matrix{Float64}}, Tuple{typeof(Trixi.entropy_timederivative)}, SVector{1, Float64}, NamedTuple{(:u_local, :u_tmp1, :x_local, :x_tmp1), Tuple{StrideArraysCore.StaticStrideArray{Float64, 3, (1, 2, 3), Tuple{Static.StaticInt{1}, Static.StaticInt{7}, Static.StaticInt{7}}, Tuple{Nothing, Nothing, Nothing}, Tuple{Static.StaticInt{1}, Static.StaticInt{1}, Static.StaticInt{1}}, 49}, StrideArraysCore.StaticStrideArray{Float64, 3, (1, 2, 3), Tuple{Static.StaticInt{1}, Static.StaticInt{7}, Static.StaticInt{4}}, Tuple{Nothing, Nothing, Nothing}, Tuple{Static.StaticInt{1}, Static.StaticInt{1}, Static.StaticInt{1}}, 28}, StrideArraysCore.StaticStrideArray{Float64, 3, (1, 2, 3), Tuple{Static.StaticInt{2}, Static.StaticInt{7}, Static.StaticInt{7}}, Tuple{Nothing, Nothing, Nothing}, Tuple{Static.StaticInt{1}, Static.StaticInt{1}, Static.StaticInt{1}}, 98}, StrideArraysCore.StaticStrideArray{Float64, 3, (1, 2, 3), Tuple{Static.StaticInt{2}, Static.StaticInt{7}, Static.StaticInt{4}}, Tuple{Nothing, Nothing, Nothing}, Tuple{Static.StaticInt{1}, Static.StaticInt{1}, Static.StaticInt{1}}, 56}}}}, typeof(Trixi.initialize!), typeof(SciMLBase.FINALIZE_DEFAULT)}, DiscreteCallback{StepsizeCallback{Float64}, StepsizeCallback{Float64}, typeof(Trixi.initialize!), typeof(SciMLBase.FINALIZE_DEFAULT)}, DiscreteCallback{SaveSolutionCallback{Int64, typeof(cons2prim)}, SaveSolutionCallback{Int64, typeof(cons2prim)}, typeof(Trixi.initialize_save_cb!), typeof(SciMLBase.FINALIZE_DEFAULT)}, DiscreteCallback{SaveRestartCallback, SaveRestartCallback, typeof(Trixi.initialize!), typeof(SciMLBase.FINALIZE_DEFAULT)}}}((), (SummaryCallback, DiscreteCallback{Trixi.var"#1362#1369"{Int64}, AnalysisCallback{Trixi.LobattoLegendreAnalyzer{Float64, 7, SVector{7, Float64}, Matrix{Float64}}, Tuple{typeof(Trixi.entropy_timederivative)}, SVector{1, Float64}, NamedTuple{(:u_local, :u_tmp1, :x_local, :x_tmp1), Tuple{StrideArraysCore.StaticStrideArray{Float64, 3, (1, 2, 3), Tuple{Static.StaticInt{1}, Static.StaticInt{7}, Static.StaticInt{7}}, Tuple{Nothing, Nothing, Nothing}, Tuple{Static.StaticInt{1}, Static.StaticInt{1}, Static.StaticInt{1}}, 49}, StrideArraysCore.StaticStrideArray{Float64, 3, (1, 2, 3), Tuple{Static.StaticInt{1}, Static.StaticInt{7}, Static.StaticInt{4}}, Tuple{Nothing, Nothing, Nothing}, Tuple{Static.StaticInt{1}, Static.StaticInt{1}, Static.StaticInt{1}}, 28}, StrideArraysCore.StaticStrideArray{Float64, 3, (1, 2, 3), Tuple{Static.StaticInt{2}, Static.StaticInt{7}, Static.StaticInt{7}}, Tuple{Nothing, Nothing, Nothing}, Tuple{Static.StaticInt{1}, Static.StaticInt{1}, Static.StaticInt{1}}, 98}, StrideArraysCore.StaticStrideArray{Float64, 3, (1, 2, 3), Tuple{Static.StaticInt{2}, Static.StaticInt{7}, Static.StaticInt{4}}, Tuple{Nothing, Nothing, Nothing}, Tuple{Static.StaticInt{1}, Static.StaticInt{1}, Static.StaticInt{1}}, 56}}}}, typeof(Trixi.initialize!), typeof(SciMLBase.FINALIZE_DEFAULT)}(Trixi.var"#1362#1369"{Int64}(5), AnalysisCallback{Trixi.LobattoLegendreAnalyzer{Float64, 7, SVector{7, Float64}, Matrix{Float64}}, Tuple{typeof(Trixi.entropy_timederivative)}, SVector{1, Float64}, NamedTuple{(:u_local, :u_tmp1, :x_local, :x_tmp1), Tuple{StrideArraysCore.StaticStrideArray{Float64, 3, (1, 2, 3), Tuple{Static.StaticInt{1}, Static.StaticInt{7}, Static.StaticInt{7}}, Tuple{Nothing, Nothing, Nothing}, Tuple{Static.StaticInt{1}, Static.StaticInt{1}, Static.StaticInt{1}}, 49}, StrideArraysCore.StaticStrideArray{Float64, 3, (1, 2, 3), Tuple{Static.StaticInt{1}, Static.StaticInt{7}, Static.StaticInt{4}}, Tuple{Nothing, Nothing, Nothing}, Tuple{Static.StaticInt{1}, Static.StaticInt{1}, Static.StaticInt{1}}, 28}, StrideArraysCore.StaticStrideArray{Float64, 3, (1, 2, 3), Tuple{Static.StaticInt{2}, Static.StaticInt{7}, Static.StaticInt{7}}, Tuple{Nothing, Nothing, Nothing}, Tuple{Static.StaticInt{1}, Static.StaticInt{1}, Static.StaticInt{1}}, 98}, StrideArraysCore.StaticStrideArray{Float64, 3, (1, 2, 3), Tuple{Static.StaticInt{2}, Static.StaticInt{7}, Static.StaticInt{4}}, Tuple{Nothing, Nothing, Nothing}, Tuple{Static.StaticInt{1}, Static.StaticInt{1}, Static.StaticInt{1}}, 56}}}}(0.0, 0.0, 0, 0.0, 5, false, "out", "analysis.dat", LobattoLegendreAnalyzer{Float64}(polydeg=6), [:l2_error, :linf_error], (Trixi.entropy_timederivative,), [0.0], (u_local = [0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0], u_tmp1 = [0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0], x_local = [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], x_tmp1 = [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0])), Trixi.initialize!, SciMLBase.FINALIZE_DEFAULT, Bool[0, 0]), StepsizeCallback(cfl_number=1.6), SaveSolutionCallback(interval=5), SaveRestartCallback(interval=100)))

The last step is to choose the time integration method. OrdinaryDiffEq.jl defines a wide range of ODE solvers, e.g. CarpenterKennedy2N54(williamson_condition = false). We will pass the ODE problem, the ODE solver and the callbacks to the solve function. Also, to use StepsizeCallback, we must explicitly specify the initial trial time step dt, the selected value is not important, because it will be overwritten by the StepsizeCallback. And there is no need to save every step of the solution, we are only interested in the final result.

sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), dt = 1.0,
            save_everystep = false, callback = callbacks);

████████╗██████╗ ██╗██╗  ██╗██╗
╚══██╔══╝██╔══██╗██║╚██╗██╔╝██║
   ██║   ██████╔╝██║ ╚███╔╝ ██║
   ██║   ██╔══██╗██║ ██╔██╗ ██║
   ██║   ██║  ██║██║██╔╝ ██╗██║
   ╚═╝   ╚═╝  ╚═╝╚═╝╚═╝  ╚═╝╚═╝

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SemidiscretizationHyperbolic                                                                     │
│ ════════════════════════════                                                                     │
│ #spatial dimensions: ………………………… 2                                                                │
│ mesh: ………………………………………………………………… TreeMesh{2, Trixi.SerialTree{2}} with length 341                 │
│ equations: …………………………………………………… LinearScalarAdvectionEquation2D                                  │
│ initial condition: ……………………………… initial_condition_sinpi                                          │
│ boundary conditions: ………………………… Trixi.BoundaryConditionPeriodic                                  │
│ source terms: …………………………………………… source_term_exp_sinpi                                            │
│ solver: …………………………………………………………… DG                                                               │
│ total #DOFs per field: …………………… 4096                                                             │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ TreeMesh{2, Trixi.SerialTree{2}}                                                                 │
│ ════════════════════════════════                                                                 │
│ center: …………………………………………………………… [0.0, 0.0]                                                       │
│ length: …………………………………………………………… 2.0                                                              │
│ periodicity: ……………………………………………… (true, true)                                                     │
│ current #cells: ……………………………………… 341                                                              │
│ #leaf-cells: ……………………………………………… 256                                                              │
│ maximum #cells: ……………………………………… 30000                                                            │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ LinearScalarAdvectionEquation2D                                                                  │
│ ═══════════════════════════════                                                                  │
│ #variables: ………………………………………………… 1                                                                │
│ │ variable 1: …………………………………………… scalar                                                           │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ DG{Float64}                                                                                      │
│ ═══════════                                                                                      │
│ basis: ……………………………………………………………… LobattoLegendreBasis{Float64}(polydeg=3)                         │
│ mortar: …………………………………………………………… LobattoLegendreMortarL2{Float64}(polydeg=3)                      │
│ surface integral: ………………………………… SurfaceIntegralWeakForm                                          │
│ │ surface flux: ……………………………………… flux_central                                                     │
│ volume integral: …………………………………… VolumeIntegralWeakForm                                           │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ AnalysisCallback                                                                                 │
│ ════════════════                                                                                 │
│ interval: ……………………………………………………… 5                                                                │
│ analyzer: ……………………………………………………… LobattoLegendreAnalyzer{Float64}(polydeg=6)                      │
│ │ error 1: …………………………………………………… l2_error                                                         │
│ │ error 2: …………………………………………………… linf_error                                                       │
│ │ integral 1: …………………………………………… entropy_timederivative                                           │
│ save analysis to file: …………………… no                                                               │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ StepsizeCallback                                                                                 │
│ ════════════════                                                                                 │
│ CFL number: ………………………………………………… 1.6                                                              │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SaveSolutionCallback                                                                             │
│ ════════════════════                                                                             │
│ interval: ……………………………………………………… 5                                                                │
│ solution variables: …………………………… cons2prim                                                        │
│ save initial solution: …………………… yes                                                              │
│ save final solution: ………………………… yes                                                              │
│ output directory: ………………………………… /home/runner/work/Trixi.jl/Trixi…build/tutorials/first_steps/out │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SaveRestartCallback                                                                              │
│ ═══════════════════                                                                              │
│ interval: ……………………………………………………… 100                                                              │
│ save final solution: ………………………… yes                                                              │
│ output directory: ………………………………… /home/runner/work/Trixi.jl/Trixi…build/tutorials/first_steps/out │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Time integration                                                                                 │
│ ════════════════                                                                                 │
│ Start time: ………………………………………………… 0.0                                                              │
│ Final time: ………………………………………………… 1.0                                                              │
│ time integrator: …………………………………… CarpenterKennedy2N54                                             │
│ adaptive: ……………………………………………………… false                                                            │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Environment information                                                                          │
│ ═══════════════════════                                                                          │
│ #threads: ……………………………………………………… 1                                                                │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

────────────────────────────────────────────────────────────────────────────────────────────────────
 Simulation running 'LinearScalarAdvectionEquation2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
 #timesteps:                  0                run time:       1.22300000e-06 s
 Δt:             1.00000000e+00                └── GC time:    0.00000000e+00 s (0.000%)
 sim. time:      0.00000000e+00 (0.000%)       time/DOF/rhs!:         NaN s
                                               PID:                   Inf s
 #DOFs per field:          4096                alloc'd memory:       3805.493 MiB
 #elements:                 256

 Variable:       scalar
 L2 error:       5.57371016e-06
 Linf error:     2.37926999e-05
 ∑∂S/∂U ⋅ Uₜ :   2.40654766e-17
────────────────────────────────────────────────────────────────────────────────────────────────────


────────────────────────────────────────────────────────────────────────────────────────────────────
 Simulation running 'LinearScalarAdvectionEquation2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
 #timesteps:                  5                run time:       6.36953800e-03 s
 Δt:             5.55555556e-02                └── GC time:    0.00000000e+00 s (0.000%)
 sim. time:      2.77777778e-01 (27.778%)      time/DOF/rhs!:  3.92735460e-08 s
                                               PID:            5.30710637e-08 s
 #DOFs per field:          4096                alloc'd memory:       3805.653 MiB
 #elements:                 256

 Variable:       scalar
 L2 error:       3.53168531e-01
 Linf error:     8.67376253e-01
 ∑∂S/∂U ⋅ Uₜ :   3.19973014e-02
────────────────────────────────────────────────────────────────────────────────────────────────────


────────────────────────────────────────────────────────────────────────────────────────────────────
 Simulation running 'LinearScalarAdvectionEquation2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
 #timesteps:                 10                run time:       1.15899670e-02 s
 Δt:             5.55555556e-02                └── GC time:    0.00000000e+00 s (0.000%)
 sim. time:      5.55555556e-01 (55.556%)      time/DOF/rhs!:  3.90553260e-08 s
                                               PID:            4.41451270e-08 s
 #DOFs per field:          4096                alloc'd memory:       3805.805 MiB
 #elements:                 256

 Variable:       scalar
 L2 error:       6.17502929e-01
 Linf error:     1.33304028e+00
 ∑∂S/∂U ⋅ Uₜ :  -3.54326400e-03
────────────────────────────────────────────────────────────────────────────────────────────────────


────────────────────────────────────────────────────────────────────────────────────────────────────
 Simulation running 'LinearScalarAdvectionEquation2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
 #timesteps:                 15                run time:       1.67614560e-02 s
 Δt:             5.55555556e-02                └── GC time:    0.00000000e+00 s (0.000%)
 sim. time:      8.33333333e-01 (83.333%)      time/DOF/rhs!:  3.89381479e-08 s
                                               PID:            4.38904590e-08 s
 #DOFs per field:          4096                alloc'd memory:       3805.924 MiB
 #elements:                 256

 Variable:       scalar
 L2 error:       7.94341216e-01
 Linf error:     1.66424456e+00
 ∑∂S/∂U ⋅ Uₜ :  -3.28577696e-02
────────────────────────────────────────────────────────────────────────────────────────────────────


────────────────────────────────────────────────────────────────────────────────────────────────────
 Simulation running 'LinearScalarAdvectionEquation2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
 #timesteps:                 18                run time:       2.02872070e-02 s
 Δt:             5.55555556e-02                └── GC time:    0.00000000e+00 s (0.000%)
 sim. time:      1.00000000e+00 (100.000%)     time/DOF/rhs!:  3.87635803e-08 s
                                               PID:            4.64538737e-08 s
 #DOFs per field:          4096                alloc'd memory:       3806.042 MiB
 #elements:                 256

 Variable:       scalar
 L2 error:       8.65869257e-01
 Linf error:     1.75854798e+00
 ∑∂S/∂U ⋅ Uₜ :  -1.70992020e-02
────────────────────────────────────────────────────────────────────────────────────────────────────

Finally, we print the timer summary.

summary_callback()
 ────────────────────────────────────────────────────────────────────────────────
            Trixi.jl                    Time                    Allocations
                               ───────────────────────   ────────────────────────
       Tot / % measured:           47.3ms /  45.1%           2.76MiB /  25.3%

 Section               ncalls     time    %tot     avg     alloc    %tot      avg
 ────────────────────────────────────────────────────────────────────────────────
 rhs!                      91   14.5ms   68.1%   160μs   9.33KiB    1.3%     105B
   source terms            91   8.86ms   41.5%  97.3μs     0.00B    0.0%    0.00B
   volume integral         91   3.32ms   15.5%  36.5μs     0.00B    0.0%    0.00B
   interface flux          91   1.39ms    6.5%  15.3μs     0.00B    0.0%    0.00B
   surface integral        91    434μs    2.0%  4.77μs     0.00B    0.0%    0.00B
   prolong2interfaces      91    335μs    1.6%  3.68μs     0.00B    0.0%    0.00B
   ~rhs!~                  91   90.6μs    0.4%   995ns   9.33KiB    1.3%     105B
   Jacobian                91   49.4μs    0.2%   543ns     0.00B    0.0%    0.00B
   reset ∂u/∂t             91   40.9μs    0.2%   449ns     0.00B    0.0%    0.00B
   prolong2boundaries      91   4.44μs    0.0%  48.8ns     0.00B    0.0%    0.00B
   prolong2mortars         91   4.10μs    0.0%  45.1ns     0.00B    0.0%    0.00B
   mortar flux             91   3.76μs    0.0%  41.3ns     0.00B    0.0%    0.00B
   boundary flux           91   2.68μs    0.0%  29.4ns     0.00B    0.0%    0.00B
 analyze solution           5   3.44ms   16.1%   687μs    132KiB   18.4%  26.4KiB
 I/O                        8   3.38ms   15.8%   422μs    575KiB   80.3%  71.9KiB
   save solution            5   2.26ms   10.6%   453μs    505KiB   70.5%   101KiB
   ~I/O~                    8   1.11ms    5.2%   139μs   70.4KiB    9.8%  8.80KiB
   get element vari...      5    560ns    0.0%   112ns     0.00B    0.0%    0.00B
   save mesh                5    181ns    0.0%  36.2ns     0.00B    0.0%    0.00B
   get node variables       5    161ns    0.0%  32.2ns     0.00B    0.0%    0.00B
 calculate dt              19   5.11μs    0.0%   269ns     0.00B    0.0%    0.00B
 ────────────────────────────────────────────────────────────────────────────────

Now you can plot the solution as shown below, analyze it and improve the stability, accuracy or efficiency of your setup.

Visualize the solution

In the previous part of the tutorial, we calculated the final solution of the given problem, now we want to visualize it. A more detailed explanation of visualization methods can be found in the section Visualization.

Using Plots.jl

The first option is to use the Plots.jl package directly after calculations, when the solution is saved in the sol variable. We load the package and use the plot function.

using Plots
plot(sol)
Example block output

To show the mesh on the plot, we need to extract the visualization data from the solution as a PlotData2D object. Mesh extraction is possible using the getmesh function. Plots.jl has the plot! function that allows you to modify an already built graph.

pd = PlotData2D(sol)
plot!(getmesh(pd))
Example block output

Using Trixi2Vtk.jl

Another way to visualize a solution is to extract it from a saved HDF5 file. After we used the solve function with SaveSolutionCallback there is a file with the final solution. It is located in the out folder and is named as follows: solution_index.h5. The index is the final time step of the solution that is padded to 6 digits with zeros from the beginning. With Trixi2Vtk you can convert the HDF5 output file generated by Trixi.jl into a VTK file. This can be used in visualization tools such as ParaView or VisIt to plot the solution. The important thing is that currently Trixi2Vtk.jl supports conversion only for solutions in 2D and 3D spatial domains.

If you haven't added Trixi2Vtk.jl to your project yet, you can add it as follows.

import Pkg
Pkg.add(["Trixi2Vtk"])

Now we load the Trixi2Vtk.jl package and convert the file out/solution_000018.h5 with the final solution using the trixi2vtk function saving the resulting file in the out folder.

using Trixi2Vtk
trixi2vtk(joinpath("out", "solution_000018.h5"), output_directory="out")
 ────────────────────────────────────────────────────────────────────────────────
                                        Time                    Allocations
                               ───────────────────────   ────────────────────────
       Tot / % measured:            981ms /  45.3%           57.5MiB /  63.9%

 Section               ncalls     time    %tot     avg     alloc    %tot      avg
 ────────────────────────────────────────────────────────────────────────────────
 add data to VTK file       1    317ms   71.1%   317ms   9.34MiB   25.4%  9.34MiB
   scalar                   1    163ms   36.6%   163ms   5.70MiB   15.5%  5.70MiB
   add data to VTK ...      1    513μs    0.1%   513μs    118KiB    0.3%   118KiB
     cell_ids               1    275μs    0.1%   275μs   35.6KiB    0.1%  35.6KiB
     element_ids            1    199μs    0.0%   199μs   37.8KiB    0.1%  37.8KiB
     levels                 1   30.1μs    0.0%  30.1μs   37.8KiB    0.1%  37.8KiB
 read mesh                  1   62.9ms   14.1%  62.9ms   6.17MiB   16.8%  6.17MiB
 build VTK grid (no...      1   54.3ms   12.2%  54.3ms   15.4MiB   41.9%  15.4MiB
 prepare VTK cells ...      1   8.79ms    2.0%  8.79ms   5.04MiB   13.7%  5.04MiB
 build VTK grid (ce...      1   1.07ms    0.2%  1.07ms    353KiB    0.9%   353KiB
 save VTK file              2    469μs    0.1%   235μs   4.97KiB    0.0%  2.48KiB
 interpolate data           1    400μs    0.1%   400μs    297KiB    0.8%   297KiB
 read data                  1    314μs    0.1%   314μs   68.8KiB    0.2%  68.8KiB
 prepare VTK cells ...      1    233μs    0.1%   233μs    100KiB    0.3%   100KiB
 ────────────────────────────────────────────────────────────────────────────────

Now two files solution_000018.vtu and solution_000018_celldata.vtu have been generated in the out folder. The first one contains all the information for visualizing the solution, the second one contains all the cell-based or discretization-based information.

Now let's visualize the solution from the generated files in ParaView. Follow this short instruction to get the visualization.

  • Download, install and open ParaView.
  • Press Ctrl+O and select the generated files solution_000018.vtu and solution_000018_celldata.vtu from the out folder.
  • In the upper-left corner in the Pipeline Browser window, left-click on the eye-icon near solution_000018.vtu.
  • In the lower-left corner in the Properties window, change the Coloring from Solid Color to scalar. This already generates the visualization of the final solution.
  • Now let's add the mesh to the visualization. In the upper-left corner in the Pipeline Browser window, left-click on the eye-icon near solution_000018_celldata.vtu.
  • In the lower-left corner in the Properties window, change the Representation from Surface to Wireframe. Then a white grid should appear on the visualization.

Now, if you followed the instructions exactly, you should get a similar image as shown in the section Using Plots.jl:

paraview_trixi2vtk_example

After completing this tutorial you are able to set up your own simulations with Trixi.jl. If you have an interest in contributing to Trixi.jl as a developer, refer to the third part of the introduction titled Changing Trixi.jl itself.


This page was generated using Literate.jl.