Gauss-Legendre DGSEM Solver: 2D P4estMesh periodic

examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl

@@ -0,0 +1,49 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs,
+               basis_type = GaussLegendreBasis)
+
+# Deformed rectangle that looks like a waving flag,
+# lower and upper faces are sinus curves, left and right are vertical lines.
+f1(s) = SVector(-1.0, s - 1.0)
+f2(s) = SVector(1.0, s + 1.0)
+f3(s) = SVector(s, -1.0 + sinpi(0.5 * s))
+f4(s) = SVector(s, 1.0 + sinpi(0.5 * s))
+
+# Create P4estMesh with 3 x 2 trees and 6 x 4 elements,
+# approximate the geometry with a smaller polydeg for testing.
+trees_per_dimension = (3, 2)
+mesh = P4estMesh(trees_per_dimension, polydeg = 3,
+                 faces = (f1, f2, f3, f4),
+                 periodicity = true,
+                 initial_refinement_level = 1)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,
+                                    solver;
+                                    boundary_conditions = boundary_condition_periodic)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+ode = semidiscretize(semi, (0.0, 0.2))
+
+summary_callback = SummaryCallback()
+analysis_callback = AnalysisCallback(semi, interval = 100)
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback, analysis_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
+            ode_default_options()..., callback = callbacks);

src/meshes/p4est_mesh_view.jl

@@ -97,6 +97,10 @@ function extract_interfaces(mesh::P4estMeshView, interfaces_parent)
 
     # Copy relevant entries from parent mesh
     @views interfaces.u .= interfaces_parent.u[.., mask]
+    if interfaces.normal_directions !== nothing # Needed for Gauss-Legendre basis
+        @views interfaces.normal_directions .= interfaces_parent.normal_directions[..,
+                                                                                   mask]
+    end
     @views interfaces.node_indices .= interfaces_parent.node_indices[.., mask]
     @views neighbor_ids = interfaces_parent.neighbor_ids[.., mask]
 

src/solvers/dgsem_p4est/containers.jl

@@ -190,11 +190,8 @@ function Adapt.adapt_structure(to,
                                                size(elements.surface_flux_values))
 
     new_type_params = (NDIMS,
-                       RealT,
-                       uEltype,
-                       NDIMS + 1,
-                       NDIMS + 2,
-                       NDIMS + 3,
+                       RealT, uEltype,
+                       NDIMS + 1, NDIMS + 2, NDIMS + 3,
                        typeof(inverse_jacobian),     # ArrayRealTNDIMSP1
                        typeof(node_coordinates),     # ArrayRealTNDIMSP2
                        typeof(jacobian_matrix),      # ArrayRealTNDIMSP3
@@ -213,38 +210,48 @@ function Adapt.adapt_structure(to,
                                                      _surface_flux_values)
 end
 
-mutable struct P4estInterfaceContainer{NDIMS, uEltype <: Real, NDIMSP2,
+mutable struct P4estInterfaceContainer{NDIMS, RealT <: Real, uEltype <: Real,
+                                       NDIMSP1, NDIMSP2,
                                        uArray <: DenseArray{uEltype, NDIMSP2},
+                                       NormalArray <:
+                                       Union{DenseArray{RealT, NDIMSP1}, Nothing},
                                        IdsMatrix <: DenseMatrix{Int},
                                        IndicesMatrix <:
                                        DenseMatrix{NTuple{NDIMS, Symbol}},
                                        uVector <: DenseVector{uEltype},
+                                       NormalVector <:
+                                       Union{DenseVector{RealT}, Nothing},
                                        IdsVector <: DenseVector{Int},
                                        IndicesVector <:
                                        DenseVector{NTuple{NDIMS, Symbol}}} <:
                AbstractInterfaceContainer
-    u::uArray                   # [primary/secondary, variable, i, j, interface]
-    neighbor_ids::IdsMatrix     # [primary/secondary, interface]
-    node_indices::IndicesMatrix # [primary/secondary, interface]
+    u::uArray                      # [primary/secondary, variable, i, j, interface]
+    normal_directions::NormalArray # [dimension, i, j, interface]
+    neighbor_ids::IdsMatrix        # [primary/secondary, interface]
+    node_indices::IndicesMatrix    # [primary/secondary, interface]
 
     # internal `resize!`able storage
     _u::uVector
+    _normal_directions::NormalVector
     _neighbor_ids::IdsVector
     _node_indices::IndicesVector
 end
 
+# `trivial` means that the interface normals can be taken from
+# the outer/surface element nodes
+@inline trivial_interface_normals(::LobattoLegendreBasis) = true
+# For Gauss-Legendre basis, the interface normals need to be interpolated,
+# analogous to the surface values.
+@inline trivial_interface_normals(::GaussLegendreBasis) = false
+
 @inline function ninterfaces(interfaces::P4estInterfaceContainer)
     return size(interfaces.neighbor_ids, 2)
 end
 @inline Base.ndims(::P4estInterfaceContainer{NDIMS}) where {NDIMS} = NDIMS
-@inline function Base.eltype(::P4estInterfaceContainer{NDIMS, uEltype}) where {NDIMS,
-                                                                               uEltype}
-    return uEltype
-end
 
 # See explanation of Base.resize! for the element container
 function Base.resize!(interfaces::P4estInterfaceContainer, capacity)
-    @unpack _u, _neighbor_ids, _node_indices = interfaces
+    @unpack _u, _normal_directions, _neighbor_ids, _node_indices = interfaces
 
     n_dims = ndims(interfaces)
     n_nodes = size(interfaces.u, 3)
@@ -256,6 +263,18 @@ function Base.resize!(interfaces::P4estInterfaceContainer, capacity)
                                (2, n_variables, ntuple(_ -> n_nodes, n_dims - 1)...,
                                 capacity))
 
+    if _normal_directions === nothing # trivial interface normals, e.g. for Lobatto-Legendre basis
+        interfaces.normal_directions = nothing
+    else # non-trivial interface normals, e.g. for Gauss-Legendre basis
+        resize!(_normal_directions, n_dims * n_nodes^(n_dims - 1) * capacity)
+        interfaces.normal_directions = unsafe_wrap(ArrayType,
+                                                   pointer(_normal_directions),
+                                                   (n_dims,
+                                                    ntuple(_ -> n_nodes,
+                                                           n_dims - 1)...,
+                                                    capacity))
+    end
+
     resize!(_neighbor_ids, 2 * capacity)
     interfaces.neighbor_ids = unsafe_wrap(ArrayType, pointer(_neighbor_ids),
                                           (2, capacity))
@@ -271,6 +290,7 @@ end
 function init_interfaces(mesh::Union{P4estMesh, P4estMeshView, T8codeMesh}, equations,
                          basis, elements)
     NDIMS = ndims(elements)
+    RealT = real(mesh)
     uEltype = eltype(elements)
 
     # Initialize container
@@ -283,59 +303,94 @@ function init_interfaces(mesh::Union{P4estMesh, P4estMeshView, T8codeMesh}, equa
                     (2, nvariables(equations), ntuple(_ -> nnodes(basis), NDIMS - 1)...,
                      n_interfaces))
 
+    if !trivial_interface_normals(basis)
+        _normal_directions = Vector{RealT}(undef,
+                                           NDIMS * nnodes(basis)^(NDIMS - 1) *
+                                           n_interfaces)
+        normal_directions = unsafe_wrap(Array, pointer(_normal_directions),
+                                        (NDIMS,
+                                         ntuple(_ -> nnodes(basis), NDIMS - 1)...,
+                                         n_interfaces))
+    else
+        _normal_directions = nothing
+        normal_directions = nothing
+    end
+
     _neighbor_ids = Vector{Int}(undef, 2 * n_interfaces)
     neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids), (2, n_interfaces))
 
     _node_indices = Vector{NTuple{NDIMS, Symbol}}(undef, 2 * n_interfaces)
     node_indices = unsafe_wrap(Array, pointer(_node_indices), (2, n_interfaces))
 
-    interfaces = P4estInterfaceContainer{NDIMS, uEltype, NDIMS + 2,
-                                         typeof(u), typeof(neighbor_ids),
-                                         typeof(node_indices), typeof(_u),
+    interfaces = P4estInterfaceContainer{NDIMS, RealT, uEltype,
+                                         NDIMS + 1, NDIMS + 2,
+                                         typeof(u), typeof(normal_directions),
+                                         typeof(neighbor_ids), typeof(node_indices),
+                                         typeof(_u), typeof(_normal_directions),
                                          typeof(_neighbor_ids), typeof(_node_indices)}(u,
+                                                                                       normal_directions,
                                                                                        neighbor_ids,
                                                                                        node_indices,
                                                                                        _u,
+                                                                                       _normal_directions,
                                                                                        _neighbor_ids,
                                                                                        _node_indices)
 
-    init_interfaces!(interfaces, mesh)
+    init_interfaces!(interfaces, elements, mesh, basis)
 
     return interfaces
 end
 
-function init_interfaces!(interfaces, mesh::Union{P4estMesh, P4estMeshView})
+function init_interfaces!(interfaces, elements,
+                          mesh::Union{P4estMesh, P4estMeshView}, basis)
     init_surfaces!(interfaces, nothing, nothing, mesh)
+    init_normal_directions!(interfaces, basis, elements)
 
     return interfaces
 end
 
 function Adapt.parent_type(::Type{<:P4estInterfaceContainer{<:Any, <:Any, <:Any,
+                                                            <:Any, <:Any,
                                                             ArrayT}}) where {ArrayT}
     return ArrayT
 end
 
 # Manual adapt_structure since we have aliasing memory
-function Adapt.adapt_structure(to, interfaces::P4estInterfaceContainer)
+function Adapt.adapt_structure(to,
+                               interfaces::P4estInterfaceContainer{NDIMS,
+                                                                   RealT,
+                                                                   uEltype}) where {
+                                                                                    NDIMS,
+                                                                                    uEltype,
+                                                                                    RealT
+                                                                                    }
     # Adapt underlying storage
     _u = adapt(to, interfaces._u)
+    _normal_directions = interfaces._normal_directions === nothing ? nothing :
+                         adapt(to, interfaces._normal_directions)
     _neighbor_ids = adapt(to, interfaces._neighbor_ids)
     _node_indices = adapt(to, interfaces._node_indices)
     # Wrap arrays again
     u = unsafe_wrap_or_alloc(to, _u, size(interfaces.u))
+    normal_directions = _normal_directions === nothing ? nothing : # Check if interface normals are trivial
+                        unsafe_wrap_or_alloc(to, _normal_directions,
+                                             size(interfaces.normal_directions))
     neighbor_ids = unsafe_wrap_or_alloc(to, _neighbor_ids,
                                         size(interfaces.neighbor_ids))
     node_indices = unsafe_wrap_or_alloc(to, _node_indices,
                                         size(interfaces.node_indices))
 
-    NDIMS = ndims(interfaces)
     new_type_params = (NDIMS,
-                       eltype(_u),
-                       NDIMS + 2,
-                       typeof(u), typeof(neighbor_ids), typeof(node_indices),
-                       typeof(_u), typeof(_neighbor_ids), typeof(_node_indices))
-    return P4estInterfaceContainer{new_type_params...}(u, neighbor_ids, node_indices,
-                                                       _u, _neighbor_ids, _node_indices)
+                       RealT, eltype(_u),
+                       NDIMS + 1, NDIMS + 2,
+                       typeof(u), typeof(normal_directions),
+                       typeof(neighbor_ids), typeof(node_indices),
+                       typeof(_u), typeof(_normal_directions),
+                       typeof(_neighbor_ids), typeof(_node_indices))
+    return P4estInterfaceContainer{new_type_params...}(u, normal_directions,
+                                                       neighbor_ids, node_indices,
+                                                       _u, _normal_directions,
+                                                       _neighbor_ids, _node_indices)
 end
 
 mutable struct P4estBoundaryContainer{NDIMS, uEltype <: Real, NDIMSP1,
@@ -618,8 +673,7 @@ function Adapt.adapt_structure(to, mortars::P4estMortarContainer)
     NDIMS = ndims(mortars)
     new_type_params = (NDIMS,
                        eltype(_u),
-                       NDIMS + 1,
-                       NDIMS + 3,
+                       NDIMS + 1, NDIMS + 3,
                        typeof(u), typeof(neighbor_ids), typeof(node_indices),
                        typeof(_u), typeof(_neighbor_ids), typeof(_node_indices))
     return P4estMortarContainer{new_type_params...}(u, neighbor_ids, node_indices,
@@ -655,7 +709,12 @@ function reinitialize_containers!(mesh::P4estMesh, equations, dg::DGSEM, cache)
 
     # re-initialize containers together to reduce
     # the number of iterations over the mesh in `p4est`
-    return init_surfaces!(interfaces, mortars, boundaries, mesh)
+    init_surfaces!(interfaces, mortars, boundaries, mesh)
+
+    # init_normal_directions! requires that `node_indices` have been initialized
+    init_normal_directions!(interfaces, dg.basis, elements)
+
+    return nothing
 end
 
 # A helper struct used in initialization methods below