Skip to content
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 Laplace ∇⁻²! via kernel in ∇²! #142

Merged
merged 6 commits into from
Sep 14, 2022
Merged

Inverse Laplace ∇⁻²! via kernel in ∇²! #142

merged 6 commits into from
Sep 14, 2022

Conversation

milankl
Copy link
Member

@milankl milankl commented Sep 9, 2022

fixes #141

@milankl milankl added the transform ⬅️ ➡️ Our spherical harmonic transform, grid to spectral, spectral to grid label Sep 9, 2022
@milankl milankl mentioned this pull request Sep 9, 2022
Copy link
Member

@maximilian-gelbrecht maximilian-gelbrecht left a comment

Choose a reason for hiding this comment

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

Looks good to me!

@@ -182,8 +181,6 @@ function SpectralTransform( ::Type{NF}, # Number format NF

# EIGENVALUES (on unit sphere, hence 1/radius²-scaling is omitted)
eigenvalues = [-l*(l+1) for l in 0:lmax+1]
eigenvalues⁻¹ = inv.(eigenvalues)
eigenvalues⁻¹[1] = 0 # set the integration constant to 0
Copy link
Member Author

Choose a reason for hiding this comment

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

I just realised that we shouldn't replace eigenvalues⁻¹ by 1/eigenvalues as that makes the integration constant infinity. We currently don't test for this (and it's not used by the dynamical core atm) hence the tests still pass

julia> alms = LowerTriangularMatrix(randn(ComplexF32,33,32));
julia> alms[:,1] .= real.(alms[:,1]);
julia> ∇⁻²alms = copy(alms)
julia> S = SpectralTransform(alms);
julia> SpeedyWeather.∇²!(∇⁻²alms,alms,S,inverse=true)
33×32 LowerTriangularMatrix{ComplexF32}:
   -Inf+NaN*im                0.0+0.0im                    0.0+0.0im
    -0.369622-0.0im      0.199208-0.408254im                0.0+0.0im
...

@maximilian-gelbrecht
Copy link
Member

I am on vacation for a week. I can look into it when I am back, if you don’t do it until then.

@milankl
Copy link
Member Author

milankl commented Sep 12, 2022

Don't worry, I'll resolve that, just don't merge it yet. Enjoy your vacation 🌴

@milankl milankl added this to the v0.4 milestone Sep 12, 2022
@milankl milankl changed the title inverse Laplace as operator via kernel in _∇²! Inverse Laplace ∇⁻²! via kernel in ∇²! Sep 14, 2022
@milankl
Copy link
Member Author

milankl commented Sep 14, 2022

Also @maximilian-gelbrecht, a little change I made in ∇²! and ∇⁻²! is that previously it excluded the the last row = highest degree l. This is because scalar quantities have lmax==mmax, but vector quantities use one more degree due to the recursion in the meridional gradient. But all LowerTriangularMatrices are of size lmax+1,mmax (as if they were all vector quantities). I changed ∇²! and ∇⁻²! to include all degrees and orders, to have a generally valid (inverse) Laplace operator (the higest degree l otherwise wouldn't be touched, and values may carry over which I don't like). I believe this does not have any unforeseen consequences but in case it does, let me know.

@milankl milankl merged commit 52c96a0 into main Sep 14, 2022
@milankl milankl deleted the mk/laplace branch September 14, 2022 13:53
@maximilian-gelbrecht
Copy link
Member

Great, thanks @milankl

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
transform ⬅️ ➡️ Our spherical harmonic transform, grid to spectral, spectral to grid
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Inverse Laplacian
2 participants