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

Fix uninitialized min_rand variable in Thompson MP when using SPP #892

Merged
merged 11 commits into from
Apr 20, 2022

Conversation

JeffBeck-NOAA
Copy link
Contributor

@JeffBeck-NOAA JeffBeck-NOAA commented Apr 2, 2022

Description

This suite of PRs fixes the uninitialized min_rand variable in Thompson MP when using SPP and is explained in detail in the ufs-weather-model PR here.

Issue(s) addressed

Fixes NCAR/ccpp-physics Issue #891.

Testing

The RT suite on Hera passed for all normal and debug tests, except for the regional_spp_sppt_shum_skeb RT test, which will now have a new baseline result.

@jwolff-ncar, @willmayfield, @bluefinweiwei, @michelleharrold, @judithberner, @gsketefian

@JeffBeck-NOAA JeffBeck-NOAA changed the title Fix to original min_rand variable. Fix uninitialized min_rand variable in Thompson MP when using SPP Apr 2, 2022
Copy link
Collaborator

@gthompsnWRF gthompsnWRF left a comment

Choose a reason for hiding this comment

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

Kinda important to catch the minimum value of rand_pert. Approved.

Copy link
Collaborator

@SamuelTrahanNOAA SamuelTrahanNOAA left a comment

Choose a reason for hiding this comment

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

This PR will not have consistent results between thread counts due to this line in module_mp_thompson.F90

abs_min_rand = ABS(MINVAL(rand_pert(:,1)))

@SamuelTrahanNOAA
Copy link
Collaborator

@grantfirl

Kinda important to catch the minimum value of rand_pert. Approved.

It doesn't solve the problem correctly though, since it is only getting the minimum of the current block, not the whole domain. The results will vary depending on threading and blocking settings.

@climbfuji
Copy link
Collaborator

climbfuji commented Apr 4, 2022 via email

@SamuelTrahanNOAA
Copy link
Collaborator

The random fields are generally normalized somehow, right? Perhaps [-1,1] or [0,1]. If you know how it was normalized, then you can use that value (-1 or 0 in those examples) instead of computing them.

@gthompsnWRF
Copy link
Collaborator

The random fields are generally normalized somehow, right? Perhaps [-1,1] or [0,1]. If you know how it was normalized, then you can use that value (-1 or 0 in those examples) instead of computing them.

I was going to ask about thread safe, but I just don't know how CCPP code is structured as well as I do WRF where I had implemented this in a thread-safe way.

Yes, actually I believe the min value should always be set to -1.0, so in this case, I think it is actually best to set abs_min_rand as a universal constant of 1.0. I realize now how this could be a flaw in my original design because instead of always asking for the lowest possible negative number to add to the original number, I was doing the abs of lowest number at the time, so in the case of slow spin-up, the min_rand number is hardly ever as low as -1.0. So I now believe this can be changed to a simple constant of 1.0 (which would be abs(-1.0) basically.

@JeffBeck-NOAA
Copy link
Contributor Author

JeffBeck-NOAA commented Apr 4, 2022

All blocks are consolidated into a single 1D array in the stochastic physics wrapper (stochastic_physics_wrapper.F90) before spp_wts_mp is passed to module_mp_thompson.F90 as rand_pert:

do nb=1,Atm_block%nblks
      GFS_Data(nb)%Coupling%spp_wts_mp(:,:) = spp_wts(nb,1:GFS_Control%blksz(nb),:,n)
end do

@SamuelTrahanNOAA
Copy link
Collaborator

It will still depend on MPI decomposition (layout_x and layout_y) since different ranks will have different parts of the domain.

@JeffBeck-NOAA
Copy link
Contributor Author

JeffBeck-NOAA commented Apr 4, 2022

The random fields are generally normalized somehow, right? Perhaps [-1,1] or [0,1]. If you know how it was normalized, then you can use that value (-1 or 0 in those examples) instead of computing them.

The random fields aren't normalized, but are capped based on namelist settings (magnitude*cutoff), so they will vary from run to run.

@JeffBeck-NOAA
Copy link
Contributor Author

The random fields are generally normalized somehow, right? Perhaps [-1,1] or [0,1]. If you know how it was normalized, then you can use that value (-1 or 0 in those examples) instead of computing them.

The ramdon fields aren't normalized, but are capped based on a namelist setting, so they will vary from run to run.

It will still depend on MPI decomposition (layout_x and layout_y) since different ranks will have different parts of the domain.

OK, that makes sense now. It looks like only a constant will work in this case.

@SamuelTrahanNOAA
Copy link
Collaborator

The only way I see of doing this is an MPI_Allreduce at the model level, and pass the result down to CCPP. It is likely the stochastic physics component already does that calculation, but it may not send the minimum to the model.

@gthompsnWRF
Copy link
Collaborator

Please waste no more energy on this and consider setting this as constant since I think it is the right thing to do. If Jeff has other thoughts, please speak out.

@climbfuji
Copy link
Collaborator

climbfuji commented Apr 4, 2022 via email

@JeffBeck-NOAA
Copy link
Contributor Author

JeffBeck-NOAA commented Apr 4, 2022

Are the following variables (rand1 and rand2) non-reproducable as well? It seems like if rand_pert(:,1) won't work, rand_pert(i,1) won't work either?

if (rand_perturb_on .ne. 0) then
     if (MOD(rand_perturb_on,2) .ne. 0) rand1 = rand_pert(i,1)
     m = RSHIFT(ABS(rand_perturb_on),1)
     if (MOD(m,2) .ne. 0) rand2 = rand_pert(i,1)*2.
     m = RSHIFT(ABS(rand_perturb_on),2)
     if (MOD(m,2) .ne. 0) rand3 = 0.25*(rand_pert(i,1)+abs_min_rand)
     m = RSHIFT(ABS(rand_perturb_on),3)
endif

@SamuelTrahanNOAA
Copy link
Collaborator

Are the following variables (rand1 and rand2) non-reproducable as well?

No. Those are in an i,j loop, so rand1 is a temporary variable containing data at that i,j point. (Likewise with rand2 and rand3).

@SamuelTrahanNOAA
Copy link
Collaborator

@gthompsnWRF If you want to have a meeting, please contact your desired attendees by email, and we'll arrange something. I'd rather not discuss our schedules and locations on github. Joe Olson should be there too, since his PR is the one that sparked all this.

@gthompsnWRF
Copy link
Collaborator

Jeff: yes, this is the best possible solution! I forgot that the random field is not intended to be [-1,1]. This solution is the proper one for certain.

@JeffBeck-NOAA
Copy link
Contributor Author

Jeff: yes, this is the best possible solution! I forgot that the random field is not intended to be [-1,1]. This solution is the proper one for certain.

Thanks, Greg! It doesn't look like it will be too difficult to bring these namelist entries into CCPP (two of the three are already available), so I should have it ready in a day or so. Thanks to @grantfirl for his help here.

@ligiabernardet ligiabernardet added the CCPP v6 Needed for CCPP v6 public release label Apr 5, 2022
@@ -1027,7 +1028,8 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
re_cloud, re_ice, re_snow
INTEGER, INTENT(IN) :: rand_perturb_on, kme_stoch
REAL, DIMENSION(:,:), INTENT(IN) :: &
rand_pert
rand_pert,spp_prt_list,spp_stddev_cutoff,n_var_spp, &
Copy link
Collaborator

Choose a reason for hiding this comment

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

I believe rand_pert is a 2D var, but how are these namelist variables dimensioned as (:,:) instead of 1-D? Shouldn't these additions be 1-D or at least separated from the (i,j) assumption of the rand_pert variable? Also, I do not find prt_list as intuitive. My first guess of this variable name would be something like "print_list" so is there anything more intuitive?

Copy link
Contributor Author

@JeffBeck-NOAA JeffBeck-NOAA Apr 6, 2022

Choose a reason for hiding this comment

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

@gthompsnWRF, none of this code is really ready to be reviewed yet, as I was just pushing things to my fork to then clone ufs-weather-model and run the stochastic physics RT. I probably should have kept all changes local until I got done with this work. Sorry for the confusion. That said, thanks for finding an error in the declarations here, those fields are 1D, and shouldn't go in this 2D line.

@@ -1101,7 +1103,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
REAL, DIMENSION(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic
REAL:: dt, pptrain, pptsnow, pptgraul, pptice
REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max
REAL:: rand1, rand2, rand3, min_rand
REAL:: rand1, rand2, rand3, spp_mp_mag_times_cutoff
Copy link
Collaborator

Choose a reason for hiding this comment

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

Couldn't this be made a simpler (and shorter) and intuitive name like: rand_pert_max?

Copy link
Contributor Author

@JeffBeck-NOAA JeffBeck-NOAA Apr 6, 2022

Choose a reason for hiding this comment

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

Yes, I was trying to keep some provenance for the source of this field, but it can be found in the routine itself (spp_prt_list(k)*spp_stddev_cutoff(k)), and rand_pert_max is a better idea.

do k =1,n_var_spp
select case (spp_var_list(k))
case('mp')
spp_mp_mag_times_cutoff = spp_prt_list(k)*spp_stddev_cutoff(k)
Copy link
Collaborator

Choose a reason for hiding this comment

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

There is no else clause here, so is it safe to have this variable without being initialized somewhere? Can it safely be set to zero (or perhaps did I miss it getting a default value)?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes, I still need to initialize the variable. I think setting it to zero will be fine (it's not used in that case anyway).

enddo
endif

print*, ' spp_mp_mag_times_cutoff is = ', spp_mp_mag_times_cutoff
Copy link
Collaborator

Choose a reason for hiding this comment

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

A print statement like this should be reserved for debugging only not in production code.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes, this was just for my testing purposes with this intermediate code. I will remove it once I confirm things are getting passed into the routine and applied correctly.

@JeffBeck-NOAA
Copy link
Contributor Author

@gthompsnWRF, @SamuelTrahanNOAA, please see the latest changes to this PR, which I believe will resolve your concerns. The code is now working (ran the RT test; baselines change) and calculates the maximum random perturbation value based on namelist settings brought in from stochastic_physics.

Copy link
Collaborator

@SamuelTrahanNOAA SamuelTrahanNOAA left a comment

Choose a reason for hiding this comment

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

Three dimension names are incorrect in the meta file. If you correct this, and update submodules, it compiles and runs on hera.intel. I also ran a -DDEBUG=ON test, and it ran as far as expected. (It failed in libpost due to enabling a post output variable it shouldn't.)

physics/mp_thompson.meta Outdated Show resolved Hide resolved
@JeffBeck-NOAA
Copy link
Contributor Author

Thanks, @SamuelTrahanNOAA. I missed this file in my latest commit. It's now fixed.

@SamuelTrahanNOAA SamuelTrahanNOAA self-requested a review April 7, 2022 21:40
@@ -644,7 +650,9 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, &
diagflag=diagflag, do_radar_ref=do_radar_ref_mp, &
has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, &
rand_perturb_on=spp_mp_opt, kme_stoch=kme_stoch, &
rand_pert=spp_wts_mp, &
rand_pert=spp_wts_mp, spp_var_list=spp_var_list, &
Copy link
Collaborator

Choose a reason for hiding this comment

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

It is sorta stupid request, but couldn't the new added arguments be listed in the same order in this block as in the block 40 lines lower? Just think it makes it easier. The order of arguments differs a little.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Sure, I can make that change.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Fixed.

dimensions = ()
type = integer
intent = in
[spp_prt_list]
Copy link
Collaborator

Choose a reason for hiding this comment

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

Clearly I am still being thick, but I don't see how the name "prt_list" gives me any clue of a magnitude of spp perturbation.

Copy link
Contributor Author

@JeffBeck-NOAA JeffBeck-NOAA Apr 8, 2022

Choose a reason for hiding this comment

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

Yes, it could be better named. We followed the same naming convention for SPP as what was already in place for the land perturbation scheme in stochastic_physics. For that scheme, "lndp_prt_list" is used for the perturbation array. At some point, changing to something like *_mag_list would be more descriptive.

Copy link
Collaborator

Choose a reason for hiding this comment

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

You could use a different name in mp_thompson than GFS_Typedefs so the variable will make sense in both contexts.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Since spp_prt_list is used across multiple physics suites and UFS submodules, I'd prefer keeping the name the same everywhere. We can certainly address this in a future PR, though.

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

Successfully merging this pull request may close these issues.

9 participants