Evaluate time dependent tendencies at correct RK4 stage time value

[sbrus89]

This PR corrects a longstanding issue with RK4 and tendency terms which have an explicit time dependence (see E3SM-Project#5364 and MPAS-Dev/MPAS-Model#375). Big thanks to @gcapodag and @jeremy-lilly for providing the fix in their LTS development branch linked in E3SM-Project#5364.

Here, I have migrated these changes from their MPAS-Dev branch to E3SM and added the modifications to correct the convergence issues for the manufactured solution test case: E3SM-Project/polaris#72

[sbrus89]

@gcapodag and @jeremy-lilly, I wasn't able to add you as reviewers, I'd appreciate your review here if you have time to make sure I accurately migrated the time varying forcing changes from your branch. Thanks for making these changes!

[cbegeman]

@sbrus89 By visual inspection it looks good. It seems like we need a bit of testing in cases that use tidal, land ice, and atmospheric forcing (any other kinds of forcing needed as well?). Let me know when and how you would like me to help with testing. Maybe after @gcapodag and @jeremy-lilly have a chance to look it over?

[mark-petersen]

@sbrus89 this looks great. Can you provide more detail here on what you changed, and the result that confirmed it is working?

[gcapodag]

Hi @sbrus89 I am planning to take a look by the end of this week. One thing that @jeremy-lilly and I noticed back then was that some intermediate dts caused weird behavior for the forcing time increment, for instance if they were irrational or had too many decimal digits. We did not have time to investigate the cause then but please do make sure that it is not still happening before you merge these changes.

[jeremy-lilly]

Hey @sbrus89, thanks for putting this together!

Above, @gcapodag mentions that

...some intermediate dts caused weird behavior for the forcing time increment, for instance if they were irrational or had too many decimal digits...

The details are a little fuzzy since this was a while ago, but here is what I remember: The problem should be in mpas_advance_forcing_clock(). This takes in a dt in seconds represented as a decimal number and converts it to a rational number S_n / S_d. Then, S_n and S_d are passed to mpas_set_timeInterval() which returns an appropriate timeStep to be passed to mpas_advance_clock().

This was all done because, for some reason I can't remember, mpas_set_timeInterval(timeStep, dt=dt) on line 1258 in the previous version returned a timeStep representing a whole number of seconds, effectively taking the floor of dt.

So, either there is a better way to do this entirely or the problem is almost definitely in how S_n and S_d are determined. A possible, cheap solution could be to increase the bound powerOfTen < 11 on line 1276 to something larger throwing an exception if it breaks the imposed bound.

[sbrus89]

@mark-petersen, this change adds the $0, h/2, h/2, h$ factors that are added to $t_n$ and used to evaluate the tendency terms at each RK stage:

$$ k_1 = f(t_n, y_n) $$

$$ k_2 = f\left(t_n + \frac{h}{2}, y_n + h\frac{k_1}{2}\right) $$

$$ k_3 = f\left(t_n +\frac{h}{2}, y_n + h\frac{k_2}{2}\right) $$

$$ k_4 = f\left(t_n + h, y_n + hk_3\right) $$

These had been missing since all previous tendencies didn't have explicit time dependencies. However, this change is required to correctly evaluate the time varying atmospheric/land ice forcing, tidal potential, and manufactured solution tendencies, which all explicitly depend on $t$.

[sbrus89]

With this change the manufactured solution test case converges at second order accuracy:

See E3SM-Project#5725 for the sub-optimal convergence result prior to this fix.

The Icos7 tides case in compass is essentially unchanged with the following constituent errors:

Constituent	Baseline	This PR
K1	0.03181129	0.03181336
M2	0.10802920	0.10859477
N2	0.02325495	0.02336509
O1	0.01966013	0.01967789
S2	0.05344151	0.05367388

I still need to test the hurricane case...

[sbrus89]

@gcapodag and @jeremy-lilly, thanks for the background on mpas_advance_forcing_clock(). I wasn't 100% clear on those changes at first, but that explains why they were necessary. I'll look into it and do some testing.

[gcapodag]

Hi @sbrus89 why is the order of convergence 2 and not 4?

[gcapodag]

It might be worth comparing to a numerical solution obtained on the same mesh to check if you actually get 4. I guess the 2 might be due to the spatial discretization error if you are using an actual analytical solution.

[sbrus89]

@gcapodag, yes the second order convergence is due to the spatial error since we're refining in space+time and comparing to an analytical solution in this case. I'll do a dt refinement to verify 4th order in time.

[gcapodag]

Thanks @sbrus89 !

[sbrus89]

@gcapodag, For the manufactured solution case, I ran the 200km mesh with a series of refined timesteps (20, 10, 5, 2.5 min) and compared the error against a 15 second reference solution. The results confirm fourth order convergence in time:

I'll plan on doing the same type of study for the hurricane case to confirm things are working as expected with time varying forcing. I believe you and Jeremy had done that previously for the LTS work, is that correct?

[gcapodag]

That is great, thanks @sbrus89 ! Yes, when we were working on our latest LTS paper (the one we are all in) @jeremy-lilly and I had the changes on the time varying forcing in the branch we used for the numerical results. I still need to look at your changes for this PR but back then we were modifying something called forcingTimeIncrement directly in the LTS time stepping routine (as well as in RK4). That LTS scheme is not in the E3SM master branch. In what is being pushed here E3SM-Project#5460, LTS has an operator splitting where the time varying forcing is only computed at the beginning of the time step so intermediate evaluations are not necessary.

[sbrus89]

@gcapodag and @jeremy-lilly -- I found a commit in the framework E3SM-Project@40f8ce7 that I think may address the whole-second representation that you were seeing with mpas_set_timeInterval(timeStep, dt=dt). It looks like this change wasn't in the branch you were working with at the time, which required the rational number solution. My hope is that this commit fixes that issue and that the changes in mpas_advance_forcing_clock() aren't required any longer.

[gcapodag]

I wonder if now that we have the capability to compute the forcing correctly, we should allow the possibility of not doing it. For instance, instead of having a condition on the type of integrator here, we could have something like config_update_forcing_at_intermediate_timesteps that is also used in the RK4 timestepping code and sets all the forcingTimeIncrementRK4 to zero if false. In this way, for coupled runs, one can decide to have a lighter coupling, but if one wants to do convergence tests and use RK4 as a reference solution to validate other time stepping schemes it is still possible turning that flag on.

[xylar]

Suggested change

	if (timeIntegratorChoice /= 4) then
	if (timeIntegratorChoice /= timeIntRK4) then

[gcapodag]

Hi @sbrus89, I looked at the commit carefully and it looks like you caught all the changes that @jeremy-lilly and I did in our branch. Thanks for integrating them in the E3SM code base. Good job also on finding a possible work around for the changes we added to mpas_advance_forcing_clock(). Once you give those a go, I would recommend to repeat the convergence test you did before including dts in seconds that have a few decimals digits. In that case it might be inconvenient to use a reference solution, so maybe just use three subsequent RK4 numerical solutions on the same mesh. Unfortunately, with RK4 the forcingTimeIncrements do not have divisions by 3 like SSPRK3 so we are spared a few nastier situations to test.

[gcapodag]

Actually, maybe instead of setting the forcingTimeIncrementRK4 to zero in terms of having a light coupling we should not advance the forcing clock here. Please check if this statement is true. Thanks!

[gcapodag]

Hi @sbrus89 what does the last commit do?

[xylar]

@sbrus89, I have a few suggestions related to an enum that should be made public and used rather than hard-coded integer values.

Other than that, I would be curious about the response to @gcapodag's suggestions to have a flag to allow for lighter weight coupling even with RK4. My feeling would be that we really don't intend to use RK4 for coupled production runs to we don't need it to be cheap. It is basically always going to be used for high-fidelity runs that we can treat as a benchmark.

[xylar]

Suggested change

	if (timeIntegratorChoice /= 4) then
	if (timeIntegratorChoice /= timeIntRK4) then

[xylar]

@sbrus89, once my small changes are addressed, I think this should migrate over to E3SM, @mark-petersen can review there and @cbegeman can re-review.

[cbegeman]

Approving based on testing with the manufactured solution test case prior to addressing Xylar's suggestions. I'll do additional testing when this PR is migrated to E3SM. Thanks for your work on this @sbrus89!

[sbrus89]

@xylar, I tend to agree with your take on the lighter RK4 coupling. In addition to RK4 not being performance-critical for production runs, I think the extra expense of doing the correct stage t evaluation is negligible anyway.

[mark-petersen]

@sbrus89 thanks for your work on this long-standing problem, and your convergence testing above. Approving based on the testing and discussions above.

[xylar]

@mark-petersen, just so you know, this has moved to E3SM already:
E3SM-Project#6695

I'm going to close this to avoid confusion.