Flux differencing on simplicial elements

.github/workflows/ci.yml

@@ -63,8 +63,7 @@ jobs:
           - 2d_part1
           - 2d_part2
           - 2d_part3
-          - 2d_simplices
-          - 3d_simplices          
+          - dgmulti
           - 2d_mpi
           - 2d_threaded
           - 3d_part1

Project.toml

@@ -22,6 +22,7 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+StartUpDG = "472ebc20-7c99-4d4b-9470-8fde4e9faa0f"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StrideArrays = "d1fa6d79-ef01-42a6-86c9-f7c551f8593b"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
@@ -47,6 +48,7 @@ RecipesBase = "1.1"
 Reexport = "1.0"
 Requires = "1.1"
 StaticArrays = "1.0"
+StartUpDG = "0.10.2"
 StrideArrays = "0.1.11"
 StructArrays = "0.6"
 SummationByPartsOperators = "0.5.4"

docs/src/meshes/mesh_data_meshes.md

@@ -113,14 +113,20 @@ approximation on a triangle. There are also several parameters which can be twea
   [StartUpDG.jl docs](https://jlchan.github.io/StartUpDG.jl/dev/RefElemData/#RefElemData-based-on-SBP-finite-differences).
   Trixi will also specialize certain parts of the solver based on the `SBP` approximation type.
 
-## Trixi elixirs on simplicial meshes
+## Trixi elixirs on simplicial and tensor product element meshes
 
-Example elixirs on simplicial meshes can be found in the `examples/simplicial_mesh` folder.
-Some key elixirs to look at:
+Example elixirs with triangular, quadrilateral, and tetrahedral meshes can be found in
+the `examples/dgmulti_2d` and `examples/dgmulti_3d` folders. Some key elixirs to look at:
 
-* `elixir_euler_triangular_mesh.jl`: basic weak form DG discretization on a uniform triangular mesh.
-* `elixir_euler_periodic_triangular_mesh.jl`: same as above, but enforces periodicity in the ``x,y`` directions.
-* `elixir_euler_triangulate_pkg_mesh.jl`: uses a `TriangulateIO` unstructured mesh generated by
+* `examples/dgmulti_2d/elixir_euler_weakform.jl`: basic weak form DG discretization on a uniform triangular mesh.
+  Changing `element_type = Quad()` or `approximation_type = SBP()` will switch to a quadrilateral mesh
+  or an SBP-type discretization. Changing `surface_integral = SurfaceIntegralWeakForm(flux_ec)` and
+  `volume_integral = VolumeIntegralFluxDifferencing(flux_ec)` for some entropy conservative flux
+  (e.g., `flux_chandrashekar` or `flux_ranocha`) will switch to an entropy conservative formulation.
+* `examples/dgmulti_2d/elixir_euler_triangulate_pkg_mesh.jl`: uses a `TriangulateIO` unstructured mesh generated by
   [Triangulate.jl](https://github.com/JuliaGeometry/Triangulate.jl).
-* `elixir_ape_sbp_triangular_mesh.jl`: uses a multi-dimensional SBP discretization in weak form.
-* `elixir_euler_tet_mesh.jl`: basic weak form DG discretization on a uniform tet mesh.
+* `examples/dgmulti_3d/elixir_euler_weakform.jl`: basic weak form DG discretization on a uniform tetrahedral mesh.
+  Changing `element_type = Hex()` will switch to a hexahedral mesh. Changing
+  `surface_integral = SurfaceIntegralWeakForm(flux_ec)` and
+  `volume_integral = VolumeIntegralFluxDifferencing(flux_ec)` for some entropy conservative flux
+  (e.g., `flux_chandrashekar` or `flux_ranocha`) will switch to an entropy conservative formulation.

examples/simplicial_2d_dg/elixir_euler_triangulate_pkg_mesh.jl → examples/dgmulti_2d/elixir_euler_triangulate_pkg_mesh.jl

@@ -1,23 +1,21 @@
 # !!! warning "Experimental features"
 
 using Triangulate
-using StartUpDG, StructArrays
 using Trixi, OrdinaryDiffEq
 
-polydeg = 3
-rd = RefElemData(Tri(), polydeg)
-dg = DG(rd, nothing #= mortar =#, 
-        SurfaceIntegralWeakForm(FluxHLL()), VolumeIntegralWeakForm())
+dg = DGMulti(polydeg = 3, element_type = Tri(),
+             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             volume_integral = VolumeIntegralWeakForm())
 
 equations = CompressibleEulerEquations2D(1.4)
 initial_condition = initial_condition_convergence_test
 source_terms = source_terms_convergence_test
 
 meshIO = StartUpDG.square_hole_domain(.25) # pre-defined Triangulate geometry in StartUpDG
 
-# the pre-defined Triangulate geometry in StartUpDG has integer boundary tags. this routine 
+# the pre-defined Triangulate geometry in StartUpDG has integer boundary tags. this routine
 # assigns boundary faces based on these integer boundary tags.
-mesh = VertexMappedMesh(meshIO, rd, Dict(:bottom=>1, :right=>2, :top=>3, :left=>4))
+mesh = VertexMappedMesh(meshIO, dg, Dict(:bottom=>1, :right=>2, :top=>3, :left=>4))
 
 boundary_condition_convergence_test = BoundaryConditionDirichlet(initial_condition)
 boundary_conditions = (; :bottom => boundary_condition_convergence_test,
@@ -26,10 +24,10 @@ boundary_conditions = (; :bottom => boundary_condition_convergence_test,
                          :left => boundary_condition_convergence_test)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
-                                    source_terms = source_terms, 
-                                    boundary_conditions = boundary_conditions) 
+                                    source_terms = source_terms,
+                                    boundary_conditions = boundary_conditions)
 
-tspan = (0.0, 0.1)
+tspan = (0.0, 0.2)
 ode = semidiscretize(semi, tspan)
 
 summary_callback = SummaryCallback()
@@ -41,9 +39,6 @@ callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
 ###############################################################################
 # run the simulation
 
-dt0 = StartUpDG.estimate_h(rd,mesh.md) / StartUpDG.inverse_trace_constant(rd)
 sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
-            dt = 0.5*dt0, save_everystep=false, callback=callbacks);
+            dt = 0.5 * estimate_dt(mesh, dg), save_everystep=false, callback=callbacks);
 summary_callback() # print the timer summary
-
-l2,linf = analysis_callback(sol)

examples/simplicial_2d_dg/elixir_euler_quadrilateral_mesh.jl → examples/dgmulti_2d/elixir_euler_weakform.jl

@@ -1,12 +1,10 @@
 # !!! warning "Experimental features"
 
-using StartUpDG, StructArrays
 using Trixi, OrdinaryDiffEq
 
-polydeg = 3
-rd = RefElemData(Quad(), polydeg)
-dg = DG(rd, nothing #= mortar =#,
-        SurfaceIntegralWeakForm(FluxHLL()), VolumeIntegralWeakForm())
+dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(),
+             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             volume_integral = VolumeIntegralWeakForm())
 
 equations = CompressibleEulerEquations2D(1.4)
 initial_condition = initial_condition_convergence_test
@@ -16,8 +14,10 @@ source_terms = source_terms_convergence_test
 top_boundary(x,y,tol=50*eps()) = abs(y-1)<tol
 rest_of_boundary(x,y,tol=50*eps()) = !top_boundary(x,y,tol)
 is_on_boundary = Dict(:top => top_boundary, :rest => rest_of_boundary)
-vertex_coordinates_x, vertex_coordinates_y, EToV = StartUpDG.uniform_mesh(rd.elementType, 8)
-mesh = VertexMappedMesh(vertex_coordinates_x, vertex_coordinates_y, EToV, rd, is_on_boundary = is_on_boundary)
+
+cells_per_dimension = (8,8) # detected by `extract_initial_resolution` for convergence tests
+vertex_coordinates_x, vertex_coordinates_y, EToV = StartUpDG.uniform_mesh(dg.basis.elementType, cells_per_dimension...)
+mesh = VertexMappedMesh(vertex_coordinates_x, vertex_coordinates_y, EToV, dg, is_on_boundary = is_on_boundary)
 
 boundary_condition_convergence_test = BoundaryConditionDirichlet(initial_condition)
 boundary_conditions = (; :top => boundary_condition_convergence_test,
@@ -27,7 +27,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
                                     source_terms = source_terms,
                                     boundary_conditions = boundary_conditions)
 
-tspan = (0.0, 0.1)
+tspan = (0.0, 0.4)
 ode = semidiscretize(semi, tspan)
 
 summary_callback = SummaryCallback()
@@ -39,9 +39,6 @@ callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
 ###############################################################################
 # run the simulation
 
-dt0 = StartUpDG.estimate_h(rd,mesh.md) / StartUpDG.inverse_trace_constant(rd)
 sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
-            dt = 0.5*dt0, save_everystep=false, callback=callbacks);
+            dt = 0.5 * estimate_dt(mesh, dg), save_everystep=false, callback=callbacks);
 summary_callback() # print the timer summary
-
-l2,linf = analysis_callback(sol)

examples/simplicial_2d_dg/elixir_euler_periodic_triangular_mesh.jl → examples/dgmulti_2d/elixir_euler_weakform_periodic.jl

@@ -1,23 +1,21 @@
 # !!! warning "Experimental features"
 
-using StartUpDG, StructArrays
 using Trixi, OrdinaryDiffEq
 
-polydeg = 3
-rd = RefElemData(Tri(), polydeg) # equivalent to a "basis"
-dg = DG(rd, nothing #= mortar =#,
-        SurfaceIntegralWeakForm(FluxHLL()), VolumeIntegralWeakForm())
+dg = DGMulti(polydeg = 3, element_type = Tri(),
+             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             volume_integral = VolumeIntegralWeakForm())
 
 equations = CompressibleEulerEquations2D(1.4)
 initial_condition = initial_condition_convergence_test
 source_terms = source_terms_convergence_test
 
-vertex_coordinates_x, vertex_coordinates_y, EToV = StartUpDG.uniform_mesh(rd.elementType, 8)
-mesh = VertexMappedMesh(vertex_coordinates_x, vertex_coordinates_y, EToV, rd, is_periodic=(true,true))
+vertex_coordinates_x, vertex_coordinates_y, EToV = StartUpDG.uniform_mesh(dg.basis.elementType, 4)
+mesh = VertexMappedMesh(vertex_coordinates_x, vertex_coordinates_y, EToV, dg, is_periodic=(true,true))
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
                                     source_terms = source_terms)
 
-tspan = (0.0, 0.1)
+tspan = (0.0, 0.4)
 ode = semidiscretize(semi, tspan)
 
 summary_callback = SummaryCallback()
@@ -29,9 +27,6 @@ callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
 ###############################################################################
 # run the simulation
 
-dt0 = StartUpDG.estimate_h(rd,mesh.md) / StartUpDG.inverse_trace_constant(rd)
 sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
-            dt = 0.5*dt0, save_everystep=false, callback=callbacks);
+            dt = 0.5 * estimate_dt(mesh, dg), save_everystep=false, callback=callbacks);
 summary_callback() # print the timer summary
-
-l2,linf = analysis_callback(sol)

examples/simplicial_3d_dg/elixir_euler_tet_mesh.jl → examples/dgmulti_3d/elixir_euler_weakform.jl

@@ -1,31 +1,30 @@
 # !!! warning "Experimental features"
 
-using StartUpDG, StructArrays
 using Trixi, OrdinaryDiffEq
 
-polydeg = 3
-rd = RefElemData(Tet(), polydeg)
-dg = DG(rd, nothing #= mortar =#, 
-        SurfaceIntegralWeakForm(FluxHLL()), VolumeIntegralWeakForm())
+dg = DGMulti(polydeg = 3, element_type = Tet(),
+             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             volume_integral = VolumeIntegralWeakForm())
 
 equations = CompressibleEulerEquations3D(1.4)
 initial_condition = initial_condition_convergence_test
 source_terms = source_terms_convergence_test
 
 # example where we tag two separate boundary segments of the mesh
-top_boundary(x, y, z, tol=50*eps()) = abs(y - 1) < tol 
+top_boundary(x, y, z, tol=50*eps()) = abs(y - 1) < tol
 rest_of_boundary(x, y, z, tol=50*eps()) = !top_boundary(x, y, z, tol)
 is_on_boundary = Dict(:top => top_boundary, :rest => rest_of_boundary)
-VX, VY, VZ, EToV = StartUpDG.uniform_mesh(Tet(), 4)
-mesh = VertexMappedMesh(VX, VY, VZ, EToV, rd, is_on_boundary = is_on_boundary)
+vertex_coordinates_x, vertex_coordinates_y, vertex_coordinates_z, EToV = StartUpDG.uniform_mesh(dg.basis.elementType, 4)
+mesh = VertexMappedMesh((vertex_coordinates_x, vertex_coordinates_y, vertex_coordinates_z), EToV,
+                         dg, is_on_boundary = is_on_boundary)
 
 boundary_condition_convergence_test = BoundaryConditionDirichlet(initial_condition)
 boundary_conditions = (; :top => boundary_condition_convergence_test,
                         :rest => boundary_condition_convergence_test)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
-                                    source_terms = source_terms, 
-                                    boundary_conditions = boundary_conditions) 
+                                    source_terms = source_terms,
+                                    boundary_conditions = boundary_conditions)
 
 tspan = (0.0, 0.1)
 ode = semidiscretize(semi, tspan)
@@ -39,9 +38,6 @@ callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
 ###############################################################################
 # run the simulation
 
-dt0 = StartUpDG.estimate_h(rd,mesh.md) / StartUpDG.inverse_trace_constant(rd)
 sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
-            dt = 0.5*dt0, save_everystep=false, callback=callbacks);
+            dt = 0.5 * estimate_dt(mesh, dg), save_everystep=false, callback=callbacks);
 summary_callback() # print the timer summary
-
-l2,linf = analysis_callback(sol)

examples/dgmulti_3d/elixir_euler_weakform_periodic.jl

@@ -0,0 +1,34 @@
+# !!! warning "Experimental features"
+
+using Trixi, OrdinaryDiffEq
+
+dg = DGMulti(polydeg = 3, element_type = Tet(),
+             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             volume_integral = VolumeIntegralWeakForm())
+
+equations = CompressibleEulerEquations3D(1.4)
+initial_condition = initial_condition_convergence_test
+source_terms = source_terms_convergence_test
+
+vertex_coordinates_x, vertex_coordinates_y, vertex_coordinates_z, EToV = StartUpDG.uniform_mesh(dg.basis.elementType, 4)
+mesh = VertexMappedMesh((vertex_coordinates_x, vertex_coordinates_y, vertex_coordinates_z), EToV, dg,
+                        is_periodic = (true, true, true))
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
+                                    source_terms = source_terms)
+
+tspan = (0.0, 0.1)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval=10)
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
+callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt = 0.5 * estimate_dt(mesh, dg), save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary