Merge branch 'semidiscretization-doc' of https://github.com/ArseniyKh…

…olod/Trixi.jl into semidiscretization-doc

.github/workflows/FormatCheck.yml

@@ -30,7 +30,7 @@ jobs:
         #       format(".")
         run: |
           julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter"))'
-          julia  -e 'using JuliaFormatter; format(["benchmark", "ext", "src", "utils"])'
+          julia  -e 'using JuliaFormatter; format(["benchmark", "ext", "src", "test", "utils"])'
       - name: Format check
         run: |
           julia -e '

NEWS.md

@@ -12,7 +12,9 @@ for human readability.
 - Non-uniform `TreeMesh` available for hyperbolic-parabolic equations.
 - Capability to set truly discontinuous initial conditions in 1D.
 - Wetting and drying feature and examples for 1D and 2D shallow water equations
+- Implementation of the polytropic Euler equations in 2D
 - Subcell positivity limiting support for conservative variables in 2D for `TreeMesh`
+- AMR for hyperbolic-parabolic equations on 2D/3D `TreeMesh`
 
 #### Changed
 

Project.toml

@@ -1,7 +1,7 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor Gassner <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrew R. Winters <[email protected]>", "Jesse Chan <[email protected]>"]
-version = "0.5.45-pre"
+version = "0.5.48-pre"
 
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"

docs/src/callbacks.md

@@ -30,6 +30,11 @@ An example elixir using AMR can be found at [`examples/tree_2d_dgsem/elixir_adve
 The [`AnalysisCallback`](@ref) can be used to analyze the numerical solution, e.g. calculate
 errors or user-specified integrals, and print the results to the screen. The results can also be
 saved in a file. An example can be found at [`examples/tree_2d_dgsem/elixir_euler_vortex.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_euler_vortex.jl).
+Note that the errors (e.g. `L2 error` or `Linf error`) are computed with respect to the initial condition.
+The percentage of the simulation time refers to the ratio of the current time and the final time, i.e. it does
+not consider the maximal number of iterations. So the simulation could finish before 100% are reached.
+Note that, e.g., due to AMR or smaller time step sizes, the simulation can actually take longer than
+the percentage indicates.
 In [Performance metrics of the `AnalysisCallback`](@ref performance-metrics) you can find a detailed
 description of the different performance metrics the `AnalysisCallback` computes.
 

docs/src/parallelization.md

@@ -70,7 +70,32 @@ the same for P4est.jl and T8code.jl. This could e.g. be `libp4est.so` that usual
 in `lib/` or `local/lib/` in the installation directory of `t8code`.
 In total, in your active Julia project you should have a LocalPreferences.toml file with sections
 `[MPIPreferences]`, `[T8code]` and `[P4est]` as well as an entry `MPIPreferences` in your
-Project.toml to use a custom MPI installation.
+Project.toml to use a custom MPI installation. A `LocalPreferences.toml` file 
+created as described above might look something like the following:
+```toml
+[HDF5]
+libhdf5 = "/usr/lib/x86_64-linux-gnu/hdf5/openmpi/libhdf5.so"
+libhdf5_hl = "/usr/lib/x86_64-linux-gnu/hdf5/openmpi/libhdf5_hl.so"
+
+[MPIPreferences]
+__clear__ = ["preloads_env_switch"]
+_format = "1.0"
+abi = "OpenMPI"
+binary = "system"
+cclibs = []
+libmpi = "/lib/x86_64-linux-gnu/libmpi.so"
+mpiexec = "mpiexec"
+preloads = []
+
+[P4est]
+libp4est = "/home/mschlott/hackathon/libtrixi/t8code/install/lib/libp4est.so"
+libsc = "/home/mschlott/hackathon/libtrixi/t8code/install/lib/libsc.so"
+
+[T8code]
+libp4est = "/home/mschlott/hackathon/libtrixi/t8code/install/lib/libp4est.so"
+libsc = "/home/mschlott/hackathon/libtrixi/t8code/install/lib/libsc.so"
+libt8 = "/home/mschlott/hackathon/libtrixi/t8code/install/lib/libt8.so"
+```
 
 
 ### [Usage](@id parallel_usage)

examples/p4est_2d_dgsem/elixir_euler_double_mach_amr.jl

@@ -27,15 +27,15 @@ See Section IV c on the paper below for details.
     phi = pi / 6
     sin_phi, cos_phi = sincos(phi)
 
-    rho =  8
+    rho =  8.0
     v1  =  8.25 * cos_phi
     v2  = -8.25 * sin_phi
     p   =  116.5
   else
     rho = 1.4
-    v1  = 0
-    v2  = 0
-    p   = 1
+    v1  = 0.0
+    v2  = 0.0
+    p   = 1.0
   end
 
   prim = SVector(rho, v1, v2, p)

examples/structured_2d_dgsem/elixir_eulerpolytropic_convergence.jl

@@ -0,0 +1,63 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the polytropic Euler equations
+
+gamma = 1.4
+kappa = 0.5     # Scaling factor for the pressure.
+equations = PolytropicEulerEquations2D(gamma, kappa)
+
+initial_condition = initial_condition_convergence_test
+
+volume_flux = flux_winters_etal
+solver = DGSEM(polydeg = 3, surface_flux = flux_hll,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (1.0, 1.0)
+
+cells_per_dimension = (4, 4)
+
+mesh = StructuredMesh(cells_per_dimension,
+                      coordinates_min,
+                      coordinates_max)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.1)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_errors = (:l2_error_primitive,
+                                                              :linf_error_primitive))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.1)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+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);
+summary_callback() # print the timer summary

examples/structured_2d_dgsem/elixir_eulerpolytropic_ec.jl

@@ -0,0 +1,80 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the polytropic Euler equations
+
+gamma = 1.4
+kappa = 0.5     # Scaling factor for the pressure.
+equations = PolytropicEulerEquations2D(gamma, kappa)
+
+initial_condition = initial_condition_weak_blast_wave
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = flux_winters_etal
+solver = DGSEM(polydeg=4, surface_flux=flux_winters_etal,
+               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# Get the curved quad mesh from a mapping function
+
+# Mapping as described in https://arxiv.org/abs/2012.12040, but reduced to 2D
+function mapping(xi_, eta_)
+  # Transform input variables between -1 and 1 onto [0,3]
+  xi = 1.5 * xi_ + 1.5
+  eta = 1.5 * eta_ + 1.5
+
+  y = eta + 3/8 * (cos(1.5 * pi * (2 * xi - 3)/3) *
+                   cos(0.5 * pi * (2 * eta - 3)/3))
+
+  x = xi + 3/8 * (cos(0.5 * pi * (2 * xi - 3)/3) *
+                  cos(2 * pi * (2 * y - 3)/3))
+
+  return SVector(x, y)
+end
+
+cells_per_dimension = (16, 16)
+
+# Create curved mesh with 16 x 16 elements
+mesh = StructuredMesh(cells_per_dimension, mapping)
+
+###############################################################################
+# create the semi discretization object
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true)
+
+stepsize_callback = StepsizeCallback(cfl=1.0)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+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);
+summary_callback() # print the timer summary

examples/structured_2d_dgsem/elixir_eulerpolytropic_isothermal_wave.jl

@@ -0,0 +1,83 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the polytropic Euler equations
+
+gamma = 1.0     # With gamma = 1 the system is isothermal.
+kappa = 1.0     # Scaling factor for the pressure.
+equations = PolytropicEulerEquations2D(gamma, kappa)
+
+# Linear pressure wave in the negative x-direction.
+function initial_condition_wave(x, t, equations::PolytropicEulerEquations2D)
+    rho = 1.0
+    v1 = 0.0
+    if x[1] > 0.0
+        rho = ((1.0 + 0.01 * sin(x[1] * 2 * pi)) / equations.kappa)^(1 / equations.gamma)
+        v1 = ((0.01 * sin((x[1] - 1 / 2) * 2 * pi)) / equations.kappa)
+    end
+    v2 = 0.0
+
+    return prim2cons(SVector(rho, v1, v2), equations)
+end
+initial_condition = initial_condition_wave
+
+volume_flux = flux_winters_etal
+solver = DGSEM(polydeg = 2, surface_flux = flux_hll,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (-2.0, -1.0)
+coordinates_max = (2.0, 1.0)
+
+cells_per_dimension = (64, 32)
+
+boundary_conditions = (x_neg = boundary_condition_periodic,
+                       x_pos = boundary_condition_periodic,
+                       y_neg = boundary_condition_periodic,
+                       y_pos = boundary_condition_periodic)
+
+mesh = StructuredMesh(cells_per_dimension,
+                      coordinates_min,
+                      coordinates_max)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 200
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 50,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 1.7)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (1.0e-4, 1.0e-4),
+                                                     variables = (Trixi.density, pressure))
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, 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()