Weak form modal DG methods on simplicial (and quad/hex meshes)

.github/workflows/ci.yml

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

.gitignore

@@ -21,3 +21,5 @@ coverage_report/
 **/solution_*.txt
 **/restart_*.txt
 .vscode/
+
+.DS_Store

NEWS.md

@@ -22,6 +22,7 @@ for human readability.
 - New unstructured, curvilinear, adaptive (non-conforming) mesh type `P4estMesh` in 2D and 3D (experimental)
 - Experimental support for finite difference (FD) summation-by-parts (SBP) methods via
   [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl)
+- New support for modal DG and SBP-DG methods on triangular and tetrahedral meshes via [StartUpDG.jl](https://github.com/jlchan/StartUpDG.jl) (experimental). 
 
 #### Changed
 

Project.toml

@@ -13,6 +13,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
+Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 P4est = "7d669430-f675-4ae7-b43e-fab78ec5a902"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
@@ -36,6 +37,7 @@ LinearMaps = "2.7, 3.0"
 LoopVectorization = "0.12.35"
 MPI = "0.16, 0.17, 0.18"
 OffsetArrays = "1.3"
+Octavian = "0.2.20"
 P4est = "0.2.1"
 Polyester = "0.3"
 RecipesBase = "1.1"

docs/make.jl

@@ -54,6 +54,7 @@ makedocs(
                 "Tree mesh" => joinpath("meshes", "tree_mesh.md"),
                 "Structured mesh" => joinpath("meshes", "structured_mesh.md"),
                 "Unstructured mesh" => joinpath("meshes", "unstructured_quad_mesh.md"),
+                "Simplicial mesh" => joinpath("meshes", "mesh_data_meshes.md"),
             ],
             "Time integration" => "time_integration.md",
             "Callbacks" => "callbacks.md",

docs/src/meshes/mesh_data_meshes.md

@@ -0,0 +1,123 @@
+# Unstructured simplicial meshes
+
+Support for simplicial meshes uses the [StartUpDG.jl](https://github.com/jlchan/StartUpDG.jl) 
+package. These meshes are constructed by specifying a list of vertex coordinates `VX`, `VY`, `VZ` 
+and a connectivity matrix `EToV` where `EToV[e,:]` gives the vertices which correspond to element `e`. 
+
+We make a few simplifying assumptions about simplicial meshes:
+* meshes consist of a single type of element
+* meshes are _conforming_ (e.g., each face of an element is shared with at most one other element). 
+* the geometric mapping from reference to physical elements is polynomial (currently, only affine 
+mappings are supported).
+
+## `AbstractMeshData` wrapper type
+
+Simplicial meshes in Trixi are represented using an `AbstractMeshData{Dim, ElemType}` wrapper. 
+For example, `VertexMappedMesh{Dim, Tri} <: AbstractMeshData{Dim, Tri}` describes a mesh whose
+reference-to-physical mapping can be constructed using only the vertex positions (e.g., an affine 
+simplicial mesh or a bilinear quadrilateral mesh). These quantities are generated by most simplicial 
+mesh generators, and `StartUpDG.jl` includes both simple uniform meshes via `uniform_mesh`, as 
+well as support for triangular meshes constructed using 
+[Triangulate.jl](https://github.com/JuliaGeometry/Triangulate.jl), a wrapper around Jonathan Shewchuk's 
+[Triangle](https://www.cs.cmu.edu/~quake/triangle.html) package.
+
+The only fields are `md::MeshData`, which contains geometric terms derived from the mapping between 
+the reference and physical elements, and `boundary_faces`, which contains a `Dict` of boundary 
+segment names (symbols) and list of faces which lie on that boundary segment. 
+
+## Variable naming conventions
+
+We use the convention that coordinates on the reference element are ``r`` in 1D, ``r, s`` in 2D, 
+or ``r, s, t`` in 3D. Physical coordinates use the standard conventions ``x``, ``x, y``, and 
+``x, y, z`` in 1D, 2D, and 3D. 
+
+!["Ref-to-physical mapping"](https://jlchan.github.io/StartUpDG.jl/dev/assets/mapping_diagram.png)
+
+Derivatives of reference coordinates with respect to physical coordinates are abbreviated, e.g., 
+``\frac{\partial r}{\partial x} = r_x``. Additionally, ``J`` is used to denote the determinant of 
+the Jacobian of the reference-to-physical mapping. 
+
+## Variable meanings and conventions in `StartUpDG.jl`
+
+`StartUpDG.jl` exports structs `RefElemData{Dim, ElemShape, ...}` (which contains data associated 
+with the reference element, such as interpolation points, quadrature rules, face nodes, normals, 
+and differentiation/interpolation/projection matrices) and `MeshData{Dim}` (which contains geometric 
+data associated with a mesh). These are currently used for evaluating DG formulations in a matrix-free 
+fashion. These structs contain fields similar (but not identical) to those in 
+`Globals1D, Globals2D, Globals3D` in the Matlab codes from "Nodal Discontinuous Galerkin Methods" 
+by Hesthaven and Warburton (2007). 
+
+In general, we use the following code conventions:
+* variables such as `r, s,...` and `x, y,...` correspond to values at nodal interpolation points. 
+* variables ending in `q` (e.g., `rq, sq,...` and `xq, yq,...`) correspond to values at volume 
+quadrature points. 
+* variables ending in `f` (e.g., `rf, sf,...` and `xf, yf,...`) correspond to values at face 
+quadrature points. 
+* variables ending in `p` (e.g., `rp, sp,...`) correspond to "plotting" points, which are usually 
+a fine grid of equispaced points.
+* `V` matrices correspond to interpolation matrices from nodal interpolation points, e.g., `Vq` 
+interpolates to volume quadrature points, `Vf` interpolates to face quadrature points. 
+* geometric quantities in `MeshData` are stored as matrices of dimension 
+``\text{number of points per element } \times \text{number of elements}``. 
+
+Quantities in `rd::RefElemData`: 
+* `rd.Np, rd.Nq, rd.Nf`: the number of nodal interpolation points, volume quadrature points, and 
+face quadrature points on the reference element, respectively. 
+* `rd.Vq`: interpolation matrices from values at nodal interpolation points to volume quadrature points
+* `rd.wq`: volume quadrature weights on the reference element
+* `rd.Vf`: interpolation matrices from values at nodal interpolation points to face quadrature points
+* `rd.wf`: a vector containing face quadrature weights on the reference element
+* `rd.M`: the quadrature-based mass matrix, computed via `rd.Vq' * diagm(rd.wq) * rd.Vq`.
+* `rd.Pq`: a quadrature-based ``L^2`` projection matrix `rd.Pq = rd.M \ rd.Vq' * diagm(rd.wq)` 
+which maps between values at quadrature points and values at nodal points. 
+* `Dr, Ds, Dt` matrices are nodal differentiation matrices with respect to the ``r,s,t`` coordinates, 
+e.g., `Dr*f.(r,s)` approximates the derivative of ``f(r,s)`` at nodal points. 
+
+Quantities in `md::MeshData`: 
+* `md.xyz` is a tuple of matrices `md.x`, `md.y`, `md.z`, where column `e` contains coordinates of 
+physical interpolation points. 
+* `md.xyzq` is a tuple of matrices `md.xq`, `md.yq`, `md.zq`, where column `e` contains coordinates 
+of physical quadrature points. 
+* `md.rxJ, md.sxJ, ...` are matrices where column `e` contains values of 
+``J\frac{\partial r}{\partial x}``, ``J\frac{\partial s}{\partial x}``, etc. at nodal interpolation 
+points on the element `e`.
+* `md.J` is a matrix where column `e` contains values of the Jacobian ``J`` at nodal interpolation points.
+* `md.Jf` is a matrix where column `e` contains values of the face Jacobian (e.g., determinant of 
+the geometric mapping between a physical face and a reference face) at face quadrature points.
+* `md.nxJ, md.nyJ, ...` are matrices where column `e` contains values of components of the unit 
+normal scaled by the face Jacobian `md.Jf` at face quadrature points.
+
+For more details, please see the [StartUpDG.jl docs](https://jlchan.github.io/StartUpDG.jl/dev/). 
+## Special options
+
+Trixi solvers on simplicial meshes use `DG` solver types with the basis field set as a `RefElemData` 
+type. By default, `RefElemData(Tri(), polydeg)` will generate data for a degree `polydeg` polynomial 
+approximation on a triangle. There are also several parameters which can be tweaked:
+
+* `RefElemData(Tri(), polydeg, quad_rule_vol = quad_nodes(Tri(), Nq))` will substitute in a volume 
+quadrature rule of degree `Nq` instead of the default (which is a quadrature rule of degree `polydeg`).
+Here, a degree `Nq` rule will be exact for at least degree `2*Nq` integrands (such that the mass 
+matrix is integrated exactly). Quadrature rules of which exactly integrate degree `Nq` integrands 
+may also be specified by setting `quad_rule_vol = StartUpDG.quad_nodes_tri(Nq)`. 
+* `RefElemData(Tri(), polydeg, quad_rule_face = quad_nodes(Line(), Nq))` will use a face quadrature rule 
+of degree `Nq` rather than the default. This rule is also exact for at least degree `2*Nq` integrands. 
+* `RefElemData(Tri(), SBP(), polydeg)` will generate a reference element corresponding to a 
+multi-dimensional SBP discretization. Types of SBP discretizations available include 
+`SBP{Kubatko{LobattoFaceNodes}}()` (the default choice), `SBP{Kubatko{LegendreFaceNodes}}()`, and 
+[`SBP{Hicken}()`](https://doi.org/10.1007/s10915-020-01154-8). For `polydeg = 1, ..., 4`, the 
+`SBP{Kubatko{LegendreFaceNodes}}()` SBP nodes are identical to the SBP nodes of 
+[Chen and Shu](https://doi.org/10.1016/j.jcp.2017.05.025). 
+More detailed descriptions of each SBP node set can be found in the 
+[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
+
+Example elixirs on simplicial meshes can be found in the `examples/simplicial_mesh` folder. 
+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 `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/simplicial_2d_dg/elixir_ape_sbp_triangular_mesh.jl

@@ -0,0 +1,47 @@
+# !!! warning "Experimental features"
+
+using StartUpDG, StructArrays
+using Trixi, OrdinaryDiffEq
+
+polydeg = 3
+rd = RefElemData(Tri(), SBP(), polydeg)
+dg = DG(rd, nothing #= mortar =#, 
+        SurfaceIntegralWeakForm(FluxLaxFriedrichs()), VolumeIntegralWeakForm())
+
+v_mean_global = (0.25, 0.25)
+c_mean_global = 1.0
+rho_mean_global = 1.0
+equations = AcousticPerturbationEquations2D(v_mean_global, c_mean_global, rho_mean_global)
+
+initial_condition = initial_condition_convergence_test
+source_terms = source_terms_convergence_test
+
+# example where we tag two separate boundary segments of the mesh
+VX, VY, EToV = StartUpDG.uniform_mesh(Tri(), 8)
+mesh = VertexMappedMesh(VX, VY, EToV, rd)
+
+boundary_condition_convergence_test = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = (; :entire_boundary => boundary_condition_convergence_test)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
+                                    source_terms = source_terms, 
+                                    boundary_conditions = boundary_conditions) 
+
+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
+
+dt0 = StartUpDG.estimate_h(rd,mesh.md) / StartUpDG.inverse_trace_constant(rd)
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt = .5*dt0, 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

@@ -0,0 +1,37 @@
+# !!! 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())
+
+equations = CompressibleEulerEquations2D(1.4)
+initial_condition = initial_condition_convergence_test
+source_terms = source_terms_convergence_test
+
+VX, VY, EToV = StartUpDG.uniform_mesh(rd.elementType, 8)
+mesh = VertexMappedMesh(VX, VY, EToV, rd, is_periodic=(true,true))
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
+                                    source_terms = source_terms) 
+
+tspan = (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
+
+dt0 = StartUpDG.estimate_h(rd,mesh.md) / StartUpDG.inverse_trace_constant(rd)
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt = .5*dt0, 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

@@ -0,0 +1,47 @@
+# !!! warning "Experimental features"
+
+using StartUpDG, StructArrays
+using Trixi, OrdinaryDiffEq
+
+polydeg = 3
+rd = RefElemData(Quad(), polydeg)
+dg = DG(rd, nothing #= mortar =#, 
+        SurfaceIntegralWeakForm(FluxHLL()), VolumeIntegralWeakForm())
+
+equations = CompressibleEulerEquations2D(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,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)
+VX, VY, EToV = StartUpDG.uniform_mesh(rd.elementType, 8)
+mesh = VertexMappedMesh(VX, VY, EToV, rd, 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) 
+
+tspan = (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
+
+dt0 = StartUpDG.estimate_h(rd,mesh.md) / StartUpDG.inverse_trace_constant(rd)
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt = .5*dt0, save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary
+
+l2,linf = analysis_callback(sol)

examples/simplicial_2d_dg/elixir_euler_triangular_mesh.jl

@@ -0,0 +1,47 @@
+# !!! warning "Experimental features"
+
+using StartUpDG, StructArrays
+using Trixi, OrdinaryDiffEq
+
+polydeg = 3
+rd = RefElemData(Tri(), polydeg)
+dg = DG(rd, nothing #= mortar =#, 
+        SurfaceIntegralWeakForm(FluxHLL()), VolumeIntegralWeakForm())
+
+equations = CompressibleEulerEquations2D(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,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)
+VX, VY, EToV = StartUpDG.uniform_mesh(Tri(), 8)
+mesh = VertexMappedMesh(VX, VY, EToV, rd, 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) 
+
+tspan = (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
+
+dt0 = StartUpDG.estimate_h(rd,mesh.md) / StartUpDG.inverse_trace_constant(rd)
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt = .5*dt0, save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary
+
+l2,linf = analysis_callback(sol)