Initial support for P4estMesh

.gitignore

@@ -8,6 +8,7 @@
 *.avi
 *.ogv
 *.mesh
+*.inp
 **/Manifest.toml
 out*/
 docs/build

Project.toml

@@ -14,6 +14,7 @@ LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+P4est = "7d669430-f675-4ae7-b43e-fab78ec5a902"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
@@ -34,6 +35,7 @@ LinearMaps = "2.7, 3.0"
 LoopVectorization = "0.12.22"
 MPI = "0.16, 0.17, 0.18"
 OffsetArrays = "1.3"
+P4est = "0.2.1"
 Polyester = "0.3"
 RecipesBase = "1.1"
 Reexport = "1.0"

docs/src/development.md

@@ -17,7 +17,7 @@ Julia can be avoided by using the [Revise.jl](https://github.com/timholy/Revise.
 package, which tracks changed files and re-loads them automatically. Therefore,
 it is *highly recommended* to first install Revise with the following command in Julia:
 To enter the package REPL mode, press `]` in the standard Julia REPL mode. Then, execute
-```julia
+```julia-repl
 (@v1.6) pkg> add Revise
 ```
 Now you are able to run Trixi from the REPL, change Trixi code between runs,
@@ -27,7 +27,7 @@ Revise regularly, please be aware of some of the [Pitfalls when using Revise](@r
 Another recommended package for working from the REPL is
 [OhMyREPL.jl](https://github.com/KristofferC/OhMyREPL.jl). It can be installed
 by running
-```julia
+```julia-repl
 (@v1.6) pkg> add OhMyREPL
 ```
 and adds syntax highlighting, bracket highlighting, and other helpful
@@ -42,11 +42,11 @@ begin by executing:
 julia
 ```
 This will start the Julia REPL. Then, run
-```julia
+```julia-repl
 julia> using Revise; using Trixi
 ```
 You can run a simulation by executing
-```julia
+```julia-repl
 julia> trixi_include(default_example())
 ```
 Together, all of these commands can take some time, roughly half a minute on a
@@ -67,7 +67,7 @@ For example, to run Trixi this way, you need to start the REPL with
 julia --project=path/to/Trixi.jl/
 ```
 and execute
-```julia
+```julia-repl
 julia> using Revise; using Trixi
 ```
 to load Revise and Trixi. You can then proceed with the usual commands and run Trixi as in
@@ -129,7 +129,7 @@ than can increase your productivity in the Julia REPL.
 
 - Use the [REPL help mode](https://docs.julialang.org/en/v1/stdlib/REPL/#Help-mode)
   entered by typing `?`.
-  ```julia
+  ```julia-repl
   julia> using Trixi
 
   help?> trixi_include
@@ -154,25 +154,25 @@ than can increase your productivity in the Julia REPL.
 - You can copy and paste REPL history including `julia>` prompts into the REPL.
 - Use tab completion in the REPL, both for names of functions/types/variables and
   for function arguments.
-  ```julia
+  ```julia-repl
   julia> flux_ranocha( # and TAB
   flux_ranocha(u_ll, u_rr, orientation, equations::CompressibleEulerEquations1D) in Trixi at ~/.julia/dev/Trixi/src/equations/1d/compressible_euler.jl:416
   flux_ranocha(u_ll, u_rr, orientation, equations::CompressibleEulerEquations2D) in Trixi at ~/.julia/dev/Trixi/src/equations/2d/compressible_euler.jl:865
   flux_ranocha(u_ll, u_rr, orientation, equations::CompressibleEulerEquations3D) in Trixi at ~/.julia/dev/Trixi/src/equations/3d/compressible_euler.jl:710
   ```
 - Use `methodswith` to discover methods associated to a given type etc.
-  ```julia
+  ```julia-repl
   julia> methodswith(CompressibleEulerEquations2D)
   [1] initial_condition_convergence_test(x, t, equations::CompressibleEulerEquations2D) in Trixi at ~/.julia/dev/Trixi/src/equations/2d/compressible_euler.jl:38
   [...]
   ```
 - Use `@which` (or `@edit`) for method calls.
-  ```julia
+  ```julia-repl
   julia> @which trixi_include(default_example())
   trixi_include(elixir::AbstractString; kwargs...) in Trixi at ~/.julia/dev/Trixi/src/run.jl:72
   ```
 - Use `apropos` to search through the documentation and docstrings.
-  ```julia
+  ```julia-repl
   julia> apropos("MHD")
   Trixi.initial_condition_constant
   Trixi.initial_condition_rotor
@@ -268,3 +268,16 @@ the new documentation are generated at `https://trixi-framework.github.io/Trixi.
 where `XXX` is the number of the PR.
 This does not work for PRs from forks for security reasons (since anyone could otherwise push
 arbitrary stuff to the Trixi website, including malicious code).
+
+
+
+## Developing Trixi2Vtk (@id trixi2vtk-dev)
+
+Trixi2Vtk has Trixi as dependency and uses Trixi's implementation to, e.g., load mesh files.
+When developing Trixi2Vtk, one may want to change functions in Trixi to allow them to be reused
+in Trixi2Vtk.
+To use a locally modified Trixi clone instead of a Trixi release, one can tell Pkg
+to use the local source code of Trixi instead of a registered version by running
+```julia-repl
+(@v1.6) pkg> develop path/to/Trixi.jl
+```

docs/src/visualization.md

@@ -179,7 +179,7 @@ In a very similar fashion to [`PlotData2D`](@ref), you can customize your plot:
 * `plot(pd["rho", "p"])` only plots specific variables. In this case `rho` and `p`.
 * `plot!(getmesh(pd))` adds mesh lines after creating a plot.
 * Any attributes from [Plots](https://docs.juliaplots.org/latest/) can be used, e.g., `plot(pd, yguide=:temperature)`.
-* `pd = PlotData1D(adapt_to_mesh_level(sol, 4)...)` adapts the mesh before plotting 
+* `pd = PlotData1D(adapt_to_mesh_level(sol, 4)...)` adapts the mesh before plotting
   (in this example to a mesh with refinement level 4).
 
 You can also customize the [`PlotData1D`](@ref) object itself by passing attributes
@@ -317,6 +317,9 @@ uses the `time` attribute in solution/restart files to inform ParaView about the
 solution time. A comprehensive list of all possible arguments for
 `trixi2vtk` can be found in the [Trixi2Vtk.jl API](@ref).
 
+Further information regarding the development of Trixi2Vtk can be found in the
+[development section](@id trixi2vtk-dev).
+
 
 ## Trixi2Img
 !!! note "Trixi2Img is deprecated"

examples/2d/elixir_advection_basic_p4est.jl

@@ -0,0 +1,61 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advectionvelocity = (1.0, 1.0)
+equations = LinearScalarAdvectionEquation2D(advectionvelocity)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+
+coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max = ( 1.0,  1.0) # maximum coordinates (max(x), max(y))
+
+trees_per_dimension = (8, 8)
+
+# Create P4estMesh with 8 x 8 trees and 16 x 16 elements
+mesh = P4estMesh(trees_per_dimension, polydeg=3,
+                 coordinates_min=coordinates_min, coordinates_max=coordinates_max,
+                 initial_refinement_level=1)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+ode = semidiscretize(semi, (0.0, 1.0));
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval=100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval=100,
+                                     solution_variables=cons2prim)
+
+# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl=1.6)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks);
+
+# Print the timer summary
+summary_callback()

examples/2d/elixir_advection_extended_p4est.jl

@@ -0,0 +1,88 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advectionvelocity = (1.0, 1.0)
+# advectionvelocity = (0.2, -0.3)
+equations = LinearScalarAdvectionEquation2D(advectionvelocity)
+
+initial_condition = initial_condition_convergence_test
+
+# you can either use a single function to impose the BCs weakly in all
+# 1*ndims == 2 directions or you can pass a tuple containing BCs for
+# each direction
+# Note: "boundary_condition_periodic" indicates that it is a periodic boundary and can be omitted on
+#       fully periodic domains. Here, however, it is included to allow easy override during testing
+boundary_conditions = boundary_condition_periodic
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+
+# The initial condition is 2-periodic
+coordinates_min = (-1.5, 1.3) # minimum coordinates (min(x), min(y))
+coordinates_max = ( 0.5, 5.3) # maximum coordinates (max(x), max(y))
+
+trees_per_dimension = (19, 37)
+
+# Create curved mesh with 19 x 37 elements
+mesh = P4estMesh(trees_per_dimension, polydeg=3,
+                 coordinates_min=coordinates_min, coordinates_max=coordinates_max)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions=boundary_conditions)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan);
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
+                                     extra_analysis_integrals=(entropy, energy_total))
+
+# The AliveCallback prints short status information in regular intervals
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+# The SaveRestartCallback allows to save a file from which a Trixi simulation can be restarted
+save_restart = SaveRestartCallback(interval=100,
+                                   save_final_restart=true)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
+cfl=1.6
+stepsize_callback = StepsizeCallback(cfl=cfl)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_restart, save_solution,
+                        stepsize_callback)
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks, maxiters=1e5);
+
+# Print the timer summary
+summary_callback()

examples/2d/elixir_advection_p4est_non_conforming.jl

@@ -0,0 +1,82 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advectionvelocity = (1.0, 1.0)
+equations = LinearScalarAdvectionEquation2D(advectionvelocity)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+
+# Deformed rectangle that looks like a waving flag,
+# lower and upper faces are sinus curves, left and right are vertical lines.
+f1(s) = SVector(-1.0, s - 1.0)
+f2(s) = SVector( 1.0, s + 1.0)
+f3(s) = SVector(s, -1.0 + sin(0.5 * pi * s))
+f4(s) = SVector(s,  1.0 + sin(0.5 * pi * s))
+
+trees_per_dimension = (3, 2)
+
+# Create P4estMesh with 8 x 8 trees and 16 x 16 elements
+mesh = P4estMesh(trees_per_dimension, polydeg=3,
+                 faces=(f1, f2, f3, f4),
+                 initial_refinement_level=1)
+
+# Refine bottom left quadrant of each forest to level 4
+function refine_fn(p4est, which_tree, quadrant)
+  if quadrant.x == 0 && quadrant.y == 0 && quadrant.level < 4
+    # return true (refine)
+    return Cint(1)
+  else
+    # return false (don't refine)
+    return Cint(0)
+  end
+end
+
+# Refine recursively until each bottom left quadrant of a forest has level 4
+# The mesh will be rebalanced before the simulation starts
+refine_fn_c = @cfunction(refine_fn, Cint, (Ptr{Trixi.p4est_t}, Ptr{Trixi.p4est_topidx_t}, Ptr{Trixi.p4est_quadrant_t}))
+Trixi.p4est_refine(mesh.p4est, true, refine_fn_c, C_NULL)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+ode = semidiscretize(semi, (0.0, 0.2));
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval=100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval=100,
+                                     solution_variables=cons2prim)
+
+# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl=1.6)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks);
+
+# Print the timer summary
+summary_callback()