Implement 3D support for P4estMesh (#637)

* Implement structured 3D simulations with P4estMesh * Fix loading and saving of P4estMesh in 3D * Fix unstructured P4estMesh simulation in 3D * Revise 3D orientations * Implement P4estMesh simulation for non-conforming 3D meshes * Add non-conforming example * Clean up * Remove unused destroy functions * Clean up more * Add restart example with P4estMesh in 3D * Implement AMR with P4estMesh in 3D * Add AMR example with P4estMesh in 3D * First version of fast calc_contravariant_vectors! * Make calc_contravariant_vectors! more readable * Improve documentation of calc_contravariant_vectors! * Calculate both summands of contravariant vectors in one loop * Move @turbo to inner loop * Only use one loop over n * Optimize calc_jacobian_matrix! * Optimize calc_inverse_jacobian! * Use tmp arrays for multiply_dimensionwise! * Update Euler FSP errors after optimizing code * Remove duplicate code * Dispatch by function type * Compute contravariant vectors in one loop for each summand * Add more 3D `P4estMesh` examples * Make curved example actually curved * Make calc_node_coordinates! consistent with 2D * Remove sign Jacobian from normal vector p4est doesn't support left-handed coordinates anyway. * Add FSP example with P4estMesh * Add P4estMesh examples to tests * Clean up * Update examples/3d/elixir_advection_amr_p4est.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update examples/3d/elixir_advection_basic_p4est.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/solvers/dg_p4est/containers_2d.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Only use Base functions for sc_array load and wrap * Implement suggestions * Implement suggestions * Add convergence test for 3D `P4estMesh` * Implement suggestions Co-authored-by: Hendrik Ranocha <[email protected]>
trixi-framework · Jun 24, 2021 · e6aa4bd · e6aa4bd
1 parent b2d0d98
commit e6aa4bd
Show file tree

Hide file tree

Showing 38 changed files with 2,620 additions and 708 deletions.
diff --git a/examples/2d/elixir_advection_amr_p4est_unstructured_flag.jl b/examples/2d/elixir_advection_amr_p4est_unstructured_flag.jl
@@ -37,9 +37,9 @@ 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)
+mesh = P4estMesh{2}(mesh_file, polydeg=3,
+                    mapping=mapping_flag,
+                    initial_refinement_level=1)
 
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,

diff --git a/examples/2d/elixir_advection_p4est_non_conforming_flag.jl b/examples/2d/elixir_advection_p4est_non_conforming_flag.jl
@@ -40,7 +40,7 @@ end
 # Refine recursively until each bottom left quadrant of a tree 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)
+Trixi.refine_p4est!(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)
@@ -49,7 +49,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergen
 ###############################################################################
 # ODE solvers, callbacks etc.
 
-# Create ODE problem with time span from 0.0 to 1.0
+# Create ODE problem with time span from 0.0 to 0.2
 ode = semidiscretize(semi, (0.0, 0.2));
 
 # At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup

diff --git a/examples/2d/elixir_advection_p4est_non_conforming_flag_unstructured.jl b/examples/2d/elixir_advection_p4est_non_conforming_flag_unstructured.jl
@@ -36,9 +36,9 @@ 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=0)
+mesh = P4estMesh{2}(mesh_file, polydeg=3,
+                    mapping=mapping_flag,
+                    initial_refinement_level=0)
 
 # Refine bottom left quadrant of each tree to level 3
 function refine_fn(p4est, which_tree, quadrant)
@@ -54,7 +54,7 @@ end
 # Refine recursively until each bottom left quadrant of a tree 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)
+Trixi.refine_p4est!(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, solver, boundary_conditions=boundary_conditions)
@@ -63,7 +63,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
 ###############################################################################
 # ODE solvers, callbacks etc.
 
-# Create ODE problem with time span from 0.0 to 1.0
+# Create ODE problem with time span from 0.0 to 0.2
 ode = semidiscretize(semi, (0.0, 0.2));
 
 # At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup

diff --git a/examples/2d/elixir_euler_free_stream_p4est.jl b/examples/2d/elixir_euler_free_stream_p4est.jl
@@ -36,7 +36,7 @@ isfile(mesh_file) || download("https://gist.githubusercontent.com/efaulhaber/a07
                               mesh_file)
 
 # Map the unstructured mesh with the mapping above
-mesh = P4estMesh(mesh_file, polydeg=3, mapping=mapping, initial_refinement_level=1)
+mesh = P4estMesh{2}(mesh_file, polydeg=3, mapping=mapping, initial_refinement_level=1)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                     boundary_conditions=Dict(

diff --git a/examples/2d/elixir_euler_nonperiodic_p4est_unstructured.jl b/examples/2d/elixir_euler_nonperiodic_p4est_unstructured.jl
@@ -28,7 +28,7 @@ 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)
+mesh = P4estMesh{2}(mesh_file, initial_refinement_level=0)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                     source_terms=source_terms,

diff --git a/examples/3d/elixir_advection_amr_p4est.jl b/examples/3d/elixir_advection_amr_p4est.jl
@@ -0,0 +1,72 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advectionvelocity = (0.2, -0.7, 0.5)
+equations = LinearScalarAdvectionEquation3D(advectionvelocity)
+
+initial_condition = initial_condition_gauss
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+
+coordinates_min = (-5, -5, -5)
+coordinates_max = ( 5,  5,  5)
+trees_per_dimension = (1, 1, 1)
+mesh = P4estMesh(trees_per_dimension, polydeg=1,
+                 coordinates_min=coordinates_min, coordinates_max=coordinates_max,
+                 initial_refinement_level=4)
+
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.3)
+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=4,
+                                      med_level=5, med_threshold=0.1,
+                                      max_level=6, 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=1.2)
+
+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.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
diff --git a/examples/3d/elixir_advection_amr_p4est_unstructured_curved.jl b/examples/3d/elixir_advection_amr_p4est_unstructured_curved.jl
@@ -0,0 +1,105 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advectionvelocity = (0.2, -0.7, 0.5)
+equations = LinearScalarAdvectionEquation3D(advectionvelocity)
+
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+
+initial_condition = initial_condition_gauss
+boundary_condition = BoundaryConditionDirichlet(initial_condition)
+
+boundary_conditions = Dict(
+  :all => boundary_condition
+)
+
+# Mapping as described in https://arxiv.org/abs/2012.12040, but with less warping.
+# The original mapping applied to this unstructured mesh creates extreme angles,
+# which require a high resolution for proper results.
+function mapping(xi, eta, zeta)
+  # Don't transform input variables between -1 and 1 onto [0,3] to obtain curved boundaries
+  # xi = 1.5 * xi_ + 1.5
+  # eta = 1.5 * eta_ + 1.5
+  # zeta = 1.5 * zeta_ + 1.5
+
+  y = eta + 1/4 * (cos(1.5 * pi * (2 * xi - 3)/3) *
+                   cos(0.5 * pi * (2 * eta - 3)/3) *
+                   cos(0.5 * pi * (2 * zeta - 3)/3))
+
+  x = xi + 1/4 * (cos(0.5 * pi * (2 * xi - 3)/3) *
+                  cos(2 * pi * (2 * y - 3)/3) *
+                  cos(0.5 * pi * (2 * zeta - 3)/3))
+
+  z = zeta + 1/4 * (cos(0.5 * pi * (2 * x - 3)/3) *
+                    cos(pi * (2 * y - 3)/3) *
+                    cos(0.5 * pi * (2 * zeta - 3)/3))
+
+  # Transform the weird deformed cube to be approximately the size of [-5,5]^3 to match IC
+  return SVector(5 * x, 5 * y, 5 * z)
+end
+
+# Unstructured mesh with 48 cells of the cube domain [-1, 1]^3
+mesh_file = joinpath(@__DIR__, "cube_unstructured_2.inp")
+isfile(mesh_file) || download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp",
+                              mesh_file)
+
+mesh = P4estMesh{3}(mesh_file, polydeg=3,
+                    mapping=mapping,
+                    initial_refinement_level=1)
+
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions=boundary_conditions)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 8.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=1.2)
+
+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.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
diff --git a/examples/3d/elixir_advection_basic_p4est.jl b/examples/3d/elixir_advection_basic_p4est.jl
@@ -0,0 +1,63 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advectionvelocity = (1.0, 1.0, 1.0)
+equations = LinearScalarAdvectionEquation3D(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.5, -0.9, 0.0) # minimum coordinates (min(x), min(y), min(z))
+coordinates_max = ( 0.5,  1.1, 4.0) # maximum coordinates (max(x), max(y), max(z))
+
+# Create P4estMesh with 8 x 10 x 16 elements (note `refinement_level=1`)
+trees_per_dimension = (4, 5, 8)
+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 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,
+                                     solution_variables=cons2prim)
+
+# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl=1.2)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_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);
+
+# Print the timer summary
+summary_callback()
diff --git a/examples/3d/elixir_advection_p4est_non_conforming.jl b/examples/3d/elixir_advection_p4est_non_conforming.jl
@@ -0,0 +1,76 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advectionvelocity = (1.0, 1.0, 1.0)
+equations = LinearScalarAdvectionEquation3D(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, -1.0) # minimum coordinates (min(x), min(y), min(z))
+coordinates_max = ( 1.0,  1.0,  1.0) # maximum coordinates (max(x), max(y), max(z))
+
+trees_per_dimension = (1, 1, 1)
+mesh = P4estMesh(trees_per_dimension, polydeg=3,
+                 coordinates_min=coordinates_min, coordinates_max=coordinates_max,
+                 initial_refinement_level=2)
+
+# Refine bottom left quadrant of each tree to level 4
+function refine_fn(p8est, which_tree, quadrant)
+  if quadrant.x == 0 && quadrant.y == 0 && quadrant.z == 0 && quadrant.level < 3
+    # 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 tree has level 2
+# The mesh will be rebalanced before the simulation starts
+refine_fn_c = @cfunction(refine_fn, Cint, (Ptr{Trixi.p8est_t}, Ptr{Trixi.p4est_topidx_t}, Ptr{Trixi.p8est_quadrant_t}))
+Trixi.refine_p4est!(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, 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()