From 20bdbf0c69e27f221182b7bb1a6cffae46fcb5a9 Mon Sep 17 00:00:00 2001 From: kmdeck Date: Tue, 2 Jul 2024 11:19:54 -0700 Subject: [PATCH] use method in Domains only --- .../standalone/Soil/richards_runoff.jl | 6 +- src/shared_utilities/Domains.jl | 32 ++++++ src/standalone/Soil/Runoff/Runoff.jl | 11 +- src/standalone/Soil/boundary_conditions.jl | 101 +++--------------- src/standalone/Soil/energy_hydrology.jl | 10 +- 5 files changed, 57 insertions(+), 103 deletions(-) diff --git a/experiments/standalone/Soil/richards_runoff.jl b/experiments/standalone/Soil/richards_runoff.jl index 70821caae8..f2a3d46672 100644 --- a/experiments/standalone/Soil/richards_runoff.jl +++ b/experiments/standalone/Soil/richards_runoff.jl @@ -321,9 +321,8 @@ if context.device isa ClimaComms.CPUSingleThreaded θ_sfc_end = ClimaCore.Remapping.interpolate( remapper, - ClimaLand.Soil.get_top_surface_field( + ClimaLand.Domains.top_center_to_surface( oceans_to_zero.(field_to_error.(sol.u[end].soil.ϑ_l), mask), - surface_space, ), ) @@ -340,12 +339,11 @@ if context.device isa ClimaComms.CPUSingleThreaded Δθ_sfc = ClimaCore.Remapping.interpolate( remapper, - ClimaLand.Soil.get_top_surface_field( + ClimaLand.Domains.top_center_to_surface( oceans_to_zero.( field_to_error.(sol.u[end].soil.ϑ_l .- sol.u[1].soil.ϑ_l), mask, ), - surface_space, ), ) ax2 = Axis(fig[1, 3], xlabel = "Longitude", title = "θ_sfc Δ") diff --git a/src/shared_utilities/Domains.jl b/src/shared_utilities/Domains.jl index e9f023f64b..7343843f98 100644 --- a/src/shared_utilities/Domains.jl +++ b/src/shared_utilities/Domains.jl @@ -717,6 +717,37 @@ function linear_interpolation_to_surface!(sfc_field, center_field, z, Δz_top) @. (f1 - f2) / (z1 - z2) * (Δz_top + z1 - z2) + f2 end +""" + bottom_center_to_surface(center_field::ClimaCore.Fields.Field) + +Creates and returns a ClimaCore.Fields.Field defined on the space +corresponding to the surface of the space on which `center_field` +is defined, with values equal to the those at the level of the bottom +center. + +For example, given a `center_field` defined on 1D center finite difference space, +this would return a field defined on the Point space of the surface of +the column. The value would be the value of the oroginal `center_field` +at the bottommost location. Given a `center_field` defined on a 3D +extruded center finite difference space, this would return a 2D field +corresponding to the surface, with values equal to the bottommost level. +""" +function bottom_center_to_surface(center_field::ClimaCore.Fields.Field) + center_space = axes(center_field) + surface_space = obtain_surface_space(center_space) + return ClimaCore.Fields.Field( + ClimaCore.Fields.field_values(ClimaCore.Fields.level(center_field, 1)), + surface_space, + ) +end + +""" + bottom_center_to_surface(val) + +When `val` is a scalar (e.g. a single float or struct), returns `val`. +""" +bottom_center_to_surface(val) = val + """ get_Δz(z::ClimaCore.Fields.Field) @@ -792,6 +823,7 @@ export coordinates, obtain_face_space, obtain_surface_space, top_center_to_surface, + bottom_center_to_surface, top_face_to_surface, obtain_surface_domain, linear_interpolation_to_surface!, diff --git a/src/standalone/Soil/Runoff/Runoff.jl b/src/standalone/Soil/Runoff/Runoff.jl index d57d06d042..69f3a1f064 100644 --- a/src/standalone/Soil/Runoff/Runoff.jl +++ b/src/standalone/Soil/Runoff/Runoff.jl @@ -230,12 +230,8 @@ Computes the soil infiltration capacity on the surface space Currently approximates i_c = -K_sat at the surface. """ function soil_infiltration_capacity(model::RichardsModel, Y, p) - surface_space = model.domain.space.surface @. p.soil.subsfc_scratch = -1 * model.parameters.K_sat - return ClimaLand.Soil.get_top_surface_field( - p.soil.subsfc_scratch, - surface_space, - ) + return ClimaLand.Domains.top_center_to_surface(p.soil.subsfc_scratch) end """ @@ -258,10 +254,7 @@ function soil_infiltration_capacity(model::EnergyHydrology, Y, p) Ω, ) * ClimaLand.Soil.viscosity_factor(p.soil.T, γ, γT_ref) - return ClimaLand.Soil.get_top_surface_field( - p.soil.subsfc_scratch, - surface_space, - ) + return ClimaLand.Domains.top_center_to_surface(p.soil.subsfc_scratch) end diff --git a/src/standalone/Soil/boundary_conditions.jl b/src/standalone/Soil/boundary_conditions.jl index 155cbc2bb2..31524346c5 100644 --- a/src/standalone/Soil/boundary_conditions.jl +++ b/src/standalone/Soil/boundary_conditions.jl @@ -16,71 +16,6 @@ export TemperatureStateBC, RichardsAtmosDrivenFluxBC, WaterHeatBC -# Helper functions -""" - get_top_surface_field( - center_field::ClimaCore.Fields.Field, - surface_space, - ) - -A helper function to extract the top level of a center field and -cast it onto the surface face space. -""" -function get_top_surface_field( - center_field::ClimaCore.Fields.Field, - surface_space, -) - nz = Spaces.nlevels(axes(center_field)) - return Fields.Field( - Fields.field_values(Fields.level(center_field, nz)), - surface_space, - ) -end - -""" - get_top_surface_field( - center_val, - _, - ) - -A helper function for the case where we use `get_top_surface_field` on -a parameter that is a scalar rather than a field. Returns the scalar value. -""" -function get_top_surface_field(center_val, _) - return center_val -end - -""" - get_bottom_surface_field( - center_field::ClimaCore.Fields.Field, - bottom_space, - ) - -A helper function to extract the bottom level of a center field and -cast it onto the bottom face space. -""" -function get_bottom_surface_field( - center_field::ClimaCore.Fields.Field, - bottom_space, -) - return Fields.Field( - Fields.field_values(Fields.level(center_field, 1)), - bottom_space, - ) -end - -""" - get_bottom_surface_field( - center_val, - _, - ) - -A helper function for the case where we use `get_bottom_surface_field` on -a parameter that is a scalar rather than a field. Returns the scalar value. -""" -function get_bottom_surface_field(center_val, _) - return center_val -end # New BC type for Richards Equation (AbstractWaterBC) """ @@ -339,19 +274,18 @@ function boundary_flux( t, )::ClimaCore.Fields.Field FT = eltype(Δz) - surface_space = axes(Δz) # First extract the value of the top layer of the pressure head # and cast onto the face space - ψ_c = get_top_surface_field(p.soil.ψ, surface_space) + ψ_c = ClimaLand.Domains.top_center_to_surface(p.soil.ψ) # Calculate pressure head using boundary condition on ϑ_l = θ_bc # We first need to extract the parameters of the soil in the top layer # Again, we need to cast them onto the face space (; hydrology_cm, θ_r, ν, S_s) = model.parameters - hcm_bc = get_top_surface_field(hydrology_cm, surface_space) - θ_r_bc = get_top_surface_field(θ_r, surface_space) - ν_bc = get_top_surface_field(ν, surface_space) - S_s_bc = get_top_surface_field(S_s, surface_space) + hcm_bc = ClimaLand.Domains.top_center_to_surface(hydrology_cm) + θ_r_bc = ClimaLand.Domains.top_center_to_surface(θ_r) + ν_bc = ClimaLand.Domains.top_center_to_surface(ν) + S_s_bc = ClimaLand.Domains.top_center_to_surface(S_s) θ_bc = FT.(rre_bc.bc(p, t)) ψ_bc = @. pressure_head(hcm_bc, θ_r_bc, θ_bc, ν_bc, S_s_bc) @@ -361,7 +295,7 @@ function boundary_flux( # currently we approximate this as equal to the center value at the top layer (K_c) # More accurate would be to compute the mean between K_c and K evaluated at the boundary # condition. - K_eff = get_top_surface_field(p.soil.K, surface_space) + K_eff = ClimaLand.Domains.top_center_to_surface(p.soil.K) # Pass in (ψ_bc .+ Δz) as to account for contribution of gravity (∂(ψ+z)/∂z return ClimaLand.diffusive_flux(K_eff, ψ_bc .+ Δz, ψ_c, Δz) @@ -390,20 +324,19 @@ function boundary_flux( t, )::ClimaCore.Fields.Field FT = eltype(Δz) - surface_space = axes(Δz) # First extract the value of the bottom layer of the pressure head # and cast onto the face space - ψ_c = get_bottom_surface_field(p.soil.ψ, surface_space) + ψ_c = ClimaLand.Domains.bottom_center_to_surface(p.soil.ψ) # Calculate pressure head using boundary condition on ϑ_l = θ_bc # We first need to extract the parameters of the soil in the bottom layer # Again, we need to cast them onto the face space (; hydrology_cm, θ_r, ν, S_s) = model.parameters - hcm_bc = get_bottom_surface_field(hydrology_cm, surface_space) - θ_r_bc = get_bottom_surface_field(θ_r, surface_space) - ν_bc = get_bottom_surface_field(ν, surface_space) - S_s_bc = get_bottom_surface_field(S_s, surface_space) + hcm_bc = ClimaLand.Domains.bottom_center_to_surface(hydrology_cm) + θ_r_bc = ClimaLand.Domains.bottom_center_to_surface(θ_r) + ν_bc = ClimaLand.Domains.bottom_center_to_surface(ν) + S_s_bc = ClimaLand.Domains.bottom_center_to_surface(S_s) θ_bc = FT.(rre_bc.bc(p, t)) ψ_bc = @. pressure_head(hcm_bc, θ_r_bc, θ_bc, ν_bc, S_s_bc) @@ -413,7 +346,7 @@ function boundary_flux( # currently we approximate this as equal to the center value at the top layer (K_c) # More accurate would be to compute the mean between K_c and K evaluated at the boundary # condition. - K_eff = get_bottom_surface_field(p.soil.K, surface_space) + K_eff = ClimaLand.Domains.bottom_center_to_surface(p.soil.K) # At the bottom boundary, ψ_c is at larger z than ψ_bc # so we swap their order in the derivative calc @@ -596,10 +529,9 @@ function boundary_flux( )::ClimaCore.Fields.Field FT = eltype(Δz) # Approximate κ_bc ≈ κ_c (center closest to the boundary) - surface_space = axes(Δz) # We need to project the center values onto the face space. - T_c = get_top_surface_space(p.soil.T, surface_space) - κ_c = get_top_surface_space(p.soil.κ, surface_space) + T_c = ClimaLand.Domains.top_center_to_surface(p.soil.T) + κ_c = ClimaLand.Domains.top_center_to_surface(p.soil.κ) T_bc = FT.(heat_bc.bc(p, t)) return ClimaLand.diffusive_flux(κ_c, T_bc, T_c, Δz) @@ -629,10 +561,9 @@ function boundary_flux( )::ClimaCore.Fields.Field FT = eltype(Δz) # Approximate κ_bc ≈ κ_c (center closest to the boundary) - surface_space = axes(Δz) # We need to project the center values onto the face space. - T_c = get_bottom_surface_space(p.soil.T, surface_space) - κ_c = get_bottom_surface_space(p.soil.κ, surface_space) + T_c = ClimaLand.Domains.bottom_center_to_surface(p.soil.T) + κ_c = ClimaLand.Domains.bottom_center_to_surface(p.soil.κ) T_bc = FT.(heat_bc.bc(p, t)) return ClimaLand.diffusive_flux(κ_c, T_c, T_bc, Δz) end diff --git a/src/standalone/Soil/energy_hydrology.jl b/src/standalone/Soil/energy_hydrology.jl index f5cf74d331..4ab48c8e85 100644 --- a/src/standalone/Soil/energy_hydrology.jl +++ b/src/standalone/Soil/energy_hydrology.jl @@ -837,11 +837,11 @@ function ClimaLand.surface_specific_humidity( Thermodynamics.Liquid(), ) .* (@. exp(g * ψ_sfc * M_w / (R * T_sfc))) ice_frac = p.soil.ice_frac - surface_space = model.domain.space.surface - ν_sfc = get_top_surface_field(model.parameters.ν, surface_space) - θ_r_sfc = get_top_surface_field(model.parameters.θ_r, surface_space) - S_c_sfc = - get_top_surface_field(model.parameters.hydrology_cm.S_c, surface_space) + ν_sfc = ClimaLand.Domains.top_center_to_surface(model.parameters.ν) + θ_r_sfc = ClimaLand.Domains.top_center_to_surface(model.parameters.θ_r) + S_c_sfc = ClimaLand.Domains.top_center_to_surface( + model.parameters.hydrology_cm.S_c, + ) θ_l_sfc = ClimaLand.Domains.top_center_to_surface(p.soil.θ_l) θ_i_sfc = ClimaLand.Domains.top_center_to_surface(Y.soil.θ_i) @. ice_frac = ice_fraction(θ_l_sfc, θ_i_sfc, ν_sfc, θ_r_sfc)