search

gridapapps / GridapGeosciences.jl

Gridap drivers for geoscience applications
MIT License
35 stars 3 forks source link

Shallow water explicit #15

Closed davelee2804 closed 3 years ago

davelee2804 commented 3 years ago

This pull request is to add a second order explicit shallow water integrator and associated Williamson2 analytic test case into GridapGeosciences.

It includes additions to the src/ directory as follows: 1) additional diagnostics for the shallow water equations (mass, vorticity, energy, power) 2) a reference domain grad_perp operator (low level) 3) diagnostic equation to compute the vorticity (low level) 4) diagnostic equation to compute the potential vorticity (mixed low level/high level) 5) explicit integrator (nonlinearly stable, exactly balances energy exchanges in space and time) 6) integration driver function, takes in an integrator as a function pointer and cycles over this for a specified number of time steps

It also includes an analytic test (williamson2), run at multiple spatial and temporal resolutions in order to determine the L2 errors in the velocity and the fluid depth. L2 errors at 3 resolutions using RT1 and supplied in order to test against and verify second order accuracy

codecov-commenter commented 3 years ago

Codecov Report

Merging #15 (8dc7a4c) into master (68281e2) will increase coverage by 17.99%. The diff coverage is 96.12%.

Impacted file tree graph

@@             Coverage Diff             @@
##           master      #15       +/-   ##
===========================================
+ Coverage   44.72%   62.72%   +17.99%     
===========================================
  Files          10       11        +1     
  Lines         313      440      +127     
===========================================
+ Hits          140      276      +136     
+ Misses        173      164        -9
Impacted Files Coverage Δ
src/ShallowWaterExplicit.jl 95.76% <95.76%> (ø)
src/DiagnosticTools.jl 58.33% <100.00%> (+50.64%) :arrow_up:
src/CoordinateTransformations.jl 76.19% <0.00%> (+21.64%) :arrow_up:

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 68281e2...8dc7a4c. Read the comment docs.

davelee2804 commented 3 years ago

Hi @amartinhuertas , I've tried to address your helpful comments as best I can. The code looks a lot better now thanks! (results are the same as before for the test, I'm kinda surprised that there is no difference in the L2 errors given the solve() has been replaced with the ldiv!() ...)

The only thing I have not addressed (of which I'm aware) is the removal of the grad_perp_ref_domain(), and the rest of the routines in the Weak2DCurlOps.jl file. If you feel strongly that these should be removed and replaced with high level physical space terms (and an AD Jacobian) I will do so. My only concern is that without a formal proof that the two are equivalent vorticity will not be conserved and the benefits of compatible FEM will not be gaurunteed.

Let me know what you advise, and thanks again for your review.

amartinhuertas commented 3 years ago

Hi @davelee2804,

Hi @amartinhuertas , I've tried to address your helpful comments as best I can. The code looks a lot better now thanks! (results are the same as before for the test, I'm kinda surprised that there is no difference in the L2 errors given the solve() has been replaced with the ldiv!() ...)

Cool. Thanks for your work! The code looks pretty! I have marked as resolved all those comments that were already addressed.

were the following comments addressed?

In general, in PRs, if we decide not to go ahead with a recommendation, we also typically write some words on why not on the corresponing comment text box. The code reviewer might not have the whole picture of the work done, and may provide a missleading recommendation due to that. Like the rebuttal letter to the paper referees :-)

The only thing I have not addressed (of which I'm aware) is the removal of the grad_perp_ref_domain(), and the rest of the routines in the Weak2DCurlOps.jl file. If you feel strongly that these should be removed and replaced with high level physical space terms (and an AD Jacobian) I will do so. My only concern is that without a formal proof that the two are equivalent vorticity will not be conserved and the benefits of compatible FEM will not be gaurunteed.

First of all, just in case it is not clear, we can either use isoparametric or analytical mapping independently of what we want to do with the right hand side of the weak div perp operator. Second, I would go for the high level API term. Numerically, they seem to be equivalent (up to a minus sign). There is no mathematical proof yet they are equivalent, that is right.

At this point, I you are still concerned, what I suggest you to do is the following.

  1. Let us create a branch of the current branch with the current implementation, shallow_water_explicit_grad_perp_ref_vorticity.
  2. Let us transform the current implementation in shallow_water_explicit using the high level API terms. You will see that the code will look even prettier, more general, and abstract.
  3. Let us compare systematically with a benchmark problem that you are particular worried about the results of the codes in branches 1 and 2. This would provide more confidence on their equivalence.
  4. If it turns out that there are differences among 2. and 3., let us think how can we abstract the grad_perp_ref_domain term using the high level API. Hopefully they turn to be equivalent and we just have to merge 2. into main.

What do you think?

davelee2804 commented 3 years ago

Hey @amartinhuertas , sounds like a good plan. I will make this secondary branch for the grad_perp_ref_domain, and then revert the existing branch to the high level physical space version as you suggest. Since we are computing vorticity conservation as a diagnostic we can monitor this for any bias with the physical space implementation and revisit if necessary.

Sorry for not commenting on the changes I had not made, I'll be more careful with this in future....

davelee2804 commented 3 years ago

Hi @amartinhuertas , I've made the grad_perp branch as you suggested, then removed the low level reference operators from the shallow_water_explicit branch and replaced these with the high level physical space terms as you suggested. The errors are similar but not the same. Here are the conservation and l2 errors for the final 10 time steps for both configurations using the C16 mesh and RT1 elements, first the high level physical space curl:

step | mass conservation | vorticity conservation | kinetic energy | potential energy | energy conservation | power
70  -1.9114360408403353e-15 5.960464477539063e-8    6.633793145418441e20    1.477262020091013e22    -6.469974972230993e-10  2.9468087973284225e9
71  -1.9114360408403353e-15 1.7881393432617188e-7   6.633793153306743e20    1.4772620200009785e22   -6.542217757846671e-10  5.013208853950043e9
72  -1.9114360408403353e-15 0.0 6.633793164808813e20    1.4772620198748801e22   -6.613983670972954e-10  7.013473197860233e9
73  -1.9114360408403353e-15 -1.7881393432617188e-7  6.63379317773161e20 1.4772620197346446e22   -6.685294449387136e-10  7.772140080217796e9
74  -1.9114360408403353e-15 0.0 6.633793188754459e20    1.4772620196134932e22   -6.756057707535717e-10  6.685436270421591e9
75  -1.9114360408403353e-15 5.960464477539063e-8    6.633793195407092e20    1.4772620195361106e22   -6.82638756874949e-10   4.243648649419777e9
76  -1.9114360408403353e-15 0.0 6.633793197822413e20    1.4772620195011748e22   -6.896241916084949e-10  1.8942620096129227e9
77  -1.9114360408403353e-15 -3.5762786865234375e-7  6.633793199121693e20    1.47726201947747e22 -6.965637052875065e-10  1.280609604765152e9
78  -2.1238178231559282e-15 -5.960464477539063e-8   6.633793204058822e20    1.4772620194174609e22   -7.034553958564704e-10  3.2850473646575394e9
79  -1.9114360408403353e-15 5.364418029785156e-7    6.633793216577788e20    1.477262019281706e22    -7.102999426209271e-10  7.441219225807159e9
80  -2.1238178231559282e-15 1.1920928955078125e-7   6.633793237602512e20    1.4772620190609624e22   -7.170997910808223e-10  1.2075524820265934e10
n=16,   err_u: 0.0007415233680146857,   err_h: 0.00036813020391681306

...and here are the low level reference space grad-perp conservation and l2 errors:

step | mass conservation | vorticity conservation | kinetic energy | potential energy | energy conservation | power
70  -1.9114360408403353e-15 0.0 6.633793145418448e20    1.4772620200910134e22   -6.46997089639775e-10   2.9468087698963394e9
71  -1.9114360408403353e-15 1.7881393432617188e-7   6.633793153306746e20    1.4772620200009787e22   -6.542216399235591e-10  5.013208836234894e9
72  -1.9114360408403353e-15 1.1920928955078125e-7   6.633793164808811e20    1.4772620198748807e22   -6.613978236528631e-10  7.0134731802047e9
73  -1.9114360408403353e-15 1.1920928955078125e-7   6.633793177731617e20    1.4772620197346446e22   -6.685294449387136e-10  7.772140071719622e9
74  -1.9114360408403353e-15 1.7881393432617188e-7   6.633793188754468e20    1.4772620196134934e22   -6.756053631702474e-10  6.685436286510313e9
75  -1.9114360408403353e-15 2.384185791015625e-7    6.633793195407095e20    1.4772620195361104e22   -6.826388927360571e-10  4.2436486332578735e9
76  -2.1238178231559282e-15 2.980232238769531e-7    6.633793197822415e20    1.4772620195011752e22   -6.896237840251706e-10  1.894261969870346e9
77  -1.9114360408403353e-15 2.384185791015625e-7    6.633793199121695e20    1.4772620194774694e22   -6.965641128708307e-10  1.2806096133585129e9
78  -2.336199605471521e-15  0.0 6.633793204058827e20    1.477262019417461e22    -7.034552599953623e-10  3.2850473504682846e9
79  -1.9114360408403353e-15 2.384185791015625e-7    6.633793216577798e20    1.477262019281705e22    -7.103002143431434e-10  7.441219187366409e9
80  -2.1238178231559282e-15 1.7881393432617188e-7   6.63379323760251e20 1.4772620190609633e22   -7.1709924763639e-10    1.2075524782441132e10
n=16,   err_u: 0.0007415233680147283,   err_h: 0.0003681302039168138

Since there is no obvious bias in the vorticity conservation errors, and they are the same magnitude in both cases, we are probably OK, though I will continue to monitor these vorticity conservation errors for more advanced tests...

Are there any more outstanding issues I should address?

Thanks again.

amartinhuertas commented 3 years ago

Hi @davelee2804,

nice work!

If I understand correctly the results in the tables (and as you say), it seems both quantitatively and qualitatively that both approaches achieve the desired outcome right? It is surprising that the vorticity becomes exactly 0.0 in some steps (happens in both cases).

If you still feel that a single experiment is still not enough, we may want to systematize experiments with the help of DrWatson.jl as well. I expect this work to be also necessary for the future, to compare systematically different algorithms and parameter values (mesh resolution, time step, FE order, geometrical order, etc.). I could help with this.

In regards to the PR, I will start now a new round of revision, with additional comments.

Thanks for your work!

amartinhuertas commented 3 years ago

Hi @davelee2804 ... I am going out for a walk. I will continue later adding more comments to my revision.

amartinhuertas commented 3 years ago

Comment out temporarily Manifest.toml from .gitignore and commit and push Manifest.toml/.gitignore. Note that we are temporarily pointing to a particular branch in Gridap.jl github repo. (while these devs are not released, currently under revision). Thus, this is needed to let the tests pass in Github CI actions. By the way, I would stop using [ci skip] for a while in this PR so that we double check that the tests pass in Github CI actions. We want to merge if and only if the tests are passing there.

amartinhuertas commented 3 years ago

@davelee2804 ... FYI ... tests are failing in this branch at Github Actions. I am positive it is related to this:

https://github.com/gridapapps/GridapGeosciences.jl/pull/15#issuecomment-913085099

Also take into account the modified signatures for assemble_vector! and assemble_matrix_and_vector! signatures.

davelee2804 commented 3 years ago

Hey @amartinhuertas , hopefully I've addressed all the issues you raised (except for the ones I have already commented on). The CI test is failing at: https://github.com/gridapapps/GridapGeosciences.jl/blob/shallow_water_explicit/src/ShallowWaterExplicit.jl#L27 but presumably this is because the CI testing is referencing the older version of assemble_vector!() with the different argument ordering. Please let me know if I've missed anything, and thanks again for the super helpful and educational review!

amartinhuertas commented 3 years ago

Hey @amartinhuertas , hopefully I've addressed all the issues you raised (except for the ones I have already commented on). The CI test is failing at: https://github.com/gridapapps/GridapGeosciences.jl/blob/shallow_water_explicit/src/ShallowWaterExplicit.jl#L27 but presumably this is because the CI testing is referencing the older version of assemble_vector!() with the different argument ordering. Please let me know if I've missed anything, and thanks again for the super helpful and educational review!

I will fix the failing test, no worries.

On the other hand, have you addressed this?

https://github.com/gridapapps/GridapGeosciences.jl/pull/15#discussion_r700965121

Finally, I would give a shot to splitting your long functions into subfunctions. I think having functions that only do one thing is essential. If you do not find sufficiently confident how to do it, I can do it. Let me know.

davelee2804 commented 3 years ago

Hey @amartinhuertas , I don't know how to split up the explicit function. I know it is long, and as a two stage function is make sense to split it. However each of the two stages is different and relies on a different amount of data (the first stage relies on a single time level of data, the second stage on two).

If you have any idea on how to split it I'm happy to try, but there is no obvious way to me (It is not like this is straight Runge-Kutta where you do the same thing on different data at each stage)

amartinhuertas commented 3 years ago

Hey @amartinhuertas , I don't know how to split up the explicit function. I know it is long, and as a two stage function is make sense to split it. However each of the two stages is different and relies on a different amount of data (the first stage relies on a single time level of data, the second stage on two). If you have any idea on how to split it I'm happy to try, but there is no obvious way to me (It is not like this is straight Runge-Kutta where you do the same thing on different data at each stage)

You can see the result of the work here: https://github.com/gridapapps/GridapGeosciences.jl/blob/60313734c8d369d57fa2803681e88cff492ef83d/src/ShallowWaterExplicit.jl#L68

Let me know if that makes sense.

Please note that the function has been split into many different reusable functions which only do one thing. I guess that these functions might be used for other solvers as well. Also note that it was impossible to reuse code from the previous function without replicating it.

davelee2804 commented 3 years ago

Thanks for all your help here @amartinhuertas , I've learnt a lot! I made one final change before merging, the L2 errors were different in the last couple of decimals (perhaps the order of operations changed on something). Nothing big enough to trigger a test failure (set at 1.0e-12), and perhaps we are not bit-4-bit on different machines anyway, but just for consistency.

davelee2804 commented 3 years ago

Hey @amartinhuertas , the merge just failed. Is this related to the commit you just made?

amartinhuertas commented 3 years ago

Hey @amartinhuertas , the merge just failed. Is this related to the commit you just made?

I do not see that it failed (still on the way). Can you send the link to Github Actions of the failed job?

amartinhuertas commented 3 years ago

I made one final change before merging, the L2 errors were different in the last couple of decimals (perhaps the order of operations changed on something). Nothing big enough to trigger a test failure (set at 1.0e-12), and perhaps we are not bit-4-bit on different machines anyway, but just for consistency.

FYI ... Julia offers the built-in isapprox function to facilitate the implementation of this kind of floating point equality tests.

davelee2804 commented 3 years ago

Hmm, I can't see it in the actions, and the error message is no longer here, I'm confused. I'll wait til your CI tests have run and try again...

davelee2804 commented 3 years ago

Merged now :)

amartinhuertas commented 3 years ago

Note that the order of the arguments has changed for some functions ... When you pull from master into the other brnaches this will affect you.

davelee2804 commented 3 years ago

OK, thanks I'll try and keep that in mind