Increase coverage with P4estMesh (#635)

* Add test for non-periodic structured P4estMesh

* Add Euler gravity test with P4estMesh

* Remove nboundaries(UnstructuredQuadSortedBoundaryTypes)

* Move p4est Euler gravity tests to other p4est tests

* Add news entry for P4estMesh

NEWS.md

@@ -19,6 +19,7 @@ for human readability.
   wave speed estimates were added in [#493](https://github.com/trixi-framework/Trixi.jl/pull/493)
 - New structured, curvilinear, conforming mesh type `CurvedMesh` (experimental)
 - New unstructured, curvilinear, conforming mesh type `UnstructuredQuadMesh` in 2D (experimental)
+- New unstructured, curvilinear, adaptive (non-conforming) mesh type `P4estMesh` in 2D (experimental)
 
 #### Changed
 

examples/2d/elixir_euler_nonperiodic_p4est.jl

@@ -15,20 +15,25 @@ source_terms = source_terms_convergence_test
 # BCs must be passed as Dict
 boundary_condition = BoundaryConditionDirichlet(initial_condition)
 boundary_conditions = Dict(
-  :all => boundary_condition
+  :x_neg => boundary_condition,
+  :x_pos => boundary_condition,
+  :y_neg => boundary_condition,
+  :y_pos => boundary_condition
 )
 
 solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
 
 ###############################################################################
 # Get the uncurved mesh from a file (downloads the file if not available locally)
 
-# 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)
+coordinates_min = (0.0, 0.0)
+coordinates_max = (2.0, 2.0)
 
-mesh = P4estMesh(mesh_file, initial_refinement_level=0)
+trees_per_dimension = (4, 4)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg=1, initial_refinement_level=2,
+                 coordinates_min=coordinates_min, coordinates_max=coordinates_max,
+                 periodicity=false)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                     source_terms=source_terms,
@@ -38,7 +43,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
 ###############################################################################
 # ODE solvers, callbacks etc.
 
-tspan = (0.0, 1.0)
+tspan = (0.0, 2.0)
 ode = semidiscretize(semi, tspan)
 
 summary_callback = SummaryCallback()

examples/2d/elixir_euler_nonperiodic_p4est_unstructured.jl

@@ -0,0 +1,71 @@
+
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+initial_condition = initial_condition_convergence_test
+
+source_terms = source_terms_convergence_test
+
+# BCs must be passed as Dict
+boundary_condition = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = Dict(
+  :all => boundary_condition
+)
+
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+
+###############################################################################
+# Get the uncurved mesh from a file (downloads the file if not available locally)
+
+# 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, initial_refinement_level=0)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms=source_terms,
+                                    boundary_conditions=boundary_conditions)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.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_restart = SaveRestartCallback(interval=100,
+                                   save_final_restart=true)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl=1.0)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_restart, 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/2d/elixir_eulergravity_eoc_test_p4est.jl

@@ -0,0 +1,87 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+
+initial_condition = initial_condition_eoc_test_coupled_euler_gravity
+
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+gamma = 2.0
+equations_euler = CompressibleEulerEquations2D(gamma)
+
+polydeg = 3
+solver_euler = DGSEM(polydeg, flux_hll)
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (2.0, 2.0)
+
+trees_per_dimension = (1, 1)
+mesh = P4estMesh(trees_per_dimension, polydeg=1,
+                 coordinates_min=coordinates_min, coordinates_max=coordinates_max,
+                 initial_refinement_level=2)
+
+semi_euler = SemidiscretizationHyperbolic(mesh, equations_euler, initial_condition, solver_euler,
+                                         source_terms=source_terms_eoc_test_coupled_euler_gravity)
+
+
+###############################################################################
+# semidiscretization of the hyperbolic diffusion equations
+equations_gravity = HyperbolicDiffusionEquations2D()
+
+solver_gravity = DGSEM(polydeg, flux_lax_friedrichs)
+
+semi_gravity = SemidiscretizationHyperbolic(mesh, equations_gravity, initial_condition, solver_gravity,
+                                            source_terms=source_terms_harmonic)
+
+
+###############################################################################
+# combining both semidiscretizations for Euler + self-gravity
+parameters = ParametersEulerGravity(background_density=2.0, # aka rho0
+                                    # rho0 is (ab)used to add a "+8π" term to the source terms
+                                    # for the manufactured solution
+                                    gravitational_constant=1.0, # aka G
+                                    cfl=1.1,
+                                    resid_tol=1.0e-10,
+                                    n_iterations_max=1000,
+                                    timestep_gravity=timestep_gravity_erk52_3Sstar!)
+
+semi = SemidiscretizationEulerGravity(semi_euler, semi_gravity, parameters)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan);
+
+summary_callback = SummaryCallback()
+
+stepsize_callback = StepsizeCallback(cfl=0.8)
+
+save_solution = SaveSolutionCallback(interval=10,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+save_restart = SaveRestartCallback(interval=100,
+                                   save_final_restart=true)
+
+analysis_interval = 100
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+analysis_callback = AnalysisCallback(semi_euler, interval=analysis_interval,
+                                     save_analysis=true)
+
+callbacks = CallbackSet(summary_callback, stepsize_callback,
+                        save_restart, save_solution,
+                        analysis_callback, alive_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
+println("Number of gravity subcycles: ", semi.gravity_counter.ncalls_since_readout)

src/solvers/dg_unstructured_quad/sort_boundary_conditions.jl

@@ -15,8 +15,6 @@ mutable struct UnstructuredQuadSortedBoundaryTypes{N, BCs<:NTuple{N, Any}}
   boundary_dictionary::Dict{Symbol, Any} # boundary conditions as set by the user in the elixir file
 end
 
-@inline nboundaries(container::UnstructuredQuadSortedBoundaryTypes) = sum(length, container.boundary_indices)
-
 
 # constructor that "eats" the original boundary condition dictionary and sorts the information
 # from the `UnstructuredBoundaryContainer2D` in cache.boundaries according to the boundary types

test/test_examples_2d_p4est.jl

@@ -47,6 +47,12 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d")
 
   @testset "elixir_euler_nonperiodic_p4est.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_nonperiodic_p4est.jl"),
+      l2   = [2.3653424742684444e-6, 2.1388875095440695e-6, 2.1388875095548492e-6, 6.010896863397195e-6],
+      linf = [1.4080465931654018e-5, 1.7579850587257084e-5, 1.7579850592586155e-5, 5.956893531156027e-5])
+  end
+
+  @testset "elixir_euler_nonperiodic_p4est_unstructured.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_nonperiodic_p4est_unstructured.jl"),
       l2   = [0.005689496253354319, 0.004522481295923261, 0.004522481295922983, 0.009971628336802528],
       linf = [0.05125433503504517, 0.05343803272241532, 0.053438032722404216, 0.09032097668196482])
   end
@@ -58,6 +64,13 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d")
       atol = 5e-13, # required to make CI tests pass on macOS
     )
   end
+
+  @testset "elixir_eulergravity_eoc_test_p4est.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_eoc_test_p4est.jl"),
+      l2   = [0.0002456507718735638, 0.0003374236900407321, 0.00033742369004074064, 0.0007198398465056179],
+      linf = [0.0013653498212873316, 0.0018996803803537077, 0.0018996803803517093, 0.004384818684607161],
+      tspan = (0.0, 0.1))
+  end
 end
 
 end # module