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

replace global sum with full-precision reproducible sum #6456

Merged
merged 2 commits into from
Jul 10, 2024

Conversation

maltrud
Copy link
Contributor

@maltrud maltrud commented May 31, 2024

The global sum in the surface salinity restoring routine was forced to be BFB using a truncation calculation that resulted in reduced precision and some unexplained behavior. This PR replaces this calculation with an mpas standard reproducible sum with full precision.

Testing was done using a G-case with the IcoswISC30E3r5 grid for one year, then calculating the global average SSS restoring tendency term, which should be effectively zero. Also confirmed that the new code resulted in BFB solution when processor layout was changed from 10 nodes on chrysalis to 15 nodes.

Fixes #6446
[non-BFB] for C- and G-cases with salinity restoring

@maltrud maltrud added mpas-ocean non-BFB PR makes roundoff changes to answers. labels May 31, 2024
Copy link

github-actions bot commented May 31, 2024

PR Preview Action v1.4.7
🚀 Deployed preview to https://E3SM-Project.github.io/E3SM/pr-preview/pr-6456/
on branch gh-pages at 2024-06-20 21:03 UTC

@maltrud
Copy link
Contributor Author

maltrud commented May 31, 2024

Here's a time series of the log10 of the absolute value of the globally averaged SSS tendency (monthly averages) for 4 simulations.

Screenshot 2024-05-30 at 5 26 39 PM

The black line is the old global sum with real4 output, and the green is the old sum with real8 output. The red line is the new reproducible sum with r4 output, and the blue is the new sum with real8 output.

@@ -414,11 +422,14 @@ subroutine ocn_get_surfaceSalinityData( streamManager, &

endif

indexForReproSum(1) = 1
indexForReproSum(2) = 0
Copy link
Contributor

Choose a reason for hiding this comment

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

This is probably just do to the construction of the repro sum, but I'm a bit confused by the two length array of indexForReproSum. It seems position 1 is always just 1 and position 2 is the total number of entries in the sum. Is that right? What is position 1 of this index array used for? It never appears to be modified and then just get used in the globalSum call.

Copy link
Contributor

Choose a reason for hiding this comment

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

@vanroekel The reproducible sums in MPAS allow sums over full arrays and over restricted address ranges in each dimension. For full generality, the range must be specified in min/max pairs for each array dimension. In this case, the array is getting packed in advance so the min is always 1 but the max can vary. But both need to be provided for the interface.

Copy link
Contributor

@alicebarthel alicebarthel Jun 3, 2024

Choose a reason for hiding this comment

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

Thanks @philipwjones for the explanation. It didn't make sense to me either. Can we add a comment to explain that, if it's cleaner to retain the indexForReproSum with dimensions nSums and indexForReproSum(1) = min = 1 ?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@vanroekel this has always been part of the algorithm (essentially the "if deltaS .ne. 0.0" part) to account for possible cases where we explicitly don't want to apply restoring to some cells (eg, under land ice maybe?). then we don't actually want to include every cell in these sums. does that make sense?

Copy link
Contributor

@alicebarthel alicebarthel Jun 3, 2024

Choose a reason for hiding this comment

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

Suggested change
indexForReproSum(2) = 0
indexForReproSum(2) = 0
! For reproducible sums like mpas_global_sum_mod(), the range must be specified in min/max pairs for each array dimension.
! In this case, the array is getting packed in advance so the min is always 1 but the max can vary.
! Both need to be provided for the interface.

Copy link
Contributor

Choose a reason for hiding this comment

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

Thanks @philipwjones for the help understanding this. I'd like to include the comment from @alicebarthel in the PR for clarity. If it is too long for a comment, perhaps we can think of a shorter version.

real (kind=RKIND), dimension(nSums) :: reductions, sums
real (kind=RKIND), dimension(nSums) :: reductions

integer, dimension(nSums) :: indexForReproSum
Copy link
Contributor

@alicebarthel alicebarthel Jun 3, 2024

Choose a reason for hiding this comment

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

@philipwjones It's a minor detail, but I'm curious. Can you confirm if I understood this correctly: in this case, it's a coincidence that indexForReproSum and reductions have the same dimensions of nSums?
Technically, indexForReproSum should always have 2 (the min, max you mentioned), while nSums could be anything given the needs of the routine? (here, 2 since we need the area sum and the area-weighted salinity sum).
Does it make sense to use nSums to define indexForReproSum then? It seems misleading to me.

Copy link
Contributor

Choose a reason for hiding this comment

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

@alicebarthel Good catch. The index array should be dimensioned 2*(number of dimensions being summed over) - in this case just 2 since the arrays are just 1-d. So nSums is coincidentally correct but in incorrect in general.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@alicebarthel I wouldn't say it is "a coincidence" that nSums=2; this routine isn't meant to be general which is why I hardwired it to be correct for this case.

Copy link
Contributor

@alicebarthel alicebarthel Jun 4, 2024

Choose a reason for hiding this comment

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

@maltrud sure, I understand it's not meant to be general. I don't mind hard-coding the dimensions of indexForReproSum to 2. I mind hard-coding it to nSums if it's not supposed to have nSums as a dimension in general.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

ok, maybe I still don't understand the question. we have 2 sums to do: area and areadeltaS, so nSums = 2. so we need the same number of outputs from the 2 global sum calls: global_sum(area) and global_sum(areadeltaS) which are put into reductions(1) and reductions(2), respectively. so reductions should always be dimensioned nSums. I guess I could have just put sumAreaGlobal and sumAreaDeltaSGlobal directly in the global sum calls and not needed reductions(), but I wanted to keep the new version as similar to the current version as possible. does this make sense?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

ah, ok, now I understand. you're right, that is misleading. should I just place it with "2"?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

and a comment? thanks for clarifying.

Copy link
Contributor

Choose a reason for hiding this comment

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

I think you should create a variable ThroatWarblerMangrove that's pronounced "two"...

Copy link
Contributor

Choose a reason for hiding this comment

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

@maltrud Sounds good to me, thanks!

Copy link
Contributor

Choose a reason for hiding this comment

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

Could I suggest Rhabarberbarbarabarbarbarenbartbarbierbier as a suitable alternative?

@@ -250,7 +251,11 @@ subroutine ocn_get_surfaceSalinityData( streamManager, &
integer :: ierr

integer, parameter :: nSums = 2
real (kind=RKIND), dimension(nSums) :: reductions, sums
real (kind=RKIND), dimension(nSums) :: reductions
Copy link
Contributor

Choose a reason for hiding this comment

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

@maltrud just curious: is there any reason you went with reductions as a name for these variables? Sums is more intuitive to me, but maybe I'm not familiar with the jargon?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@alicebarthel I wanted to make as few code changes as possible so I retained the same variable names that already existed whenever possible. in the new code, "sums" isn't needed anymore, while "reductions" serves the same purpose in both the current and new versions. probably should have done things differently many years ago, but that's the rationale.

thanks for your careful reviewing.

Copy link
Contributor

Choose a reason for hiding this comment

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

@maltrud I totally understand the wisdom of making minimal code changes. As a newcomer to these routines, I allow myself to make minor comments as to readability to new users. Take it or leave it :-)

Copy link
Contributor

@alicebarthel alicebarthel left a comment

Choose a reason for hiding this comment

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

Approving based on visual inspection, provided that change to L256 is made.

@mark-petersen mark-petersen self-requested a review June 20, 2024 17:24
Copy link
Contributor

@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.

I compiled this with gnu (debug and optimized) on chicoma and intel on chrysalis and tested with the nightly test suite. That tests the compiling, but not the function of this code.

Thanks for running the G case tests. That is convincing. Since you are able to show reproducible sums on different processor counts, I'm happy to have this merged in.

@maltrud
Copy link
Contributor Author

maltrud commented Jun 20, 2024

I've made the recommended changes and push changes to the branch.

@jonbob
Copy link
Contributor

jonbob commented Jun 24, 2024

@vanroekel -- do you approve now with the changes?

Copy link
Contributor

@vanroekel vanroekel left a comment

Choose a reason for hiding this comment

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

Approving based on visual inspection and testing from @maltrud and @alicebarthel

jonbob added a commit that referenced this pull request Jul 9, 2024
Replace global sum with full-precision reproducible sum

The global sum in the surface salinity restoring routine was forced to
be BFB using a truncation calculation that resulted in reduced precision
and some unexplained behavior. This PR replaces this calculation with an
mpas standard reproducible sum with full precision.

Testing was done using a G-case with the IcoswISC30E3r5 grid for one
year, then calculating the global average SSS restoring tendency term,
which should be effectively zero. Also confirmed that the new code
resulted in BFB solution when processor layout was changed from 10 nodes
on chrysalis to 15 nodes.

Fixes #6446
[non-BFB] for C- and G-cases with salinity restoring
@jonbob
Copy link
Contributor

jonbob commented Jul 9, 2024

As expected, passes:

  • ERP_Ld3.ne30pg2_r05_IcoswISC30E3r5.WCYCL1850.chrysalis_intel.allactive-pioroot1
  • PET_Ln5.T62_oQU240wLI.GMPAS-DIB-IAF-DISMF.chrysalis_intel (G-case but no restoring file)

and DIFFs for:

  • SMS_D_Ld1.TL319_IcoswISC30E3r5.GMPAS-JRA1p5-DIB-PISMF.chrysalis_intel.mpaso-jra_1958 (G-case with restoring file)

merged to next

@jonbob jonbob merged commit 429dc36 into master Jul 10, 2024
13 checks passed
@jonbob jonbob deleted the maltrud/ocean/sss-restoring-reprosum branch July 10, 2024 15:51
@jonbob
Copy link
Contributor

jonbob commented Jul 10, 2024

merged to master and expected DIFFs blessed

xylar added a commit to xylar/compass that referenced this pull request Aug 8, 2024
This merge updates the E3SM-Project submodule from [c7d7998](https://github.com/E3SM-Project/E3SM/tree/c7d7998) to [727ad81](https://github.com/E3SM-Project/E3SM/tree/727ad81).

This update includes the following MPAS-Ocean and MPAS-Frameworks PRs (check mark indicates bit-for-bit with previous PR in the list):
- [ ]  (ocn) E3SM-Project/E3SM#6456
- [ ]  (ocn) E3SM-Project/E3SM#6510
- [ ]  (ocn) E3SM-Project/E3SM#6535
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug fix PR mpas-ocean non-BFB PR makes roundoff changes to answers.
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Global sum BFB code in the surface salinity restoring routine reduces precision.
7 participants