Microphysics graupel

[OngChia]

Port of the ICON microphysics graupel scheme. The graupel scheme code has been refactored in favor of a more readable style. Currently, we verify it by using the serialized data from Weisman-Klemp experiment performed on the torus grid. Because the ported graupel scheme can only be run on embedded backend, which is slow, a manual slicing is performed in the test. So, only individual stencils are tested in this PR. Saturation adjustment may be called right after graupel scheme. Currently, I put a an option in the configuration to do that without taking into account the tendencies resulted from the graupel scheme. This needs to be changed too when slow performance is resolved.

The graupel scheme is split into three parts:

1. the main algorithm for computing the temperature and tracer tendencies due to cloud microphyisal processes. It is wrapped in a scan operator because the original graupel scheme contains implicit computation for the sedimentation fluxes. A number of local variables are also output for determination of rain, snow, graupel, and ice sedimentation fluxes in the following steps.
2. computation of rain, snow, graupel, and ice sedimentation fluxes above the ground.
3. computation of rain, snow, graupel, and ice sedimentation fluxes on the surface.
   All the pre-computed coefficients in gscp_data.f90 are now stored in the FrozenNameSpace class SingleMomentSixClassIconGraupelParams. They can be directly used in the scan operator.

TODO:

1. SingleMomentSixClassIconGraupelParams contains some constants in constants.py and saturation_adjustment.py. Refactor them when gt4py allows direct import use for compile-time constants.
2. Better renaming of the local variables to be done later.
3. And making it to conform to the physics interface protocol of icon4py should be done in a follow-up PR.
4. Remove the slicing when gtfn backend is supported.

[OngChia]

cscs-ci run default

[OngChia]

cscs-ci run default

[OngChia]

cscs-ci run default

[OngChia]

cscs-ci run default

[halungge]

You have the exact convergence criteria twice. Wouldn't it make more sense to turn this into a while type loop:

num_iter = 0
while(not converged and num_iter < self.config.max_iter):
    
    # perform updates
    num_iter = num_iter + 1

where converged = xp.any() is short for your convergence criteria
or if you want to have to raise the ConvergenceError

num_iter = 0
while(not converged):
    if num_iter >= self.config.max_iter:
          raise ConvergenceError(f"No convergence reached after maximal number of iterations: {self.config.max_iter}")
    # perform updates
    (...)
    num_iter = num_iter + 1

The raise will already break your control flow, so you don't need the extra break

[OngChia]

Thanks for the nice suggestion. Indeed my original code is not good. I would adopt the your second idea that allows me to raise convergence error.

[halungge]

you might want to rename this to test_single_moment_six_class_gscp_graupel.py

[halungge]

Out of curiosity: what does the name mean. I know the gscp comes from ICON but I don't know what it stands for, what is the "six_class" ?

[OngChia]

I also do not know what gscp stands for. I could not find its meaning in the source code and any ICON documentation. six_class refers to the number of hydrometeors you solve. In ICON graupel scehem, it treats moisture (qv), cloud droplets (qc), rain droplets (qr), ice crystals (qi), snow crystals (qs), and graupel (qg) as prognostic variables. There are six of them. So, we call it six-class microphysics. The single moment means that it only solves the mass of these hydrometeors while the number of particles are diagnosed from the mass and other factors. For example, double moment means that the number of particles are also treated as prognostic variables. So, ### moment ### class microphysics is kind of a standard name which lets scientists immediately know what kind of microphysics scheme we are using.
I retain the gscp_graupel just so that ICON community knows easily it is the graupel scheme in ICON.

[halungge]

For example, once using the component interface, you would not have to setup this structs which do not match you needs (except for the tracers) but you would just get the rho, temperature and pressure from the state.

[halungge]

Can you convert these option into enum values?

class LiquidAutoConversion(enum.Enum):
    KESSLER = 1
    SEIFERT_BEHENG = 2

Suggested change

	liquid_autoconversion_option: gtx.int32 = gtx.int32(1)
	liquid_autoconversion_option: LiquidAutoConversion = LiquidAutoConversion.KESSLER

Actually, I would do this for all "namelist switches". That is configuration parameters that have a couple of fixed valid entries. These are essentially the ones that are integer type in ICON namelist.
If you need the int value (for example inside a stencil) you can always do LiquidAutoConvertions.KESSLER.value

[OngChia]

Okay. Thanks.
I noticed that if liquid_autoconversion_option == LiquidAutoConversion.KESSLER.value: does not work in a field_operator. I add a TODO whenever this happens because I think it is better to be able to directly use the value in LiquidAutoConversion instead of if liquid_autoconversion_option == 0 which may lead to an out-of-sync problem in case that LiquidAutoConversion.KESSLER is set to other values than 0 in the future.

[halungge]

could you use the general DBL_EPS value that we define in constants.py: which again is system dependent. It simply DBL_EPS = sys.float_info.epsilon or do you need that precise value?

[OngChia]

From my understanding, I can simply use DBL_EPS for eps because it is only a small number that prevent division-by-zero that should be system dependent. For qmin, I would need that precise value. A different value will greatly alter the result based on previous my experience with different microphysics. Users may need to tune the parameters again for better predictability if qmin changes.

[halungge]

Generally, I would give the naming of the parameters (config and the constants) a second thought. For me, they can still be more explicit (and longer) but editors and IDE s come with auto-completion so short is not an asset. As an extreme example:
rain_v_sedi_min does not gain anything compared to v_sedi_rain_min then I would rather use the ICON version to keep consistency, or

[OngChia]

From you suggestion of a new naming for ice_sedi_density_factor_exp, I came up with a better naming minimum_rain_fall_speed for v_sedi_rain_min as I wrote in my comment above. Following your suggestion, I will retain remaining parameters which I still do not have a good idea of a better name.

[halungge]

As for the stencils (field_operators and scan_operator) I cannot really judge , I assume you do it right. And you figured out a way to surf around the gt4py deficiencies for physics...

[OngChia]

Hannes and Till once looked at the older version of this graupel scheme. I only sporadically communicate with them about the structure of the scan_operator and field_operator. Perhaps it is better for me to also ask Till to briefly look at my PR for comments on how I use the scan_operator and field_operator in the graupel scheme.

[OngChia]

cscs-ci run default

[OngChia]

launch jenkins spack

[OngChia]

launch jenkins spack

[github-actions]

Mandatory Tests

Please make sure you run these tests via comment before you merge!

• cscs-ci run default
• launch jenkins spack

Optional Tests

To run benchmarks you can use:

• cscs-ci run benchmark

To run tests and benchmarks with the DaCe backend you can use:

• cscs-ci run dace

In case your change might affect downstream icon-exclaim, please consider running

• launch jenkins icon

For more detailed information please look at CI in the EXCLAIM universe.