diff --git a/docs/literate/src/files/innards_of_the_basic_setup.jl b/docs/literate/src/files/innards_of_the_basic_setup.jl new file mode 100644 index 00000000000..5193b8800dc --- /dev/null +++ b/docs/literate/src/files/innards_of_the_basic_setup.jl @@ -0,0 +1,244 @@ +#src # Core aspects of the basic setup + +# This tutorial aims to provide insights into the Trixi.jl framework, with a focus on commonly +# used complex functions. These functions may lack clarity in their names and brief documentation +# descriptions. + +# We will guide you through a simplified setup, emphasizing functionalities that may not be evident +# to an average user. This setup is based on stable parts of Trixi.jl that are unlikely to undergo +# significant changes in the near future. Nevertheless, it will clarify fundamental concepts that +# continue to be applied in more recent and flexible configurations. + + +# ## Basic setup + +# Import essential libraries and specify an equation. + +using Trixi, OrdinaryDiffEq +equations=LinearScalarAdvectionEquation2D((-0.2,0.7)) + +# Generate a geometrical instruction for spatial discretization using [`TreeMesh`](@ref) with +# pre-coarsened cell. + +coordinates_min = (-2.0, -2.0) +coordinates_max = (2.0, 2.0) + +coarsening_patches = ( # Coarsen cell in the lower-right quarter + (type="box", coordinates_min=[0.0, -2.0], coordinates_max=[2.0, 0.0]), +) + + mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level=2, + n_cells_max=30_000, + coarsening_patches=coarsening_patches) + +# Created `TreeMesh` looks like: + +# ![TreeMesh_example](https://github.com/ArseniyKholod/Test/assets/119304909/df675af6-9eba-4e18-9e68-454a401b1292) + +# Instantiate a [`DGSEM`](@ref) solver with a user-specified polynomial degree. The solver's +# objective is to define the `polydeg + 1` Gauss-Lobatto nodes, and their associated weights within +# the reference interval [-1, 1]. These nodes will be subsequently used to approximate solutions +# on each leaf cell of the `TreeMesh`. + +solver = DGSEM(polydeg=3) + +# Gauss-Lobatto nodes with `N=4`: + +# ![Gauss-Lobatto_nodes_example](https://github.com/ArseniyKholod/Test/assets/119304909/ce7ecdc8-23b6-47b6-9f3d-3ff58b4e5677) + + +# ## Overview of the [`SemidiscretizationHyperbolic`](@ref) function + +# At this stage, all the necessary components for configuring domain discretization are in place. +# The remaining task is to combine these components into a single structure that will be used +# throughout the entire solving process. This is where [`SemidiscretizationHyperbolic`](@ref) comes +# into play. + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver) + +# The `SemidiscretizationHyperbolic` function calls numerous sub-functions to perform the necessary +# steps. A brief description of the key sub-functions is provided below. + +# - `init_elements(leaf_cell_ids, mesh, equations, dg.basis, RealT, uEltype)` + +# As was already mentioned, the fundamental elements for approximating a solution are the leaf +# cells. This implies that on each leaf cell, the solution is treated as a polynomial of the +# degree specified in the `DGSEM` solver and evaluated at the Gauss-Lobatto nodes, which were +# previously illustrated. To provide this, the `init_elements` function extracts these leaf cells +# from the `TreeMesh`, assigns them the label "elements", records their coordinates, and maps the +# Gauss-Lobatto nodes from the 1D interval [-1, 1] onto each axis of every element. + +# ![elements_example](https://github.com/ArseniyKholod/Test/assets/119304909/637a6efb-9e04-47bd-9dd0-dc4b8ee48dd1) + +# The visualization of elements with nodes includes spaces between elements, which do not exist +# in reality. This spacing is included only for illustrative purposes to underscore the +# separation of elements and the independent projection of nodes onto each element. + +# - `init_interfaces(leaf_cell_ids, mesh, elements)` + +# At this point, the elements with nodes have been defined; however, they lack the necessary +# communication functionality. This is crucial because the solutions on the elements are not +# independent of each other. Furthermore, nodes on the boundary of adjacent elements share +# the same spatial location, requiring a method to combine their solutions. + +# As demonstrated earlier, the elements can have varying sizes. Let us initially consider +# neighbors with equal size. For these elements, the `init_interfaces` function generates +# interfaces that store information about adjacent elements, their relative positions, and +# allocate containers for sharing solutions between neighbors during the solving process. + +# In our visualization, these interfaces would conceptually resemble tubes connecting the +# corresponding elements. + +# ![interfaces_example](https://github.com/ArseniyKholod/Test/assets/119304909/4f6c7835-c98e-4e46-a025-2df81b28883b) + +# - `init_mortars(leaf_cell_ids, mesh, elements, dg.mortar)` + +# Returning to the consideration of different sizes among adjacent elements, within the +# `TreeMesh`, adjacent leaf cells can vary in size by a maximum factor of two. This implies +# that a coarsened element on each side, excluding the domain boundaries, has one neighbor of +# equal size with a connection through an interface, or two neighbors with sizes two times +# smaller, requiring a connection through so called "mortars". + +# Mortars store information about the connected elements, their relative positions, and allocate +# containers for sharing solutions between these elements along their boundaries. + +# In our visualization, mortars are represented as branched tubes. + +# ![mortars_example](https://github.com/ArseniyKholod/Test/assets/119304909/1c603899-c02c-4923-85b9-75d0e9e0426e) + +# - `init_boundaries(leaf_cell_ids, mesh, elements)` + +# In order to apply boundary conditions, it is necessary to identify the locations of the +# boundaries. Therefore, we initialize a "boundaries" object, which records the elements that +# contain boundaries, specifies which side of an element is a boundary, stores the coordinates +# of boundary nodes, and allocates containers for managing solutions at these boundaries. + +# In our visualization, boundaries and their corresponding nodes are highlighted with green, +# semi-transparent lines. + +# ![boundaries_example](https://github.com/ArseniyKholod/Test/assets/119304909/db0c6d62-e520-49a1-925f-049e3a378841) + +# All the structures mentioned earlier are packed as a cache of type `Tuple`. Subsequently, an +# object of type `SemidiscretizationHyperbolic` is initialized using this cache, initial and +# boundary conditions, equations, mesh and solver. + +# In conclusion, `Semidiscretization`'s primary function is to collect equations, the geometric +# representation of the domain, and approximation instructions, creating specialized structures to +# interconnect these components in a manner that enables their utilization for the numerical +# solution of partial differential equations (PDEs). + +# As evident from the earlier description of `SemidiscretizationHyperbolic`, it comprises numerous +# functions called recursively. Without delving into details, the structure of the primary calls +# can be illustrated as follows: + +# ![SemidiscretizationHyperbolic_structure](https://github.com/ArseniyKholod/Test/assets/119304909/ac95eb2e-d963-4b16-be1c-65dd731368bd) + + +# ## Overview of the [`semidiscretize`](@ref) function + +# At this stage, we have defined the equations and configured the domain's discretization. The +# final step before solving is to select a suitable time span and apply the corresponding initial +# conditions, which are already stored in the initialized `SemidiscretizationHyperbolic` object. + +# The purpose of the `semidiscretize` function is to wrap the semi-discretization as an +# `ODEProblem` within the specified time interval, while also applying the initial conditions at +# the left boundary of this interval. This `ODEProblem` can be subsequently passed to the `solve` +# function from the OrdinaryDiffEq.jl package. + +ode = semidiscretize(semi, (0.0, 1.0)); + +# The `semidiscretize` function involves a deep tree of recursive calls, with the primary ones +# explained below. + +# - `allocate_coefficients(mesh, equations, solver, cache)` + +# To apply initial conditions, a container needs to be generated to store the initial values of +# the target variables for each node within each element. The `allocate_coefficients` function +# initializes `u_ode` as a 1D zero-vector with a length that depends on the number of variables, +# elements, nodes, and dimensions. The use of a 1D vector format is consistent with the +# requirements of the `solve` function from OrdinaryDiffEq.jl. Therefore, Trixi.jl follows this +# format to be able to utilize the functionalities of OrdinaryDiffEq.jl. + +# - `wrap_array(u_ode, semi)` + +# As previously noted, `u_ode` is constructed as a 1D vector to ensure compatibility with +# OrdinaryDiffEq.jl. However, for internal use within Trixi.jl, identification which part of the +# vector relates to specific variables, elements and nodes can be challenging. + +# This is why the `u_ode` vector is wrapped by the `wrap_array` function to create a +# multidimensional array `u`, with each dimension representing variables, nodes and elements. +# Consequently, navigation within this multidimensional array becomes noticeably easier. + +# "Wrapping" in this context involves the creation of a reference to the same storage location +# but with an alternative structural representation. This approach enables the use of both +# instances `u` and `u_ode` as needed, so that changes are simultaneously reflected in both. +# This is possible because, from a storage perspective, they share the same stored data, while +# access to this data is provided in different ways. + +# - `compute_coefficients!(u, initial_conditions, t, mesh::DG, equations, solver, cache)` + +# Now the variable `u`, intended to store solutions, has been allocated and wrapped, it is time +# to apply the initial conditions. The `compute_coefficients!` function calculates the initial +# conditions for each variable at every node within each element and properly stores them in the +# `u` array. + +# At this stage, the `semidiscretize` function has all the necessary components to initialize and +# return an `ODEProblem` object, which will be used by the `solve` function to compute the +# solution. + +# In summary, the internal workings of `semidiscretize` with brief descriptions can be presented +# as follows. + +# ![semidiscretize_structure](https://github.com/ArseniyKholod/Test/assets/119304909/fa775789-7b4e-4909-9a44-e43b5b3a4a75) + + +# ## Functions `solve` and `rhs!` + +# Once the `ODEProblem` object is initialized, the `solve` function from the OrdinaryDiffEq.jl +# package can be utilized to compute an approximate solution using the instructions contained in +# the `ODEProblem` object. + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), dt=0.01, save_everystep=false); + +# Since the `solve` function is defined in another package without knowledge of how to handle +# discretizations performed in Trixi.jl, it is necessary to define the right-hand-side function, +# `rhs!`, within Trixi.jl. + +# Trixi.jl includes a set of `rhs!` functions designed to compute `du` according to the structure +# of the setup. These `rhs!` functions calculate interface, mortars, and boundary fluxes, in +# addition to surface and volume integrals, in order to construct the `du` vector. This `du` vector +# is then used by the time integration method to derive solutions in the subsequent time step. +# The `rhs!` function is called by time integration methods in each iteration of the solve-loop +# within OrdinaryDiffEq.jl, with arguments `du`, `u`, `Semidiscretization`, and the time step. + +# The problem is that `rhs!` functions within Trixi.jl are specialized for specific solver and mesh +# types. However, the types of arguments passed to `rhs!` by time integration methods do not +# explicitly provide this information. Consequently, Trixi.jl uses a two-levels approach for `rhs!` +# functions. The first level is limited to a single function for each `semidiscretization` type, +# and its role is to redirect data to the target `rhs!`. It performs this by extracting the +# necessary data from the integrator and passing them, along with the originally received +# arguments, to the specialized for solver and mesh types `rhs!` function, which is +# responsible for calculating `du`. + +# Path from the `solve` function call to the appropriate `rhs!` function call: + +# ![rhs_structure](https://github.com/ArseniyKholod/Test/assets/119304909/6603223c-8287-41f9-b609-1dcff7a5a8eb) + +# Computed solution: + +using Plots +plot(sol) +pd = PlotData2D(sol) +plot!(getmesh(pd)) + + +# ## Package versions + +# These results were obtained using the following versions. + +using InteractiveUtils +versioninfo() + +using Pkg +Pkg.status(["Trixi", "OrdinaryDiffEq", "Plots", "ForwardDiff"], + mode=PKGMODE_MANIFEST) diff --git a/docs/make.jl b/docs/make.jl index df8ac04be12..3ffce0a1666 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -68,7 +68,8 @@ files = [ # Topic: other stuff "Explicit time stepping" => "time_stepping.jl", "Differentiable programming" => "differentiable_programming.jl", - "Custom semidiscretizations" => "custom_semidiscretization.jl" + "Custom semidiscretizations" => "custom_semidiscretization.jl", + "Core aspects of the basic setup" => "innards_of_the_basic_setup.jl", ] tutorials = create_tutorials(files)