Skip to content

Conversation

@glwagner
Copy link
Member

I am expecting the GPU+Float32 to fail here, but not sure.

I am expecting the GPU+Float32 to fail here, but not sure.
@giordano
Copy link
Collaborator

Someone broke isapprox on Fields, I wonder who could have been 😩

@glwagner
Copy link
Member Author

Someone broke isapprox on Fields, I wonder who could have been 😩

ah no its not your fault! it looks like it is not defined for array/scalar comparison, eg

julia> a = ones(10, 10)
10×10 Matrix{Float64}:
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0

julia> a  1
ERROR: MethodError: no method matching isapprox(::Matrix{Float64}, ::Int64)
The function `isapprox` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  isapprox(::Missing, ::Any; kwargs...)
   @ Base missing.jl:90
  isapprox(::Integer, ::Integer; atol, rtol, nans, norm)
   @ Base floatfuncs.jl:229
  isapprox(::Number, ::Number; atol, rtol, nans, norm)
   @ Base floatfuncs.jl:220
  ...

Stacktrace:
 [1] top-level scope
   @ REPL[37]:1

@glwagner
Copy link
Member Author

We do have an abstraction called ConstantField(c) that might work here, I will try it.

@giordano
Copy link
Collaborator

all(isapprox(1), c)? But it requires scalar indexing.

@glwagner
Copy link
Member Author

all(isapprox(1), c)? But it requires scalar indexing.

We have an abstraction that represents a constant broadcast to an array (but without allocating memory for it). I guess this doesn't exist in julia base (or maybe it does somehow) but we can use it here, its called ConstantField

@glwagner
Copy link
Member Author

We also need scalar indexing for ordinary approx right?

Or maybe I am confused. Do we hit the same codepath into GPUArrays with Field magically?

@glwagner
Copy link
Member Author

Well, this works

julia> grid = RectilinearGrid(CPU(); size=(8, 8, 8), x=(0, 1_000), y=(0, 1_000), z=(0, 1_000))
8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── Periodic x  [0.0, 1000.0) regularly spaced with Δx=125.0
├── Periodic y  [0.0, 1000.0) regularly spaced with Δy=125.0
└── Bounded  z  [0.0, 1000.0] regularly spaced with Δz=125.0

julia> c = CenterField(grid)
8×8×8 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
└── data: 14×14×14 OffsetArray(::Array{Float64, 3}, -2:11, -2:11, -2:11) with eltype Float64 with indices -2:11×-2:11×-2:11
    └── max=0.0, min=0.0, mean=0.0

julia> c  Field{Nothing, Nothing, Nothing}(grid)
true

@glwagner
Copy link
Member Author

ConstantField does not work because it does not have a grid and therefore no notion of size. That is pretty easily changed though

@giordano
Copy link
Collaborator

We also need scalar indexing for ordinary approx right?

Not after CliMA/Oceananigans.jl#4924 🙂 (well, there's a call to norm which internally does a scalar indexing, but it's on a one-element array after a reduction, not on the whole field)

@glwagner
Copy link
Member Author

We also need scalar indexing for ordinary approx right?

Not after CliMA/Oceananigans.jl#4924 🙂 (well, there's a call to norm which internally does a scalar indexing, but it's on a one-element array after a reduction, not on the whole field)

ok so we also need to remove all the allowscalar!

@giordano
Copy link
Collaborator

ok so we also need to remove all the allowscalar!

How was it? Great minds...? 5a8233d

@glwagner
Copy link
Member Author

ok so we also need to remove all the allowscalar!

How was it? Great minds...? 5a8233d

😂

@navidcy navidcy changed the title Add tests for set! Add tests for set! Nov 13, 2025
Comment on lines 20 to 21
abstract type AbstractEquilibrium end

Copy link
Collaborator

@giordano giordano Nov 29, 2025

Choose a reason for hiding this comment

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

I'm assuming that this wasn't meant to be removed: it disappeared in a merge commit, which feels unlikely to be done on purpose? Edit: tentatively changed in 53c7a56

Comment on lines 133 to 146
g = thermo.gravitational_acceleration
z = znode(i, j, k, grid, c, c, c)
Copy link
Collaborator

@giordano giordano Nov 29, 2025

Choose a reason for hiding this comment

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

These lines also got removed in another merged commit: 4e642d1. Was that done on purpose? Edit: tentatively changed in 53c7a56

Comment on lines 188 to 202
@inline function adjust_state(𝒰₀, T, thermo, equilibrium)
@inline function adjust_state(𝒰₀::AbstractThermodynamicState{FT}, T::FT,
thermo::ThermodynamicConstants{FT},
equilibrium::MixedPhaseEquilibrium{FT}) where FT
Copy link
Collaborator

Choose a reason for hiding this comment

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

This stricter signature caused

MethodError: no method matching adjust_state(::Breeze.Thermodynamics.PotentialTemperatureState{Float64}, ::Float64, ::ThermodynamicConstants{Float64, CondensedPhase{Float64}, CondensedPhase{Float64}}, ::WarmPhaseEquilibrium)
The function `adjust_state` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  adjust_state(::Breeze.Thermodynamics.AbstractThermodynamicState{FT}, ::FT, ::ThermodynamicConstants{FT}, !Matched::MixedPhaseEquilibrium{FT}) where FT
   @ Breeze ~/work/Breeze.jl/Breeze.jl/src/Microphysics/saturation_adjustment.jl:200

Do we need to define a different method for WarmPhaseEquilibrium? Or the last argument can be AbstractEquilibrium?

Copy link
Collaborator

Choose a reason for hiding this comment

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

diff --git a/src/Microphysics/saturation_adjustment.jl b/src/Microphysics/saturation_adjustment.jl
index 54bd973..393258f 100644
--- a/src/Microphysics/saturation_adjustment.jl
+++ b/src/Microphysics/saturation_adjustment.jl
@@ -199,7 +199,7 @@ end
 
 @inline function adjust_state(𝒰₀::AbstractThermodynamicState{FT}, T::FT,
                               thermo::ThermodynamicConstants{FT},
-                              equilibrium::MixedPhaseEquilibrium{FT}) where FT
+                              equilibrium::AbstractEquilibrium) where FT
 
     pᵣ = 𝒰₀.reference_pressure
     qᵗ = total_specific_moisture(𝒰₀)
@@ -210,7 +210,7 @@ end
 
 @inline function saturation_adjustment_residual(T::FT, 𝒰₀::AbstractThermodynamicState{FT},
                                                 thermo::ThermodynamicConstants{FT},
-                                                equilibrium::MixedPhaseEquilibrium{FT}) where FT
+                                                equilibrium::AbstractEquilibrium) where FT
 
     𝒰₁ = adjust_state(𝒰₀, T, thermo, equilibrium)
     T₁ = temperature(𝒰₁, thermo)

would fix that issue and make all the saturation_adjustment tests pass, but I'm not sure this is the intended design, so I'm leaving that to you.

Copy link
Collaborator

Choose a reason for hiding this comment

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

I tentatively pushed 705c0b9 to implement the above, let me know if you didn't want it.

"""
struct WarmPhaseEquilibrium <: AbstractEquilibrium end
@inline equilibrated_surface(::WarmPhaseEquilibrium, T) = PlanarLiquidSurface()
convert_equilibrium(FT, ::WarmPhaseEquilibrium) = WarmPhaseEquilibrium()
Copy link
Collaborator

Choose a reason for hiding this comment

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

Is this method necessary? I presume this is already covered by

convert_equilibrium(FT, equil) = equil # fallback
?

@giordano
Copy link
Collaborator

giordano commented Dec 3, 2025

I should have resolved all merge conflicts, but now there are some tests failures in atmosphere_model_unit.

This PR has had a few conflicts with other PRs over the last few days, it'd be good to get this in soon if it's still desired, otherwise it'll keep conflicting and breaking forever. Edit: the failures in atmosphere_model_unit were already in aa22c7b, at least I didn't break anything today!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants