Function space SBP operators

.gitignore

@@ -9,3 +9,6 @@ docs/build
 docs/src/code_of_conduct.md
 docs/src/contributing.md
 docs/src/license.md
+
+run
+run/*

Project.toml

@@ -27,12 +27,14 @@ Unrolled = "9602ed7d-8fef-5bc8-8597-8f21381861e8"
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
 
 [extensions]
 SummationByPartsOperatorsBandedMatricesExt = "BandedMatrices"
 SummationByPartsOperatorsDiffEqCallbacksExt = "DiffEqCallbacks"
 SummationByPartsOperatorsForwardDiffExt = "ForwardDiff"
+SummationByPartsOperatorsOptimExt = "Optim"
 SummationByPartsOperatorsStructArraysExt = "StructArrays"
 
 [compat]
@@ -46,6 +48,7 @@ InteractiveUtils = "1"
 LinearAlgebra = "1"
 LoopVectorization = "0.12.22"
 MuladdMacro = "0.2"
+Optim = "1"
 PolynomialBases = "0.4.15"
 PrecompileTools = "1.0.1"
 RecursiveArrayTools = "2.11, 3"
@@ -64,4 +67,5 @@ julia = "1.6"
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"

ext/SummationByPartsOperatorsOptimExt.jl

@@ -0,0 +1,19 @@
+module SummationByPartsOperatorsOptimExt
+
+if isdefined(Base, :get_extension)
+    using Optim: LBFGS, optimize
+else
+    using ..Optim: LBFGS, optimize
+end
+
+using SummationByPartsOperators: SummationByPartsOperators, GlaubitzNordströmÖffner2023
+
+function SummationByPartsOperators.construct_function_space_operator(basis, x_min, x_max, nodes, ::GlaubitzNordströmÖffner2023)
+    # TODO: implement the construction of the function space operator by solving an optimization problem
+    N = length(nodes)
+    weights = ones(N)
+    D = zeros(N, N)
+    return weights, D
+end
+
+end # module

src/SummationByPartsOperators.jl

@@ -100,6 +100,7 @@ function __init__()
     @require BandedMatrices="aae01518-5342-5314-be14-df237901396f" include("../ext/SummationByPartsOperatorsBandedMatricesExt.jl")
     @require DiffEqCallbacks="459566f4-90b8-5000-8ac3-15dfb0a30def" include("../ext/SummationByPartsOperatorsDiffEqCallbacksExt.jl")
     @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" include("../ext/SummationByPartsOperatorsForwardDiffExt.jl")
+    @require Optim="429524aa-4258-5aef-a3af-852621145aeb" include("../ext/SummationByPartsOperatorsOptimExt.jl")
     @require StructArrays="09ab397b-f2b6-538f-b94a-2f83cf4a842a" include("../ext/SummationByPartsOperatorsStructArraysExt.jl")
   end
 end
@@ -109,6 +110,7 @@ include("fourier_operators.jl")
 include("fourier_operators_2d.jl")
 include("legendre_operators.jl")
 include("matrix_operators.jl")
+include("function_space_operators.jl")
 include("upwind_operators.jl")
 include("SBP_coefficients/MattssonNordström2004.jl")
 include("SBP_coefficients/MattssonSvärdNordström2004.jl")
@@ -139,7 +141,7 @@ export PeriodicDerivativeOperator, PeriodicDissipationOperator,
        FourierPolynomialDerivativeOperator, FourierRationalDerivativeOperator,
        FourierDerivativeOperator2D,
        LegendreDerivativeOperator, LegendreSecondDerivativeOperator,
-       MatrixDerivativeOperator,
+       MatrixDerivativeOperator, FunctionSpaceOperator,
        UpwindOperators, PeriodicUpwindOperators
 export FilterCallback, ConstantFilter, ExponentialFilter
 export SafeMode, FastMode, ThreadedMode
@@ -169,6 +171,7 @@ export MattssonNordström2004, MattssonSvärdNordström2004, MattssonSvärdShoey
        SharanBradyLivescu2022
 export Tadmor1989, MadayTadmor1989, Tadmor1993,
        TadmorWaagan2012Standard, TadmorWaagan2012Convergent
+export GlaubitzNordströmÖffner2023
 
 export BurgersPeriodicSemidiscretization, BurgersNonperiodicSemidiscretization,
        CubicPeriodicSemidiscretization, CubicNonperiodicSemidiscretization,

src/function_space_operators.jl

@@ -0,0 +1,49 @@
+
+"""
+    GlaubitzNordströmÖffner2023()
+
+Function space SBP (FSBP) operators given in
+- Glaubitz, Nordström, Öffner (2023)
+  Summation-by-parts operators for general function spaces.
+  SIAM Journal on Numerical Analysis 61, 2, pp. 733-754.
+
+See also
+- Glaubitz, Nordström, Öffner (2024)
+  An optimization-based construction procedure for function space based
+  summation-by-parts operators on arbitrary grids.
+  arXiv, arXiv:2405.08770v1.
+"""
+struct GlaubitzNordströmÖffner2023 <: SourceOfCoefficients end
+
+function Base.show(io::IO, source::GlaubitzNordströmÖffner2023)
+  if get(io, :compact, false)
+    summary(io, source)
+  else
+      print(io,
+          "Glaubitz, Nordström, Öffner (2023) \n",
+          "  Summation-by-parts operators for general function spaces \n",
+          "  SIAM Journal on Numerical Analysis 61, 2, pp. 733-754. \n",
+          "See also \n",
+          "  Glaubitz, Nordström, Öffner (2024) \n",
+          "  An optimization-based construction procedure for function \n",
+          "    space based summation-by-parts operators on arbitrary grids \n",
+          "  arXiv, arXiv:2405.08770v1.")
+  end
+end
+
+"""
+    FunctionSpaceOperator(basis, x_min, x_max, nodes, source)
+
+Construct a `MatrixDerivativeOperator` that represents a derivative operator in a function space.
+"""
+function FunctionSpaceOperator(basis::Vector{BasisFunctions},
+                               x_min::T, x_max::T, nodes::Vector{T},
+                               source::SourceOfCoefficients) where {BasisFunctions, T, SourceOfCoefficients}
+
+  weights, D = construct_function_space_operator(basis, x_min, x_max, nodes, source)
+  accuracy_order = length(basis) # TODO: this might needs to be adjusted
+  return MatrixDerivativeOperator(x_min, x_max, nodes, weights, D, accuracy_order, source)
+end
+
+# This function is extended in the package extension SummationByPartsOperatorsOptimExt
+function construct_function_space_operator end