-
Notifications
You must be signed in to change notification settings - Fork 25
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Inverse Laplacian #141
Comments
Yes, we can make that functional again. I don't think it's necessary to use precalculated eigenvalues, as they are so simple, so I'd probably drop
|
I know the documentation is far from up to date, but at least that is already described https://milankl.github.io/SpeedyWeather.jl/dev/spectral_transform/#U,V-from-vorticity-and-divergence |
it should also work to just extend function ∇²!( ∇²alms::LowerTriangularMatrix{Complex{NF}}, # Output: Laplacian of alms
alms::LowerTriangularMatrix{Complex{NF}}, # Input: spectral coefficients
S::SpectralTransform{NF}; # precomputed eigenvalues
add::Bool=false, # add to output array or overwrite
flipsign::Bool=false, # -∇² or ∇²
inverse::Bool=false # ∇⁻² or ∇²
) where {NF<:AbstractFloat}
@boundscheck size(alms) == size(∇²alms) || throw(BoundsError)
lmax,mmax = size(alms) .- (2,1) # degree l, order m of the Legendre polynomials
@unpack eigenvalues = S
eigenvalues = inverse ? inv.(eigenvalues) : eigenvalues
@boundscheck length(eigenvalues) >= lmax+1 || throw(BoundsError)
kernel(o,a) = flipsign ? (add ? (o-a) : -a) :
(add ? (o+a) : a)
lm = 0
@inbounds for m in 1:mmax+1 # order m = 0:mmax but 1-based
for l in m:lmax+1 # degree l = m:lmax but 1-based
lm += 1
∇²alms[lm] = kernel(∇²alms[lm],alms[lm]*eigenvalues[l])
end
lm += 1
end
return ∇²alms
end I didn't test it, though |
Yes, good idea, I didn't realise I had already written ∇²! with a kernel. With #142 I decided to externalise the kernel such that we can also define ∇⁻² which might be nice to avoid confusion? And I didn't like the |
Oh we could also leave the kernel inside and define |
Internalised kernel is committed now, interestingly for 3x ?-clauses in a function you need to @inline kernel(o,a,b) = inverse ? (flipsign ? (add ? (o - a/b) : -a/b) :
(add ? (o + a/b) : a/b)) :
(flipsign ? (add ? (o - a*b) : -a*b) :
(add ? (o + a*b) : a*b)) |
For the record, writing a kernel like this # use eigenvalues⁻¹/eigenvalues for ∇⁻²/∇² based but name both eigenvalues
eigenvalues = inverse ? S.eigenvalues⁻¹ : S.eigenvalues
@boundscheck length(eigenvalues) >= lmax+1 || throw(BoundsError)
@inline kernel(o,a) = flipsign ? (add ? (o-a) : -a) :
(add ? (o+a) : a)
lm = 0
@inbounds for m in 1:mmax+1 # order m = 0:mmax but 1-based
for l in m:lmax+1 # degree l = m:lmax but 1-based
lm += 1
∇²alms[lm] = kernel(∇²alms[lm],alms[lm]*eigenvalues[l])
end
end Gives the same speed as hardcoding |
I’d need the inverse laplacian to convert vorticity to streamfunction. I’ve seen that this is available in the code, but commented out in
spectral_gradients!
as∇⁻²!
. Why is that so? Could we make it functional again?SpectralTransform
still has theeigenvalues⁻¹
field that I also only used by∇⁻²!
UV_from_vor!
also somehow inverts the laplacian but I don't understand how. Anyway though, I need to compute the stream function not the wind field.The text was updated successfully, but these errors were encountered: