Implement non-periodic boundaries with P4estMesh (#598)

* Replace @timeit_debug with @timed * Extract orientation from face info * Reuse code from UnstructuredQuadMesh to implement non-periodic boundaries * Fix save/restart * Add nonperiodic example * Add assertion * Remove periodicity from P4estMesh * Fix loading meshes from different paths * Add constructor to build P4estMeshes from ABAQUS files * Allow initial_refinement_level > 1 * Add mapping as additional parameter to P4estMesh from file * Rename calc_node_coordinates! to avoid ambiguity * Add curved p4est Euler FSP example * Add documentation to new constructor * Check BCs for integrity * Add EOC test for non-periodic Euler on unstructured mesh * Implement suggestions * Implement suggestions * Fix tests * Implement suggestions * Fix spacing
trixi-framework · May 25, 2021 · d1af12b · d1af12b
1 parent ade432c
commit d1af12b
Show file tree

Hide file tree

Showing 13 changed files with 676 additions and 71 deletions.
diff --git a/.gitignore b/.gitignore
@@ -8,6 +8,7 @@
 *.avi
 *.ogv
 *.mesh
+*.inp
 **/Manifest.toml
 out*/
 docs/build

diff --git a/examples/2d/elixir_euler_free_stream_p4est.jl b/examples/2d/elixir_euler_free_stream_p4est.jl
@@ -0,0 +1,77 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+initial_condition = initial_condition_constant
+
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+
+# 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
+
+###############################################################################
+# Get the uncurved mesh from a file (downloads the file if not available locally)
+
+# Unstructured mesh with 48 cells of the square domain [-1, 1]^n
+mesh_file = joinpath(@__DIR__, "square_unstructured_1.inp")
+isfile(mesh_file) || download("https://gist.githubusercontent.com/efaulhaber/a075f8ec39a67fa9fad8f6f84342cbca/raw/a7206a02ed3a5d3cadacd8d9694ac154f9151db7/square_unstructured_1.inp",
+                              mesh_file)
+
+# Map the unstructured mesh with the mapping above
+mesh = P4estMesh(mesh_file, polydeg=3, mapping=mapping, initial_refinement_level=1)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions=Dict(
+                                      :all => BoundaryConditionDirichlet(initial_condition)
+                                    ))
+
+
+###############################################################################
+# 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_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl=2.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
diff --git a/examples/2d/elixir_euler_nonperiodic_p4est.jl b/examples/2d/elixir_euler_nonperiodic_p4est.jl
@@ -0,0 +1,70 @@
+
+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
diff --git a/src/mesh/mesh_io.jl b/src/mesh/mesh_io.jl
@@ -141,21 +141,22 @@ function save_mesh_file(mesh::P4estMesh, output_directory)
   mkpath(output_directory)
 
   filename = joinpath(output_directory, "mesh.h5")
-  p4est_filename = joinpath(output_directory, "p4est_data")
+  p4est_filename = "p4est_data"
+  p4est_file = joinpath(output_directory, p4est_filename)
 
   # Save the complete connectivity/p4est data to disk.
-  p4est_save(p4est_filename, mesh.p4est, false)
+  p4est_save(p4est_file, mesh.p4est, false)
 
   # Open file (clobber existing content)
   h5open(filename, "w") do file
     # Add context information as attributes
     attributes(file)["mesh_type"] = get_name(mesh)
     attributes(file)["ndims"] = ndims(mesh)
     attributes(file)["p4est_file"] = p4est_filename
-    attributes(file)["periodicity"] = collect(mesh.periodicity)
 
     file["tree_node_coordinates"] = mesh.tree_node_coordinates
     file["nodes"] = mesh.nodes
+    file["boundary_names"] = mesh.boundary_names .|> String
   end
 
   return filename
@@ -234,22 +235,27 @@ function load_mesh_serial(mesh_file::AbstractString; n_cells_max, RealT)
     end
     mesh = UnstructuredQuadMesh(mesh_filename; RealT=RealT, periodicity=periodicity_, unsaved_changes=false)
   elseif mesh_type == "P4estMesh"
-    p4est_file, periodicity_, tree_node_coordinates, nodes = h5open(mesh_file, "r") do file
+    p4est_filename, tree_node_coordinates,
+        nodes, boundary_names_ = h5open(mesh_file, "r") do file
       return read(attributes(file)["p4est_file"]),
-             read(attributes(file)["periodicity"]),
              read(file["tree_node_coordinates"]),
-             read(file["nodes"])
+             read(file["nodes"]),
+             read(file["boundary_names"])
     end
 
-    periodicity = Tuple(periodicity_)
+    boundary_names = boundary_names_ .|> Symbol
+
+    p4est_file = joinpath(dirname(mesh_file), p4est_filename)
+    # Prevent Julia crashes when p4est can't find the file
+    @assert isfile(p4est_file)
 
     conn_vec = Vector{Ptr{p4est_connectivity_t}}(undef, 1)
     p4est = p4est_load(p4est_file, C_NULL, 0, false, C_NULL, pointer(conn_vec))
     ghost = p4est_ghost_new(p4est, P4EST_CONNECT_FACE)
     p4est_mesh = p4est_mesh_new(p4est, ghost, P4EST_CONNECT_FACE)
 
     mesh = P4estMesh{ndims, RealT, ndims+2}(p4est, p4est_mesh, tree_node_coordinates,
-                                            nodes, periodicity, "", false)
+                                            nodes, boundary_names, "", false)
   else
     error("Unknown mesh type!")
   end