boutproject / BOUT-dev

BOUT++: Plasma fluid finite-difference simulation code in curvilinear coordinate systems
http://boutproject.github.io/
GNU Lesser General Public License v3.0
184 stars 95 forks source link

MMS testing W3 derivative #1049

Open johnomotani opened 6 years ago

johnomotani commented 6 years ago

The W3 scheme for a centred derivative https://github.com/boutproject/BOUT-dev/blob/8389878169b4e4824c2d486dad7d9d65fa0db8c5/src/mesh/index_derivs.cxx#L171-L210 does not converge at 3rd order for the test in tests/MMS/derivatives3, even though the solution is smooth. At 64->128 grid points the apparent order of convergence is ~2.3 while at high resolution (256->512) it is ~7. Increasing resolution further the order remains above 3 until 2048 before the error reaches a minimum at 4096 where presumably it hits double-precision round-off error.

Is there a better way to test whether W3 is performing correctly?

dschwoerer commented 6 years ago

tests/MMS/advection/ for weno3 also checks for 2nd order convergence ...

dschwoerer commented 6 years ago

Some plots of the errors The function tested was 2+sin(x) although the offset did not help. The W3 seems strange, as the behaviour changes as we approach w3 128 w3_128 w3 256 w3_256 w3 512 w3_512 w2 128 w2_128 w2 256 w2_256 w2 1024 w2_1024 c4 1024 (D2D2) c4_1024

They seem to have issues if the derivative is flat. I changed WENO_SMALL to 1e-4, and got this result with L_inf:

DD: 2+sin(%s) - W2 - x - CENTRE -> CENTRE  is not working. Expected 2.000000 got 3.246411 (64 -> 128)
DD: 2+sin(%s) - W2 - x - CENTRE -> CENTRE  is not working. Expected 2.000000 got 3.175009 (128 -> 256)
DD: 2+sin(%s) - W2 - x - CENTRE -> CENTRE  is not working. Expected 2.000000 got 3.034516 (256 -> 512)
DD: 2+sin(%s) - W2 - x - CENTRE -> CENTRE  is not working. Expected 2.000000 got 2.713054 (512 -> 1024)
DD: 2+sin(%s) - W3 - x - CENTRE -> CENTRE  is not working. Expected 3.000000 got 6.335953 (64 -> 128)
DD: 2+sin(%s) - W3 - x - CENTRE -> CENTRE  is not working. Expected 3.000000 got 4.435499 (128 -> 256)
DD: 2+sin(%s) - W3 - x - CENTRE -> CENTRE  is not working. Expected 3.000000 got 4.032094 (256 -> 512)
DD: 2+sin(%s) - W3 - x - CENTRE -> CENTRE  is not working. Expected 3.000000 got 4.000035 (512 -> 1024)

and with L_2:

DD: 2+sin(%s) - W2 - x - CENTRE -> CENTRE  is not working. Expected 2.000000 got 2.734049 (64 -> 128)
DD: 2+sin(%s) - W2 - x - CENTRE -> CENTRE  is not working. Expected 2.000000 got 2.685778 (128 -> 256)
DD: 2+sin(%s) - W2 - x - CENTRE -> CENTRE  is not working. Expected 2.000000 got 2.720727 (256 -> 512)
DD: 2+sin(%s) - W2 - x - CENTRE -> CENTRE  is not working. Expected 2.000000 got 2.588341 (512 -> 1024)
DD: 2+sin(%s) - W3 - x - CENTRE -> CENTRE  is not working. Expected 3.000000 got 5.959308 (64 -> 128)
DD: 2+sin(%s) - W3 - x - CENTRE -> CENTRE  is not working. Expected 3.000000 got 3.994329 (128 -> 256)
DD: 2+sin(%s) - W3 - x - CENTRE -> CENTRE  is not working. Expected 3.000000 got 3.999927 (256 -> 512)
DD: 2+sin(%s) - W3 - x - CENTRE -> CENTRE  is not working. Expected 3.000000 got 3.999984 (512 -> 1024)

Where the convergence is 4, the error is again a smooth sin/cos, as expected for a smooth scheeme. Smaller values have higher periodicy contributions: w2_256_-4

ZedThree commented 6 years ago

A brief update on this -- the solution for f in the MMS test:

[f]
solution = sin(t)*sin(3.0*x + 2.0*z) + cos(4.0*x^2 + z)

is badly aliased for nx < 128, so I'm not sure it's much of a surprise that there's a problem around there. However, I'm still not convinced it is working as expected, but that's partly because I'm not sure we know exactly how it (or any of the flux limiting schemes) should behave for any given input. I feel like if the "order" of the scheme is allowed to change through a scan in resolution, the order may not be well defined.

ZedThree commented 5 years ago

See #1499 for more discussion -- possibly some problem with the boundaries. Maybe make them periodic and test boundaries separately?

dschwoerer commented 5 years ago

Note that we don't need boundaries here. The input-function is purely the analytical one. The derivative is calculated once, and compared to the expected solution, excluding boundaries. Thus the discussion about boundaries does not help here.