diff --git a/src/SummationByPartsOperators.jl b/src/SummationByPartsOperators.jl index 3e3011c8..fdc57fc3 100644 --- a/src/SummationByPartsOperators.jl +++ b/src/SummationByPartsOperators.jl @@ -69,7 +69,7 @@ include("second_order_eqs/wave_eq.jl") # exports export PeriodicDerivativeOperator, PeriodicDissipationOperator, - PeriodicRationalDerivativeOperator, + PeriodicRationalDerivativeOperator, PeriodicDerivativeOperatorQuotient, DerivativeOperator, DissipationOperator, VarCoefDerivativeOperator, SourceOfCoefficients, FourierDerivativeOperator, FourierConstantViscosity, diff --git a/src/fourier_operators.jl b/src/fourier_operators.jl index 1cd01699..2e1745f2 100644 --- a/src/fourier_operators.jl +++ b/src/fourier_operators.jl @@ -82,7 +82,7 @@ function mul!(dest::AbstractVector{T}, D::FourierDerivativeOperator, u::Abstract if iseven(N) @inbounds tmp[end] = zero(eltype(tmp)) else - @inbounds tmp[end] *= (N-1)*im * jac + @inbounds tmp[end] *= (length(tmp)-1)*im * jac end mul!(dest, brfft_plan, tmp) end @@ -103,7 +103,7 @@ end # if iseven(N) # @inbounds tmp[end] = zero(eltype(tmp)) # else -# @inbounds tmp[end] *= (N-1)*im * jac +# @inbounds tmp[end] *= (length(tmp)-1)*im * jac # end # mul!(dest, brfft_plan, tmp, α, β) # end @@ -149,6 +149,15 @@ function fourier_derivative_matrix(N, xmin::Real=0.0, xmax::Real=2π) end +function _coef_end(coef, D1::FourierDerivativeOperator{T}) where {T} + if isodd(size(D1, 1)) + return coef + else + return ntuple(i -> isodd(i) ? coef[i] : zero(T), length(coef)) + end +end +# fallback for periodic FD operators etc. +_coef_end(coef, D1) = coef struct FourierPolynomialDerivativeOperator{T<:Real, Grid, RFFT, BRFFT, N} <: AbstractPeriodicDerivativeOperator{T} D1::FourierDerivativeOperator{T,Grid,RFFT,BRFFT} @@ -156,10 +165,7 @@ struct FourierPolynomialDerivativeOperator{T<:Real, Grid, RFFT, BRFFT, N} <: Abs coef_end::NTuple{N,T} function FourierPolynomialDerivativeOperator(D1::FourierDerivativeOperator{T,Grid,RFFT,BRFFT}, coef::NTuple{N,T}) where {T<:Real, Grid, RFFT, BRFFT, N} - coef_end = ntuple(i -> isodd(i) ? coef[i] : zero(T), length(coef)) - if isodd(size(D1, 1)) - coef_end = coef - end + coef_end = _coef_end(coef, D1) new{T,Grid,RFFT,BRFFT,N}(D1, coef, coef_end) end end @@ -360,7 +366,7 @@ function mul!(dest::AbstractVector{T}, poly::FourierPolynomialDerivativeOperator tmp[j] *= evalpoly((j-1)*im * jac*N, coef) / N end # see e.g. Steven G. Johnson (2011) Notes on FFT based differentiation - @inbounds tmp[end] *= evalpoly((N-1)*im * jac*N, coef_end) / N + @inbounds tmp[end] *= evalpoly((length(tmp)-1)*im * jac*N, coef_end) / N mul!(dest, brfft_plan, tmp) end @@ -379,7 +385,7 @@ function LinearAlgebra.ldiv!(dest::AbstractVector{T}, rat::FourierPolynomialDeri tmp[j] /= (evalpoly((j-1)*im * jac*N, coef) * N) end # see e.g. Steven G. Johnson (2011) Notes on FFT based differentiation - @inbounds tmp[end] /= (evalpoly((N-1)*im * jac*N, coef_end) * N) + @inbounds tmp[end] /= (evalpoly((length(tmp)-1)*im * jac*N, coef_end) * N) mul!(dest, brfft_plan, tmp) end @@ -393,12 +399,8 @@ struct FourierRationalDerivativeOperator{T<:Real, Grid, RFFT, BRFFT, Nnum, Nden} den_coef_end::NTuple{Nden,T} function FourierRationalDerivativeOperator(D1::FourierDerivativeOperator{T,Grid,RFFT,BRFFT}, num_coef::NTuple{Nnum,T}, den_coef::NTuple{Nden,T}) where {T<:Real, Grid, RFFT, BRFFT, Nnum, Nden} - num_coef_end = ntuple(i -> isodd(i) ? num_coef[i] : zero(T), length(num_coef)) - den_coef_end = ntuple(i -> isodd(i) ? den_coef[i] : zero(T), length(den_coef)) - if isodd(size(D1, 1)) - num_coef_end = num_coef - den_coef_end = den_coef - end + num_coef_end = _coef_end(num_coef, D1) + den_coef_end = _coef_end(den_coef, D1) new{T,Grid,RFFT,BRFFT,Nnum,Nden}(D1, num_coef, num_coef_end, den_coef, den_coef_end) end end @@ -450,10 +452,6 @@ function Base.:/(num::FourierPolynomialDerivativeOperator, den::FourierDerivativ num / FourierPolynomialDerivativeOperator(den) end -function Base.://(num::FourierPolynomialDerivativeOperator, den::FourierPolynomialDerivativeOperator) - num / den -end - function Base.inv(rat::FourierRationalDerivativeOperator) FourierRationalDerivativeOperator(rat.D1, rat.den_coef, rat.num_coef) end @@ -547,18 +545,6 @@ function Base.:/(rat1::Union{FourierDerivativeOperator,FourierPolynomialDerivati FourierRationalDerivativeOperator(rat1) / rat2 end -function Base.://(rat1::FourierRationalDerivativeOperator, rat2::FourierRationalDerivativeOperator) - rat1 / rat2 -end - -function Base.://(rat1::FourierRationalDerivativeOperator, rat2::Union{FourierDerivativeOperator,FourierPolynomialDerivativeOperator}) - rat1 / rat2 -end - -function Base.://(rat1::Union{FourierDerivativeOperator,FourierPolynomialDerivativeOperator}, rat2::FourierRationalDerivativeOperator) - rat1 / rat2 -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 @@ -575,7 +561,7 @@ function mul!(dest::AbstractVector{T}, rat::FourierRationalDerivativeOperator, u tmp[j] *= evalpoly((j-1)*im * jac*N, num_coef) / (N * evalpoly((j-1)*im * jac*N, den_coef)) end # see e.g. Steven G. Johnson (2011) Notes on FFT based differentiation - @inbounds tmp[end] *= evalpoly((N-1)*im * jac*N, num_coef_end) / (N * evalpoly((N-1)*im * jac*N, den_coef_end)) + @inbounds tmp[end] *= evalpoly((length(tmp)-1)*im * jac*N, num_coef_end) / (N * evalpoly((length(tmp)-1)*im * jac*N, den_coef_end)) mul!(dest, brfft_plan, tmp) end @@ -606,6 +592,227 @@ end +struct PeriodicDerivativeOperatorQuotient{T<:Real, numDtype<:Union{PeriodicDerivativeOperator{T},FourierDerivativeOperator{T}}, denDtype<:Union{PeriodicDerivativeOperator{T},FourierDerivativeOperator{T}}, Nnum, Nden, RFFT, IRFFT} <: AbstractPeriodicDerivativeOperator{T} + num_D::numDtype + den_D::denDtype + num_coef::NTuple{Nnum,T} + den_coef::NTuple{Nden,T} + num_coef_end::NTuple{Nnum,T} + den_coef_end::NTuple{Nden,T} + tmp::Vector{Complex{T}} + num_eigval::Vector{Complex{T}} + den_eigval::Vector{Complex{T}} + rfft_plan::RFFT + irfft_plan::IRFFT + + function PeriodicDerivativeOperatorQuotient(num_D::numDtype, den_D::denDtype, num_coef::NTuple{Nnum,T}, den_coef::NTuple{Nden,T}, tmp::Vector{Complex{T}}, num_eigval::Vector{Complex{T}}, den_eigval::Vector{Complex{T}}, rfft_plan::RFFT, irfft_plan::IRFFT) where {T<:Real, numDtype<:Union{PeriodicDerivativeOperator{T},FourierDerivativeOperator{T}}, denDtype<:Union{PeriodicDerivativeOperator{T},FourierDerivativeOperator{T}}, Nnum, Nden, RFFT, IRFFT} + @argcheck length(irfft_plan) == length(tmp) DimensionMismatch + @argcheck length(irfft_plan) == length(num_eigval) DimensionMismatch + @argcheck length(irfft_plan) == length(den_eigval) DimensionMismatch + @argcheck length(irfft_plan) == (length(rfft_plan)÷2)+1 DimensionMismatch + @argcheck length(grid(num_D)) == length(rfft_plan) DimensionMismatch + @argcheck length(grid(den_D)) == length(rfft_plan) DimensionMismatch + + num_coef_end = _coef_end(num_coef, num_D) + den_coef_end = _coef_end(den_coef, den_D) + + new{T,numDtype,denDtype,Nnum,Nden,RFFT,IRFFT}(num_D, den_D, num_coef, den_coef, num_coef_end, den_coef_end, tmp, num_eigval, den_eigval, rfft_plan, irfft_plan) + end +end + +Base.size(quot::PeriodicDerivativeOperatorQuotient) = size(quot.num_D) +grid(quot::PeriodicDerivativeOperatorQuotient) = grid(quot.num_D) + +function Base.show(io::IO, quot::PeriodicDerivativeOperatorQuotient) + print(io, "Quotient of the polynomial with coefficients\n") + print(io, quot.num_coef) + print(io, "\nof the operator:\n") + print(io, quot.num_D) + print(io, "\nand the polynomial with coefficients\n") + print(io, quot.den_coef) + print(io, "\nof the operator:\n") + print(io, quot.den_D) +end + + +function _eigvals!(eigval::Vector{Complex{T}}, D::FourierDerivativeOperator{T}) where {T<:Real} + N = length(grid(D)) + jac = D.jac * N + @inbounds for idx in 1:(length(eigval)) + eigval[idx] = (idx-1)*im * jac + end + + eigval +end + +function Base.://(num::FourierDerivativeOperator, den::PeriodicRationalDerivativeOperator) + T = eltype(num) + @argcheck T == eltype(den) ArgumentError + @argcheck grid(num) == grid(den) ArgumentError + @argcheck den.den_coef == (one(T),) ArgumentError + + @unpack tmp, rfft_plan = num + irfft_plan = plan_irfft(tmp, length(grid(num))) + + num_eigval = similar(tmp) + _eigvals!(num_eigval, num) + + den_eigval = similar(tmp) + _eigvals!(den_eigval, den.D) + + PeriodicDerivativeOperatorQuotient(num, den.D, (zero(T), one(T)), den.num_coef, tmp, num_eigval, den_eigval, rfft_plan, irfft_plan) +end + +function Base.://(num::FourierPolynomialDerivativeOperator, den::PeriodicRationalDerivativeOperator) + T = eltype(num) + @argcheck T == eltype(den) ArgumentError + @argcheck grid(num) == grid(den) ArgumentError + @argcheck den.den_coef == (one(T),) ArgumentError + + @unpack tmp, rfft_plan = num.D1 + irfft_plan = plan_irfft(tmp, length(grid(num))) + + num_eigval = similar(tmp) + _eigvals!(num_eigval, num.D1) + + den_eigval = similar(tmp) + _eigvals!(den_eigval, den.D) + + PeriodicDerivativeOperatorQuotient(num.D1, den.D, num.coef, den.num_coef, tmp, num_eigval, den_eigval, rfft_plan, irfft_plan) +end + +function Base.://(num::PeriodicDerivativeOperator, den::PeriodicRationalDerivativeOperator) + T = eltype(num) + @argcheck T == eltype(den) ArgumentError + @argcheck grid(num) == grid(den) ArgumentError + @argcheck den.den_coef == (one(T),) ArgumentError + + x = grid(num) + u = zero.(x) + rfft_plan = plan_rfft(u) + tmp = rfft_plan * u + irfft_plan = plan_irfft(tmp, length(grid(num))) + + num_eigval = similar(tmp) + _eigvals!(num_eigval, num) + + den_eigval = similar(tmp) + _eigvals!(den_eigval, den.D) + + PeriodicDerivativeOperatorQuotient(num, den.D, (zero(T), one(T)), den.num_coef, tmp, num_eigval, den_eigval, rfft_plan, irfft_plan) +end + +function Base.://(num::PeriodicRationalDerivativeOperator, den::PeriodicRationalDerivativeOperator) + T = eltype(num) + @argcheck T == eltype(den) ArgumentError + @argcheck grid(num) == grid(den) ArgumentError + @argcheck num.den_coef == (one(T),) ArgumentError + @argcheck den.den_coef == (one(T),) ArgumentError + + x = grid(num) + u = zero.(x) + rfft_plan = plan_rfft(u) + tmp = rfft_plan * u + irfft_plan = plan_irfft(tmp, length(grid(num))) + + num_eigval = similar(tmp) + _eigvals!(num_eigval, num.D) + + den_eigval = similar(tmp) + _eigvals!(den_eigval, den.D) + + PeriodicDerivativeOperatorQuotient(num.D, den.D, num.num_coef, den.num_coef, tmp, num_eigval, den_eigval, rfft_plan, irfft_plan) +end + +function Base.://(num::FourierDerivativeOperator, den::FourierPolynomialDerivativeOperator) + T = eltype(num) + @argcheck T == eltype(den) ArgumentError + @argcheck grid(num) == grid(den) ArgumentError + + @unpack tmp, rfft_plan = den.D1 + irfft_plan = plan_irfft(tmp, length(grid(num))) + + num_eigval = similar(tmp) + _eigvals!(num_eigval, num) + + den_eigval = similar(tmp) + _eigvals!(den_eigval, den.D1) + + PeriodicDerivativeOperatorQuotient(num, den.D1, (zero(T), one(T)), den.coef, tmp, num_eigval, den_eigval, rfft_plan, irfft_plan) +end + +function Base.://(num::FourierPolynomialDerivativeOperator, den::FourierPolynomialDerivativeOperator) + T = eltype(num) + @argcheck T == eltype(den) ArgumentError + @argcheck grid(num) == grid(den) ArgumentError + + @unpack tmp, rfft_plan = den.D1 + irfft_plan = plan_irfft(tmp, length(grid(num))) + + num_eigval = similar(tmp) + _eigvals!(num_eigval, num.D1) + + den_eigval = similar(tmp) + _eigvals!(den_eigval, den.D1) + + PeriodicDerivativeOperatorQuotient(num.D1, den.D1, num.coef, den.coef, tmp, num_eigval, den_eigval, rfft_plan, irfft_plan) +end + +function Base.://(num::PeriodicDerivativeOperator, den::FourierPolynomialDerivativeOperator) + T = eltype(num) + @argcheck T == eltype(den) ArgumentError + @argcheck grid(num) == grid(den) ArgumentError + + @unpack tmp, rfft_plan = den.D1 + irfft_plan = plan_irfft(tmp, length(grid(num))) + + num_eigval = similar(tmp) + _eigvals!(num_eigval, num) + + den_eigval = similar(tmp) + _eigvals!(den_eigval, den.D1) + + PeriodicDerivativeOperatorQuotient(num, den.D1, (zero(T), one(T)), den.coef, tmp, num_eigval, den_eigval, rfft_plan, irfft_plan) +end + +function Base.://(num::PeriodicRationalDerivativeOperator, den::FourierPolynomialDerivativeOperator) + T = eltype(num) + @argcheck T == eltype(den) ArgumentError + @argcheck grid(num) == grid(den) ArgumentError + @argcheck num.den_coef == (one(T),) ArgumentError + + @unpack tmp, rfft_plan = den.D1 + irfft_plan = plan_irfft(tmp, length(grid(num))) + + num_eigval = similar(tmp) + _eigvals!(num_eigval, num.D) + + den_eigval = similar(tmp) + _eigvals!(den_eigval, den.D1) + + PeriodicDerivativeOperatorQuotient(num.D, den.D1, num.num_coef, den.coef, tmp, num_eigval, den_eigval, rfft_plan, irfft_plan) +end + + +function mul!(dest::AbstractVector{T}, quot::PeriodicDerivativeOperatorQuotient, u::AbstractVector{T}) where {T} + @unpack num_D, num_coef, den_coef, num_coef_end, den_coef_end, tmp, num_eigval, den_eigval, rfft_plan, irfft_plan = quot + N, _ = size(num_D) + @boundscheck begin + @argcheck N == length(u) + @argcheck N == length(dest) + end + + mul!(tmp, rfft_plan, u) + @inbounds @simd for j in Base.OneTo(length(tmp)-1) + tmp[j] *= evalpoly(num_eigval[j], num_coef) / evalpoly(den_eigval[j], den_coef) + end + tmp[end] *= evalpoly(num_eigval[end], num_coef_end) / evalpoly(den_eigval[end], den_coef_end) + mul!(dest, irfft_plan, tmp) +end + + + + """ ConstantFilter(D::FourierDerivativeOperator, filter, TmpEltype=T) diff --git a/src/periodic_operators.jl b/src/periodic_operators.jl index a6fb78b5..69e58cdb 100644 --- a/src/periodic_operators.jl +++ b/src/periodic_operators.jl @@ -865,19 +865,13 @@ struct PeriodicRationalDerivativeOperator{T<:Real, Dtype<:PeriodicDerivativeOper end end -function PeriodicRationalDerivativeOperator(D::PeriodicDerivativeOperator{T}) where {T<:Real} - x = grid(D) - u = zero.(x) - rfft_plan = plan_rfft(u) - tmp = rfft_plan * u - irfft_plan = plan_irfft(tmp, length(x)) - eigval = similar(tmp) +function _eigvals!(eigval::Vector{Complex{T}}, D::PeriodicDerivativeOperator{T}) where {T<:Real} @unpack factor = D @unpack lower_coef, central_coef, upper_coef = D.coefficients - N = length(x) + N = length(grid(D)) - for idx in 1:length(eigval) + @inbounds for idx in 1:length(eigval) tmp_eigval = zero(T) for j in Base.OneTo(length(lower_coef)) tmp_eigval += factor * lower_coef[j] * (exp(2π * im * (idx-1) * (-j) / N)) @@ -889,6 +883,19 @@ function PeriodicRationalDerivativeOperator(D::PeriodicDerivativeOperator{T}) wh eigval[idx] = tmp_eigval end + eigval +end + +function PeriodicRationalDerivativeOperator(D::PeriodicDerivativeOperator{T}) where {T<:Real} + x = grid(D) + u = zero.(x) + rfft_plan = plan_rfft(u) + tmp = rfft_plan * u + irfft_plan = plan_irfft(tmp, length(x)) + + eigval = similar(tmp) + _eigvals!(eigval, D) + PeriodicRationalDerivativeOperator(D, (zero(T), one(T)), (one(T),), tmp, eigval, rfft_plan, irfft_plan) end @@ -914,14 +921,6 @@ function Base.show(io::IO, rat::PeriodicRationalDerivativeOperator) print(io, rat.D) end -function LinearAlgebra.issymmetric(rat::PeriodicRationalDerivativeOperator) - @unpack num_coef, den_coef = rat - num_is_even = all(iszero, num_coef[idx] for idx in eachindex(num_coef) if iseven(idx)) - den_is_even = all(iszero, den_coef[idx] for idx in eachindex(den_coef) if iseven(idx)) - - num_is_even == den_is_even -end - function mul_poly(coef1, coef2) T = promote_type(eltype(coef1), eltype(coef2)) @@ -1123,18 +1122,6 @@ function Base.:/(rat1::PeriodicDerivativeOperator, rat2::PeriodicRationalDerivat PeriodicRationalDerivativeOperator(rat1) / rat2 end -function Base.://(rat1::PeriodicRationalDerivativeOperator, rat2::PeriodicRationalDerivativeOperator) - rat1 / rat2 -end - -function Base.://(rat1::PeriodicRationalDerivativeOperator, rat2::PeriodicDerivativeOperator) - rat1 / rat2 -end - -function Base.://(rat1::PeriodicDerivativeOperator, rat2::PeriodicRationalDerivativeOperator) - rat1 / rat2 -end - function mul!(dest::AbstractVector{T}, rat::PeriodicRationalDerivativeOperator, u::AbstractVector{T}) where {T} @unpack D, num_coef, den_coef, tmp, eigval, rfft_plan, irfft_plan = rat diff --git a/test/periodic_operators_test.jl b/test/periodic_operators_test.jl index c203533e..e9e32eaa 100644 --- a/test/periodic_operators_test.jl +++ b/test/periodic_operators_test.jl @@ -801,8 +801,8 @@ end for T in (Float32, Float64) xmin = -one(T) xmax = one(T) - for N in (8, 9), derivative_order in (1, 2), accuracy_order in (2, 3, 4) - D = periodic_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N) + for N in (8, 9), der_order in (1, 2), acc_order in (2, 3, 4) + D = periodic_derivative_operator(der_order, acc_order, xmin, xmax, N) x = grid(D) u = @. sinpi(x) - cospi(x)^2 + exp(sinpi(x)) println(devnull, D) @@ -810,9 +810,9 @@ for T in (Float32, Float64) @test Matrix(D^2) ≈ Matrix(D)^2 @test Matrix(D^2) ≈ Matrix(D * D) ≈ Matrix((I * D) * D) ≈ Matrix(D * (D * I)) - @test issymmetric(I - D^2) - @test !issymmetric(I + D) - @test !issymmetric(D - I) + @test Matrix(I - D^2) ≈ I - Matrix(D^2) + @test Matrix(I + D) ≈ I + Matrix(D) + @test Matrix(D - I) ≈ Matrix(D) - I rat = (I + 2D + 5*D^2) * (2I * D - D^3 * 5I) * (D*2 - D^2 * 5) @test rat.num_coef == (0.0, 0.0, 4.0, -2.0, -10.0, -45.0, 0.0, 125.0) @@ -841,4 +841,28 @@ for T in (Float32, Float64) @test integrate(u, D) ≈ sum(mass_matrix(D) * u) @test integrate(u->u^2, u, D) ≈ dot(u, mass_matrix(D), u) end + + for N in (8, 9), acc_order in (2, 3, 4) + D1a = periodic_derivative_operator(1, acc_order, xmin, xmax, N+1) + D2a = D1a^2 + D1b = fourier_derivative_operator(xmin, xmax, N) + D2b = D1b^2 + D2c = periodic_derivative_operator(2, acc_order, xmin, xmax, N+1) + + @test Matrix(D1a // (I - D2a)) ≈ Matrix(D1a) / (I - Matrix(D2a)) + @test Matrix(D1b // (I - D2a)) ≈ Matrix(D1b) / (I - Matrix(D2a)) + @test Matrix(D1a // (I - D2b)) ≈ Matrix(D1a) / (I - Matrix(D2b)) + @test Matrix(D1b // (I - D2b)) ≈ Matrix(D1b) / (I - Matrix(D2b)) + @test Matrix(D1b // (I - D2b)) ≈ Matrix(D1b / (I - D2b)) + @test Matrix(D1a // (I - D2c)) ≈ Matrix(D1a) / (I - Matrix(D2c)) + @test Matrix(D1b // (I - D2c)) ≈ Matrix(D1b) / (I - Matrix(D2c)) + + @test Matrix((I + D1a) // (I - D2a)) ≈ (I + Matrix(D1a)) / (I - Matrix(D2a)) + @test Matrix((I + D1b) // (I - D2a)) ≈ (I + Matrix(D1b)) / (I - Matrix(D2a)) + @test Matrix((I + D1a) // (I - D2b)) ≈ (I + Matrix(D1a)) / (I - Matrix(D2b)) + @test Matrix((I + D1b) // (I - D2b)) ≈ (I + Matrix(D1b)) / (I - Matrix(D2b)) + @test Matrix((I + D1b) // (I - D2b)) ≈ Matrix((I + D1b) / (I - D2b)) + @test Matrix((I + D1a) // (I - D2c)) ≈ (I + Matrix(D1a)) / (I - Matrix(D2c)) + @test Matrix((I + D1b) // (I - D2c)) ≈ (I + Matrix(D1b)) / (I - Matrix(D2c)) + end end