Skip to content

Commit

Permalink
Merge pull request #258 from CliMA/ck/fix_warnings
Browse files Browse the repository at this point in the history
Fix some warnings
  • Loading branch information
charleskawczynski authored Feb 29, 2024
2 parents e160d77 + fd56e83 commit 31f0bc9
Show file tree
Hide file tree
Showing 3 changed files with 65 additions and 70 deletions.
6 changes: 4 additions & 2 deletions src/utilities/convergence_condition.jl
Original file line number Diff line number Diff line change
Expand Up @@ -92,8 +92,10 @@ updated_cache((; rate, order)::MinimumRateOfConvergence, cache, val, err, iter)
Checks multiple `ConvergenceCondition`s, combining their results with either
`all` or `any`.
"""
struct MultipleConditions{CC <: Union{typeof(all), typeof(any)}, C <: Tuple{Vararg{<:ConvergenceCondition}}} <:
ConvergenceCondition
struct MultipleConditions{
CC <: Union{typeof(all), typeof(any)},
C <: Tuple{Vararg{T} where {T <: ConvergenceCondition}},
} <: ConvergenceCondition
condition_combiner::CC
conditions::C
end
Expand Down
6 changes: 5 additions & 1 deletion test/integrator_utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,11 @@ data structure has an instance of any types in `t`.
"""
function any_reltype(found, obj, name, ets, pc = (); warn = true)
for pn in propertynames(obj)
prop = getproperty(obj, pn)
prop = if obj isa Base.Pairs
values(obj)
else
getproperty(obj, pn)
end
pc_full = (pc..., ".", pn)
pc_string = name * string(join(pc_full))
for et in ets
Expand Down
123 changes: 56 additions & 67 deletions test/problems.jl
Original file line number Diff line number Diff line change
@@ -1,16 +1,7 @@
using DiffEqBase, ClimaTimeSteppers, LinearAlgebra, StaticArrays
using ClimaCore
using ClimaComms
import ClimaCore.Domains as Domains
import ClimaCore.Geometry as Geometry
import ClimaCore.Meshes as Meshes
import ClimaCore.Topologies as Topologies
import ClimaCore.Spaces as Spaces
import ClimaCore.Fields as Fields
import ClimaCore.Operators as Operators

import Krylov
Krylov.ktypeof(x::Fields.FieldVector) = ClimaComms.array_type(x){eltype(parent(x)), 1}
import ClimaCore: Domains, Geometry, Meshes, Topologies, Spaces, Fields, Operators, Limiters

"""
Single variable linear ODE
Expand Down Expand Up @@ -449,6 +440,7 @@ Wfact!(W, Y, p, dtγ, t) = nothing
2D diffusion test problem. See [`2D diffusion problem`](@ref) for more details.
"""
function climacore_2Dheat_test_cts(::Type{FT}) where {FT}
context = ClimaComms.context()
dss_tendency = true

n_elem_x = 2
Expand All @@ -465,7 +457,7 @@ function climacore_2Dheat_test_cts(::Type{FT}) where {FT}
Domains.IntervalDomain(Geometry.YPoint(FT(0)), Geometry.YPoint(FT(1)), periodic = true),
)
mesh = Meshes.RectilinearMesh(domain, n_elem_x, n_elem_y)
topology = Topologies.Topology2D(mesh)
topology = Topologies.Topology2D(context, mesh)
quadrature = Spaces.Quadratures.GLL{n_poly + 1}()
space = Spaces.SpectralElementSpace2D(topology, quadrature)
(; x, y) = Fields.coordinate_field(space)
Expand Down Expand Up @@ -564,6 +556,7 @@ end
# Implemented in flux form, with an optional limiter and hyperdiffusion.
# TODO: Use this as an integration test.
function deformational_flow_test(::Type{FT}; use_limiter = true, use_hyperdiffusion = true) where {FT}
context = ClimaComms.context()
# Table III
# Note: the paper uses "a" in place of "R"
R = FT(6371220) # radius of Earth [m]
Expand Down Expand Up @@ -593,30 +586,28 @@ function deformational_flow_test(::Type{FT}; use_limiter = true, use_hyperdiffus
# hyperviscosity coefficient [m^4] (specified in the limiter paper)
D₄ = FT(6.6e14)

centers = ClimaCore.Geometry.LatLongZPoint.(rad2deg(φ_c), rad2deg.((λ_c1, λ_c2)), FT(0))
centers = Geometry.LatLongZPoint.(rad2deg(φ_c), rad2deg.((λ_c1, λ_c2)), FT(0))

# custom discretization (paper's discretization results in a very slow test)
vert_nelems = 10
horz_nelems = 4
horz_npoly = 3

vert_domain = ClimaCore.Domains.IntervalDomain(
ClimaCore.Geometry.ZPoint(FT(0)),
ClimaCore.Geometry.ZPoint(z_top);
boundary_names = (:bottom, :top),
)
vert_mesh = ClimaCore.Meshes.IntervalMesh(vert_domain, nelems = vert_nelems)
vert_cent_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(vert_mesh)
vert_domain =
Domains.IntervalDomain(Geometry.ZPoint(FT(0)), Geometry.ZPoint(z_top); boundary_names = (:bottom, :top))
vert_mesh = Meshes.IntervalMesh(vert_domain, nelems = vert_nelems)
z_topology = Topologies.IntervalTopology(context, vert_mesh)
vert_cent_space = Spaces.CenterFiniteDifferenceSpace(z_topology)

horz_domain = ClimaCore.Domains.SphereDomain(R)
horz_mesh = ClimaCore.Meshes.EquiangularCubedSphere(horz_domain, horz_nelems)
horz_topology = ClimaCore.Topologies.Topology2D(horz_mesh)
horz_quad = ClimaCore.Spaces.Quadratures.GLL{horz_npoly + 1}()
horz_space = ClimaCore.Spaces.SpectralElementSpace2D(horz_topology, horz_quad)
horz_domain = Domains.SphereDomain(R)
horz_mesh = Meshes.EquiangularCubedSphere(horz_domain, horz_nelems)
horz_topology = Topologies.Topology2D(context, horz_mesh)
horz_quad = Spaces.Quadratures.GLL{horz_npoly + 1}()
horz_space = Spaces.SpectralElementSpace2D(horz_topology, horz_quad)

cent_space = ClimaCore.Spaces.ExtrudedFiniteDifferenceSpace(horz_space, vert_cent_space)
cent_coords = ClimaCore.Fields.coordinate_field(cent_space)
face_space = ClimaCore.Spaces.FaceExtrudedFiniteDifferenceSpace(cent_space)
cent_space = Spaces.ExtrudedFiniteDifferenceSpace(horz_space, vert_cent_space)
cent_coords = Fields.coordinate_field(cent_space)
face_space = Spaces.FaceExtrudedFiniteDifferenceSpace(cent_space)

# initial density (Equation 8)
cent_ρ = @. p_0 / (R_d * T_0) * exp(-cent_coords.z / H)
Expand All @@ -627,7 +618,7 @@ function deformational_flow_test(::Type{FT}; use_limiter = true, use_hyperdiffus
φ = deg2rad(coord.lat)

ds = map(centers) do center
r = ClimaCore.Geometry.great_circle_distance(coord, center, horz_space.global_geometry)
r = Geometry.great_circle_distance(coord, center, horz_space.global_geometry)
return min(1, (r / R_t)^2 + ((z - z_c) / Z_t)^2)
end
in_slot = z > z_c && φ_c - FT(0.125) < φ < φ_c + FT(0.125)
Expand All @@ -640,11 +631,11 @@ function deformational_flow_test(::Type{FT}; use_limiter = true, use_hyperdiffus
return (; q1, q2, q3, q4, q5)
end

init_state = ClimaCore.Fields.FieldVector(; cent_ρ, cent_ρq = cent_ρ .* cent_q)
init_state = Fields.FieldVector(; cent_ρ, cent_ρq = cent_ρ .* cent_q)

# current wind vector (Equations 15--26)
current_cent_wind_vector = ClimaCore.Fields.Field(ClimaCore.Geometry.UVWVector{FT}, cent_space)
current_face_wind_vector = ClimaCore.Fields.Field(ClimaCore.Geometry.UVWVector{FT}, face_space)
current_cent_wind_vector = Fields.Field(Geometry.UVWVector{FT}, cent_space)
current_face_wind_vector = Fields.Field(Geometry.UVWVector{FT}, face_space)
function wind_vector(coord, ρ, t)
z = coord.z
φ = deg2rad(coord.lat)
Expand All @@ -667,45 +658,42 @@ function deformational_flow_test(::Type{FT}; use_limiter = true, use_hyperdiffus
ω = ω_0 * sin(λ′) * cos(φ) * cos(2 * FT(π) * t / τ) * s
w = -ω / (g * ρ)

return ClimaCore.Geometry.UVWVector(u, v, w)
return Geometry.UVWVector(u, v, w)
end

horz_div = ClimaCore.Operators.Divergence()
horz_wdiv = ClimaCore.Operators.WeakDivergence()
horz_grad = ClimaCore.Operators.Gradient()
horz_div = Operators.Divergence()
horz_wdiv = Operators.WeakDivergence()
horz_grad = Operators.Gradient()
cent_χ = similar(cent_q)
function T_lim!(tendency, state, _, t)
@. current_cent_wind_vector = wind_vector(cent_coords, state.cent_ρ, t)
@. tendency.cent_ρ = -horz_div(state.cent_ρ * current_cent_wind_vector)
@. tendency.cent_ρq = -horz_div(state.cent_ρq * current_cent_wind_vector)
use_hyperdiffusion || return nothing
@. cent_χ = horz_wdiv(horz_grad(state.cent_ρq / state.cent_ρ))
ClimaCore.Spaces.weighted_dss!(cent_χ)
Spaces.weighted_dss!(cent_χ)
@. tendency.cent_ρq += -D₄ * horz_wdiv(state.cent_ρ * horz_grad(cent_χ))
return nothing
end

limiter = ClimaCore.Limiters.QuasiMonotoneLimiter(cent_q; rtol = FT(0))
limiter = Limiters.QuasiMonotoneLimiter(cent_q; rtol = FT(0))
function lim!(state, _, t, ref_state)
use_limiter || return nothing
ClimaCore.Limiters.compute_bounds!(limiter, ref_state.cent_ρq, ref_state.cent_ρ)
ClimaCore.Limiters.apply_limiter!(state.cent_ρq, state.cent_ρ, limiter)
Limiters.compute_bounds!(limiter, ref_state.cent_ρq, ref_state.cent_ρ)
Limiters.apply_limiter!(state.cent_ρq, state.cent_ρ, limiter)
return nothing
end

vert_div = ClimaCore.Operators.DivergenceF2C()
vert_interp = ClimaCore.Operators.InterpolateC2F(
top = ClimaCore.Operators.Extrapolate(),
bottom = ClimaCore.Operators.Extrapolate(),
)
vert_div = Operators.DivergenceF2C()
vert_interp = Operators.InterpolateC2F(top = Operators.Extrapolate(), bottom = Operators.Extrapolate())
function T_exp!(tendency, state, _, t)
@. current_face_wind_vector = wind_vector(face_coords, vert_interp(state.cent_ρ), t)
@. tendency.cent_ρ = -vert_div(vert_interp(state.cent_ρ) * current_face_wind_vector)
@. tendency.cent_ρq = -vert_div(vert_interp(state.cent_ρq) * current_face_wind_vector)
end

function dss!(state, _, t)
ClimaCore.Spaces.weighted_dss!(state.q)
Spaces.weighted_dss!(state.q)
end

function analytic_sol(t)
Expand All @@ -731,6 +719,7 @@ end
# (https://gmd.copernicus.org/articles/5/887/2012/gmd-5-887-2012.pdf)
# Implemented in flux form, with an optional limiter and hyperdiffusion.
function horizontal_deformational_flow_test(::Type{FT}; use_limiter = true, use_hyperdiffusion = true) where {FT}
context = ClimaComms.context()
# constants (using the same notation as deformational_flow_test)
R = FT(6371220) # radius of Earth [m]
τ = 60 * 60 * 24 * FT(12) # period of motion (12 days) [s]
Expand All @@ -740,18 +729,18 @@ function horizontal_deformational_flow_test(::Type{FT}; use_limiter = true, use_
R_t = R / 2 # horizontal half-width of tracers [m]
D₄ = FT(6.6e14) # hyperviscosity coefficient [m^4] (specified in the limiter paper)

centers = ClimaCore.Geometry.LatLongPoint.(rad2deg(φ_c), rad2deg.((λ_c1, λ_c2)))
centers = Geometry.LatLongPoint.(rad2deg(φ_c), rad2deg.((λ_c1, λ_c2)))

# 1.5° resolution on the equator: 360° / (4 * nelems * npoly) = 1.5°
nelems = 20
npoly = 3

domain = ClimaCore.Domains.SphereDomain(R)
mesh = ClimaCore.Meshes.EquiangularCubedSphere(domain, nelems)
topology = ClimaCore.Topologies.Topology2D(mesh)
quad = ClimaCore.Spaces.Quadratures.GLL{npoly + 1}()
space = ClimaCore.Spaces.SpectralElementSpace2D(topology, quad)
coords = ClimaCore.Fields.coordinate_field(space)
domain = Domains.SphereDomain(R)
mesh = Meshes.EquiangularCubedSphere(domain, nelems)
topology = Topologies.Topology2D(context, mesh)
quad = Spaces.Quadratures.GLL{npoly + 1}()
space = Spaces.SpectralElementSpace2D(topology, quad)
coords = Fields.coordinate_field(space)

# initial conditions (Section 2.2)
ρ = ones(space)
Expand All @@ -760,15 +749,15 @@ function horizontal_deformational_flow_test(::Type{FT}; use_limiter = true, use_
λ = deg2rad(coord.long)

hs = map(centers) do center
center′ = ClimaCore.Geometry.CartesianPoint(center, space.global_geometry)
coord′ = ClimaCore.Geometry.CartesianPoint(coord, space.global_geometry)
center′ = Geometry.CartesianPoint(center, space.global_geometry)
coord′ = Geometry.CartesianPoint(coord, space.global_geometry)
dist_squared = (coord′.x1 - center′.x1)^2 + (coord′.x2 - center′.x2)^2 + (coord′.x3 - center′.x3)^2
# Note: the paper doesn't divide by R^2, which only works if R = 1
return FT(0.95) * exp(-5 * dist_squared / R^2)
end
gaussian_hills = hs[1] + hs[2]
rs = map(centers) do center
return ClimaCore.Geometry.great_circle_distance(coord, center, space.global_geometry)
return Geometry.great_circle_distance(coord, center, space.global_geometry)
end
cosine_bells = if rs[1] < R_t
FT(0.1) + FT(0.9) * (1 + cos(FT(π) * rs[1] / R_t)) / 2
Expand All @@ -791,10 +780,10 @@ function horizontal_deformational_flow_test(::Type{FT}; use_limiter = true, use_

return (; gaussian_hills, cosine_bells, slotted_cylinders)
end
init_state = ClimaCore.Fields.FieldVector(; ρ, ρq = ρ .* q)
init_state = Fields.FieldVector(; ρ, ρq = ρ .* q)

# current wind vector (Section 2.3)
current_wind_vector = ClimaCore.Fields.Field(ClimaCore.Geometry.UVVector{FT}, space)
current_wind_vector = Fields.Field(Geometry.UVVector{FT}, space)
function wind_vector(coord, t)
φ = deg2rad(coord.lat)
λ = deg2rad(coord.long)
Expand All @@ -805,35 +794,35 @@ function horizontal_deformational_flow_test(::Type{FT}; use_limiter = true, use_
u = k * sin(λ′)^2 * sin(2 * φ) * cos(FT(π) * t / τ) + 2 * FT(π) * R / τ * cos(φ)
v = k * sin(2 * λ′) * cos(φ) * cos(FT(π) * t / τ)

return ClimaCore.Geometry.UVVector(u, v)
return Geometry.UVVector(u, v)
end

div = ClimaCore.Operators.Divergence()
wdiv = ClimaCore.Operators.WeakDivergence()
grad = ClimaCore.Operators.Gradient()
div = Operators.Divergence()
wdiv = Operators.WeakDivergence()
grad = Operators.Gradient()
χ = similar(q)
function T_lim!(tendency, state, _, t)
@. current_wind_vector = wind_vector(coords, t)
@. tendency.ρ = -div(state.ρ * current_wind_vector)
@. tendency.ρq = -div(state.ρq * current_wind_vector)
use_hyperdiffusion || return nothing
@. χ = wdiv(grad(state.ρq / state.ρ))
ClimaCore.Spaces.weighted_dss!(χ)
Spaces.weighted_dss!(χ)
@. tendency.ρq += -D₄ * wdiv(state.ρ * grad(χ))
return nothing
end

limiter = ClimaCore.Limiters.QuasiMonotoneLimiter(q; rtol = FT(0))
limiter = Limiters.QuasiMonotoneLimiter(q; rtol = FT(0))
function lim!(state, _, t, ref_state)
use_limiter || return nothing
ClimaCore.Limiters.compute_bounds!(limiter, ref_state.ρq, ref_state.ρ)
ClimaCore.Limiters.apply_limiter!(state.ρq, state.ρ, limiter)
Limiters.compute_bounds!(limiter, ref_state.ρq, ref_state.ρ)
Limiters.apply_limiter!(state.ρq, state.ρ, limiter)
return nothing
end

function dss!(state, _, t)
ClimaCore.Spaces.weighted_dss!(state.ρ)
ClimaCore.Spaces.weighted_dss!(state.ρq)
Spaces.weighted_dss!(state.ρ)
Spaces.weighted_dss!(state.ρq)
end

function analytic_sol(t)
Expand Down

0 comments on commit 31f0bc9

Please sign in to comment.