Skip to content

Commit

Permalink
more Fourier operator functions (#280)
Browse files Browse the repository at this point in the history 
* more methods for multiplication and division of Fourier operators and numbers

* add tests

* operations with UniformScaling
  • Loading branch information
@ranocha
ranocha committed Jul 23, 2024
1 parent 87bfd2c commit 560c6d8
Show file tree
Hide file tree
Showing 3 changed files with 128 additions and 17 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "SummationByPartsOperators"
uuid = "9f78cca6-572e-554e-b819-917d2f1cf240"
author = ["Hendrik Ranocha"]
version = "0.5.64"
version = "0.5.65"

[deps]
ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197"
Expand Down
88 changes: 88 additions & 0 deletions src/fourier_operators.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -530,6 +530,23 @@ function Base.:+(rat1::Union{FourierDerivativeOperator,FourierPolynomialDerivati
FourierRationalDerivativeOperator(rat1) + rat2
end

function Base.:+(rat::FourierRationalDerivativeOperator, scaling::UniformScaling)
rat + scaling.λ * rat.D1^0
end

function Base.:+(scaling::UniformScaling, rat::FourierRationalDerivativeOperator)
rat + scaling
end

function Base.:-(rat::FourierRationalDerivativeOperator)
@unpack num_coef = rat
for idx in eachindex(num_coef)
num_coef = Base.setindex(num_coef, -num_coef[idx], idx)
end

FourierRationalDerivativeOperator(rat.D1, num_coef, rat.den_coef)
end

function Base.:-(rat1::FourierRationalDerivativeOperator, rat2::FourierRationalDerivativeOperator)
T = eltype(rat1)
@argcheck T == eltype(rat2) ArgumentError
Expand All @@ -551,6 +568,14 @@ function Base.:-(rat1::Union{FourierDerivativeOperator,FourierPolynomialDerivati
FourierRationalDerivativeOperator(rat1) - rat2
end

function Base.:-(rat::FourierRationalDerivativeOperator, scaling::UniformScaling)
rat - scaling.λ * rat.D1^0
end

function Base.:-(scaling::UniformScaling, rat::FourierRationalDerivativeOperator)
scaling + (-rat)
end

function Base.:*(rat1::FourierRationalDerivativeOperator, rat2::FourierRationalDerivativeOperator)
T = eltype(rat1)
@argcheck T == eltype(rat2) ArgumentError
Expand All @@ -572,6 +597,27 @@ function Base.:*(rat1::Union{FourierDerivativeOperator,FourierPolynomialDerivati
FourierRationalDerivativeOperator(rat1) * rat2
end

function Base.:*(factor::Union{Real,Integer}, rat::FourierRationalDerivativeOperator)
@unpack num_coef = rat
for idx in eachindex(num_coef)
num_coef = Base.setindex(num_coef, factor * num_coef[idx], idx)
end

FourierRationalDerivativeOperator(rat.D1, num_coef, rat.den_coef)
end

function Base.:*(poly::FourierRationalDerivativeOperator, factor::Union{Real,Integer})
factor * poly
end

function Base.:*(rat::FourierRationalDerivativeOperator, scaling::UniformScaling)
scaling.λ * rat
end

function Base.:*(scaling::UniformScaling, rat::FourierRationalDerivativeOperator)
rat * scaling
end

function Base.:/(rat1::FourierRationalDerivativeOperator, rat2::FourierRationalDerivativeOperator)
T = eltype(rat1)
@argcheck T == eltype(rat2) ArgumentError
Expand All @@ -593,6 +639,48 @@ function Base.:/(rat1::Union{FourierDerivativeOperator,FourierPolynomialDerivati
FourierRationalDerivativeOperator(rat1) / rat2
end

function Base.:/(factor::Union{Real,Integer}, rat::FourierRationalDerivativeOperator)
@unpack den_coef = rat
for idx in eachindex(den_coef)
den_coef = Base.setindex(den_coef, factor * den_coef[idx], idx)
end

FourierRationalDerivativeOperator(rat.D1, den_coef, rat.num_coef)
end

function Base.:/(rat::FourierRationalDerivativeOperator, factor::Union{Real,Integer})
@unpack num_coef = rat
for idx in eachindex(num_coef)
num_coef = Base.setindex(num_coef, num_coef[idx] / factor, idx)
end

FourierRationalDerivativeOperator(rat.D1, num_coef, rat.den_coef)
end

function Base.:\(factor::Union{Real,Integer}, rat::FourierRationalDerivativeOperator)
rat / factor
end

function Base.:\(rat::FourierRationalDerivativeOperator, factor::Union{Real,Integer})
inv(rat) * factor
end

function Base.:/(rat::FourierRationalDerivativeOperator, scaling::UniformScaling)
rat / scaling.λ
end

function Base.:/(scaling::UniformScaling, rat::FourierRationalDerivativeOperator)
scaling.λ / rat
end

function Base.:\(rat::FourierRationalDerivativeOperator, scaling::UniformScaling)
rat \ scaling.λ
end

function Base.:\(scaling::UniformScaling, rat::FourierRationalDerivativeOperator)
scaling.λ \ rat
end


function mul!(dest::AbstractVector{T}, rat::FourierRationalDerivativeOperator, u::AbstractVector{T}) where {T}
@unpack D1, num_coef, num_coef_end, den_coef, den_coef_end = rat
Expand Down
55 changes: 39 additions & 16 deletions test/fourier_operators_test.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,13 +54,13 @@ for T in (Float32, Float64)
end


# Check Fourier polynomial operators
# Check Fourier polynomial/rational operators
for T in (Float32, Float64)
xmin = -one(T)
xmax = one(T)
for N in (8, 9)
D = fourier_derivative_operator(xmin, xmax, N)
x = grid(D)
D = @inferred fourier_derivative_operator(xmin, xmax, N)
x = @inferred grid(D)
u = @. sinpi(x) - cospi(x)^2 + exp(sinpi(x))
println(devnull, D)

Expand All @@ -72,37 +72,60 @@ 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
@test @inferred(SummationByPartsOperators.xmin(D^2)) ≈ xmin
@test @inferred(SummationByPartsOperators.xmax(D^2)) ≈ xmax

poly = (I + 2D + 5*D^2) * (2I * D - D^3 * 5I) * (D*2 - D^2 * 5)
poly = @inferred (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)

@test (I + one(T)/2*D) * u ≈ (u + D*u ./ 2)
@test @inferred(I + one(T)/2*D) * u ≈ (u + D*u ./ 2)

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
@test @inferred(SummationByPartsOperators.xmin(inv(I - D^2))) ≈ xmin
@test @inferred(SummationByPartsOperators.xmax(inv(I - D^2))) ≈ xmax

v = (I - D^2) \ u
v = @inferred(I - D^2) \ u
@test D * v ≈ (D / (I - D^2)) * u

rat = (I - D^2) / (I + D^4)
rat = @inferred((I - D^2) / (I + D^4))
println(devnull, rat)
v = rat * u
@test (I - D^2) \ (v + D^4 * v) ≈ u

rat1 = (I - D^2) / (I + D^4)
rat2 = (I + D^4) / (I - D^2)
rat3 = (I - D^4) / (I - D^2)
@test (rat2 + rat3) * u ≈ 2 * ((I - D^2) \ u)
@test (rat1 * rat2) * u ≈ u
rat1 = @inferred((I - D^2) / (I + D^4))
rat2 = @inferred((I + D^4) / (I - D^2))
rat3 = @inferred((I - D^4) / (I - D^2))
@test @inferred(rat2 + rat3) * u ≈ 2 * ((I - D^2) \ u)
@test @inferred(rat1 * rat2) * u ≈ u

@test integrate(u, D) ≈ sum(mass_matrix(D) * u)
@test integrate(u->u^2, u, D) ≈ dot(u, mass_matrix(D), u)

# combine rational operators and scalars
@test @inferred(2 * rat1) * u ≈ rat1 * (2 * u)
@test @inferred(rat1 * 2) * u ≈ rat1 * (2 * u)
@test @inferred(2 / rat1) * u ≈ inv(rat1) * (2 * u)
@test @inferred(rat1 / 2) * u ≈ rat1 * (u / 2)
@test @inferred(2 \ rat1) * u ≈ rat1 * (u / 2)
@test @inferred(rat1 \ 2) * u ≈ inv(rat1) * (2 * u)

# combine rational operators and uniform scaling via *, /, \
@test @inferred((2 * I) * rat1) == 2 * rat1
@test @inferred(rat1 * (2 * I)) == 2 * rat1
@test @inferred((2 * I) / rat1) == 2 / rat1
@test @inferred(rat1 / (2 * I)) == rat1 / 2
@test @inferred((2 * I) \ rat1) == rat1 / 2
@test @inferred(rat1 \ (2 * I)) == 2 / rat1

# combine rational operators and uniform scaling via +, -
@test @inferred(I + rat1) == D^0 + rat1
@test @inferred(rat1 + I) == D^0 + rat1
@test @inferred(I - rat1) == D^0 - rat1
@test @inferred(rat1 - I) == rat1 - D^0
@test @inferred(-rat1) == @inferred(-1 * rat1)
end
end

Expand Down

2 comments on commit 560c6d8

@ranocha
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/111578

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.5.65 -m "<description of version>" 560c6d830a8dfa47676671b858024e9c552bf843
git push origin v0.5.65

Please sign in to comment.