CLUBB update and clubb_intr improvements

[huebleruwm]

There's two parts here - getting a new version of clubb in, and enabling the descending grid mode in clubb_intr. Both these goals were split up over many commits. Fixes #1411

New CLUBB

The first commits are dedicated to getting a new version of CLUBB in. Because of clubb_intr diverging between cam_development, which had significant changes to enable GPUization, and UWM's branch, which had redimensioning changes, the merging was done manually.

The first 6 commits here include changes that can be matched with a certain version of CLUBB release:
3d40d8c0e03a298ae3925564bc4db2f0df5e9437 works with clubb hash 1632cf12
e67fc4719fa21f8e449b024a0e3b6df2d0a7f8cb works with clubb hash 673beb05
187d7b536c2f36968fc7f5e1b9d1167e430ad03f works with clubb hash dc302b95
e4b71220b33aeaddb0afc68c9103555edccb59eb works with clubb hash dc302b95
703aca60ed1e0b6b24f2cd890c3a4497041d25b8 works with clubb hash d5957b30
4d9b1b8a528ca532d964c1799e1860e96e068a12 works with clubb hash d5957b30
(to use the clubb hash, go to src/physics/clubb and run the git checkout there)

These commits all have to do with just getting a new version of clubb in, so we need to ensure that at least 4d9b1b8a528ca532d964c1799e1860e96e068a12 is working correctly. The later commits have more complicated changes to clubb_intr, so if we find any problems with this branch, we should test commit 4d9b1b8a528ca532d964c1799e1860e96e068a12 as a next step.

clubb_intr improvements

The next commits after getting new clubb in are increasingly ambitious changes aimed at simplifying clubb_intr and reducing its runtime cost.

e60848b4ec4df90a3060ffd7f664fab42e847509 introduces a parameter, clubb_grid_dir, in clubb_intr that controls which direction the grid is when calling advance_clubb_core. When using -O0 to compile, and setting l_test_grid_generalization = .true. in src/clubb/src/CLUBB_core/model_flag.F90, the results are BFB in either grid direction.

dddff494966bf2bf4341fa4a7526b2f8b0f3d16e separates the flipping code from the copying code (copying cam sized arrays to clubb sizes arrays). Should all be BFB.

055e53f70741531a58d0f7da788b824c76fef087 pushes flipping code inward until it's directly around the call to advance_clubb_core, and the way of controlling the grid direction has been changed to a flag, l_ascending_grid. This should all be BFB as well, and I tested with clubb_cloudtop_cooling, clubb_rainevap_turb, and do_clubb_mf all true to ensure BFBness in either ascending or descending mode. One caveat - the clubb_mf code assumes an ascending grid, so before calling it we flip the arrays, then flip the outputs to descending.

2fb2ba9bd6b1ec5e5c60039f90d0c1020663d0c9 is a pretty safe intermediate commit, mainly moving stuff around in preparation for redimensioning.

bded8a561131e4dbbccad293f14226e5e8c0e856 is some very safe dimensionings.

fce8e1b1b5e8d3e93232c2129be7671fae79db23 is the big one that redimensions most pbuf arrays, and uses them instead of the clubb_sizes local ones. This allows us to avoid the vast majority of the the data copying, and delete the local version of the clubb arrays.

The rest are safe and easy things mainly, or commits to make GPU code work.

Testing

I plan to run an ECT test to compare the cam_development head I started with to the head of this branch. Answers are expected to change slightly due to the ghost point removal in clubb, so I think it unlikely that this passes, but if it does that would be fantastic, and might be the only testing we really need.

If that first ECT test doesn't pass, then I'll will assume that the difference is due to the ghost point removal (or other small bit changing commits in clubb over the past year), and rely on @adamrher to check if the results are good.

The biggest concern is if the answer changes are acceptable, and the only real differences expected are from the new version of CLUBB. If the answers from this branch look bad, we should go back to commit 4d9b1b8a528ca532d964c1799e1860e96e068a12 and check the answers from that, since it only includes the new version of CLUBB, and NO unnecessary clubb_intr changes. If the answers still look bad in that commit, then we have a few more we can step back through to try to figure out which commit introduces the differences. If the hypothetical problem is present in the first commit (3d40d8c0e03a298ae3925564bc4db2f0df5e9437), then the problem is harder, because that includes (pretty much only) the removal of the ghost point in clubb, which is expected to change answers, but hopefully not significantly.

Again if that first ECT test fails, then I can still run another ECT test between the version where clubb is up to date (4d9b1b8a528ca532d964c1799e1860e96e068a12) and the head of this branch. The changes between that commit and the head may be slightly bit changing without (-O0), but definitely shouldn't be answer changing. If this ECT test fails, then I've made a mistake in the changes meant to avoid the flipping/copying, and I'll have to step back through and figure out what bad thing I did.

Some other tests I've ran along the way help confirm at least some of these changes:

• I've ensured BFBness between ascending and descending modes since it was introduces in commit e60848b4ec4df90a3060ffd7f664fab42e847509
• The ERP test passes
• Commit 4d9b1b8a528ca532d964c1799e1860e96e068a12 should be working on the GPU, but later commits are untested

Next steps

I left a number of things messy for now, such as comments and gptl timers. Once we confirm these changes, I'd like to go through and make some of those nicer as a final step.

Performance

In addition to the clubb_intr changes improving performance, we should use this as an opportunity to try to flags that should be significantly faster:

• clubb_penta_solve_method = 2 uses our custom pentadiagonal matrix solvers, which should be significantly faster than lapack and should pass an ECT test
• clubb_tridiag_solve_method = 2 uses our custom tridiagonal solver, which should also be faster and pass an ECT test
• clubb_fill_holes_type = 4 uses a different hole filling algorithm, that in theory is just all around better and faster than our current one, and I suspect it will pass ECT, but I have yet to test it
• clubb_l_call_pdf_closure_twice = .false. will avoid a pdf_closure call (which is significant in cost) and reduce the memory footprint , and I think will have minimal effect (based on my visual analysis of what that affects in the code), but is the most likely to break the ECT test

[cacraigucar]

@huebleruwm - After reading through this text, I moved this PR to draft. Once it is ready for the SEs to review and process it for bringing into CAM, please move it out of draft.

[huebleruwm]

@adamrher I have done a number of ECT tests and it was pretty much as expected. I made a baseline from ESCOMP/cam_development commit 4c480c8e8c7e3138d8d5cb3118f635c71886ecf7 (which has cam6_4_128 as the last tag). This commit from this branch (4d9b1b8a528ca532d964c1799e1860e96e068a12), which is just getting in the newest version of clubb) fails the ECT test when compared to the cam_development baseline. So we will need some help confirming the changes as acceptable.

To check the more "dangerous" changes to clubb_intr, I made another baseline from commit 4d9b1b8a528ca532d964c1799e1860e96e068a12 (with new clubb), then compared that to the head of this branch, and that passed the ECT test. I also made sure to run with settings that have top_lev = 3 so that we can test the redimensioning changes, which I did by setting trop_cloud_top_press to 1.D3 in the namelist instead of 1.D2.

Other things I tested:

• When running with -O0 and l_test_grid_generalization = .true. (clubb internal parameter), then the code is BFB between the ascending and descending modes in clubb_intr (changed by setting l_ascending_grid)
• Also when running with -O0 the code is to be BFB between 4d9b1b8a528ca532d964c1799e1860e96e068a12 and the head of this branch, outside of the few output variables that differ above top_lev, see the comment from this commit for more detail.
• The ERP test passes
• The GPU code runs to completion, haven't check results with ECT test yet, but that's the only thing on the list left to do, and any further changes to make that work (if needed) won't affect CPU results.

So I think the clubb_intr changes are pretty well tested and confirmed. But the changes from new clubb don't pass the ECT test, and we'll need to confirm them as acceptable.

Note

Also, I found a number of examples where clubb_intr is interacting above top_lev, and I'm not sure exactly what we should do about that. For now I've left the errors in place for BFBness, but I've made little code comments explaining the issue (all contain a TODO). There's also a TODO with how we're setting concld_pbuf, which seems suspicious, but not a clear bug. @adamrher could you take a look at these "TODO" comments and see if you have any insight please. I added what I think the "correct" lines should be, and what the errors were that I left in for BFBness.

[huebleruwm]

I ran some ECT tests with the performance options:

• clubb_tridiag_solve_method = 2 passes the ECT test and should be significantly faster
• clubb_penta_solve_method = 2 passes the ECT test and should be significantly faster, but requires clubb_release commit bba14c86a748da7993af72c056a9f6505cf43f1b or newer, since I just fixed a bug in it
• clubb_fill_holes_type = 4 FAILS the ECT test, which was expected since it's fills holes in the same conservative way, but over completely different ranges, making it answer changing. This should also be decently faster, so it would be nice to test if possible.

[adamrher]

@huebleruwm I checked out this PR branch, and while the clubb_intr.F90 changes seem to be there, the .gitmodules are pointing to the clubb externals currently on cam_development:

[submodule "clubb"]
        path = src/physics/clubb
        url = https://github.com/larson-group/clubb_release
        fxrequired = AlwaysRequired
        fxsparse = ../.clubb_sparse_checkout
        fxtag = clubb_4ncar_20240605_73d60f6_gpufixes_posinf
        fxDONOTUSEurl = https://github.com/larson-group/clubb_release

This should be pointing to a different tag / hash containing support for descending / ascending options, no?

[huebleruwm]

This should be pointing to a different tag / hash containing support for descending / ascending options, no?

Yes. In hindsight, it would've made much more sense to update that with the right hash each commit, rather than just including it in the commit comment.

I've just been going to src/physics/clubb and running git checkout master to get the newest version. Or git checkout hash when testing the first handful of commits:
cam 3d40d8c0e03a298ae3925564bc4db2f0df5e9437 works with clubb hash 1632cf12
cam e67fc4719fa21f8e449b024a0e3b6df2d0a7f8cb works with clubb hash 673beb05
cam 187d7b536c2f36968fc7f5e1b9d1167e430ad03f works with clubb hash dc302b95
cam e4b71220b33aeaddb0afc68c9103555edccb59eb works with clubb hash dc302b95
cam 703aca60ed1e0b6b24f2cd890c3a4497041d25b8 works with clubb hash d5957b30 or newer
cam 4d9b1b8a528ca532d964c1799e1860e96e068a12 works with clubb hash d5957b30 or newer

[adamrher]

I just pushed my last commit, so it's all yours now. It changes clubb_mf to operate on descending arrays, and I've confirmed it's b4b with the ascending version of the code. We can now get rid of the ascending code entirely. I'll leave that to you, along with making a new clubb tag and putting it in .gitmodules (Externals_CAM is deprecated).

I should be able to finish up the science validation by tomorrow -- I'll let you and Vince know.

[adamrher]

Is the clubb update @HingOng's PR incorporating non-traditional coriolis terms? That would be great to include here.

[huebleruwm]

Is the clubb update @HingOng's PR incorporating non-traditional coriolis terms? That would be great to include here.

Yes, that's the one. Just waiting on the nightly tests to make sure they all pass. It should all be BFB with the new flags set to false.

[vlarson]

I just pushed my last commit, so it's all yours now. It changes clubb_mf to operate on descending arrays, and I've confirmed it's b4b with the ascending version of the code. We can now get rid of the ascending code entirely. I'll leave that to you

Which ascending code is that? Is this ascending code you're referring to specific to CLUBB-MF or does it apply more generally to CLUBB? I thought we were going to leave in some ascending code just to allow testing of the ascending version of CLUBB within CAM.

[adamrher]

Which ascending code is that? Is this ascending code you're referring to specific to CLUBB-MF or does it apply more generally to CLUBB? I thought we were going to leave in some ascending code just to allow testing of the ascending version of CLUBB within CAM.

I forgot that we did agree to keep the ascending code option for CLUBB for debugging purposes. As for CLUBB-MF, the integrate_mf call in clubb_intr.F90 operates on ascending arrays, and so here I just flipped it so it operates on descending arrays. I suppose we can either add an endrun statement if clubb_l_ascending_grid=.true. and do_clubb_mf=.true. or I can add logic so that CLUBB-MF can run in ascending mode as well. I would need to make clubb_l_ascending_grid available in clubb_mf.F90, which contains the relevant integrate_mf changes.

[huebleruwm]

I suppose we can either add an endrun statement if clubb_l_ascending_grid=.true. and do_clubb_mf=.true. or I can add logic so that CLUBB-MF can run in ascending mode as well.

Setting clubb_l_ascending_grid only flips things directly before and after the advance_clubb_core call, so it won't affect integrate_mf at all. So I don't think any further changes are needed. I suppose we could confirm this by testing clubb_l_ascending_grid=.true. with clubb_mf, since they should work together and produce the same answer, but it's probably safe to skip.

[adamrher]

Setting clubb_l_ascending_grid only flips things directly before and after the advance_clubb_core call, so it won't affect integrate_mf at all. So I don't think any further changes are needed. I suppose we could confirm this by testing clubb_l_ascending_grid=.true. with clubb_mf, since they should work together and produce the same answer, but it's probably safe to skip.

Ah, yeah, you're right. Bottom line is that I changed integrate_mf to operate on descending arrays, which allowed me to remove the array flipping before and after the call to integrate_mf in clubb_intr.F90. This reduces the CLUBB-MF footprint in clubb_intr.F90 and so I think it's a welcome change, and is independent of whether or not CLUBB operates in ascending or descending mode.

[huebleruwm]

I added the changes needed for the new clubb coriolis code and made a tag that I updated the gitmodules with (clubb_4ncar_20260109_ddf5110).

[huebleruwm]

I did some more targetting timings with the options were considering turning on by default. All speedups are measured against a baseline run with the current default settings.

• setting clubb_penta_solve_method = 2 reduces the cost of advance_clubb_core by ~10%
• setting clubb_tridiag_solve_method = 2 reduces the cost of advance_clubb_core by ~4%
• setting clubb_fill_holes_type = 4 reduces the cost of advance_clubb_core by a dissapointing ~2%

The small speedup from the new fill_holes method seems to suggest that filling is either rare in cam or being done on sparser hole patterns than what's seen in clubb_standalone. The upside is that if the new hole filler is both changing answers and barely a speedup, then we can just leave it set to 2 and we still get most of the speedup we were hoping for and no answer changes - assuming the new matrix solver options aren't also the causing answer changes.

[adamrher]

The upside is that if the new hole filler is both changing answers and barely a speedup, then we can just leave it set to 2 and we still get most of the speedup we were hoping for and no answer changes

This turns out to be the case. I ran a 20 year F-case with the new solver options and using the default hole filling option 2 (7d_pe2), and the answers aren't any different from the control. The table below shows that the area of the globe with significant answer changes is not outside of natural variability (0a*-0a) for most variables (note I don't put much weight on the T850/T950 result since the areas of significance occur where my interpolation of T to a fixed pressure level intersects the surface, which I don't treat properly for this field).

[vlarson]

The upside is that if the new hole filler is both changing answers and barely a speedup, then we can just leave it set to 2 and we still get most of the speedup we were hoping for and no answer changes

This turns out to be the case. I ran a 20 year F-case with the new solver options and using the default hole filling option 2 (7d_pe2), and the answers aren't any different from the control.

For this configuration, what speed-up do you see?

[adamrher]

@vlarson just looking at one of the 2 year segments for each case, the mean clubb_tend_cam time in 7d_pe2 is reduced by 9.5% relative to 6d. Compare that to 7d_pe3, which is reduced by 10.3% relative to 6d. Gunther's numbers would suggest combining using the two new solvers would create a 14% reduction in 7d_pe2 (provided clubb_intr is identical in all these runs, which I believe is true). There's an envelope of error due to load imbalancing, which I tend to think means longer runs have more reliable timings.

[adamrher]

@huebleruwm this is my first pass at a code review, making it only to the beginning of clubb_tend_cam. I had a few questions and clarifications but no major issues. There's a lot of code clean up and variable clarification which is awesome.

I'm not sure how much time I'll have in the next two weeks to finish this code review, with AMWG coming up two weeks from today. I'll get to it as soon as I can.

[adamrher]

@huebleruwm I made some more progress on the code review, but I'll have to finish it up another day. I really like how you've added the _pbuf modifier to all the pbuf vars. It will make this code a lot easier to understand going forward!

[adamrher]

Publishing my morning review now so it doesn't get lost. Will pick up with a separate review this afternoon.

[adamrher]

If this is not the surface pressure (state_loc%pint(i,pverp)) than it needs to named something else. Perusing advance_clubb_core.F90 it is consistently referred to as surface pressure in the comments, so I'm skeptical of setting this to state_loc%pmid(i,pver) (the mid point pressure in the lowest model level).

[huebleruwm]

After talking with Vince and checking it out in the code - p_sfc is the right name and it should be set to state_loc%pint. So this is a bug.

The only place it's used in the code is pdf_closure, see here, and is used when l_call_pdf_closure_twice = .true., which it is in cam. So this will change results.

I also did a test of changing this from %pmid to %pint, and doing so breaks the ECT test for ~40 variables, many of which report >5 sigma differences. So it seems like we'll need another science evaluation if we want to change this. I'll leave it as is for now and add a little TODO comment explaining the problem.

[adamrher]

Darn. I'll bring this up in our SE mtg tomorrow and report back on next steps.

[adamrher]

Just to clarify, p_sfc was not an argument to advance_clubb_core in the clubb tag currently on cam_development. Is it a bug on cam_development currently? Or did we introduce this bug in this PR?

[huebleruwm]

Looks like the way this was handled before was by using p_in_Pa(i,1) in place of p_sfc(i).

• Before the ghost point removal we used: p_in_Pa_zm(i,1) = p_in_Pa(i,1)
• After the removal that became: p_in_Pa_zm(i,1) = p_sfc(i)

The comment in that section of pdf_closure where it's used says this:

! The pressure at thermodynamic level k = 1 has been set to be the surface
! (or model lower boundary) pressure.  Since the surface (or model lower
! boundary) is located at momentum level k = 1, the pressure there is
! p_sfc, which is p_in_Pa(1).  Thus, p_in_Pa_zm(1) = p_in_Pa(1).

So internally, clubb was assuming that the ghost point in p_in_Pa contained the "true" (momentum) surface pressure, and that IS how it was being set in clubb_standalone, using

p_in_Pa(1) = p_in_Pa_zm(1)

but cam was setting

p_in_Pa(i,k+1)         = state1%pmid(i,pver-k+1)
p_in_Pa(i,1)           = p_in_Pa(i,2)

and I think that should have been p_in_Pa(i,1) = state1%pint(i,pver)

So this was a bug that's been in cam the whole time - I tracked it back 7 years to a commit named CAM6 initial release where it first appeared, and clubb has been making the same internal assumption that p_in_Pa(1) was the momentum surface pressure the whole time.

[adamrher]

I'm preforming a science eval on this bug fix. I'll provide an update here ASAP.

[adamrher]

@huebleruwm I've finished my review of your PR. I have several comments that should be addressed. I do think you found a bug in our code for top_lev > 1 in the energy fixer and the applied tendencies, so thanks for catching that!

[adamrher]

I agree with @huebleruwm that that this is a bug in the code that would lead to spurious ptend%s tendencies above top_lev. While top_lev determines the top of the moist physics, and therefore where clouds can form, they may be advected to higher levels by the dycore. I endorse @huebleruwm's fixes (which are currently commented out to give bfb with the buggy code)

[edit -- april 9 -- I have committed this fix in b740f05c8a76e0985283057760f3a33a5a036d68]

[huebleruwm]

@adamrher I think I've addressed all the PR comments now and made the changes over the last 2 commits. They were all pretty small changes - moving things around, renaming stuff, comment tweaks, etc. I ran the usual ECT tests on CPU and GPU and both are passing.

Unless I've missed something, all that's left to do is deal with the p_sfc issue.

[cacraigucar]

Still one TODO here . Please remove or modify the comment

[huebleruwm]

This is the p_sfc issue that Adam is looking into right now, #1441 (comment) - I just figured I'd mark it to make clear that we're not 100% done yet.