fix xmax for Fourier operators (#256)

ranocha · Mar 11, 2024 · 65aa3cd · 65aa3cd · ranocha · Mar 11, 2024
1 parent 360c87f
commit 65aa3cd
Show file tree

Hide file tree

Showing 4 changed files with 22 additions and 3 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "SummationByPartsOperators"
 uuid = "9f78cca6-572e-554e-b819-917d2f1cf240"
 author = ["Hendrik Ranocha"]
-version = "0.5.55"
+version = "0.5.56"
 
 [deps]
 ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197"

diff --git a/src/fourier_operators.jl b/src/fourier_operators.jl
@@ -11,8 +11,8 @@ See also [`fourier_derivative_operator`](@ref).
 @auto_hash_equals struct FourierDerivativeOperator{T<:Real, Grid, RFFT, BRFFT} <: AbstractPeriodicDerivativeOperator{T}
     jac::T
     Δx::T
-    grid_compute::Grid   # N-1 nodes, including the left and excluding the right boundary
-    grid_evaluate::Grid #  N  nodes, including both boundaries
+    grid_compute::Grid   # N nodes, including the left and excluding the right boundary
+    grid_evaluate::Grid #  N+1 nodes, including both boundaries
     tmp::Vector{Complex{T}}
     rfft_plan::RFFT
     brfft_plan::BRFFT
@@ -84,6 +84,9 @@ function Base.show(io::IO, D::FourierDerivativeOperator)
     end
 end
 
+xmin(D::FourierDerivativeOperator) = first(D.grid_evaluate)
+xmax(D::FourierDerivativeOperator) = last(D.grid_evaluate)
+
 
 function mul!(dest::AbstractVector{T}, D::FourierDerivativeOperator, u::AbstractVector{T}) where {T}
     @unpack jac, tmp, rfft_plan, brfft_plan = D
@@ -212,6 +215,9 @@ function Base.show(io::IO, poly::FourierPolynomialDerivativeOperator)
     end
 end
 
+xmin(poly::FourierPolynomialDerivativeOperator) = xmin(poly.D1)
+xmax(poly::FourierPolynomialDerivativeOperator) = xmax(poly.D1)
+
 
 function Base.:*(D1::FourierDerivativeOperator, D2::FourierDerivativeOperator)
     T = eltype(D1)
@@ -431,6 +437,9 @@ end
 Base.size(rat::FourierRationalDerivativeOperator) = size(rat.D1)
 grid(rat::FourierRationalDerivativeOperator) = grid(rat.D1)
 
+xmin(rat::FourierRationalDerivativeOperator) = xmin(rat.D1)
+xmax(rat::FourierRationalDerivativeOperator) = xmax(rat.D1)
+
 function Base.show(io::IO, rat::FourierRationalDerivativeOperator)
     if get(io, :compact, false)
         print(io, "Rational Fourier operator")

diff --git a/test/fourier_operators_test.jl b/test/fourier_operators_test.jl
@@ -19,6 +19,8 @@ for T in (Float32, Float64)
         println(devnull, D)
         @test SummationByPartsOperators.derivative_order(D) == 1
         @test issymmetric(D) == false
+        @test SummationByPartsOperators.xmin(D) ≈ xmin
+        @test SummationByPartsOperators.xmax(D) ≈ xmax
         u = compute_coefficients(zero, D)
         res = D*u
         for k in 0:(N÷2)-1
@@ -58,6 +60,9 @@ for T in (Float32, Float64)
         @test !issymmetric(I + D)
         @test !issymmetric(D - I)
 
+        @test SummationByPartsOperators.xmin(D^2) ≈ xmin
+        @test SummationByPartsOperators.xmax(D^2) ≈ xmax
+
         poly = (I + 2D + 5*D^2) * (2I * D - D^3 * 5I) * (D*2 - D^2 * 5)
         @test poly.coef == (0.0, 0.0, 4.0, -2.0, -10.0, -45.0, 0.0, 125.0)
         println(devnull, poly)
@@ -67,6 +72,9 @@ for T in (Float32, Float64)
         v = (I - D^2) * u
         @test inv(I - D^2) * v ≈ u
 
+        @test SummationByPartsOperators.xmin(inv(I - D^2)) ≈ xmin
+        @test SummationByPartsOperators.xmax(inv(I - D^2)) ≈ xmax
+
         v = (I - D^2) \ u
         @test D * v ≈ (D / (I - D^2)) * u
 

diff --git a/test/legendre_operators_test.jl b/test/legendre_operators_test.jl
@@ -21,6 +21,8 @@ for T in (Float32, Float64)
         println(devnull, D)
         @test SummationByPartsOperators.derivative_order(D) == 1
         @test issymmetric(D) == false
+        @test SummationByPartsOperators.xmin(D) ≈ xmin
+        @test SummationByPartsOperators.xmax(D) ≈ xmax
         u = compute_coefficients(zero, D)
         res = D*u
         for k in 1:N-1