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

Poisson solver test: recovery of an analytic solution. #227

Merged
merged 2 commits into from
May 22, 2019

Conversation

ali-ramadhan
Copy link
Member

Tests whether the Poisson solver can recover an analytic solution that is a product of sines and cosines. It can!

Tests whether the Poisson solver can recover an analytic solution that 
is a product of sines and cosines.
@ali-ramadhan ali-ramadhan added numerics 🧮 So things don't blow up and boil the lobsters alive testing 🧪 Tests get priority in case of emergency evacuation labels May 22, 2019
@ali-ramadhan ali-ramadhan added this to the v1.0 milestone May 22, 2019
@ali-ramadhan ali-ramadhan self-assigned this May 22, 2019
@ali-ramadhan ali-ramadhan changed the title Poisson analytic test Poisson solver test: recovery of an analytic solution. May 22, 2019
@codecov
Copy link

codecov bot commented May 22, 2019

Codecov Report

Merging #227 into master will not change coverage.
The diff coverage is n/a.

Impacted file tree graph

@@           Coverage Diff           @@
##           master     #227   +/-   ##
=======================================
  Coverage   66.51%   66.51%           
=======================================
  Files          18       18           
  Lines         675      675           
=======================================
  Hits          449      449           
  Misses        226      226
Impacted Files Coverage Δ
src/Oceananigans.jl 66.66% <ø> (ø) ⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 6c0f252...f34afe8. Read the comment docs.

1 similar comment
@codecov
Copy link

codecov bot commented May 22, 2019

Codecov Report

Merging #227 into master will not change coverage.
The diff coverage is n/a.

Impacted file tree graph

@@           Coverage Diff           @@
##           master     #227   +/-   ##
=======================================
  Coverage   66.51%   66.51%           
=======================================
  Files          18       18           
  Lines         675      675           
=======================================
  Hits          449      449           
  Misses        226      226
Impacted Files Coverage Δ
src/Oceananigans.jl 66.66% <ø> (ø) ⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 6c0f252...f34afe8. Read the comment docs.

@ali-ramadhan ali-ramadhan merged commit 4254c69 into master May 22, 2019
Copy link
Member

@glwagner glwagner 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, main issue is choosing rtol. Would actually be nice if we could set rtol = C / N^2 (since we believe the solver is 2nd order accurate) --- but might be asking too much for right now.

by giving it the source term or right hand side (RHS), which is ``f(x, y, z) = \\nabla^2 \\Psi(x, y, z) =
-((\\pi m_z / L_z)^2 + (2\\pi m_y / L_y)^2 + (2\\pi m_x/L_x)^2) \\Psi(x, y, z)``.
"""
function poisson_ppn_recover_sine_cosine_solution(ft, Nx, Ny, Nz, Lx, Ly, Lz, mx, my, mz)
Copy link
Member

Choose a reason for hiding this comment

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

This is just semantics, but I don't like that the ordinary action of the Poisson solver is called 'recovery'... perhaps just call this poisson_solver_test_sinusoid or something?

Copy link
Member Author

Choose a reason for hiding this comment

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

That's definitely a better name. Will update.


error = norm(ϕ.data - Ψ.(xC, yC, zC)) / √(Nx*Ny*Nz)

@info "Error (ℓ²-norm), $ft, N=($Nx, $Ny, $Nz), m=($mx, $my, $mz): $error"
Copy link
Member

Choose a reason for hiding this comment

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

I applaud \ell.


@info "Error (ℓ²-norm), $ft, N=($Nx, $Ny, $Nz), m=($mx, $my, $mz): $error"

isapprox(real.(ϕ.data), Ψ.(xC, yC, zC); rtol=5e-2)
Copy link
Member

Choose a reason for hiding this comment

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

What determines rtol here? I think it's at least necessary to leave a comment on how this was decided. Honestly 5e-2 seems quite large. Also since the error depends on both the resolution and mx, my, and mz, it might be a good idea to hard code the resolution and the wavenumber into the function rather than allow them to be changed in the function signature, where they could actually break the code even if the solver still works as expected.

Copy link
Member Author

Choose a reason for hiding this comment

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

For sure. It's the rtol of the worst scenario (low resolution, high wavenumbers) but that's not a very good metric... Your C / N^2 suggestion would be much better.

I guess it's okay to be verbose in the tests so we could have the wave numbers hardcoded in.

@testset "Analytic solution reconstruction" begin
println(" Testing analytic solution reconstruction...")
for N in [32, 48, 64], m in [1, 2, 3]
@test poisson_ppn_recover_sine_cosine_solution(Float64, N, N, N, 100, 100, 100, m, m, m)
Copy link
Member

Choose a reason for hiding this comment

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

Any reason for not using just 1, 1, 1 for the length?

Does it make sense to provide multiple error tolerances if we are testing more than one resolution? Seems like we don't get much out of testing the lower resolutions here: if the higher resolutions pass, the lower ones are basically guaranteed to.

Copy link
Member

Choose a reason for hiding this comment

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

I take back my suggestion for 1, 1, 1 --- you should definitely use 2π, 2π, π for simplicity!

Copy link
Member Author

Choose a reason for hiding this comment

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

Not really, (1, 1, 1) makes more sense. I just usually default to 100, 100, 100 I guess since most of the simulations I've done were on that domain lol.

I agree it would be good to have multiple error tolerances.

The reason I did multiple resolutions and logged the error is just to see if the error decreases as expected. I'm going to start doing this more often so that we get some useful debug/diagnostic information from the test suite.

@ali-ramadhan
Copy link
Member Author

ali-ramadhan commented May 22, 2019

Ah sorry I merged too quickly! Lots of easy improvements we can make though!

I can make another PR with the improvements.

@glwagner
Copy link
Member

Yes, totally fine for this to be TODO. We have plenty on our plate.

@ali-ramadhan
Copy link
Member Author

Looks good, main issue is choosing rtol. Would actually be nice if we could set rtol = C / N^2 (since we believe the solver is 2nd order accurate) --- but might be asking too much for right now.

What's C? I assume N = Nx*Ny*Nz.

@glwagner
Copy link
Member

C is a constant that you'd have to determine.

@ali-ramadhan
Copy link
Member Author

Ah ok, just a bit of trial and error then.

@glwagner
Copy link
Member

glwagner commented May 23, 2019

The idea is just that in the 'asymptotic limit' the solver error should decrease with N^2, where N = Nx = Ny = Nz, more or less --- with a constant of proportionality C that is not known, in general. We should actually demonstrate this at some point --- it may not be easy to incorporate the idea into a test (I don't know what resolution we need to be in that asymptotic limit) but the log-log plot of solver error versus resolution is certainly a plot you want to make at some point.

@ali-ramadhan
Copy link
Member Author

Yes, it would be good to test that it converges quadratically as expected. If it's an expensive test that requires large solves we can throw it into a set of integration tests (#139) that runs less frequently.

@ali-ramadhan ali-ramadhan deleted the poisson-analytic-test branch May 24, 2019 23:47
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
numerics 🧮 So things don't blow up and boil the lobsters alive testing 🧪 Tests get priority in case of emergency evacuation
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants