From bfb9b6f3028f820031de9cdcfd1162c705b25c74 Mon Sep 17 00:00:00 2001 From: Julia Sloan Date: Tue, 21 Oct 2025 17:50:27 -0700 Subject: [PATCH] initial ClimaSeaIce infrastructure changes --- docs/src/fieldexchanger.md | 10 +- docs/src/fluxcalculator.md | 13 +- docs/src/interfacer.md | 15 ++ .../components/ocean/climaocean_helpers.jl | 129 ++++++++++++++++ .../components/ocean/oceananigans.jl | 138 +----------------- .../components/ocean/prescr_seaice.jl | 2 + src/FieldExchanger.jl | 17 ++- src/Interfacer.jl | 4 + test/field_exchanger_tests.jl | 9 +- 9 files changed, 188 insertions(+), 149 deletions(-) create mode 100644 experiments/ClimaEarth/components/ocean/climaocean_helpers.jl diff --git a/docs/src/fieldexchanger.md b/docs/src/fieldexchanger.md index e89e50c654..6be11ffa0e 100644 --- a/docs/src/fieldexchanger.md +++ b/docs/src/fieldexchanger.md @@ -12,7 +12,7 @@ The specific fields that are exchanged depend on the requirements of the compone The fields imported from the atmosphere to the coupler are specified by extending `FieldExchanger.import_atmos_fields!` The default `import_atmos_fields!` imports radiative fluxes, liquid precipitation, and snow precipitation. -The fields of a component model that get updated by the coupler are specified by extending `FieldExchanger.update_sim!` +The fields of a component model that get updated by the coupler are specified by extending `FieldExchanger.update_sim!`. The default `update_sim!` for an atmosphere model updates the direct and diffuse surface albedos, the surface temperature, and the turbulent fluxes. The default `update_sim!` for a surface model updates the air density, radiative fluxes, @@ -21,13 +21,17 @@ These updates are done via the `update_field!` function, which must be extended particular variable and component model. If an `update_field!` function is not defined for a particular component model, it will be ignored. +Note that turbulent fluxes are not updated in `update_sim!`, but rather via +`FluxCalculator.update_turbulent_fluxes!`, where fluxes are computed between +the atmosphere and each surface model. + ## FieldExchanger API ```@docs + ClimaCoupler.FieldExchanger.exchange! ClimaCoupler.FieldExchanger.update_sim! ClimaCoupler.FieldExchanger.step_model_sims! ClimaCoupler.FieldExchanger.update_surface_fractions! - ClimaCoupler.FieldExchanger.exchange! ClimaCoupler.FieldExchanger.set_caches! ``` @@ -35,4 +39,6 @@ If an `update_field!` function is not defined for a particular component model, ```@docs ClimaCoupler.FieldExchanger.combine_surfaces! + ClimaCoupler.FieldExchanger.resolve_ocean_ice_fractions! + ClimaCoupler.FieldExchanger.import_atmos_fields! ``` diff --git a/docs/src/fluxcalculator.md b/docs/src/fluxcalculator.md index 23c1a18840..eba49cc4b6 100644 --- a/docs/src/fluxcalculator.md +++ b/docs/src/fluxcalculator.md @@ -11,24 +11,23 @@ turbulent fluxes and ancillary quantities, such as the Obukhov length, using function is called at the end of each coupling step. All the quantities computed in `turbulent_fluxes!` are calculated -separately for each surface type using the +separately for each surface model using the [`FluxCalculator.compute_surface_fluxes!`](@ref) function. This function can be extended by component models if they need specific type of flux calculation, and a default is provided for models that can use the standard flux calculation. -The final result of `turbulent_fluxes!` is an area-weighted sum of -all the contributions of the various surfaces. The default method of [`FluxCalculator.compute_surface_fluxes!`](@ref), in turn, calls [`FluxCalculator.get_surface_fluxes`](@ref). This function uses a thermal state obtained by using the model surface temperature, extrapolates atmospheric density adiabatically to the surface, and with the surface humidity (if available, if not, assuming a saturation specific humidity for liquid phase). -`compute_surface_fluxes!` also updates the component internal fluxes fields with +`compute_surface_fluxes!` also updates the component internal fluxes fields via [`FluxCalculator.update_turbulent_fluxes!`](@ref), and adds the area-weighted contribution from this component model to the `CoupledSimulation` fluxes fields. -Any extension of this function for a particular surface model is also expected -to update both the component models' internal fluxes and the CoupledSimulation -object's fluxes fields. + +Any extension of `FluxCalculator.compute_surface_fluxes!` for a particular +surface model is also expected to update both the component models' internal +fluxes and the CoupledSimulation object's fluxes fields. [`FluxCalculator.compute_surface_fluxes!`](@ref) sets: - the flux of momenta, `F_turb_ρτxz`, `F_turb_ρτyz`; diff --git a/docs/src/interfacer.md b/docs/src/interfacer.md index 400d87a4ea..3b207712d4 100644 --- a/docs/src/interfacer.md +++ b/docs/src/interfacer.md @@ -225,11 +225,26 @@ for the following properties: | `surface_direct albedo` | bulk direct surface albedo | | | `surface_diffuse albedo` | bulk diffuse surface albedo | | | `surface_temperature` | surface temperature | K | +| `ice_concentration` | sea ice concentration (*sea ice models only*) | | !!! note `area_fraction` is expected to be defined on the boundary space of the simulation, while all other fields will likely be on the simulation's own space. +!!! note "Sea ice concentration vs. area fraction" +Sea ice models are expected to provide both `area_fraction` and `ice_concentration`. +This may seem redundant, but there are subtle differences between the two. +`ice_concentration` is internal to the ice model and may be determined +via a prognostic variable, prescribed data, etc. `area_fraction` is defined +at the coupler level and may follow some constraints that don't apply to `ice_concentration`. +For example, we require that surface model area fractions sum to 1 at all points; +this constraint is enforced for `area_fraction`, but not for `ice_concentration`. +Additionally, since `area_fraction` is a coupler-defined concept, it is defined on +the coupler boundary space, whereas `ice_concentration` exists on the model's space. +Generally, `ice_concentration` and `area_fraction` should largely agree, +with differences only arising from `area_fraction` corrections and +the fields existing on different spaces. + - `update_field!(::SurfaceModelSimulation, ::Val{property}, field)`: A function to update the value of property in the component model simulation, using the values in `field` passed from the coupler diff --git a/experiments/ClimaEarth/components/ocean/climaocean_helpers.jl b/experiments/ClimaEarth/components/ocean/climaocean_helpers.jl new file mode 100644 index 0000000000..5927ab5327 --- /dev/null +++ b/experiments/ClimaEarth/components/ocean/climaocean_helpers.jl @@ -0,0 +1,129 @@ +""" + to_node(pt::CC.Geometry.LatLongPoint) + +Transform `LatLongPoint` into a tuple (long, lat, 0), where the 0 is needed because we only +care about the surface. +""" +@inline to_node(pt::CC.Geometry.LatLongPoint) = pt.long, pt.lat, zero(pt.lat) +# This next one is needed if we have "LevelGrid" +@inline to_node(pt::CC.Geometry.LatLongZPoint) = pt.long, pt.lat, zero(pt.lat) + +""" + map_interpolate(points, oc_field::OC.Field) + +Interpolate the given 3D field onto the target points. + +If the underlying grid does not contain a given point, return 0 instead. + +Note: `map_interpolate` does not support interpolation from `Field`s defined on +`OrthogononalSphericalShellGrids` such as the `TripolarGrid`. + +TODO: Use a non-allocating version of this function (simply replace `map` with `map!`) +""" +function map_interpolate(points, oc_field::OC.Field) + loc = map(L -> L(), OC.Fields.location(oc_field)) + grid = oc_field.grid + data = oc_field.data + + # TODO: There has to be a better way + min_lat, max_lat = extrema(OC.φnodes(grid, OC.Center(), OC.Center(), OC.Center())) + + map(points) do pt + FT = eltype(pt) + + # The oceananigans grid does not cover the entire globe, so we should not + # interpolate outside of its latitude bounds. Instead we return 0 + min_lat < pt.lat < max_lat || return FT(0) + + fᵢ = OC.Fields.interpolate(to_node(pt), data, loc, grid) + convert(FT, fᵢ)::FT + end +end + +""" + surface_flux(f::OC.AbstractField) + +Extract the top boundary conditions for the given field. +""" +function surface_flux(f::OC.AbstractField) + top_bc = f.boundary_conditions.top + if top_bc isa OC.BoundaryCondition{<:OC.BoundaryConditions.Flux} + return top_bc.condition + else + return nothing + end +end + +function Interfacer.remap(field::OC.Field, target_space) + return map_interpolate(CC.Fields.coordinate_field(target_space), field) +end + +function Interfacer.remap(operation::OC.AbstractOperations.AbstractOperation, target_space) + evaluated_field = OC.Field(operation) + OC.compute!(evaluated_field) + return Interfacer.remap(evaluated_field, target_space) +end + +""" + set_from_extrinsic_vector!(vector, grid, u_cc, v_cc) + +Given the extrinsic vector components `u_cc` and `v_cc` as `Center, Center` +fields, rotate them onto the target grid and remap to `Face, Center` and +`Center, Face` fields, respectively. +""" +function set_from_extrinsic_vector!(vector, grid, u_cc, v_cc) + arch = OC.Architectures.architecture(grid) + + # Rotate vector components onto the grid + OC.Utils.launch!(arch, grid, :xy, _rotate_vector!, u_cc, v_cc, grid) + + # Fill halo regions with the rotated vector components so we can use them to interpolate + OC.fill_halo_regions!(u_cc) + OC.fill_halo_regions!(v_cc) + + # Interpolate the vector components to face/center and center/face respectively + OC.Utils.launch!( + arch, + grid, + :xy, + _interpolate_vector!, + vector.u, + vector.v, + grid, + u_cc, + v_cc, + ) + return nothing +end + +""" + _rotate_vector!(τx, τy, grid) + +Rotate the velocities from the extrinsic coordinate system to the intrinsic +coordinate system. +""" +@kernel function _rotate_vector!(τx, τy, grid) + # Use `k = 1` to index into the reduced Fields + i, j = @index(Global, NTuple) + # Rotate u, v from extrinsic to intrinsic coordinate system + τxr, τyr = OC.Operators.intrinsic_vector(i, j, 1, grid, τx, τy) + @inbounds begin + τx[i, j, 1] = τxr + τy[i, j, 1] = τyr + end +end + +""" + _interpolate_vector!(τx, τy, grid, τx_cc, τy_cc) + +Interpolate the input fluxes `τx_cc` and `τy_cc`, which are Center/Center +Fields to Face/Center and Center/Face coordinates, respectively. +""" +@kernel function _interpolate_vector!(τx, τy, grid, τx_cc, τy_cc) + # Use `k = 1` to index into the reduced Fields + i, j = @index(Global, NTuple) + @inbounds begin + τx[i, j, 1] = OC.Operators.ℑxᶠᵃᵃ(i, j, 1, grid, τx_cc) + τy[i, j, 1] = OC.Operators.ℑyᵃᶠᵃ(i, j, 1, grid, τy_cc) + end +end diff --git a/experiments/ClimaEarth/components/ocean/oceananigans.jl b/experiments/ClimaEarth/components/ocean/oceananigans.jl index cb43859427..b665fb981b 100644 --- a/experiments/ClimaEarth/components/ocean/oceananigans.jl +++ b/experiments/ClimaEarth/components/ocean/oceananigans.jl @@ -7,6 +7,8 @@ import Thermodynamics as TD import ClimaOcean.EN4: download_dataset using KernelAbstractions: @kernel, @index, @inbounds +include("climaocean_helpers.jl") + """ OceananigansSimulation{SIM, A, OPROP, REMAP} @@ -62,10 +64,6 @@ function OceananigansSimulation( z = OC.ExponentialDiscretization(Nz, -depth, 0; scale = 0.85 * depth) # Regular LatLong because we know how to do interpolation there - - # TODO: When moving to TripolarGrid, note that we need to be careful about - # ensuring the coordinate systems align (ie, rotate vectors on the OC grid) - underlying_grid = OC.LatitudeLongitudeGrid( arch; size = resolution_points, @@ -222,70 +220,6 @@ end Interfacer.step!(sim::OceananigansSimulation, t) = OC.time_step!(sim.ocean, float(t) - sim.ocean.model.clock.time) -# We always want the surface, so we always set zero(pt.lat) for z -""" - to_node(pt::CA.ClimaCore.Geometry.LatLongPoint) - -Transform `LatLongPoint` into a tuple (long, lat, 0), where the 0 is needed because we only -care about the surface. -""" -@inline to_node(pt::CA.ClimaCore.Geometry.LatLongPoint) = pt.long, pt.lat, zero(pt.lat) -# This next one is needed if we have "LevelGrid" -@inline to_node(pt::CA.ClimaCore.Geometry.LatLongZPoint) = pt.long, pt.lat, zero(pt.lat) - -""" - map_interpolate(points, oc_field::OC.Field) - -Interpolate the given 3D field onto the target points. - -If the underlying grid does not contain a given point, return 0 instead. - -TODO: Use a non-allocating version of this function (simply replace `map` with `map!`) -""" -function map_interpolate(points, oc_field::OC.Field) - loc = map(L -> L(), OC.Fields.location(oc_field)) - grid = oc_field.grid - data = oc_field.data - - # TODO: There has to be a better way - min_lat, max_lat = extrema(OC.φnodes(grid, OC.Center(), OC.Center(), OC.Center())) - - map(points) do pt - FT = eltype(pt) - - # The oceananigans grid does not cover the entire globe, so we should not - # interpolate outside of its latitude bounds. Instead we return 0 - min_lat < pt.lat < max_lat || return FT(0) - - fᵢ = OC.Fields.interpolate(to_node(pt), data, loc, grid) - convert(FT, fᵢ)::FT - end -end - -""" - surface_flux(f::OC.AbstractField) - -Extract the top boundary conditions for the given field. -""" -function surface_flux(f::OC.AbstractField) - top_bc = f.boundary_conditions.top - if top_bc isa OC.BoundaryCondition{<:OC.BoundaryConditions.Flux} - return top_bc.condition - else - return nothing - end -end - -function Interfacer.remap(field::OC.Field, target_space) - return map_interpolate(CC.Fields.coordinate_field(target_space), field) -end - -function Interfacer.remap(operation::OC.AbstractOperations.AbstractOperation, target_space) - evaluated_field = OC.Field(operation) - OC.compute!(evaluated_field) - return Interfacer.remap(evaluated_field, target_space) -end - Interfacer.get_field(sim::OceananigansSimulation, ::Val{:area_fraction}) = sim.area_fraction # TODO: Better values for this @@ -323,8 +257,6 @@ and Oceananigans represents moisture moving from atmosphere to ocean as a positi so a sign change is needed when we convert from moisture to salinity flux. """ function FluxCalculator.update_turbulent_fluxes!(sim::OceananigansSimulation, fields) - # Only LatitudeLongitudeGrid are supported because otherwise we have to rotate the vectors - (; F_lh, F_sh, F_turb_ρτxz, F_turb_ρτyz, F_turb_moisture) = fields grid = sim.ocean.model.grid @@ -349,7 +281,7 @@ function FluxCalculator.update_turbulent_fluxes!(sim::OceananigansSimulation, fi # Set the momentum flux BCs at the correct locations using the remapped scratch fields oc_flux_u = surface_flux(sim.ocean.model.velocities.u) oc_flux_v = surface_flux(sim.ocean.model.velocities.v) - set_from_extrinsic_vectors!( + set_from_extrinsic_vector!( (; u = oc_flux_u, v = oc_flux_v), grid, F_turb_ρτxz_cc, @@ -390,70 +322,6 @@ function FluxCalculator.update_turbulent_fluxes!(sim::OceananigansSimulation, fi return nothing end -""" - set_from_extrinsic_vectors!(vectors, grid, u_cc, v_cc) - -Given the extrinsic vector components `u_cc` and `v_cc` as `Center, Center` -fields, rotate them onto the target grid and remap to `Face, Center` and -`Center, Face` fields, respectively. -""" -function set_from_extrinsic_vectors!(vectors, grid, u_cc, v_cc) - arch = grid.architecture - - # Rotate vectors onto the grid - OC.Utils.launch!(arch, grid, :xy, _rotate_velocities!, u_cc, v_cc, grid) - - # Fill halo regions with the rotated vectors so we can use them to interpolate - OC.fill_halo_regions!(u_cc) - OC.fill_halo_regions!(v_cc) - - # Interpolate the vectors to face/center and center/face respectively - OC.Utils.launch!( - arch, - grid, - :xy, - _interpolate_velocities!, - vectors.u, - vectors.v, - grid, - u_cc, - v_cc, - ) - return nothing -end - -""" - _rotate_velocities!(u, v, grid) - -Rotate the velocities from the extrinsic coordinate system to the intrinsic -coordinate system. -""" -@kernel function _rotate_velocities!(u, v, grid) - # Use `k = 1` to index into the reduced Fields - i, j = @index(Global, NTuple) - # Rotate u, v from extrinsic to intrinsic coordinate system - ur, vr = OC.Operators.intrinsic_vector(i, j, 1, grid, u, v) - @inbounds begin - u[i, j, 1] = ur - v[i, j, 1] = vr - end -end - -""" - _interpolate_velocities!(u, v, grid, u_cc, v_cc) - -Interpolate the input velocities `u_cc` and `v_cc`, which are Center/Center -Fields to Face/Center and Center/Face coordinates, respectively. -""" -@kernel function _interpolate_velocities!(u, v, grid, u_cc, v_cc) - # Use `k = 1` to index into the reduced Fields - i, j = @index(Global, NTuple) - @inbounds begin - u[i, j, 1] = OC.Operators.ℑxyᶠᶜᵃ(i, j, 1, grid, u_cc) - v[i, j, 1] = OC.Operators.ℑxyᶜᶠᵃ(i, j, 1, grid, v_cc) - end -end - function Interfacer.update_field!(sim::OceananigansSimulation, ::Val{:area_fraction}, field) sim.area_fraction .= field return nothing diff --git a/experiments/ClimaEarth/components/ocean/prescr_seaice.jl b/experiments/ClimaEarth/components/ocean/prescr_seaice.jl index d45d12fb14..7ba8a8aad9 100644 --- a/experiments/ClimaEarth/components/ocean/prescr_seaice.jl +++ b/experiments/ClimaEarth/components/ocean/prescr_seaice.jl @@ -167,6 +167,8 @@ end # extensions required by Interfacer Interfacer.get_field(sim::PrescribedIceSimulation, ::Val{:area_fraction}) = sim.integrator.p.area_fraction +Interfacer.get_field(sim::PrescribedIceSimulation, ::Val{:ice_concentration}) = + sim.integrator.p.area_fraction Interfacer.get_field(sim::PrescribedIceSimulation, ::Val{:emissivity}) = sim.integrator.p.params.ϵ Interfacer.get_field(sim::PrescribedIceSimulation, ::Val{:roughness_buoyancy}) = diff --git a/src/FieldExchanger.jl b/src/FieldExchanger.jl index f8a9508933..3e16043030 100644 --- a/src/FieldExchanger.jl +++ b/src/FieldExchanger.jl @@ -10,7 +10,13 @@ import ClimaCore as CC import ..Interfacer, ..FluxCalculator, ..Utilities -export update_sim!, update_model_sims!, step_model_sims!, exchange!, set_caches! +export update_sim!, + update_model_sims!, + step_model_sims!, + exchange!, + set_caches!, + update_surface_fractions!, + resolve_ocean_ice_fractions! """ update_surface_fractions!(cs::Interfacer.CoupledSimulation) @@ -39,12 +45,14 @@ function update_surface_fractions!(cs::Interfacer.CoupledSimulation) # ice and ocean fractions are dynamic if haskey(cs.model_sims, :ice_sim) ice_sim = cs.model_sims.ice_sim - ice_fraction_before = Interfacer.get_field(ice_sim, Val(:area_fraction)) + Interfacer.get_field!(cs.fields.scalar_temp1, ice_sim, Val(:ice_concentration)) + ice_concentration = cs.fields.scalar_temp1 + # max needed to avoid Float32 errors (see issue #271; Heisenbug on HPC) Interfacer.update_field!( ice_sim, Val(:area_fraction), - max.(min.(ice_fraction_before, FT(1) .- land_fraction), FT(0)), + max.(min.(ice_concentration, FT(1) .- land_fraction), FT(0)), ) ice_fraction = Interfacer.get_field(ice_sim, Val(:area_fraction)) else @@ -198,7 +206,8 @@ end """ update_sim!(sim::SurfaceModelSimulation, csf) -Updates the surface component model cache with the current coupler fields besides turbulent fluxes. +Updates the surface component model cache with the current coupler fields +*besides turbulent fluxes*, which are updated in `update_turbulent_fluxes`. # Arguments - `sim`: [Interfacer.SurfaceModelSimulation] containing a surface model simulation object. diff --git a/src/Interfacer.jl b/src/Interfacer.jl index a30eda85d5..d7049fa8eb 100644 --- a/src/Interfacer.jl +++ b/src/Interfacer.jl @@ -221,6 +221,10 @@ get_field( }, ) = get_field_error(sim, val) +# Sea ice models need to provide ice concentration +get_field(sim::SeaIceModelSimulation, val::Val{:ice_concentration}) = + get_field_error(sim, val) + """ get_field(sim::ComponentModelSimulation, val::Val) diff --git a/test/field_exchanger_tests.jl b/test/field_exchanger_tests.jl index 7a7896490d..1ab0fceaf2 100644 --- a/test/field_exchanger_tests.jl +++ b/test/field_exchanger_tests.jl @@ -64,6 +64,7 @@ struct DummyStub{C} <: Interfacer.SurfaceModelSimulation cache::C end Interfacer.get_field(sim::DummyStub, ::Val{:area_fraction}) = sim.cache.area_fraction +Interfacer.get_field(sim::DummyStub, ::Val{:ice_concentration}) = sim.cache.area_fraction function Interfacer.update_field!( sim::DummyStub, ::Val{:area_fraction}, @@ -195,10 +196,16 @@ for FT in (Float32, Float64) # Construct ice fraction of 0s. During the update, the first half should be set to 0.5 ocean_d = CC.Fields.zeros(test_space) + # Fill in only the necessary parts of the simulation + coupler_fields = Interfacer.init_coupler_fields( + FT, + Interfacer.default_coupler_fields(), + test_space, + ) cs = Interfacer.CoupledSimulation{FT}( nothing, # dates - ones(test_space), # fields + coupler_fields, # fields nothing, # conservation_checks (Int(0), Int(1000)), # tspan Int(200), # Δt_cpl