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

Changes redi tapering and slopes #65

Closed

Conversation

vanroekel
Copy link
Collaborator

In current MPAS slopes are altered directly, e.g. when S > Scrit it is set to zero. Griffies et al 1998 (Appendix C) equates this to slope clipping. The appropriate setup is to taper the redi kappa coefficient for each individual flux instead. The cross terms ("term2" and "term 3") end up being identical as they are d/dz(kappa S grad phi) so tapering kappa or S are identical. However for the vertical it is d/dz(kappa |S|^2 d phi / dz) so tapering S is different from tapering kappa. Also term 1 (laplacian diffusion) was not being tapered with slopes creating a cross isopycnal flux in steeply sloping regions. These issues are fixed here.

In addition a new slope calculation is introduced to prevent the slopes from getting very large when grad phi or d/dz(phi) get small. This function limits the slope to be between -1 and 1.

This changes tapering on the redi terms to be fully on the slopes
themselves and fixes a bug where k33 slopes are tapered by the square
and not linearly.  A 'safe slope' calculation is also introduced to
bound slopes from -1 to 1.
@vanroekel
Copy link
Collaborator Author

@jonbob @mark-petersen @maltrud @alicebarthel @golaz -- This is the first of two PRs I mentioned on the coupled call today. I have yet to test the build of the splitting and will report back once I have but wanted to get this posted for initial examination.

Note that all my testing of this feature includes #53

@vanroekel
Copy link
Collaborator Author

@dengwirda just a note that I haven't included your alternate slope calculation in this PR. Happy to have that discussion here and involve @mark-petersen as to whether or not we should make that change.

@vanroekel
Copy link
Collaborator Author

This should be ready for testing, I was able to build the code.

the redi n2 based limiter options are removed from the namelist and the
term1 limiter name is changed to reference tapering not just the n2
based tapering for clarity
@mark-petersen
Copy link

This does not compile standalone with OPENMP=true on chrysalis with gnu. With OPENMP=false it compiles with both gnu and intel on chrysals.

@dengwirda
Copy link
Collaborator

@vanroekel @mark-petersen @maltrud @jonbob @alicebarthel and all, to follow-up on the slope calc. aspects here, there currently appears to be a (seemingly unreasonable) dependence on the slope limit (config_Redi_maximum_slope) in the parabolic bowl test case --- this problem should setup simple linearaly sloping isopycnals via an uncontroversial stratification, and the slope limit should not be necessary to maintain stability. It should therefore be possible to set a larger value of config_Redi_maximum_slope than the actual slope defined in the test case and obtain consistent output. (e.g. if the defined isopycnal slope is 0.005, setting config_Redi_maximum_slope = 0.01, 0.05, 0.10, 0.20, etc should give the same answers).

Currently though, the slope limit appears to immediately become active, and both it and the monotone limiter in #53 are needed to keep the simulation on the rails. To me this suggests revisiting the slope triad formulations re: instabilities.

There are two things that I've tried:

  • The safe_slope function included in this PR, which bounds the slopes onto [-1,+1] and prevents div-by-zero errors while removing the need for a small fp-tolerance, which in itself seems to be responsible weird behaviour around 0.
  • Evaluating the slope triads at the top and bottom layer interfaces rather than the layer midplanes, as per: https://github.com/dengwirda/E3SM/blob/6c9aa01fdd11981f699b0d48c54665389cef9d66/components/mpas-ocean/src/shared/mpas_ocn_gm.F#L332. This averaging does appear to significantly improve the instability in the parabolic bowl case (allowing to run stably with a larger config_Redi_maximum_slope than the isopycnals themselves). It's unclear whether this change is necessarily compatible with the slope triad definitions themselves though.

Overall, this behaviour makes me unconfident re: the slope formulation in general, and suggests a further look.

@vanroekel
Copy link
Collaborator Author

@dengwirda The specification of a isopycnal slope and still having a dependence on max slope was indeed troubling. As I mentioned to you, Griffies et al 98 suggests the triad formulation as implemented is higher order accurate, but offer no proof of that claim. The formulation in that paper does not make much intuitive sense either (at least to me) as we are combining a vertical gradient and horizontal gradient that exist in different locations.

@dengwirda have you ever run a full G-case with your option 2 above? I had not yet, but it may be worth it. Your improvements in the bowl case are suggesting a further look at that approach I think.

@dengwirda
Copy link
Collaborator

@vanroekel I did run a few years of a G-case with these changes, and didn't see much change in AMOC, but that run didn't have your tapering updates, and also included a few other K_33 options that we've been going back forth on. Overall it could be good to re-test with a clean branch I think.

For me, understanding if / why there's instability in the triad slope formulation just in the parabolic bowl case would be good to get to the bottom of perhaps as a first step, since that should be a simple enough config. that's currently not behaving as expected. My colocated iso-slope branch seems to iron-out some of these issues, but as you say understanding whether this is consistent with Griffies et al re: triads is unclear...

@vanroekel
Copy link
Collaborator Author

@dengwirda just a quick update on your comment about slope sensitivity. I've been having a closer look and running some parallel tests. It has become pretty clear that the problem is at bathymetry. If you look at how the solutions diverge for max_slope = 0.1 and max_slope = 0.5 -- the differences emerge right along the bottom. It appears. you get a difference right on the boundary and then an opposite sign difference one layer above. e.g. --
image
These diffs then propagate throughout the domain and just get bigger.

I'll go back and have another look at the bottom boundary logic and update here if I learn anything.

Copy link

@mark-petersen mark-petersen left a comment

Choose a reason for hiding this comment

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

Based on the above discussion and testing, this is now ready to move to E3SM-Project/E3SM repo.

@mark-petersen
Copy link

Moved to E3SM-Project#5947

@xylar
Copy link
Collaborator

xylar commented Feb 3, 2024

@mark-petersen, can we delete this branch on Ocean-Discussion?

@xylar xylar deleted the vanroekel/ocean/change-redi-slopes-and-limiter branch February 3, 2024 13:45
@xylar
Copy link
Collaborator

xylar commented Feb 3, 2024

I'm going to go ahead and delete it...

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants