First step towards FD methods and surface integrals/terms

NEWS.md

@@ -20,6 +20,8 @@ for human readability.
 - New structured, curvilinear, conforming mesh type `CurvedMesh` (experimental)
 - New unstructured, curvilinear, conforming mesh type `UnstructuredQuadMesh` in 2D (experimental)
 - New unstructured, curvilinear, adaptive (non-conforming) mesh type `P4estMesh` in 2D (experimental)
+- Experimental support for finite difference (FD) summation-by-parts (SBP) methods via
+  [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl)
 
 #### Changed
 

Project.toml

@@ -22,6 +22,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StrideArrays = "d1fa6d79-ef01-42a6-86c9-f7c551f8593b"
+SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
@@ -42,6 +43,7 @@ Reexport = "1.0"
 Requires = "1.1"
 StaticArrays = "1.0"
 StrideArrays = "0.1.8"
+SummationByPartsOperators = "0.5.4"
 TimerOutputs = "0.5"
 UnPack = "1.0"
 julia = "1.6"

examples/tree_2d_fdsbp/elixir_advection_extended.jl

@@ -0,0 +1,57 @@
+# TODO: FD
+# !!! warning "Experimental feature"
+#     This is an experimental feature and may change in any future releases.
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advectionvelocity = (0.2, -0.3)
+equations = LinearScalarAdvectionEquation2D(advectionvelocity)
+
+initial_condition = initial_condition_convergence_test
+
+D_SBP = derivative_operator(SummationByPartsOperators.MattssonNordström2004(),
+                            derivative_order=1, accuracy_order=4,
+                            xmin=0.0, xmax=1.0, N=100)
+solver = DG(D_SBP, nothing #= mortar =#,
+            SurfaceIntegralStrongForm(flux_lax_friedrichs),
+            VolumeIntegralStrongForm())
+
+coordinates_min = (-1.0, -1.0)
+coordinates_max = ( 1.0,  1.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=1,
+                n_cells_max=30_000,
+                periodicity=true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+
+###############################################################################
+# 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,
+                                     extra_analysis_integrals=(energy_total,))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback)
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, RDPK3SpFSAL49(), abstol=1.0e-8, reltol=1.0e-8,
+            save_everystep=false, callback=callbacks, maxiters=1e5)
+summary_callback()

src/Trixi.jl

@@ -18,7 +18,7 @@ module Trixi
 # Include other packages that are used in Trixi
 # (standard library packages first, other packages next, all of them sorted alphabetically)
 
-using LinearAlgebra: dot, mul!, norm, cross, normalize
+using LinearAlgebra: diag, dot, mul!, norm, cross, normalize
 using Printf: @printf, @sprintf, println
 
 # import @reexport now to make it available for further imports/exports
@@ -47,6 +47,11 @@ using TimerOutputs: TimerOutputs, @notimeit, TimerOutput, print_timer, reset_tim
 @reexport using UnPack: @unpack
 using UnPack: @pack!
 
+# finite difference SBP operators
+using SummationByPartsOperators: AbstractDerivativeOperator, DerivativeOperator, grid
+import SummationByPartsOperators: integrate, left_boundary_weight, right_boundary_weight
+@reexport using SummationByPartsOperators:
+  SummationByPartsOperators, derivative_operator
 
 # Define the entry points of our type hierarchy, e.g.
 #     AbstractEquations, AbstractSemidiscretization etc.
@@ -139,9 +144,11 @@ export TreeMesh, CurvedMesh, UnstructuredQuadMesh, P4estMesh
 
 export DG,
        DGSEM, LobattoLegendreBasis,
-       VolumeIntegralWeakForm, VolumeIntegralFluxDifferencing,
+       VolumeIntegralWeakForm, VolumeIntegralStrongForm,
+       VolumeIntegralFluxDifferencing,
        VolumeIntegralPureLGLFiniteVolume,
        VolumeIntegralShockCapturingHG, IndicatorHennemannGassner,
+       SurfaceIntegralWeakForm, SurfaceIntegralStrongForm,
        MortarL2
 
 export nelements, nnodes, nvariables,

src/auxiliary/precompile.jl

@@ -351,24 +351,24 @@ function _precompile_manual_()
     mortar_type = mortar_type_dgsem(RealT, nnodes_)
 
     # 1D, serial
-    @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries),Array{Int,1},TreeMesh{1,Trixi.SerialTree{1}},Trixi.ElementContainer1D{RealT,uEltype,1,polydeg}})
-    @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces),Array{Int,1},TreeMesh{1,Trixi.SerialTree{1}},Trixi.ElementContainer1D{RealT,uEltype,1,polydeg}})
+    @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries),Array{Int,1},TreeMesh{1,Trixi.SerialTree{1}},Trixi.ElementContainer1D{RealT,uEltype}})
+    @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces),Array{Int,1},TreeMesh{1,Trixi.SerialTree{1}},Trixi.ElementContainer1D{RealT,uEltype}})
 
     # 2D, serial
-    @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries),Array{Int,1},TreeMesh{2,Trixi.SerialTree{2}},Trixi.ElementContainer2D{RealT,uEltype,1,polydeg}})
-    @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces),Array{Int,1},TreeMesh{2,Trixi.SerialTree{2}},Trixi.ElementContainer2D{RealT,uEltype,1,polydeg}})
-    @assert Base.precompile(Tuple{typeof(Trixi.init_mortars),Array{Int,1},TreeMesh{2,Trixi.SerialTree{2}},Trixi.ElementContainer2D{RealT,uEltype,1,polydeg},mortar_type})
+    @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries),Array{Int,1},TreeMesh{2,Trixi.SerialTree{2}},Trixi.ElementContainer2D{RealT,uEltype}})
+    @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces),Array{Int,1},TreeMesh{2,Trixi.SerialTree{2}},Trixi.ElementContainer2D{RealT,uEltype}})
+    @assert Base.precompile(Tuple{typeof(Trixi.init_mortars),Array{Int,1},TreeMesh{2,Trixi.SerialTree{2}},Trixi.ElementContainer2D{RealT,uEltype},mortar_type})
 
     # 2D, parallel
-    @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries),Array{Int,1},TreeMesh{2,Trixi.ParallelTree{2}},Trixi.ElementContainer2D{RealT,uEltype,1,polydeg}})
-    @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces),Array{Int,1},TreeMesh{2,Trixi.ParallelTree{2}},Trixi.ElementContainer2D{RealT,uEltype,1,polydeg}})
-    @assert Base.precompile(Tuple{typeof(Trixi.init_mortars),Array{Int,1},TreeMesh{2,Trixi.ParallelTree{2}},Trixi.ElementContainer2D{RealT,uEltype,1,polydeg},mortar_type})
-    @assert Base.precompile(Tuple{typeof(Trixi.init_mpi_interfaces),Array{Int,1},TreeMesh{2,Trixi.ParallelTree{2}},Trixi.ElementContainer2D{RealT,uEltype,1,polydeg}})
+    @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries),Array{Int,1},TreeMesh{2,Trixi.ParallelTree{2}},Trixi.ElementContainer2D{RealT,uEltype}})
+    @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces),Array{Int,1},TreeMesh{2,Trixi.ParallelTree{2}},Trixi.ElementContainer2D{RealT,uEltype}})
+    @assert Base.precompile(Tuple{typeof(Trixi.init_mortars),Array{Int,1},TreeMesh{2,Trixi.ParallelTree{2}},Trixi.ElementContainer2D{RealT,uEltype},mortar_type})
+    @assert Base.precompile(Tuple{typeof(Trixi.init_mpi_interfaces),Array{Int,1},TreeMesh{2,Trixi.ParallelTree{2}},Trixi.ElementContainer2D{RealT,uEltype}})
 
     # 3D, serial
-    @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries),Array{Int,1},TreeMesh{3,Trixi.SerialTree{3}},Trixi.ElementContainer3D{RealT,uEltype,1,polydeg}})
-    @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces),Array{Int,1},TreeMesh{3,Trixi.SerialTree{3}},Trixi.ElementContainer3D{RealT,uEltype,1,polydeg}})
-    @assert Base.precompile(Tuple{typeof(Trixi.init_mortars),Array{Int,1},TreeMesh{3,Trixi.SerialTree{3}},Trixi.ElementContainer3D{RealT,uEltype,1,polydeg},mortar_type})
+    @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries),Array{Int,1},TreeMesh{3,Trixi.SerialTree{3}},Trixi.ElementContainer3D{RealT,uEltype}})
+    @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces),Array{Int,1},TreeMesh{3,Trixi.SerialTree{3}},Trixi.ElementContainer3D{RealT,uEltype}})
+    @assert Base.precompile(Tuple{typeof(Trixi.init_mortars),Array{Int,1},TreeMesh{3,Trixi.SerialTree{3}},Trixi.ElementContainer3D{RealT,uEltype},mortar_type})
   end
 
   # various `show` methods

src/callbacks_step/amr_dg1d.jl

@@ -23,7 +23,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1},
     # re-initialize elements container
     @unpack elements = cache
     resize!(elements, length(leaf_cell_ids))
-    init_elements!(elements, leaf_cell_ids, mesh, dg.basis.nodes)
+    init_elements!(elements, leaf_cell_ids, mesh, dg.basis)
     @assert nelements(dg, cache) > old_n_elements
 
     resize!(u_ode, nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache))
@@ -137,7 +137,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1},
     # re-initialize elements container
     @unpack elements = cache
     resize!(elements, length(leaf_cell_ids))
-    init_elements!(elements, leaf_cell_ids, mesh, dg.basis.nodes)
+    init_elements!(elements, leaf_cell_ids, mesh, dg.basis)
     @assert nelements(dg, cache) < old_n_elements
 
     resize!(u_ode, nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache))

src/callbacks_step/amr_dg2d.jl

@@ -69,7 +69,7 @@ function reinitialize_containers!(mesh::TreeMesh{2}, equations, dg::DGSEM, cache
   # re-initialize elements container
   @unpack elements = cache
   resize!(elements, length(leaf_cell_ids))
-  init_elements!(elements, leaf_cell_ids, mesh, dg.basis.nodes)
+  init_elements!(elements, leaf_cell_ids, mesh, dg.basis)
 
   # re-initialize interfaces container
   @unpack interfaces = cache

src/callbacks_step/amr_dg3d.jl

@@ -23,7 +23,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{3},
     # re-initialize elements container
     @unpack elements = cache
     resize!(elements, length(leaf_cell_ids))
-    init_elements!(elements, leaf_cell_ids, mesh, dg.basis.nodes)
+    init_elements!(elements, leaf_cell_ids, mesh, dg.basis)
     @assert nelements(dg, cache) > old_n_elements
 
     resize!(u_ode, nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache))
@@ -177,7 +177,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{3},
     # re-initialize elements container
     @unpack elements = cache
     resize!(elements, length(leaf_cell_ids))
-    init_elements!(elements, leaf_cell_ids, mesh, dg.basis.nodes)
+    init_elements!(elements, leaf_cell_ids, mesh, dg.basis)
     @assert nelements(dg, cache) < old_n_elements
 
     resize!(u_ode, nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache))

src/callbacks_step/analysis.jl

@@ -20,7 +20,7 @@ solution and integrated over the computational domain.
 See `Trixi.analyze`, `Trixi.pretty_form_utf`, `Trixi.pretty_form_ascii` for further
 information on how to create custom analysis quantities.
 """
-mutable struct AnalysisCallback{Analyzer<:SolutionAnalyzer, AnalysisIntegrals, InitialStateIntegrals, Cache}
+mutable struct AnalysisCallback{Analyzer, AnalysisIntegrals, InitialStateIntegrals, Cache}
   start_time::Float64
   interval::Int
   save_analysis::Bool
@@ -178,7 +178,7 @@ function (analysis_callback::AnalysisCallback)(integrator)
     mpi_println()
     mpi_println("─"^100)
     # TODO: Taal refactor, polydeg is specific to DGSEM
-    mpi_println(" Simulation running '", get_name(equations), "' with polydeg = ", polydeg(solver))
+    mpi_println(" Simulation running '", get_name(equations), "' with ", summary(solver))
     mpi_println("─"^100)
     mpi_println(" #timesteps:     " * @sprintf("% 14d", iter) *
                 "               " *

src/callbacks_step/analysis_dg2d.jl

@@ -173,7 +173,7 @@ end
 
 function integrate(func::Func, u,
                    mesh::Union{TreeMesh{2}, CurvedMesh{2}, UnstructuredQuadMesh, P4estMesh{2}},
-                   equations, dg::DGSEM, cache; normalize=true) where {Func}
+                   equations, dg::DG, cache; normalize=true) where {Func}
   integrate_via_indices(u, mesh, equations, dg, cache; normalize=normalize) do u, i, j, element, equations, dg
     u_local = get_node_vars(u, equations, dg, i, j, element)
     return func(u_local, equations)

src/equations/equations.jl

@@ -53,6 +53,9 @@ end
 
 @inline Base.ndims(::AbstractEquations{NDIMS}) where NDIMS = NDIMS
 
+# equations act like scalars in broadcasting
+Base.broadcastable(equations::AbstractEquations) = Ref(equations)
+
 
 """
     flux(u, orientation_or_normal, equations)

src/solvers/dg_common.jl

@@ -5,7 +5,7 @@ function allocate_coefficients(mesh::AbstractMesh, equations, dg::DG, cache)
   zeros(eltype(cache.elements), nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache))
 end
 
-@inline function wrap_array(u_ode::AbstractVector, mesh::AbstractMesh, equations, dg::DG, cache)
+@inline function wrap_array(u_ode::AbstractVector, mesh::AbstractMesh, equations, dg::DGSEM, cache)
   @boundscheck begin
     @assert length(u_ode) == nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)
   end
@@ -46,6 +46,11 @@ end
   end
 end
 
+# General fallback
+@inline function wrap_array(u_ode::AbstractVector, mesh::AbstractMesh, equations, dg::DG, cache)
+  wrap_array_native(u_ode, mesh, equations, dg, cache)
+end
+
 # Like `wrap_array`, but guarantees to return a plain `Array`, which can be better
 # for interfacing with external C libraries (MPI, HDF5, visualization),
 # writing solution files etc.

src/solvers/dg_curved/dg.jl

@@ -17,30 +17,32 @@ end
 
 @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation,
                                                   boundary_condition::BoundaryConditionPeriodic,
-                                                  mesh::CurvedMesh, equations, dg::DG, cache,
+                                                  mesh::CurvedMesh, equations,
+                                                  surface_integral, dg::DG, cache,
                                                   direction, node_indices, surface_node_indices, element)
   @assert isperiodic(mesh, orientation)
 end
 
 
 @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation,
                                                   boundary_condition,
-                                                  mesh::CurvedMesh, equations,dg::DG, cache,
+                                                  mesh::CurvedMesh, equations,
+                                                  surface_integral, dg::DG, cache,
                                                   direction, node_indices, surface_node_indices, element)
   @unpack node_coordinates, contravariant_vectors, inverse_jacobian = cache.elements
-  @unpack surface_flux = dg
+  @unpack surface_flux = surface_integral
 
   u_inner = get_node_vars(u, equations, dg, node_indices..., element)
   x = get_node_coords(node_coordinates, equations, dg, node_indices..., element)
 
-  # If the mapping is orientation-reversing, the contravariant vectors' orientation 
-  # is reversed as well. The normal vector must be oriented in the direction 
+  # If the mapping is orientation-reversing, the contravariant vectors' orientation
+  # is reversed as well. The normal vector must be oriented in the direction
   # from `left_element` to `right_element`, or the numerical flux will be computed
   # incorrectly (downwind direction).
   sign_jacobian = sign(inverse_jacobian[node_indices..., element])
 
   # Contravariant vector Ja^i is the normal vector
-  normal = sign_jacobian * get_contravariant_vector(orientation, contravariant_vectors, 
+  normal = sign_jacobian * get_contravariant_vector(orientation, contravariant_vectors,
                                                     node_indices..., element)
 
   # If the mapping is orientation-reversing, the normal vector will be reversed (see above).

src/solvers/dg_curved/dg_1d.jl

@@ -13,15 +13,15 @@ function rhs!(du, u, t,
 
   # Calculate interface and boundary fluxes
   @timed timer() "interface flux" calc_interface_flux!(
-    cache, u, mesh, equations, dg)
+    cache, u, mesh, equations, dg.surface_integral, dg)
 
   # Calculate boundary fluxes
   @timed timer() "boundary flux" calc_boundary_flux!(
-    cache, u, t, boundary_conditions, mesh, equations, dg)
+    cache, u, t, boundary_conditions, mesh, equations, dg.surface_integral, dg)
 
   # Calculate surface integrals
   @timed timer() "surface integral" calc_surface_integral!(
-    du, mesh, equations, dg, cache)
+    du, u, mesh, equations, dg.surface_integral, dg, cache)
 
   # Apply Jacobian from mapping to reference element
   @timed timer() "Jacobian" apply_jacobian!(
@@ -36,8 +36,8 @@ end
 
 
 function calc_interface_flux!(cache, u, mesh::CurvedMesh{1},
-                              equations, dg::DG)
-  @unpack surface_flux = dg
+                              equations, surface_integral, dg::DG)
+  @unpack surface_flux = surface_integral
 
   @threaded for element in eachelement(dg, cache)
     left_element = cache.elements.left_neighbors[1, element]
@@ -46,11 +46,11 @@ function calc_interface_flux!(cache, u, mesh::CurvedMesh{1},
       u_ll = get_node_vars(u, equations, dg, nnodes(dg), left_element)
       u_rr = get_node_vars(u, equations, dg, 1,          element)
 
-      flux = surface_flux(u_ll, u_rr, 1, equations)
+      f1 = surface_flux(u_ll, u_rr, 1, equations)
 
       for v in eachvariable(equations)
-        cache.elements.surface_flux_values[v, 2, left_element] = flux[v]
-        cache.elements.surface_flux_values[v, 1, element] = flux[v]
+        cache.elements.surface_flux_values[v, 2, left_element] = f1[v]
+        cache.elements.surface_flux_values[v, 1, element] = f1[v]
       end
     end
   end
@@ -61,21 +61,21 @@ end
 
 # TODO: Taal dimension agnostic
 function calc_boundary_flux!(cache, u, t, boundary_condition::BoundaryConditionPeriodic,
-                             mesh::CurvedMesh{1}, equations, dg::DG)
+                             mesh::CurvedMesh{1}, equations, surface_integral, dg::DG)
   @assert isperiodic(mesh)
 end
 
 
 function calc_boundary_flux!(cache, u, t, boundary_condition,
-                             mesh::CurvedMesh{1}, equations, dg::DG)
+                             mesh::CurvedMesh{1}, equations, surface_integral, dg::DG)
   calc_boundary_flux!(cache, u, t, (boundary_condition, boundary_condition),
-                      mesh, equations, dg)
+                      mesh, equations, surface_integral, dg)
 end
 
 
 function calc_boundary_flux!(cache, u, t, boundary_conditions::Union{NamedTuple,Tuple},
-                             mesh::CurvedMesh{1}, equations, dg::DG)
-  @unpack surface_flux = dg
+                             mesh::CurvedMesh{1}, equations, surface_integral, dg::DG)
+  @unpack surface_flux = surface_integral
   @unpack surface_flux_values, node_coordinates = cache.elements
 
   orientation = 1