Implement PredictedParticlePropertyMicrophysics

[glwagner]

This PR takes a stab at implementing the most recent version of the "Predicted Particle Property" (P3) Microphysics in Breeze as contained in https://github.com/P3-microphysics/P3-microphysics and described by https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2024MS004404.

Some references for P3 (please let me know if there are others I should add):

• Morrison and Milbrandt 2015
• Morrison et al 2015
• Milbrandt and Morrison 2016
• Cholette et al 2019
• Paukert et al 2019
• Jouan et al 2020
• Milbrandt et al 2021
• Cholette et al 2023
• Morrison, Milbrandt, Cholette 2025
• Milbrandt, Morrison, Cholette 2025

github repo: https://github.com/P3-microphysics/P3-microphysics, which includes interfaces for CM1

[codecov]

Codecov Report

❌ Patch coverage is 48.15175% with 533 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
...ysics/PredictedParticleProperties/process_rates.jl	0.00%	271 Missing ⚠️
...hysics/PredictedParticleProperties/p3_interface.jl	3.44%	140 Missing ⚠️
...ysics/PredictedParticleProperties/lambda_solver.jl	59.45%	90 Missing ⚠️
...ophysics/PredictedParticleProperties/quadrature.jl	95.45%	8 Missing ⚠️
...cs/PredictedParticleProperties/cloud_properties.jl	0.00%	7 Missing ⚠️
...s/PredictedParticleProperties/size_distribution.jl	62.50%	6 Missing ⚠️
...ctedParticleProperties/cloud_droplet_properties.jl	87.50%	1 Missing ⚠️
...PredictedParticleProperties/ice_bulk_properties.jl	85.71%	1 Missing ⚠️
...sics/PredictedParticleProperties/ice_collection.jl	85.71%	1 Missing ⚠️
...sics/PredictedParticleProperties/ice_deposition.jl	85.71%	1 Missing ⚠️
... and 7 more

📢 Thoughts on this report? Let us know!

[giordano]

Here, and in many other places, all in this file:

┌ Error: failed to run `@example` block in src/microphysics/p3_size_distribution.md:204-226
│ ```@example p3_psd
│ fig = Figure(size=(600, 400))
│ ax = Axis(fig[1, 1],
│     xlabel = "Diameter D [mm]",
│     ylabel = "N'(D) [m⁻⁴]",
│     yscale = log10,
│     title = "Effect of Riming on Size Distribution\n(L = 10⁻⁴ kg/m³, N = 10⁵ m⁻³)")
│
│ L_ice = 1e-4
│ N_ice = 1e5
│
│ for (Ff, label, color) in [(0.0, "Fᶠ = 0 (unrimed)", :blue),
│                             (0.3, "Fᶠ = 0.3", :green),
│                             (0.6, "Fᶠ = 0.6", :orange)]
│     params = distribution_parameters(L_ice, N_ice, Ff, 500.0)
│     N_D = @. params.intercept * D_m^params.shape * exp(-params.slope * D_m)
│     lines!(ax, D_mm, N_D, label=label, color=color)
│ end
│
│ axislegend(ax, position=:rt)
│ ylims!(ax, 1e3, 1e12)
│ fig
│ ```
│   exception =
│    FieldError: type Breeze.Microphysics.PredictedParticleProperties.IceDistributionParameters has no field `intercept`, available fields: `N₀`, `λ`, `μ`

Not sure what this was meant to do.

[glwagner]

I suspect it is trying to verify an internal calculation with a manual calculation (I'd like the documentation to provide a sort of tutorial to the internals of the scheme. I suspect the scheme got changfed but this wasnt updated.

There are some pretty grave errors also in this PR like manual loops vs kernels. My plan is to get #394 merged and then return to this --- with #394 we can do an apples to apples comparison with the fortran and hopefully build some confidence (along with the docs) that this is doing the right thing

[giordano]

Can this be revised now?

[glwagner]

yeah, i'll work on this

[giordano]

Why this is using CPU specifically instead of default_arch?

[glwagner]

yeah thats wrong

[giordano]

The arch argument is actually unused (because calculation is always done on CPU), is that supposed to be an argument at all?

[glwagner]

I think what we have here is that we will have a parameter object PredictedParticlePropertyMicrophysics(). However, it will be slow or won't even work (?) unless we generate look-up tables. I'm trying to figure out an interface for that, but one idea is something like:

p3 = PredictedParticlePropertyMicrophysics()
tabulated_p3 = tabulate(p3, :ice_fall_speed, GPU())

we have to specify the architecture to indicate what array type to use for the tabulation.

This could also all be pushed under the hood. Basically the tabulation would only be done inside the model constructor, not by users.

But I am open to thoughts and ideas how to handle this...

(also I think this is properly just the device, since there nothing special that happens for arch::Distributed)

[giordano]

How long is the look-up table going to be? If it's under ~100 elements, a tuple may be good for all architectures (don't even need to allocate memory for it). Do you need to generate a specific table for different datatypes? That'd complicate things a little bit because you can't necessarily hardcode the table (although the conversion to a lower precision may probably be done at compile time with hardcoded coefficients for a table)

[glwagner]

ah interesting. because for a tuple, the look up is actually inlined at compile time. that is pretty cool! we should add this as a feature to TabulatedFunction in Oceananigans too.

To answer your question, I have no idea how long the lookup table needs to be, so let me figure that out.

[glwagner]

as for precision, I would guess that we can get away with low precision for the lookup since the main approximation is the linear interpolation between nodes.