Physics driver l2#1301
Conversation
msimberg
left a comment
There was a problem hiding this comment.
This is not really a review, but more questions and comments about minor and major things. Mostly for my own understanding and reminder to revisit as the PR evolves, but adding the comments in case someone else has thought about them already and has answers.
| def __call__( | ||
| self, state: dict[str, model.DataField], time_step: datetime.datetime | ||
| ) -> dict[str, model.DataField]: |
There was a problem hiding this comment.
Question for myself:I suppose this is one of the critical interfaces, though the signature could likely be tightened to something resembling state: (Something?)State, time_step: datetime.datetime? Or the whole class should inherit from a protocol.
There was a problem hiding this comment.
Minor note as well that may need broader discussion: If this is indeed the protocol, my preference is to give this a name rather than using the generic __call__ operator, but I also know we already use __call__ extensively, so may be me still not being used to __call__ being used everywhere.
There was a problem hiding this comment.
And further thinking out loud: Could this be based on a more generic TendencyAdaptor? The muphys-specific parts seem quite minimal and the interesting logic seems to be the _to_tendency parts.
There was a problem hiding this comment.
and time_step should be datetime.timedelta if we decide to use timedeltas for parametrizing (as opposed to the integer time steps)
There was a problem hiding this comment.
correction: it should be RelativeTime either from standalone_driver or time.py depending on whether #1353 is merged before this one or not
| self, | ||
| prognostic: prognostics.PrognosticState, | ||
| outputs: dict[str, fa.CellKField[ta.wpfloat]], | ||
| dt: float, |
There was a problem hiding this comment.
I guess this is in seconds? I wonder if we can add a bit more type information to this (something from datetime?) and just convert to the correct units as late as possible (i.e. when passing into gt4py programs).
There was a problem hiding this comment.
after a separate discussion with @msimberg we came up with this for the standalone driver:
https://github.com/C2SM/icon4py/pull/1304/changes#diff-1374141d85d171d3fa2c5711b656072708d65ff1080beaf0cdeedcc64f105bd7R37
should we use it here too?
Readers for the three new savepoints generated by the instrumented ICON run (icon-nwp add_exp_ape_aes): - IconTimeStepExitSavepoint: end-of-timestep prognostic state, written after integrate_nh returns (all physics applied, time levels swapped) - IconAesGraupelInit/ExitSavepoint: state + prm_tend accumulators around the mig block (cloud_mig = satad + graupel + satad), plus aggregated surface precip (rsfl, ssfl, pr, ufcs) Replace test_standalone_driver_moist_physics: instead of only checking the fields muphys does not touch against diffusion/nonhydro exits, the full prognostic state (vn, w, rho, exner, theta_v, all six tracers) is now validated against the single time-step-exit savepoint. Provisional tolerances document the known phy2dyn fidelity gaps (linearized exner update, no theta_v update, no negative-tracer clipping, full-column vs jks_cloudy) until they are ported. Add test_aes_graupel_granule: first validation of the muphys granule against Fortran ICON (previous coverage was netCDF references from the C++ reference implementation only). Compares component tendencies and surface precip against the aes-graupel-init/exit pair, restricted to the jks_cloudy..nlev levels ICON computes, and checks the granule produces (near-)zero tendencies above. Tolerances are provisional until measured on the regenerated exclaim_ape_aesPhys v05 archive (overwritten in place, no version bump). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Reader classes, provider methods, test file and function names now use the muphys nomenclature; the serialbox savepoint names (aes-graupel-init/exit) and the Fortran sources keep the aes_graupel naming so it stays clear the data interface is the icon-nwp Fortran mig block, not the icon-mpim Kokkos port. Also drop the qnc-units bug commentary obsoleted by the MuphysConfig fix (qnc now 50.0e6 m^-3, matching Fortran cloud_num) and read qnc from the config again in the datatest. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The muphys port followed the older muphys C++ reference implementation, while icon-nwp carries the newer MPIM rain-microphysics revisions. This adds the icon-nwp formulation (mo_aes_graupel.f90) as a second scheme, selected via MuphysConfig.scheme / a use_icon_nwp flag threaded through graupel_run and muphys_run: - cloud_to_rain: rho-dependent polynomial accretion kernel in log(clamped rho*qr) instead of the constant 5.25 - rain_to_vapor: exp-polynomial evaporation instead of the power law - cloud_to_snow: riming tuning factor 3.0 - snow_number/snow_lambda: operate on the snow mass density rho*qs with a lower clamp instead of the additive (qs + 2e-6) offset, and branch on rho*qs > qmin - fall speeds (vm): rain becomes a log-polynomial (with the constant term not density-corrected, faithfully reproducing the Fortran operator precedence), ice prefactor 0.80, snow prefactor 115.6 with snow_number evaluated at the clamped falling-mass density and the previous level's temperature (now carried in the scan state), and a two-sided clamp of the hydrometeor density for all species - ice_dep is computed from the vapor_x_ice deposition before the nucleation contribution is added, matching the Fortran ordering Defaults: setup_graupel/setup_muphys/MuphysComponent keep CPP_REFERENCE, so the existing netCDF-based tests are unchanged (verified: graupel_only/full_muphys/component datatests pass at rtol 1e-14 after the refactor); MuphysConfig defaults to ICON_NWP, so the standalone driver and the aes-graupel serialbox datatest use the ported scheme matching the Fortran that generates the reference data. Also fix the fused muphys_run to return the temperature updated by the final saturation adjustment (as the Fortran does), consistent with the adjusted qv/qc; the separate-programs path already did this. New stencil tests validate every ported formula against independent numpy references (pass on gtfn_cpu); an execution smoke run confirms the ICON_NWP path compiles, stays finite, and differs from CPP_REFERENCE by physically plausible magnitudes. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Resolve conflicts from PR #1347 (APE experiment) landing on both branches plus main's #1324 (analytical geometry), #1303, #1345: - Took main's canonical/refined versions of the shared #1347 files (constants, thermodynamic_functions, pressure [incl. the surface_pressure_k cleanup], driver_io, jablonowski_williamson, utils, initial_condition/config). - config.py: kept both the muphys config wiring and main's GeometryConfig. - test_initial_conditions.py: kept our qv atol=1e-13 (default atol=0 fails on qv). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Scheme names now say what each formulation is: KOKKOS_MUPHYS is the muphys C++/Kokkos reference implementation, AES_GRAUPEL is the icon-nwp mo_aes_graupel.f90 formulation. All _icon_nwp suffixes (field operators, programs, constants class, stencil tests, the use_icon_nwp flag) become _aes_graupel accordingly. The Python-level scheme dispatch in setup_graupel/setup_muphys uses match/case with an explicit case per scheme. Inside the GT4Py operators the branches stay scalar-bool if/else with a labeled else: the DSL has no match statement, and an if/elif pair on a bool is rejected by its definedness analysis (it cannot prove the two conditions exhaustive, so names assigned in the branches count as possibly undefined). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> # Conflicts: # model/standalone_driver/tests/standalone_driver/integration_tests/test_standalone_driver.py
There was a problem hiding this comment.
todo after you fix the merge:
I had not noticed this change from t to te in the other PR.
I think these two should be handled with if use_icon_nwp like everything else, i.e.:
if use_icon_nwp:
# return the temperature updated by the final saturation adjustment (as the
# Fortran does), consistent with the adjusted qv/qc
return te, Q(v=qve, c=qce, r=q.r, s=q.s, i=q.i, g=q.g), pflx, pr, ps, pi, pg, pre
else:
return t, Q(v=qve, c=qce, r=q.r, s=q.s, i=q.i, g=q.g), pflx, pr, ps, pi, pg, pre
just with the corrected use_icon_nwp or match case consistent with the other changes.
| def __call__( | ||
| self, state: dict[str, model.DataField], time_step: datetime.datetime | ||
| ) -> dict[str, model.DataField]: |
There was a problem hiding this comment.
correction: it should be RelativeTime either from standalone_driver or time.py depending on whether #1353 is merged before this one or not
jcanton
left a comment
There was a problem hiding this comment.
more comments but haven't finished reviewing yet...
There was a problem hiding this comment.
| if use_aes_graupel: | |
| # return the temperature updated by the final saturation adjustment (as the | |
| # Fortran does), consistent with the adjusted qv/qc | |
| return te, Q(v=qve, c=qce, r=q.r, s=q.s, i=q.i, g=q.g), pflx, pr, ps, pi, pg, pre | |
| else: | |
| return t, Q(v=qve, c=qce, r=q.r, s=q.s, i=q.i, g=q.g), pflx, pr, ps, pi, pg, pre |
msimberg
left a comment
There was a problem hiding this comment.
These are starting to be nitpicky comments, so this is definitely getting close.
The Component/State discussion (i.e. the protocols not being ready) is my biggest concern, but we leave that out from this PR.
jcanton
left a comment
There was a problem hiding this comment.
more comments I had forgotten about
| def _require(field: fa.CellKField[ta.wpfloat] | None, name: str) -> fa.CellKField[ta.wpfloat]: | ||
| """Return ``field``, or raise if it is inactive (``None``). | ||
|
|
||
| @gtx.field_operator | ||
| def _add_scaled_tendency( | ||
| field: fa.CellKField[ta.wpfloat], | ||
| tendency: fa.CellKField[ta.wpfloat], | ||
| dt: ta.wpfloat, | ||
| ) -> fa.CellKField[ta.wpfloat]: | ||
| return field + tendency * dt | ||
|
|
||
|
|
||
| @gtx.program(grid_type=gtx.GridType.UNSTRUCTURED) | ||
| def apply_tendency( # noqa: PLR0917 # stencil params referenced in domain specs must stay positional | ||
| field: fa.CellKField[ta.wpfloat], | ||
| tendency: fa.CellKField[ta.wpfloat], | ||
| dt: ta.wpfloat, | ||
| result: fa.CellKField[ta.wpfloat], | ||
| horizontal_start: gtx.int32, | ||
| horizontal_end: gtx.int32, | ||
| vertical_start: gtx.int32, | ||
| vertical_end: gtx.int32, | ||
| ) -> None: | ||
| """``result = field + tendency * dt`` (Lie-Trotter split-update). | ||
|
|
||
| Pass ``result=field`` for an in-place update, or a distinct buffer to keep | ||
| the original (e.g. computing the new temperature without clobbering ``te``). | ||
| muphys needs all six moisture species; ``TracerState`` fields are optional | ||
| (a tracer may be inactive per ``TracerConfig``), so we fail loudly here rather | ||
| than feed ``None`` into the microphysics. | ||
| """ | ||
| _add_scaled_tendency( | ||
| field, | ||
| tendency, | ||
| dt, | ||
| out=result, | ||
| domain={ | ||
| dims.CellDim: (horizontal_start, horizontal_end), | ||
| dims.KDim: (vertical_start, vertical_end), | ||
| }, | ||
| ) | ||
| if field is None: | ||
| raise ValueError(f"muphys requires tracer '{name}' to be active in the TracerState") | ||
| return field |
There was a problem hiding this comment.
what I mean is that instead of calling _require 15-20 times in the code, you can just verify once at the setup of the component, or the first time it is called, that all tracers necessary are active by doing something like all(SPECIES in TracerState.active_fields()) (with proper syntax), makes sense?
| prepare the input fields for muphys from the prognostic state. This includes: | ||
| - binding the references for rho and q (the muphys input fields that are stored prognostically in the dycore state) | ||
| - diagnosing the muphys input fields that aren't stored prognostically (te, p) from the prognostic state. | ||
| """ |
There was a problem hiding this comment.
Instead of diagnosing all those diagnostic variables used in microphysics or each process, we need a better strategy. We need to check whether the process wants to use the most up-to-date variables or just use the values at the beginning of the time loop (after dycore).
For reference, in NWP, fast physics always uses the variables updated by preceding physics process, while slow physics uses the values after all fast physics processes have been called.
On the other hand, AES physics always uses the values at the beginning of each time loop and returns tendencies or change rate.
There was a problem hiding this comment.
we need to think about a better strategy here, in the second pr
| experiment: utils.MuphysExperiment, | ||
| ) -> None: | ||
| allocator = model_backends.get_allocator(backend_like) | ||
| inp = common.GraupelInput.load(filename=experiment.input_file, allocator=allocator) |
There was a problem hiding this comment.
Nothing wrong, It is just very unfortunate that inp stands for ice-nucleating particles in microphysics community. So, this variable name becomes confusing for microphysics researchers if they happen to look at the code.
| "horizontal_start": gtx.int32(0), | ||
| "horizontal_end": gtx.int32(self._ncells), | ||
| } | ||
| vertical_sizes = {"vertical_start": gtx.int32(0), "vertical_end": gtx.int32(self._nlev)} |
There was a problem hiding this comment.
usually in ICON or other models, there is k level above which microphysics is turned off (usually for manual performance tuning, may be for stability issues in some rare cases). In ICON, it is kstart_moist. How about passing in vertical params? Could you check whether we can remove this manual performance tuning so that we always run the microphysics for all levels?
There was a problem hiding this comment.
I will check this
| class MuphysConfig: | ||
| """Configuration for the muphys microphysics component.""" | ||
|
|
||
| qnc: float = 50.0e6 # cloud droplet number concentration [m^-3], matches Fortran cloud_num |
There was a problem hiding this comment.
Ignore my comment here. This is just for myself.
I am surprised that AES uses a fixed qnc. Check why.
| # diagnose theta_v = Tv/exner, mirroring ICON's phy2dyn coupling | ||
| # (mo_interface_iconam_aes.f90). rho is unchanged by the fast physics, so the | ||
| # exner/rho/theta_v trio stays EOS-consistent. | ||
| self._update_exner_and_theta_v( |
There was a problem hiding this comment.
Maybe not in this PR. We need to make this optional if we are following aes strategy. I remember that aes updates all prognostic variables at the very end of all physics calls.
|
When developing, you can test your changes on CSCS CI before merge with the You can pass options to override pipeline variables, for example:
Available options are:
For each option, See The Merging Once your PR is approved and ready for merging, add it to the merge queue. The Optional Tests To run benchmarks you can use:
For more detailed information please look at CI in the EXCLAIM universe. |
No description provided.