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_amr_p4est_unstructured_flag.jl

@@ -0,0 +1,95 @@
+
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advectionvelocity = (0.2, -0.3)
+equations = LinearScalarAdvectionEquation2D(advectionvelocity)
+
+initial_condition = initial_condition_gauss
+
+boundary_condition = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = Dict(
+  :all => boundary_condition
+)
+
+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(-5.0, 5 * s - 5.0)
+f2(s) = SVector( 5.0, 5 * s + 5.0)
+f3(s) = SVector(5 * s, -5.0 + 5 * sin(0.5 * pi * s))
+f4(s) = SVector(5 * s,  5.0 + 5 * sin(0.5 * pi * s))
+faces = (f1, f2, f3, f4)
+
+# This creates a mapping that transforms [-1, 1]^2 to the domain with the faces defined above.
+# It generally doesn't work for meshes loaded from mesh files because these can be meshes
+# of arbitrary domains, but the mesh below is specifically built on the domain [-1, 1]^2.
+Trixi.validate_faces(faces)
+mapping_flag = Trixi.transfinite_mapping(faces)
+
+# Unstructured mesh with 24 cells of the square domain [-1, 1]^n
+mesh_file = joinpath(@__DIR__, "square_unstructured_2.inp")
+isfile(mesh_file) || download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp",
+                              mesh_file)
+
+mesh = P4estMesh(mesh_file, polydeg=3,
+                 mapping=mapping_flag,
+                 initial_refinement_level=1)
+
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions=boundary_conditions)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 10.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
+                                     extra_analysis_integrals=(entropy,))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_restart = SaveRestartCallback(interval=100,
+                                   save_final_restart=true)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first),
+                                      base_level=1,
+                                      med_level=2, med_threshold=0.1,
+                                      max_level=3, max_threshold=0.6)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval=5,
+                           adapt_initial_condition=true,
+                           adapt_initial_condition_only_refine=true)
+
+stepsize_callback = StepsizeCallback(cfl=0.7)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_restart, save_solution,
+                        amr_callback, stepsize_callback);
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks);
+
+summary_callback() # print the timer summary

examples/2d/elixir_advection_amr_solution_independent.jl

@@ -17,8 +17,8 @@ function IndicatorSolutionIndependent(semi)
   return IndicatorSolutionIndependent{typeof(cache)}(cache)
 end
 function (indicator::IndicatorSolutionIndependent)(u::AbstractArray{<:Any,4},
-                                             equations, dg, cache;
-                                             t, kwargs...)
+                                                   equations, dg, cache;
+                                                   t, kwargs...)
 
   mesh = indicator.cache.mesh
   alpha = indicator.cache.alpha
@@ -86,8 +86,8 @@ initial_condition = initial_condition_gauss
 
 solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
 
-coordinates_min = (-5, -5)
-coordinates_max = ( 5,  5)
+coordinates_min = (-5.0, -5.0)
+coordinates_max = ( 5.0,  5.0)
 mesh = TreeMesh(coordinates_min, coordinates_max,
                 initial_refinement_level=4,
                 n_cells_max=30_000)