From 14ae29965a5352ef0ca4342bd32240779e4cd286 Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Mon, 17 Jul 2023 10:06:57 -0700 Subject: [PATCH 01/11] Fix gpu column integrals --- src/Operators/integrals.jl | 61 +- test/Fields/field.jl | 1303 ++++++++++++++++++------------------ 2 files changed, 708 insertions(+), 656 deletions(-) diff --git a/src/Operators/integrals.jl b/src/Operators/integrals.jl index c7c65d0b0f..8ef963da07 100644 --- a/src/Operators/integrals.jl +++ b/src/Operators/integrals.jl @@ -7,23 +7,74 @@ The definite vertical column integral, `col∫field`, of field `field`. """ -function column_integral_definite!(col∫field::Fields.Field, field::Fields.Field) +column_integral_definite!(col∫field::Fields.Field, field::Fields.Field) = + column_integral_definite!( + ClimaComms.device(axes(col∫field)), + col∫field, + field, + ) + +function column_integral_definite!( + ::ClimaComms.CUDADevice, + col∫field::Fields.Field, + field::Fields.Field, +) + Ni, Nj, _, _, Nh = size(Fields.field_values(col∫field)) + nthreads, nblocks = Spaces._configure_threadblock(Ni * Nj * Nh) + @cuda threads = nthreads blocks = nblocks column_integral_definite_kernel!( + col∫field, + field, + ) +end + +function column_integral_definite_kernel!( + col∫field::Fields.SpectralElementField, + field::Fields.ExtrudedFiniteDifferenceField, +) + idx = threadIdx().x + (blockIdx().x - 1) * blockDim().x + Ni, Nj, _, _, Nh = size(Fields.field_values(field)) + if idx <= Ni * Nj * Nh + i, j, h = Spaces._get_idx((Ni, Nj, Nh), idx) + _column_integral_definite!( + Spaces.column(col∫field, i, j, h), + Spaces.column(field, i, j, h), + ) + end + return nothing +end + +column_integral_definite_kernel!( + col∫field::Fields.PointField, + field::Fields.FiniteDifferenceField, +) = _column_integral_definite!(col∫field, field) + +function column_integral_definite!( + ::ClimaComms.AbstractCPUDevice, + col∫field::Fields.SpectralElementField, + field::Fields.ExtrudedFiniteDifferenceField, +) Fields.bycolumn(axes(field)) do colidx - column_integral_definite!(col∫field[colidx], field[colidx]) + _column_integral_definite!(col∫field[colidx], field[colidx]) nothing end return nothing end -function column_integral_definite!( +column_integral_definite!( + ::ClimaComms.AbstractCPUDevice, + col∫field::Fields.PointField, + field::Fields.FiniteDifferenceField, +) = _column_integral_definite!(col∫field, field) + +function _column_integral_definite!( col∫field::Fields.PointField, field::Fields.ColumnField, ) - @inbounds col∫field[] = column_integral_definite(field) + @inbounds col∫field[] = _column_integral_definite(field) return nothing end -function column_integral_definite(field::Fields.ColumnField) +function _column_integral_definite(field::Fields.ColumnField) field_data = Fields.field_values(field) Δz_data = Spaces.Δz_data(axes(field)) Nv = Spaces.nlevels(axes(field)) diff --git a/test/Fields/field.jl b/test/Fields/field.jl index cb711f0306..fd6a6f0638 100644 --- a/test/Fields/field.jl +++ b/test/Fields/field.jl @@ -17,656 +17,656 @@ include( ) import .TestUtilities as TU -function spectral_space_2D(; n1 = 1, n2 = 1, Nij = 4) - domain = Domains.RectangleDomain( - Geometry.XPoint(-3.0) .. Geometry.XPoint(5.0), - Geometry.YPoint(-2.0) .. Geometry.YPoint(8.0), - x1periodic = false, - x2periodic = false, - x1boundary = (:east, :west), - x2boundary = (:south, :north), - ) - mesh = Meshes.RectilinearMesh(domain, n1, n2) - device = ClimaComms.CPUSingleThreaded() - grid_topology = - Topologies.Topology2D(ClimaComms.SingletonCommsContext(device), mesh) - - quad = Spaces.Quadratures.GLL{Nij}() - space = Spaces.SpectralElementSpace2D(grid_topology, quad) - return space -end - -@testset "1×1 2D domain space" begin - Nij = 4 - n1 = n2 = 1 - space = spectral_space_2D(n1 = n1, n2 = n2, Nij = Nij) - - field = - Fields.Field(IJFH{ComplexF64, Nij}(ones(Nij, Nij, 2, n1 * n2)), space) - - @test sum(field) ≈ Complex(1.0, 1.0) * 8.0 * 10.0 rtol = 10eps() - @test sum(x -> 3.0, field) ≈ 3 * 8.0 * 10.0 rtol = 10eps() - @test mean(field) ≈ Complex(1.0, 1.0) rtol = 10eps() - @test mean(x -> 3.0, field) ≈ 3 rtol = 10eps() - @test norm(field) ≈ sqrt(2.0) rtol = 10eps() - @test norm(field, 1) ≈ norm(Complex(1.0, 1.0)) rtol = 10eps() - @test norm(field, Inf) ≈ norm(Complex(1.0, 1.0)) rtol = 10eps() - @test norm(field; normalize = false) ≈ sqrt(2.0 * 8.0 * 10.0) rtol = 10eps() - @test norm(field, 1; normalize = false) ≈ - norm(Complex(1.0, 1.0)) * 8.0 * 10.0 rtol = 10eps() - @test norm(field, Inf; normalize = false) ≈ norm(Complex(1.0, 1.0)) rtol = - 10eps() - - @test extrema(real, field) == (1.0, 1.0) - - @test Operators.matrix_interpolate(field, 4) ≈ - [Complex(1.0, 1.0) for i in 1:(4 * n1), j in 1:(4 * n2)] - - - field_sin = map(x -> sin((x.x) / 2), Fields.coordinate_field(space)) - M = Operators.matrix_interpolate(field_sin, 20) - @test size(M) == (20, 20) # 20 x 20 for a 1 element field - - real_field = field.re - - # test broadcasting - res = field .+ 1 - @test parent(Fields.field_values(res)) == Float64[ - f == 1 ? 2 : 1 for i in 1:Nij, j in 1:Nij, f in 1:2, h in 1:(n1 * n2) - ] - - res = field.re .+ 1 - @test parent(Fields.field_values(res)) == - Float64[2 for i in 1:Nij, j in 1:Nij, f in 1:1, h in 1:(n1 * n2)] - - # test field slab broadcasting - f1 = ones(space) - f2 = ones(space) - - for h in 1:(n1 * n2) - f1_slab = Fields.slab(f1, h) - f2_slab = Fields.slab(f2, h) - q = f1_slab .+ f2_slab - f1_slab .= q .+ f2_slab - end - @test all(parent(f1) .== 3) - - point_field = Fields.column(field, 1, 1, 1) - @test axes(point_field) isa Spaces.PointSpace -end - -# https://github.com/CliMA/ClimaCore.jl/issues/1126 -function pow_n(f) - @. f.x = f.x^2 - return nothing -end -@testset "Broadcasting with ^n" begin - FT = Float32 - device = ClimaComms.CPUSingleThreaded() # fill is broken on gpu - context = ClimaComms.SingletonCommsContext(device) - for space in TU.all_spaces(FT; context) - f = fill((; x = FT(1)), space) - pow_n(f) # Compile first - p_allocated = @allocated pow_n(f) - if space isa Spaces.SpectralElementSpace1D - @test p_allocated == 0 - else - @test p_allocated == 0 broken = (device isa ClimaComms.CUDADevice) - end - end -end - -function ifelse_broadcast_allocating(a, b, c) - FT = eltype(a) - @. a = ifelse(true || c < b * FT(1), FT(0), c) - return nothing -end - -function ifelse_broadcast_or(a, b, c) - FT = eltype(a) - val = FT(1) - @. a = ifelse(true || c < b * val, FT(0), c) - return nothing -end - -function ifelse_broadcast_simple(a, b, c) - FT = eltype(a) - @. a = ifelse(c < b * FT(1), FT(0), c) - return nothing -end - -@testset "Broadcasting ifelse" begin - FT = Float32 - device = ClimaComms.CPUSingleThreaded() # broken on gpu - context = ClimaComms.SingletonCommsContext(device) - for space in ( - TU.CenterExtrudedFiniteDifferenceSpace(FT; context), - TU.ColumnCenterFiniteDifferenceSpace(FT; context), - ) - a = Fields.level(fill(FT(0), space), 1) - b = Fields.level(fill(FT(2), space), 1) - c = Fields.level(fill(FT(3), space), 1) - - ifelse_broadcast_allocating(a, b, c) - p_allocated = @allocated ifelse_broadcast_allocating(a, b, c) - @test_broken p_allocated == 0 - - ifelse_broadcast_or(a, b, c) - p_allocated = @allocated ifelse_broadcast_or(a, b, c) - @test p_allocated == 0 - - ifelse_broadcast_simple(a, b, c) - p_allocated = @allocated ifelse_broadcast_simple(a, b, c) - @test p_allocated == 0 - end -end - -# Requires `--check-bounds=yes` -@testset "Constructing & broadcasting over empty fields" begin - FT = Float32 - for space in TU.all_spaces(FT) - f = fill((;), space) - @. f += f - end - - function test_broken_throws(f) - try - @. f += 1 - # we want to throw exception, test is broken - @test_broken false - catch - # we want to throw exception, unexpected pass - @test_broken true - end - end - empty_field(space) = fill((;), space) - - # Broadcasting over the wrong size should error - test_broken_throws(empty_field(TU.PointSpace(FT))) - test_broken_throws(empty_field(TU.SpectralElementSpace1D(FT))) - test_broken_throws(empty_field(TU.SpectralElementSpace2D(FT))) - test_broken_throws(empty_field(TU.ColumnCenterFiniteDifferenceSpace(FT))) - test_broken_throws(empty_field(TU.ColumnFaceFiniteDifferenceSpace(FT))) - test_broken_throws(empty_field(TU.SphereSpectralElementSpace(FT))) - test_broken_throws(empty_field(TU.CenterExtrudedFiniteDifferenceSpace(FT))) - test_broken_throws(empty_field(TU.FaceExtrudedFiniteDifferenceSpace(FT))) - - # TODO: performance optimization: shouldn't we do - # nothing when broadcasting over empty fields? - # This is otherwise a performance penalty if - # users regularly rely on empty fields. In particular: - # - does iterating over empty fields load data? - # - what is the overhead in iterating over empty fields? - # - what is the use case of anything useful that can be - # done by iterating over empty fields? -end - -@testset "Broadcasting interception for tuple-valued fields" begin - n1 = n2 = 1 - Nij = 4 - space = spectral_space_2D(n1 = n1, n2 = n2, Nij = Nij) - - nt_field = Fields.Field( - IJFH{NamedTuple{(:a, :b), Tuple{Float64, Float64}}, Nij}( - ones(Nij, Nij, 2, n1 * n2), - ), - space, - ) - nt_sum = sum(nt_field) - @test nt_sum isa NamedTuple{(:a, :b), Tuple{Float64, Float64}} - @test nt_sum.a ≈ 8.0 * 10.0 rtol = 10eps() - @test nt_sum.b ≈ 8.0 * 10.0 rtol = 10eps() - @test norm(nt_field) ≈ sqrt(2.0) rtol = 10eps() - - # test scalar asignment - nt_field.a .= 0.0 - @test sum(nt_field.a) == 0.0 -end - -@testset "Special case handling for broadcased norm to pass through space local geometry" begin - space = spectral_space_2D() - u = Geometry.Covariant12Vector.(ones(space), ones(space)) - @test norm.(u) ≈ hypot(4 / 8 / 2, 4 / 10 / 2) .* ones(space) -end - -@testset "FieldVector" begin - space = spectral_space_2D() - u = Geometry.Covariant12Vector.(ones(space), ones(space)) - x = Fields.coordinate_field(space) - y = [1.0, 2.0, 3.0] - z = 1.0 - Y = Fields.FieldVector(u = u, k = (x = x, y = y, z = z)) - - @test propertynames(Y) == (:u, :k) - @test propertynames(Y.k) == (:x, :y, :z) - @test Y.u === u - @test Y.k.x === x - @test Y.k.y === y - @test Y.k.z === z - - @test deepcopy(Y).u !== u - @test deepcopy(Y).k.x !== x - @test deepcopy(Y).k.y !== y - - @test getfield(deepcopy(Y).u, :space) === space - - Y1 = 2 .* Y - @test parent(Y1.u) == 2 .* parent(u) - @test parent(Y1.k.x) == 2 .* parent(x) - @test Y1.k.y == 2 .* y - @test Y1.k.z === 2 * z - - Y1 .= Y1 .+ 2 .* Y - @test parent(Y1.u) == 4 .* parent(u) - @test parent(Y1.k.x) == 4 .* parent(x) - @test Y1.k.y == 4 .* y - @test Y1.k.z === 4 * z - - Y.k.z = 3.0 - @test Y.k.z === 3.0 -end - -@testset "FieldVector array_type" begin - device = ClimaComms.device() - context = ClimaComms.SingletonCommsContext(device) - space = TU.PointSpace(Float32; context) - xcenters = Fields.coordinate_field(space).x - y = Fields.FieldVector(x = xcenters) - @test ClimaComms.array_type(y) == ClimaComms.array_type(device) - y = Fields.FieldVector(x = xcenters, y = xcenters) - @test ClimaComms.array_type(y) == ClimaComms.array_type(device) -end - -@testset "FieldVector basetype replacement and deepcopy" begin - device = ClimaComms.CPUSingleThreaded() # constructing space_vijfh is broken - context = ClimaComms.SingletonCommsContext(device) - domain_z = Domains.IntervalDomain( - Geometry.ZPoint(-1.0) .. Geometry.ZPoint(1.0), - periodic = true, - ) - mesh_z = Meshes.IntervalMesh(domain_z; nelems = 10) - topology_z = Topologies.IntervalTopology(context, mesh_z) - - domain_x = Domains.IntervalDomain( - Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0), - periodic = true, - ) - mesh_x = Meshes.IntervalMesh(domain_x; nelems = 10) - topology_x = Topologies.IntervalTopology(context, mesh_x) - - domain_xy = Domains.RectangleDomain( - Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0), - Geometry.YPoint(-1.0) .. Geometry.YPoint(1.0), - x1periodic = true, - x2periodic = true, - ) - mesh_xy = Meshes.RectilinearMesh(domain_xy, 10, 10) - topology_xy = Topologies.Topology2D(context, mesh_xy) - - quad = Spaces.Quadratures.GLL{4}() - - space_vf = Spaces.CenterFiniteDifferenceSpace(topology_z) - space_ifh = Spaces.SpectralElementSpace1D(topology_x, quad) - space_ijfh = Spaces.SpectralElementSpace2D(topology_xy, quad) - space_vifh = Spaces.ExtrudedFiniteDifferenceSpace(space_ifh, space_vf) - space_vijfh = Spaces.ExtrudedFiniteDifferenceSpace(space_ijfh, space_vf) - - space2field(space) = map( - coord -> (coord, Geometry.Covariant12Vector(1.0, 2.0)), - Fields.coordinate_field(space), - ) - - Y = Fields.FieldVector( - field_vf = space2field(space_vf), - field_if = slab(space2field(space_ifh), 1), - field_ifh = space2field(space_ifh), - field_ijf = slab(space2field(space_ijfh), 1, 1), - field_ijfh = space2field(space_ijfh), - field_vifh = space2field(space_vifh), - field_vijfh = space2field(space_vijfh), - array = [1.0, 2.0, 3.0], - scalar = 1.0, - ) - - Yf = ForwardDiff.Dual{Nothing}.(Y, 1.0) - Yf .= Yf .^ 2 .+ Y - @test all(ForwardDiff.value.(Yf) .== Y .^ 2 .+ Y) - @test all(ForwardDiff.partials.(Yf, 1) .== 2 .* Y) - - dual_field = Yf.field_vf - dual_field_original_basetype = similar(Y.field_vf, eltype(dual_field)) - @test eltype(dual_field_original_basetype) === eltype(dual_field) - @test eltype(parent(dual_field_original_basetype)) === Float64 - @test eltype(parent(dual_field)) === ForwardDiff.Dual{Nothing, Float64, 1} - - object_that_contains_Yf = (; Yf) - @test axes(deepcopy(Yf).field_vf) === space_vf - @test axes(deepcopy(object_that_contains_Yf).Yf.field_vf) === space_vf -end - -@testset "Scalar field iterator" begin - space = spectral_space_2D() - u = Geometry.Covariant12Vector.(ones(space), ones(space)) - x = Fields.coordinate_field(space) - y = [1.0, 2.0, 3.0] - z = 1.0 - Y = Fields.FieldVector(u = u, k = (x = x, y = y, z = z)) - - prop_chains = Fields.property_chains(Y) - @test length(prop_chains) == 6 - @test prop_chains[1] == (:u, :components, :data, 1) - @test prop_chains[2] == (:u, :components, :data, 2) - @test prop_chains[3] == (:k, :x, :x) - @test prop_chains[4] == (:k, :x, :y) - @test prop_chains[5] == (:k, :y) - @test prop_chains[6] == (:k, :z) - - FT = Float64 - nt = - (; x = FT(0), y = FT(0), tup = ntuple(i -> (; a = FT(1), b = FT(1)), 2)) - Y = fill(nt, space) - - prop_chains = Fields.property_chains(Y) - @test prop_chains[1] == (:x,) - @test prop_chains[2] == (:y,) - @test prop_chains[3] == (:tup, 1, :a) - @test prop_chains[4] == (:tup, 1, :b) - @test prop_chains[5] == (:tup, 2, :a) - @test prop_chains[6] == (:tup, 2, :b) - - @test Fields.single_field(Y, prop_chains[1]) === Y.x - @test Fields.single_field(Y, prop_chains[2]) === Y.y - @test Fields.single_field(Y, prop_chains[3]) === getproperty(Y.tup, 1).a - @test Fields.single_field(Y, prop_chains[4]) === getproperty(Y.tup, 1).b - @test Fields.single_field(Y, prop_chains[5]) === getproperty(Y.tup, 2).a - @test Fields.single_field(Y, prop_chains[6]) === getproperty(Y.tup, 2).b - - for (i, (var, prop_chain)) in enumerate(Fields.field_iterator(Y)) - @test prop_chains[i] == prop_chain - @test var === Fields.single_field(Y, prop_chain) - end -end - -# Test truncated field type printing: -ClimaCore.Fields.truncate_printing_field_types() = true -@testset "Truncated printing" begin - nt = (; x = Float64(0), y = Float64(0)) - Y = fill(nt, spectral_space_2D()) - @test sprint(show, typeof(Y); context = IOContext(stdout)) == - "Field{(:x, :y)} (trunc disp)" -end -ClimaCore.Fields.truncate_printing_field_types() = false - -@testset "Standard printing" begin - nt = (; x = Float64(0), y = Float64(0)) - Y = fill(nt, spectral_space_2D()) - s = sprint(show, typeof(Y)) # just make sure this doesn't break -end - -@testset "Set!" begin - space = spectral_space_2D() - FT = Float64 - nt = (; x = FT(0), y = FT(0)) - Y = fill(nt, space) - foo(local_geom) = - sin(local_geom.coordinates.x * local_geom.coordinates.y) + 3 - Fields.set!(foo, Y.x) - @test all((parent(Y.x) .> 1)) -end - -@testset "PointField" begin - device = ClimaComms.CPUSingleThreaded() # a bunch of cuda pieces are broken - context = ClimaComms.SingletonCommsContext(device) - FT = Float64 - coord = Geometry.XPoint(FT(π)) - space = Spaces.PointSpace(context, coord) - @test parent(Spaces.local_geometry_data(space)) == - FT[Geometry.component(coord, 1), 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] - field = Fields.coordinate_field(space) - @test field isa Fields.PointField - @test Fields.field_values(field)[] == coord - - if ClimaComms.device(context) isa ClimaComms.AbstractCPUDevice - @test sum(field.x) == FT(π) - elseif ClimaComms.device(context) isa ClimaComms.CUDADevice - # Not yet supported - # @test sum(field.x) == FT(π) - end - - field = ones(space) .* π - sin_field = sin.(field) - add_field = field .+ field - @test isapprox(Fields.field_values(sin_field)[], FT(0.0); atol = √eps(FT)) - @test isapprox(Fields.field_values(add_field)[], FT(2π)) -end - -@testset "Level" begin - FT = Float64 - for space in TU.all_spaces(FT) - TU.levelable(space) || continue - Y = fill((; x = FT(2)), space) - lg_space = Spaces.level(space, TU.fc_index(1, space)) - lg_field_space = axes(Fields.level(Y, TU.fc_index(1, space))) - @test all( - lg_space.local_geometry.coordinates === - lg_field_space.local_geometry.coordinates, - ) - @test all(Fields.zeros(lg_space) == Fields.zeros(lg_field_space)) - end -end - -@testset "Points from Columns" begin - FT = Float64 - for space in TU.all_spaces(FT) - if space isa Spaces.SpectralElementSpace1D - Y = fill((; x = FT(1)), space) - point_space_from_field = axes(Fields.column(Y.x, 1, 1)) - point_space = Spaces.column(space, 1, 1) - @test Fields.ones(point_space) == - Fields.ones(point_space_from_field) - end - if space isa Spaces.SpectralElementSpace2D - Y = fill((; x = FT(1)), space) - point_space_from_field = axes(Fields.column(Y.x, 1, 1, 1)) - point_space = Spaces.column(space, 1, 1, 1) - @test Fields.ones(point_space) == - Fields.ones(point_space_from_field) - end - - end -end - -@testset "(Domain/Column)-surface broadcasting" begin - FT = Float64 - function domain_surface_bc!(x, ᶜz_surf, ᶜx_surf) - @. x = x + ᶜz_surf - # exercises broadcast_shape(PointSpace, PointSpace) - @. x = x + (ᶜz_surf * ᶜx_surf) - nothing - end - function column_surface_bc!(x, ᶜz_surf, ᶜx_surf) - Fields.bycolumn(axes(x)) do colidx - @. x[colidx] = x[colidx] + ᶜz_surf[colidx] - # exercises broadcast_shape(PointSpace, PointSpace) - @. x[colidx] = x[colidx] + (ᶜz_surf[colidx] * ᶜx_surf[colidx]) - end - nothing - end - for space in TU.all_spaces(FT) - # Filter out spaces without z coordinates: - TU.has_z_coordinates(space) || continue - Y = fill((; x = FT(1)), space) - ᶜz_surf = - Spaces.level(Fields.coordinate_field(Y).z, TU.fc_index(1, space)) - ᶜx_surf = copy(Spaces.level(Y.x, TU.fc_index(1, space))) - - # Still need to define broadcast rules for surface planes with 3D domains - domain_surface_bc!(Y.x, ᶜz_surf, ᶜx_surf) - - # Skip spaces incompatible with Fields.bycolumn: - TU.bycolumnable(space) || continue - Yc = fill((; x = FT(1)), space) - column_surface_bc!(Yc.x, ᶜz_surf, ᶜx_surf) - @test Y.x == Yc.x - nothing - end - nothing -end - -@testset "Broadcasting same spaces different instances" begin - space1 = spectral_space_2D() - space2 = spectral_space_2D() - field1 = ones(space1) - field2 = 2 .* ones(space2) - @test Fields.is_diagonalized_spaces(typeof(space1), typeof(space2)) - @test_throws ErrorException( - "Broacasted spaces are the same ClimaCore.Spaces type but not the same instance", - ) field1 .= field2 - - # turn warning on - Fields.allow_mismatched_diagonalized_spaces() = true - @test_warn "Broacasted spaces are the same ClimaCore.Spaces type but not the same instance" field1 .= - field2 - @test parent(field1) == parent(field2) - Fields.allow_mismatched_diagonalized_spaces() = false -end - -struct InferenceFoo{FT} - bar::FT -end -Base.broadcastable(x::InferenceFoo) = Ref(x) -@testset "Inference failure message" begin - function ics_foo(::Type{FT}, lg, foo) where {FT} - uv = Geometry.UVVector(FT(0), FT(0)) - z = Geometry.Covariant12Vector(uv, lg) - y = foo.bingo - return (; x = FT(0) + y) - end - function ics_foo_with_field(::Type{FT}, lg, foo, f) where {FT} - uv = Geometry.UVVector(FT(0), FT(0)) - z = Geometry.Covariant12Vector(uv, lg) - ζ = f.a - y = foo.baz - return (; x = FT(0) + y - ζ) - end - function FieldFromNamedTupleBroken( - space, - ics::Function, - ::Type{FT}, - params..., - ) where {FT} - lg = Fields.local_geometry_field(space) - return ics.(FT, lg, params...) - end - FT = Float64 - foo = InferenceFoo(2.0) - device = ClimaComms.CPUSingleThreaded() # cuda fill is broken - context = ClimaComms.SingletonCommsContext(device) - for space in TU.all_spaces(FT; context) - Y = fill((; a = FT(0), b = FT(1)), space) - @test_throws ErrorException("type InferenceFoo has no field bingo") FieldFromNamedTupleBroken( - space, - ics_foo, - FT, - foo, - ) - @test_throws ErrorException("type InferenceFoo has no field baz") FieldFromNamedTupleBroken( - space, - ics_foo_with_field, - FT, - foo, - Y, - ) - end - -end - -@testset "Δz_field" begin - FT = Float64 - context = ClimaComms.SingletonCommsContext() - x = FT(1) - y = FT(2) - z = FT(3) - lat, long = FT(4), FT(5) - x1 = FT(1) - x2 = FT(2) - x3 = FT(3) - coords = [ - Geometry.ZPoint(z), - Geometry.XZPoint(x, z), - Geometry.XYZPoint(x, y, z), - Geometry.LatLongZPoint(lat, long, z), - Geometry.Cartesian3Point(x3), - Geometry.Cartesian13Point(x1, x3), - Geometry.Cartesian123Point(x1, x2, x3), - ] - all_components = [ - SMatrix{1, 1, FT}(range(1, 1)...), - SMatrix{2, 2, FT}(range(1, 4)...), - SMatrix{3, 3, FT}(range(1, 9)...), - SMatrix{3, 3, FT}(range(1, 9)...), - SMatrix{1, 1, FT}(range(1, 1)...), - SMatrix{2, 2, FT}(range(1, 4)...), - SMatrix{3, 3, FT}(range(1, 9)...), - ] - - expected_dzs = [1.0, 4.0, 9.0, 9.0, 1.0, 4.0, 9.0] - - for (components, coord, expected_dz) in - zip(all_components, coords, expected_dzs) - CoordType = typeof(coord) - AIdx = Geometry.coordinate_axis(CoordType) - at = Geometry.AxisTensor( - (Geometry.LocalAxis{AIdx}(), Geometry.CovariantAxis{AIdx}()), - components, - ) - local_geometry = Geometry.LocalGeometry(coord, FT(1.0), FT(1.0), at) - space = Spaces.PointSpace(context, local_geometry) - dz_computed = Array(parent(Fields.Δz_field(space))) - @test length(dz_computed) == 1 - @test dz_computed[1] == expected_dz - end -end - -@testset "scalar assignment" begin - device = ClimaComms.CPUSingleThreaded() # constructing space_vijfh is broken - context = ClimaComms.SingletonCommsContext(device) - domain_z = Domains.IntervalDomain( - Geometry.ZPoint(-1.0) .. Geometry.ZPoint(1.0), - periodic = true, - ) - mesh_z = Meshes.IntervalMesh(domain_z; nelems = 10) - topology_z = Topologies.IntervalTopology(context, mesh_z) - - domain_x = Domains.IntervalDomain( - Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0), - periodic = true, - ) - mesh_x = Meshes.IntervalMesh(domain_x; nelems = 10) - topology_x = Topologies.IntervalTopology(context, mesh_x) - - domain_xy = Domains.RectangleDomain( - Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0), - Geometry.YPoint(-1.0) .. Geometry.YPoint(1.0), - x1periodic = true, - x2periodic = true, - ) - mesh_xy = Meshes.RectilinearMesh(domain_xy, 10, 10) - topology_xy = Topologies.Topology2D(context, mesh_xy) - - quad = Spaces.Quadratures.GLL{4}() - - space_vf = Spaces.CenterFiniteDifferenceSpace(topology_z) - space_ifh = Spaces.SpectralElementSpace1D(topology_x, quad) - space_ijfh = Spaces.SpectralElementSpace2D(topology_xy, quad) - space_vifh = Spaces.ExtrudedFiniteDifferenceSpace(space_ifh, space_vf) - space_vijfh = Spaces.ExtrudedFiniteDifferenceSpace(space_ijfh, space_vf) - - C = map(x -> Geometry.Covariant12Vector(1.0, 1.0), zeros(space_vifh)) - @test all(==(1.0), parent(C)) - C .= Ref(zero(eltype(C))) - @test all(==(0.0), parent(C)) -end +# function spectral_space_2D(; n1 = 1, n2 = 1, Nij = 4) +# domain = Domains.RectangleDomain( +# Geometry.XPoint(-3.0) .. Geometry.XPoint(5.0), +# Geometry.YPoint(-2.0) .. Geometry.YPoint(8.0), +# x1periodic = false, +# x2periodic = false, +# x1boundary = (:east, :west), +# x2boundary = (:south, :north), +# ) +# mesh = Meshes.RectilinearMesh(domain, n1, n2) +# device = ClimaComms.CPUSingleThreaded() +# grid_topology = +# Topologies.Topology2D(ClimaComms.SingletonCommsContext(device), mesh) + +# quad = Spaces.Quadratures.GLL{Nij}() +# space = Spaces.SpectralElementSpace2D(grid_topology, quad) +# return space +# end + +# @testset "1×1 2D domain space" begin +# Nij = 4 +# n1 = n2 = 1 +# space = spectral_space_2D(n1 = n1, n2 = n2, Nij = Nij) + +# field = +# Fields.Field(IJFH{ComplexF64, Nij}(ones(Nij, Nij, 2, n1 * n2)), space) + +# @test sum(field) ≈ Complex(1.0, 1.0) * 8.0 * 10.0 rtol = 10eps() +# @test sum(x -> 3.0, field) ≈ 3 * 8.0 * 10.0 rtol = 10eps() +# @test mean(field) ≈ Complex(1.0, 1.0) rtol = 10eps() +# @test mean(x -> 3.0, field) ≈ 3 rtol = 10eps() +# @test norm(field) ≈ sqrt(2.0) rtol = 10eps() +# @test norm(field, 1) ≈ norm(Complex(1.0, 1.0)) rtol = 10eps() +# @test norm(field, Inf) ≈ norm(Complex(1.0, 1.0)) rtol = 10eps() +# @test norm(field; normalize = false) ≈ sqrt(2.0 * 8.0 * 10.0) rtol = 10eps() +# @test norm(field, 1; normalize = false) ≈ +# norm(Complex(1.0, 1.0)) * 8.0 * 10.0 rtol = 10eps() +# @test norm(field, Inf; normalize = false) ≈ norm(Complex(1.0, 1.0)) rtol = +# 10eps() + +# @test extrema(real, field) == (1.0, 1.0) + +# @test Operators.matrix_interpolate(field, 4) ≈ +# [Complex(1.0, 1.0) for i in 1:(4 * n1), j in 1:(4 * n2)] + + +# field_sin = map(x -> sin((x.x) / 2), Fields.coordinate_field(space)) +# M = Operators.matrix_interpolate(field_sin, 20) +# @test size(M) == (20, 20) # 20 x 20 for a 1 element field + +# real_field = field.re + +# # test broadcasting +# res = field .+ 1 +# @test parent(Fields.field_values(res)) == Float64[ +# f == 1 ? 2 : 1 for i in 1:Nij, j in 1:Nij, f in 1:2, h in 1:(n1 * n2) +# ] + +# res = field.re .+ 1 +# @test parent(Fields.field_values(res)) == +# Float64[2 for i in 1:Nij, j in 1:Nij, f in 1:1, h in 1:(n1 * n2)] + +# # test field slab broadcasting +# f1 = ones(space) +# f2 = ones(space) + +# for h in 1:(n1 * n2) +# f1_slab = Fields.slab(f1, h) +# f2_slab = Fields.slab(f2, h) +# q = f1_slab .+ f2_slab +# f1_slab .= q .+ f2_slab +# end +# @test all(parent(f1) .== 3) + +# point_field = Fields.column(field, 1, 1, 1) +# @test axes(point_field) isa Spaces.PointSpace +# end + +# # https://github.com/CliMA/ClimaCore.jl/issues/1126 +# function pow_n(f) +# @. f.x = f.x^2 +# return nothing +# end +# @testset "Broadcasting with ^n" begin +# FT = Float32 +# device = ClimaComms.CPUSingleThreaded() # fill is broken on gpu +# context = ClimaComms.SingletonCommsContext(device) +# for space in TU.all_spaces(FT; context) +# f = fill((; x = FT(1)), space) +# pow_n(f) # Compile first +# p_allocated = @allocated pow_n(f) +# if space isa Spaces.SpectralElementSpace1D +# @test p_allocated == 0 +# else +# @test p_allocated == 0 broken = (device isa ClimaComms.CUDADevice) +# end +# end +# end + +# function ifelse_broadcast_allocating(a, b, c) +# FT = eltype(a) +# @. a = ifelse(true || c < b * FT(1), FT(0), c) +# return nothing +# end + +# function ifelse_broadcast_or(a, b, c) +# FT = eltype(a) +# val = FT(1) +# @. a = ifelse(true || c < b * val, FT(0), c) +# return nothing +# end + +# function ifelse_broadcast_simple(a, b, c) +# FT = eltype(a) +# @. a = ifelse(c < b * FT(1), FT(0), c) +# return nothing +# end + +# @testset "Broadcasting ifelse" begin +# FT = Float32 +# device = ClimaComms.CPUSingleThreaded() # broken on gpu +# context = ClimaComms.SingletonCommsContext(device) +# for space in ( +# TU.CenterExtrudedFiniteDifferenceSpace(FT; context), +# TU.ColumnCenterFiniteDifferenceSpace(FT; context), +# ) +# a = Fields.level(fill(FT(0), space), 1) +# b = Fields.level(fill(FT(2), space), 1) +# c = Fields.level(fill(FT(3), space), 1) + +# ifelse_broadcast_allocating(a, b, c) +# p_allocated = @allocated ifelse_broadcast_allocating(a, b, c) +# @test_broken p_allocated == 0 + +# ifelse_broadcast_or(a, b, c) +# p_allocated = @allocated ifelse_broadcast_or(a, b, c) +# @test p_allocated == 0 + +# ifelse_broadcast_simple(a, b, c) +# p_allocated = @allocated ifelse_broadcast_simple(a, b, c) +# @test p_allocated == 0 +# end +# end + +# # Requires `--check-bounds=yes` +# @testset "Constructing & broadcasting over empty fields" begin +# FT = Float32 +# for space in TU.all_spaces(FT) +# f = fill((;), space) +# @. f += f +# end + +# function test_broken_throws(f) +# try +# @. f += 1 +# # we want to throw exception, test is broken +# @test_broken false +# catch +# # we want to throw exception, unexpected pass +# @test_broken true +# end +# end +# empty_field(space) = fill((;), space) + +# # Broadcasting over the wrong size should error +# test_broken_throws(empty_field(TU.PointSpace(FT))) +# test_broken_throws(empty_field(TU.SpectralElementSpace1D(FT))) +# test_broken_throws(empty_field(TU.SpectralElementSpace2D(FT))) +# test_broken_throws(empty_field(TU.ColumnCenterFiniteDifferenceSpace(FT))) +# test_broken_throws(empty_field(TU.ColumnFaceFiniteDifferenceSpace(FT))) +# test_broken_throws(empty_field(TU.SphereSpectralElementSpace(FT))) +# test_broken_throws(empty_field(TU.CenterExtrudedFiniteDifferenceSpace(FT))) +# test_broken_throws(empty_field(TU.FaceExtrudedFiniteDifferenceSpace(FT))) + +# # TODO: performance optimization: shouldn't we do +# # nothing when broadcasting over empty fields? +# # This is otherwise a performance penalty if +# # users regularly rely on empty fields. In particular: +# # - does iterating over empty fields load data? +# # - what is the overhead in iterating over empty fields? +# # - what is the use case of anything useful that can be +# # done by iterating over empty fields? +# end + +# @testset "Broadcasting interception for tuple-valued fields" begin +# n1 = n2 = 1 +# Nij = 4 +# space = spectral_space_2D(n1 = n1, n2 = n2, Nij = Nij) + +# nt_field = Fields.Field( +# IJFH{NamedTuple{(:a, :b), Tuple{Float64, Float64}}, Nij}( +# ones(Nij, Nij, 2, n1 * n2), +# ), +# space, +# ) +# nt_sum = sum(nt_field) +# @test nt_sum isa NamedTuple{(:a, :b), Tuple{Float64, Float64}} +# @test nt_sum.a ≈ 8.0 * 10.0 rtol = 10eps() +# @test nt_sum.b ≈ 8.0 * 10.0 rtol = 10eps() +# @test norm(nt_field) ≈ sqrt(2.0) rtol = 10eps() + +# # test scalar asignment +# nt_field.a .= 0.0 +# @test sum(nt_field.a) == 0.0 +# end + +# @testset "Special case handling for broadcased norm to pass through space local geometry" begin +# space = spectral_space_2D() +# u = Geometry.Covariant12Vector.(ones(space), ones(space)) +# @test norm.(u) ≈ hypot(4 / 8 / 2, 4 / 10 / 2) .* ones(space) +# end + +# @testset "FieldVector" begin +# space = spectral_space_2D() +# u = Geometry.Covariant12Vector.(ones(space), ones(space)) +# x = Fields.coordinate_field(space) +# y = [1.0, 2.0, 3.0] +# z = 1.0 +# Y = Fields.FieldVector(u = u, k = (x = x, y = y, z = z)) + +# @test propertynames(Y) == (:u, :k) +# @test propertynames(Y.k) == (:x, :y, :z) +# @test Y.u === u +# @test Y.k.x === x +# @test Y.k.y === y +# @test Y.k.z === z + +# @test deepcopy(Y).u !== u +# @test deepcopy(Y).k.x !== x +# @test deepcopy(Y).k.y !== y + +# @test getfield(deepcopy(Y).u, :space) === space + +# Y1 = 2 .* Y +# @test parent(Y1.u) == 2 .* parent(u) +# @test parent(Y1.k.x) == 2 .* parent(x) +# @test Y1.k.y == 2 .* y +# @test Y1.k.z === 2 * z + +# Y1 .= Y1 .+ 2 .* Y +# @test parent(Y1.u) == 4 .* parent(u) +# @test parent(Y1.k.x) == 4 .* parent(x) +# @test Y1.k.y == 4 .* y +# @test Y1.k.z === 4 * z + +# Y.k.z = 3.0 +# @test Y.k.z === 3.0 +# end + +# @testset "FieldVector array_type" begin +# device = ClimaComms.device() +# context = ClimaComms.SingletonCommsContext(device) +# space = TU.PointSpace(Float32; context) +# xcenters = Fields.coordinate_field(space).x +# y = Fields.FieldVector(x = xcenters) +# @test ClimaComms.array_type(y) == ClimaComms.array_type(device) +# y = Fields.FieldVector(x = xcenters, y = xcenters) +# @test ClimaComms.array_type(y) == ClimaComms.array_type(device) +# end + +# @testset "FieldVector basetype replacement and deepcopy" begin +# device = ClimaComms.CPUSingleThreaded() # constructing space_vijfh is broken +# context = ClimaComms.SingletonCommsContext(device) +# domain_z = Domains.IntervalDomain( +# Geometry.ZPoint(-1.0) .. Geometry.ZPoint(1.0), +# periodic = true, +# ) +# mesh_z = Meshes.IntervalMesh(domain_z; nelems = 10) +# topology_z = Topologies.IntervalTopology(context, mesh_z) + +# domain_x = Domains.IntervalDomain( +# Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0), +# periodic = true, +# ) +# mesh_x = Meshes.IntervalMesh(domain_x; nelems = 10) +# topology_x = Topologies.IntervalTopology(context, mesh_x) + +# domain_xy = Domains.RectangleDomain( +# Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0), +# Geometry.YPoint(-1.0) .. Geometry.YPoint(1.0), +# x1periodic = true, +# x2periodic = true, +# ) +# mesh_xy = Meshes.RectilinearMesh(domain_xy, 10, 10) +# topology_xy = Topologies.Topology2D(context, mesh_xy) + +# quad = Spaces.Quadratures.GLL{4}() + +# space_vf = Spaces.CenterFiniteDifferenceSpace(topology_z) +# space_ifh = Spaces.SpectralElementSpace1D(topology_x, quad) +# space_ijfh = Spaces.SpectralElementSpace2D(topology_xy, quad) +# space_vifh = Spaces.ExtrudedFiniteDifferenceSpace(space_ifh, space_vf) +# space_vijfh = Spaces.ExtrudedFiniteDifferenceSpace(space_ijfh, space_vf) + +# space2field(space) = map( +# coord -> (coord, Geometry.Covariant12Vector(1.0, 2.0)), +# Fields.coordinate_field(space), +# ) + +# Y = Fields.FieldVector( +# field_vf = space2field(space_vf), +# field_if = slab(space2field(space_ifh), 1), +# field_ifh = space2field(space_ifh), +# field_ijf = slab(space2field(space_ijfh), 1, 1), +# field_ijfh = space2field(space_ijfh), +# field_vifh = space2field(space_vifh), +# field_vijfh = space2field(space_vijfh), +# array = [1.0, 2.0, 3.0], +# scalar = 1.0, +# ) + +# Yf = ForwardDiff.Dual{Nothing}.(Y, 1.0) +# Yf .= Yf .^ 2 .+ Y +# @test all(ForwardDiff.value.(Yf) .== Y .^ 2 .+ Y) +# @test all(ForwardDiff.partials.(Yf, 1) .== 2 .* Y) + +# dual_field = Yf.field_vf +# dual_field_original_basetype = similar(Y.field_vf, eltype(dual_field)) +# @test eltype(dual_field_original_basetype) === eltype(dual_field) +# @test eltype(parent(dual_field_original_basetype)) === Float64 +# @test eltype(parent(dual_field)) === ForwardDiff.Dual{Nothing, Float64, 1} + +# object_that_contains_Yf = (; Yf) +# @test axes(deepcopy(Yf).field_vf) === space_vf +# @test axes(deepcopy(object_that_contains_Yf).Yf.field_vf) === space_vf +# end + +# @testset "Scalar field iterator" begin +# space = spectral_space_2D() +# u = Geometry.Covariant12Vector.(ones(space), ones(space)) +# x = Fields.coordinate_field(space) +# y = [1.0, 2.0, 3.0] +# z = 1.0 +# Y = Fields.FieldVector(u = u, k = (x = x, y = y, z = z)) + +# prop_chains = Fields.property_chains(Y) +# @test length(prop_chains) == 6 +# @test prop_chains[1] == (:u, :components, :data, 1) +# @test prop_chains[2] == (:u, :components, :data, 2) +# @test prop_chains[3] == (:k, :x, :x) +# @test prop_chains[4] == (:k, :x, :y) +# @test prop_chains[5] == (:k, :y) +# @test prop_chains[6] == (:k, :z) + +# FT = Float64 +# nt = +# (; x = FT(0), y = FT(0), tup = ntuple(i -> (; a = FT(1), b = FT(1)), 2)) +# Y = fill(nt, space) + +# prop_chains = Fields.property_chains(Y) +# @test prop_chains[1] == (:x,) +# @test prop_chains[2] == (:y,) +# @test prop_chains[3] == (:tup, 1, :a) +# @test prop_chains[4] == (:tup, 1, :b) +# @test prop_chains[5] == (:tup, 2, :a) +# @test prop_chains[6] == (:tup, 2, :b) + +# @test Fields.single_field(Y, prop_chains[1]) === Y.x +# @test Fields.single_field(Y, prop_chains[2]) === Y.y +# @test Fields.single_field(Y, prop_chains[3]) === getproperty(Y.tup, 1).a +# @test Fields.single_field(Y, prop_chains[4]) === getproperty(Y.tup, 1).b +# @test Fields.single_field(Y, prop_chains[5]) === getproperty(Y.tup, 2).a +# @test Fields.single_field(Y, prop_chains[6]) === getproperty(Y.tup, 2).b + +# for (i, (var, prop_chain)) in enumerate(Fields.field_iterator(Y)) +# @test prop_chains[i] == prop_chain +# @test var === Fields.single_field(Y, prop_chain) +# end +# end + +# # Test truncated field type printing: +# ClimaCore.Fields.truncate_printing_field_types() = true +# @testset "Truncated printing" begin +# nt = (; x = Float64(0), y = Float64(0)) +# Y = fill(nt, spectral_space_2D()) +# @test sprint(show, typeof(Y); context = IOContext(stdout)) == +# "Field{(:x, :y)} (trunc disp)" +# end +# ClimaCore.Fields.truncate_printing_field_types() = false + +# @testset "Standard printing" begin +# nt = (; x = Float64(0), y = Float64(0)) +# Y = fill(nt, spectral_space_2D()) +# s = sprint(show, typeof(Y)) # just make sure this doesn't break +# end + +# @testset "Set!" begin +# space = spectral_space_2D() +# FT = Float64 +# nt = (; x = FT(0), y = FT(0)) +# Y = fill(nt, space) +# foo(local_geom) = +# sin(local_geom.coordinates.x * local_geom.coordinates.y) + 3 +# Fields.set!(foo, Y.x) +# @test all((parent(Y.x) .> 1)) +# end + +# @testset "PointField" begin +# device = ClimaComms.CPUSingleThreaded() # a bunch of cuda pieces are broken +# context = ClimaComms.SingletonCommsContext(device) +# FT = Float64 +# coord = Geometry.XPoint(FT(π)) +# space = Spaces.PointSpace(context, coord) +# @test parent(Spaces.local_geometry_data(space)) == +# FT[Geometry.component(coord, 1), 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] +# field = Fields.coordinate_field(space) +# @test field isa Fields.PointField +# @test Fields.field_values(field)[] == coord + +# if ClimaComms.device(context) isa ClimaComms.AbstractCPUDevice +# @test sum(field.x) == FT(π) +# elseif ClimaComms.device(context) isa ClimaComms.CUDADevice +# # Not yet supported +# # @test sum(field.x) == FT(π) +# end + +# field = ones(space) .* π +# sin_field = sin.(field) +# add_field = field .+ field +# @test isapprox(Fields.field_values(sin_field)[], FT(0.0); atol = √eps(FT)) +# @test isapprox(Fields.field_values(add_field)[], FT(2π)) +# end + +# @testset "Level" begin +# FT = Float64 +# for space in TU.all_spaces(FT) +# TU.levelable(space) || continue +# Y = fill((; x = FT(2)), space) +# lg_space = Spaces.level(space, TU.fc_index(1, space)) +# lg_field_space = axes(Fields.level(Y, TU.fc_index(1, space))) +# @test all( +# lg_space.local_geometry.coordinates === +# lg_field_space.local_geometry.coordinates, +# ) +# @test all(Fields.zeros(lg_space) == Fields.zeros(lg_field_space)) +# end +# end + +# @testset "Points from Columns" begin +# FT = Float64 +# for space in TU.all_spaces(FT) +# if space isa Spaces.SpectralElementSpace1D +# Y = fill((; x = FT(1)), space) +# point_space_from_field = axes(Fields.column(Y.x, 1, 1)) +# point_space = Spaces.column(space, 1, 1) +# @test Fields.ones(point_space) == +# Fields.ones(point_space_from_field) +# end +# if space isa Spaces.SpectralElementSpace2D +# Y = fill((; x = FT(1)), space) +# point_space_from_field = axes(Fields.column(Y.x, 1, 1, 1)) +# point_space = Spaces.column(space, 1, 1, 1) +# @test Fields.ones(point_space) == +# Fields.ones(point_space_from_field) +# end + +# end +# end + +# @testset "(Domain/Column)-surface broadcasting" begin +# FT = Float64 +# function domain_surface_bc!(x, ᶜz_surf, ᶜx_surf) +# @. x = x + ᶜz_surf +# # exercises broadcast_shape(PointSpace, PointSpace) +# @. x = x + (ᶜz_surf * ᶜx_surf) +# nothing +# end +# function column_surface_bc!(x, ᶜz_surf, ᶜx_surf) +# Fields.bycolumn(axes(x)) do colidx +# @. x[colidx] = x[colidx] + ᶜz_surf[colidx] +# # exercises broadcast_shape(PointSpace, PointSpace) +# @. x[colidx] = x[colidx] + (ᶜz_surf[colidx] * ᶜx_surf[colidx]) +# end +# nothing +# end +# for space in TU.all_spaces(FT) +# # Filter out spaces without z coordinates: +# TU.has_z_coordinates(space) || continue +# Y = fill((; x = FT(1)), space) +# ᶜz_surf = +# Spaces.level(Fields.coordinate_field(Y).z, TU.fc_index(1, space)) +# ᶜx_surf = copy(Spaces.level(Y.x, TU.fc_index(1, space))) + +# # Still need to define broadcast rules for surface planes with 3D domains +# domain_surface_bc!(Y.x, ᶜz_surf, ᶜx_surf) + +# # Skip spaces incompatible with Fields.bycolumn: +# TU.bycolumnable(space) || continue +# Yc = fill((; x = FT(1)), space) +# column_surface_bc!(Yc.x, ᶜz_surf, ᶜx_surf) +# @test Y.x == Yc.x +# nothing +# end +# nothing +# end + +# @testset "Broadcasting same spaces different instances" begin +# space1 = spectral_space_2D() +# space2 = spectral_space_2D() +# field1 = ones(space1) +# field2 = 2 .* ones(space2) +# @test Fields.is_diagonalized_spaces(typeof(space1), typeof(space2)) +# @test_throws ErrorException( +# "Broacasted spaces are the same ClimaCore.Spaces type but not the same instance", +# ) field1 .= field2 + +# # turn warning on +# Fields.allow_mismatched_diagonalized_spaces() = true +# @test_warn "Broacasted spaces are the same ClimaCore.Spaces type but not the same instance" field1 .= +# field2 +# @test parent(field1) == parent(field2) +# Fields.allow_mismatched_diagonalized_spaces() = false +# end + +# struct InferenceFoo{FT} +# bar::FT +# end +# Base.broadcastable(x::InferenceFoo) = Ref(x) +# @testset "Inference failure message" begin +# function ics_foo(::Type{FT}, lg, foo) where {FT} +# uv = Geometry.UVVector(FT(0), FT(0)) +# z = Geometry.Covariant12Vector(uv, lg) +# y = foo.bingo +# return (; x = FT(0) + y) +# end +# function ics_foo_with_field(::Type{FT}, lg, foo, f) where {FT} +# uv = Geometry.UVVector(FT(0), FT(0)) +# z = Geometry.Covariant12Vector(uv, lg) +# ζ = f.a +# y = foo.baz +# return (; x = FT(0) + y - ζ) +# end +# function FieldFromNamedTupleBroken( +# space, +# ics::Function, +# ::Type{FT}, +# params..., +# ) where {FT} +# lg = Fields.local_geometry_field(space) +# return ics.(FT, lg, params...) +# end +# FT = Float64 +# foo = InferenceFoo(2.0) +# device = ClimaComms.CPUSingleThreaded() # cuda fill is broken +# context = ClimaComms.SingletonCommsContext(device) +# for space in TU.all_spaces(FT; context) +# Y = fill((; a = FT(0), b = FT(1)), space) +# @test_throws ErrorException("type InferenceFoo has no field bingo") FieldFromNamedTupleBroken( +# space, +# ics_foo, +# FT, +# foo, +# ) +# @test_throws ErrorException("type InferenceFoo has no field baz") FieldFromNamedTupleBroken( +# space, +# ics_foo_with_field, +# FT, +# foo, +# Y, +# ) +# end + +# end + +# @testset "Δz_field" begin +# FT = Float64 +# context = ClimaComms.SingletonCommsContext() +# x = FT(1) +# y = FT(2) +# z = FT(3) +# lat, long = FT(4), FT(5) +# x1 = FT(1) +# x2 = FT(2) +# x3 = FT(3) +# coords = [ +# Geometry.ZPoint(z), +# Geometry.XZPoint(x, z), +# Geometry.XYZPoint(x, y, z), +# Geometry.LatLongZPoint(lat, long, z), +# Geometry.Cartesian3Point(x3), +# Geometry.Cartesian13Point(x1, x3), +# Geometry.Cartesian123Point(x1, x2, x3), +# ] +# all_components = [ +# SMatrix{1, 1, FT}(range(1, 1)...), +# SMatrix{2, 2, FT}(range(1, 4)...), +# SMatrix{3, 3, FT}(range(1, 9)...), +# SMatrix{3, 3, FT}(range(1, 9)...), +# SMatrix{1, 1, FT}(range(1, 1)...), +# SMatrix{2, 2, FT}(range(1, 4)...), +# SMatrix{3, 3, FT}(range(1, 9)...), +# ] + +# expected_dzs = [1.0, 4.0, 9.0, 9.0, 1.0, 4.0, 9.0] + +# for (components, coord, expected_dz) in +# zip(all_components, coords, expected_dzs) +# CoordType = typeof(coord) +# AIdx = Geometry.coordinate_axis(CoordType) +# at = Geometry.AxisTensor( +# (Geometry.LocalAxis{AIdx}(), Geometry.CovariantAxis{AIdx}()), +# components, +# ) +# local_geometry = Geometry.LocalGeometry(coord, FT(1.0), FT(1.0), at) +# space = Spaces.PointSpace(context, local_geometry) +# dz_computed = Array(parent(Fields.Δz_field(space))) +# @test length(dz_computed) == 1 +# @test dz_computed[1] == expected_dz +# end +# end + +# @testset "scalar assignment" begin +# device = ClimaComms.CPUSingleThreaded() # constructing space_vijfh is broken +# context = ClimaComms.SingletonCommsContext(device) +# domain_z = Domains.IntervalDomain( +# Geometry.ZPoint(-1.0) .. Geometry.ZPoint(1.0), +# periodic = true, +# ) +# mesh_z = Meshes.IntervalMesh(domain_z; nelems = 10) +# topology_z = Topologies.IntervalTopology(context, mesh_z) + +# domain_x = Domains.IntervalDomain( +# Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0), +# periodic = true, +# ) +# mesh_x = Meshes.IntervalMesh(domain_x; nelems = 10) +# topology_x = Topologies.IntervalTopology(context, mesh_x) + +# domain_xy = Domains.RectangleDomain( +# Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0), +# Geometry.YPoint(-1.0) .. Geometry.YPoint(1.0), +# x1periodic = true, +# x2periodic = true, +# ) +# mesh_xy = Meshes.RectilinearMesh(domain_xy, 10, 10) +# topology_xy = Topologies.Topology2D(context, mesh_xy) + +# quad = Spaces.Quadratures.GLL{4}() + +# space_vf = Spaces.CenterFiniteDifferenceSpace(topology_z) +# space_ifh = Spaces.SpectralElementSpace1D(topology_x, quad) +# space_ijfh = Spaces.SpectralElementSpace2D(topology_xy, quad) +# space_vifh = Spaces.ExtrudedFiniteDifferenceSpace(space_ifh, space_vf) +# space_vijfh = Spaces.ExtrudedFiniteDifferenceSpace(space_ijfh, space_vf) + +# C = map(x -> Geometry.Covariant12Vector(1.0, 1.0), zeros(space_vifh)) +# @test all(==(1.0), parent(C)) +# C .= Ref(zero(eltype(C))) +# @test all(==(0.0), parent(C)) +# end function integrate_bycolumn!(∫y, Y) Fields.bycolumn(axes(Y.y)) do colidx Operators.column_integral_definite!(∫y[colidx], Y.y[colidx]) @@ -698,7 +698,7 @@ convergence_rate(err, Δh) = col_copy = similar(y[Fields.ColumnIndex((1, 1), 1)]) return Fields.Field(Fields.field_values(col_copy), axes(col_copy)) end - device = ClimaComms.CPUSingleThreaded() + device = ClimaComms.device() context = ClimaComms.SingletonCommsContext(device) for zelem in (2^2, 2^3, 2^4, 2^5) for space in TU.all_spaces(FT; zelem, context) @@ -706,6 +706,7 @@ convergence_rate(err, Δh) = TU.has_z_coordinates(space) || continue # Skip spaces incompatible with Fields.bycolumn: TU.bycolumnable(space) || continue + @show nameof(typeof(space)) Y = fill((; y = FT(1)), space) zcf = Fields.coordinate_field(Y.y).z From 9ac707c4432bfb8f7ea5728d63c746db36d6911d Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Wed, 19 Jul 2023 14:09:15 -0700 Subject: [PATCH 02/11] wip --- src/Operators/integrals.jl | 55 ++++++++++++++++++++++++++--------- src/Spaces/pointspace.jl | 2 +- src/Spaces/spectralelement.jl | 8 ++--- 3 files changed, 47 insertions(+), 18 deletions(-) diff --git a/src/Operators/integrals.jl b/src/Operators/integrals.jl index 8ef963da07..6725c30b73 100644 --- a/src/Operators/integrals.jl +++ b/src/Operators/integrals.jl @@ -9,7 +9,7 @@ The definite vertical column integral, `col∫field`, of field `field`. """ column_integral_definite!(col∫field::Fields.Field, field::Fields.Field) = column_integral_definite!( - ClimaComms.device(axes(col∫field)), + ClimaComms.device(axes(field)), col∫field, field, ) @@ -35,9 +35,10 @@ function column_integral_definite_kernel!( Ni, Nj, _, _, Nh = size(Fields.field_values(field)) if idx <= Ni * Nj * Nh i, j, h = Spaces._get_idx((Ni, Nj, Nh), idx) + colfield = Spaces.column(field, i, j, h) _column_integral_definite!( Spaces.column(col∫field, i, j, h), - Spaces.column(field, i, j, h), + colfield, ) end return nothing @@ -70,19 +71,47 @@ function _column_integral_definite!( col∫field::Fields.PointField, field::Fields.ColumnField, ) - @inbounds col∫field[] = _column_integral_definite(field) - return nothing -end + # field_data = Fields.field_values(field) + space = axes(field) + Δz_field = Fields.Δz_field(space) + Nv = Spaces.nlevels(space) + # d∑ = mapreduce(RecursiveApply.radd, 1:Nv) do idx + # reduction_getindex(field, idx) * reduction_getindex(Δz_field, idx) + # end -function _column_integral_definite(field::Fields.ColumnField) - field_data = Fields.field_values(field) - Δz_data = Spaces.Δz_data(axes(field)) - Nv = Spaces.nlevels(axes(field)) - ∫field = zero(eltype(field)) - @inbounds for j in 1:Nv - ∫field += field_data[j] * Δz_data[j] + col∫field .= 0 + @inbounds for idx in 1:Nv + # col∫field[] += field_data[j] * Δz_data[j] + # d∑ = getidx(space, field, Operators.Interior(), idx, hidx) * + # getidx(space, Δz_field, Operators.Interior(), idx, hidx) + # setidx!(space, col∫field, idx, hidx, d∑) + col∫field[] += reduction_getindex(field, idx) * reduction_getindex(Δz_field, idx) + # col∫field[] += field_data[j] * Δz_data[j] end - return ∫field + # reduction_setindex!(col∫field, 1, d∑) + # Operators.setidx!( + # axes(col∫field), + # col∫field, + # 1, + # (1, 1, 1), + # d∑, + # ) + return nothing end +reduction_getindex(column_field, index) = @inbounds getidx( + axes(column_field), + column_field, + Interior(), + index - 1 + left_idx(axes(column_field)), +) + +# reduction_setindex!(column_field, index, value) = @inbounds setidx!( +# axes(column_field), +# column_field, +# index - 1 + left_idx(axes(column_field)), +# (1, 1, 1), +# value, +# ) + # TODO: add support for indefinite integrals diff --git a/src/Spaces/pointspace.jl b/src/Spaces/pointspace.jl index d68ab4e13c..f346f45d89 100644 --- a/src/Spaces/pointspace.jl +++ b/src/Spaces/pointspace.jl @@ -7,7 +7,7 @@ local_geometry_data(space::AbstractPointSpace) = space.local_geometry A zero-dimensional space. """ -struct PointSpace{C <: ClimaComms.AbstractCommsContext, LG} <: +struct PointSpace{C <: Union{Nothing, ClimaComms.AbstractCommsContext}, LG} <: AbstractPointSpace context::C local_geometry::LG diff --git a/src/Spaces/spectralelement.jl b/src/Spaces/spectralelement.jl index 96ec37a70e..fcfff28649 100644 --- a/src/Spaces/spectralelement.jl +++ b/src/Spaces/spectralelement.jl @@ -644,16 +644,16 @@ Base.@propagate_inbounds slab(space::AbstractSpectralElementSpace, h) = Base.@propagate_inbounds function column(space::SpectralElementSpace1D, i, h) local_geometry = column(local_geometry_data(space), i, h) - context = ClimaComms.context(space) - PointSpace(context, local_geometry) + # context = ClimaComms.context(space) + PointSpace(nothing, local_geometry) end Base.@propagate_inbounds column(space::SpectralElementSpace1D, i, j, h) = column(space, i, h) Base.@propagate_inbounds function column(space::SpectralElementSpace2D, i, j, h) local_geometry = column(local_geometry_data(space), i, j, h) - context = ClimaComms.context(space) - PointSpace(context, local_geometry) + # context = ClimaComms.context(space) + PointSpace(nothing, local_geometry) end # XXX: this cannot take `space` as it must be constructed beforehand so From 7c490c3714b645566406b29a98be5c9839363622 Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Wed, 19 Jul 2023 14:47:15 -0700 Subject: [PATCH 03/11] remove context from PointSpace --- src/Operators/integrals.jl | 2 +- src/Spaces/finitedifference.jl | 6 ++---- src/Spaces/pointspace.jl | 10 +++++----- src/Spaces/spectralelement.jl | 6 ++---- 4 files changed, 10 insertions(+), 14 deletions(-) diff --git a/src/Operators/integrals.jl b/src/Operators/integrals.jl index 6725c30b73..55e9429fb6 100644 --- a/src/Operators/integrals.jl +++ b/src/Operators/integrals.jl @@ -79,7 +79,7 @@ function _column_integral_definite!( # reduction_getindex(field, idx) * reduction_getindex(Δz_field, idx) # end - col∫field .= 0 + col∫field[] = 0 @inbounds for idx in 1:Nv # col∫field[] += field_data[j] * Δz_data[j] # d∑ = getidx(space, field, Operators.Interior(), idx, hidx) * diff --git a/src/Spaces/finitedifference.jl b/src/Spaces/finitedifference.jl index bdffe71dcf..26844fdeb5 100644 --- a/src/Spaces/finitedifference.jl +++ b/src/Spaces/finitedifference.jl @@ -212,14 +212,12 @@ Base.@propagate_inbounds function level( v::PlusHalf, ) @inbounds local_geometry = level(local_geometry_data(space), v.i + 1) - context = ClimaComms.context(space) - PointSpace(context, local_geometry) + PointSpace(local_geometry) end Base.@propagate_inbounds function level( space::CenterFiniteDifferenceSpace, v::Int, ) local_geometry = level(local_geometry_data(space), v) - context = ClimaComms.context(space) - PointSpace(context, local_geometry) + PointSpace(local_geometry) end diff --git a/src/Spaces/pointspace.jl b/src/Spaces/pointspace.jl index f346f45d89..5dae751d51 100644 --- a/src/Spaces/pointspace.jl +++ b/src/Spaces/pointspace.jl @@ -7,19 +7,19 @@ local_geometry_data(space::AbstractPointSpace) = space.local_geometry A zero-dimensional space. """ -struct PointSpace{C <: Union{Nothing, ClimaComms.AbstractCommsContext}, LG} <: - AbstractPointSpace - context::C +struct PointSpace{LG <: DataLayouts.Data0D} <: AbstractPointSpace local_geometry::LG end ClimaComms.device(space::PointSpace) = ClimaComms.device(space.context) -ClimaComms.context(space::PointSpace) = space.context +ClimaComms.context(space::PointSpace) = ClimaComms.SingletonCommsContext(ClimaComms.CPUSingleThreaded()) +#= PointSpace(x::Geometry.LocalGeometry) = PointSpace(ClimaComms.CPUSingleThreaded(), x) PointSpace(x::Geometry.AbstractPoint) = PointSpace(ClimaComms.CPUSingleThreaded(), x) +=# function PointSpace(device::ClimaComms.AbstractDevice, x) context = ClimaComms.SingletonCommsContext(device) @@ -34,7 +34,7 @@ function PointSpace( ArrayType = ClimaComms.array_type(ClimaComms.device(context)) local_geometry_data = DataLayouts.DataF{LG}(Array{FT}) local_geometry_data[] = local_geometry - return PointSpace(context, Adapt.adapt(ArrayType, local_geometry_data)) + return PointSpace(Adapt.adapt(ArrayType, local_geometry_data)) end function PointSpace( diff --git a/src/Spaces/spectralelement.jl b/src/Spaces/spectralelement.jl index fcfff28649..81e45ca5ff 100644 --- a/src/Spaces/spectralelement.jl +++ b/src/Spaces/spectralelement.jl @@ -644,16 +644,14 @@ Base.@propagate_inbounds slab(space::AbstractSpectralElementSpace, h) = Base.@propagate_inbounds function column(space::SpectralElementSpace1D, i, h) local_geometry = column(local_geometry_data(space), i, h) - # context = ClimaComms.context(space) - PointSpace(nothing, local_geometry) + PointSpace(local_geometry) end Base.@propagate_inbounds column(space::SpectralElementSpace1D, i, j, h) = column(space, i, h) Base.@propagate_inbounds function column(space::SpectralElementSpace2D, i, j, h) local_geometry = column(local_geometry_data(space), i, j, h) - # context = ClimaComms.context(space) - PointSpace(nothing, local_geometry) + PointSpace(local_geometry) end # XXX: this cannot take `space` as it must be constructed beforehand so From 1e1e1a0b4d1d79033c2bb8541fba81c81151f1e5 Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Wed, 19 Jul 2023 15:01:53 -0700 Subject: [PATCH 04/11] formatter and cleanup --- src/Operators/integrals.jl | 39 ++++---------------------------------- src/Spaces/pointspace.jl | 3 ++- 2 files changed, 6 insertions(+), 36 deletions(-) diff --git a/src/Operators/integrals.jl b/src/Operators/integrals.jl index 55e9429fb6..9dddf35d35 100644 --- a/src/Operators/integrals.jl +++ b/src/Operators/integrals.jl @@ -8,11 +8,7 @@ The definite vertical column integral, `col∫field`, of field `field`. """ column_integral_definite!(col∫field::Fields.Field, field::Fields.Field) = - column_integral_definite!( - ClimaComms.device(axes(field)), - col∫field, - field, - ) + column_integral_definite!(ClimaComms.device(axes(field)), col∫field, field) function column_integral_definite!( ::ClimaComms.CUDADevice, @@ -36,10 +32,7 @@ function column_integral_definite_kernel!( if idx <= Ni * Nj * Nh i, j, h = Spaces._get_idx((Ni, Nj, Nh), idx) colfield = Spaces.column(field, i, j, h) - _column_integral_definite!( - Spaces.column(col∫field, i, j, h), - colfield, - ) + _column_integral_definite!(Spaces.column(col∫field, i, j, h), colfield) end return nothing end @@ -71,31 +64,15 @@ function _column_integral_definite!( col∫field::Fields.PointField, field::Fields.ColumnField, ) - # field_data = Fields.field_values(field) space = axes(field) Δz_field = Fields.Δz_field(space) Nv = Spaces.nlevels(space) - # d∑ = mapreduce(RecursiveApply.radd, 1:Nv) do idx - # reduction_getindex(field, idx) * reduction_getindex(Δz_field, idx) - # end col∫field[] = 0 @inbounds for idx in 1:Nv - # col∫field[] += field_data[j] * Δz_data[j] - # d∑ = getidx(space, field, Operators.Interior(), idx, hidx) * - # getidx(space, Δz_field, Operators.Interior(), idx, hidx) - # setidx!(space, col∫field, idx, hidx, d∑) - col∫field[] += reduction_getindex(field, idx) * reduction_getindex(Δz_field, idx) - # col∫field[] += field_data[j] * Δz_data[j] + col∫field[] += + reduction_getindex(field, idx) * reduction_getindex(Δz_field, idx) end - # reduction_setindex!(col∫field, 1, d∑) - # Operators.setidx!( - # axes(col∫field), - # col∫field, - # 1, - # (1, 1, 1), - # d∑, - # ) return nothing end @@ -106,12 +83,4 @@ reduction_getindex(column_field, index) = @inbounds getidx( index - 1 + left_idx(axes(column_field)), ) -# reduction_setindex!(column_field, index, value) = @inbounds setidx!( -# axes(column_field), -# column_field, -# index - 1 + left_idx(axes(column_field)), -# (1, 1, 1), -# value, -# ) - # TODO: add support for indefinite integrals diff --git a/src/Spaces/pointspace.jl b/src/Spaces/pointspace.jl index 5dae751d51..82b476c53d 100644 --- a/src/Spaces/pointspace.jl +++ b/src/Spaces/pointspace.jl @@ -12,7 +12,8 @@ struct PointSpace{LG <: DataLayouts.Data0D} <: AbstractPointSpace end ClimaComms.device(space::PointSpace) = ClimaComms.device(space.context) -ClimaComms.context(space::PointSpace) = ClimaComms.SingletonCommsContext(ClimaComms.CPUSingleThreaded()) +ClimaComms.context(space::PointSpace) = + ClimaComms.SingletonCommsContext(ClimaComms.CPUSingleThreaded()) #= PointSpace(x::Geometry.LocalGeometry) = From 3f637b0a88e12b85e201e4eddfc7fe45a58639d4 Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Wed, 19 Jul 2023 15:02:52 -0700 Subject: [PATCH 05/11] Revert test commenting --- test/Fields/field.jl | 1303 +++++++++++++++++++++--------------------- 1 file changed, 651 insertions(+), 652 deletions(-) diff --git a/test/Fields/field.jl b/test/Fields/field.jl index fd6a6f0638..cb711f0306 100644 --- a/test/Fields/field.jl +++ b/test/Fields/field.jl @@ -17,656 +17,656 @@ include( ) import .TestUtilities as TU -# function spectral_space_2D(; n1 = 1, n2 = 1, Nij = 4) -# domain = Domains.RectangleDomain( -# Geometry.XPoint(-3.0) .. Geometry.XPoint(5.0), -# Geometry.YPoint(-2.0) .. Geometry.YPoint(8.0), -# x1periodic = false, -# x2periodic = false, -# x1boundary = (:east, :west), -# x2boundary = (:south, :north), -# ) -# mesh = Meshes.RectilinearMesh(domain, n1, n2) -# device = ClimaComms.CPUSingleThreaded() -# grid_topology = -# Topologies.Topology2D(ClimaComms.SingletonCommsContext(device), mesh) - -# quad = Spaces.Quadratures.GLL{Nij}() -# space = Spaces.SpectralElementSpace2D(grid_topology, quad) -# return space -# end - -# @testset "1×1 2D domain space" begin -# Nij = 4 -# n1 = n2 = 1 -# space = spectral_space_2D(n1 = n1, n2 = n2, Nij = Nij) - -# field = -# Fields.Field(IJFH{ComplexF64, Nij}(ones(Nij, Nij, 2, n1 * n2)), space) - -# @test sum(field) ≈ Complex(1.0, 1.0) * 8.0 * 10.0 rtol = 10eps() -# @test sum(x -> 3.0, field) ≈ 3 * 8.0 * 10.0 rtol = 10eps() -# @test mean(field) ≈ Complex(1.0, 1.0) rtol = 10eps() -# @test mean(x -> 3.0, field) ≈ 3 rtol = 10eps() -# @test norm(field) ≈ sqrt(2.0) rtol = 10eps() -# @test norm(field, 1) ≈ norm(Complex(1.0, 1.0)) rtol = 10eps() -# @test norm(field, Inf) ≈ norm(Complex(1.0, 1.0)) rtol = 10eps() -# @test norm(field; normalize = false) ≈ sqrt(2.0 * 8.0 * 10.0) rtol = 10eps() -# @test norm(field, 1; normalize = false) ≈ -# norm(Complex(1.0, 1.0)) * 8.0 * 10.0 rtol = 10eps() -# @test norm(field, Inf; normalize = false) ≈ norm(Complex(1.0, 1.0)) rtol = -# 10eps() - -# @test extrema(real, field) == (1.0, 1.0) - -# @test Operators.matrix_interpolate(field, 4) ≈ -# [Complex(1.0, 1.0) for i in 1:(4 * n1), j in 1:(4 * n2)] - - -# field_sin = map(x -> sin((x.x) / 2), Fields.coordinate_field(space)) -# M = Operators.matrix_interpolate(field_sin, 20) -# @test size(M) == (20, 20) # 20 x 20 for a 1 element field - -# real_field = field.re - -# # test broadcasting -# res = field .+ 1 -# @test parent(Fields.field_values(res)) == Float64[ -# f == 1 ? 2 : 1 for i in 1:Nij, j in 1:Nij, f in 1:2, h in 1:(n1 * n2) -# ] - -# res = field.re .+ 1 -# @test parent(Fields.field_values(res)) == -# Float64[2 for i in 1:Nij, j in 1:Nij, f in 1:1, h in 1:(n1 * n2)] - -# # test field slab broadcasting -# f1 = ones(space) -# f2 = ones(space) - -# for h in 1:(n1 * n2) -# f1_slab = Fields.slab(f1, h) -# f2_slab = Fields.slab(f2, h) -# q = f1_slab .+ f2_slab -# f1_slab .= q .+ f2_slab -# end -# @test all(parent(f1) .== 3) - -# point_field = Fields.column(field, 1, 1, 1) -# @test axes(point_field) isa Spaces.PointSpace -# end - -# # https://github.com/CliMA/ClimaCore.jl/issues/1126 -# function pow_n(f) -# @. f.x = f.x^2 -# return nothing -# end -# @testset "Broadcasting with ^n" begin -# FT = Float32 -# device = ClimaComms.CPUSingleThreaded() # fill is broken on gpu -# context = ClimaComms.SingletonCommsContext(device) -# for space in TU.all_spaces(FT; context) -# f = fill((; x = FT(1)), space) -# pow_n(f) # Compile first -# p_allocated = @allocated pow_n(f) -# if space isa Spaces.SpectralElementSpace1D -# @test p_allocated == 0 -# else -# @test p_allocated == 0 broken = (device isa ClimaComms.CUDADevice) -# end -# end -# end - -# function ifelse_broadcast_allocating(a, b, c) -# FT = eltype(a) -# @. a = ifelse(true || c < b * FT(1), FT(0), c) -# return nothing -# end - -# function ifelse_broadcast_or(a, b, c) -# FT = eltype(a) -# val = FT(1) -# @. a = ifelse(true || c < b * val, FT(0), c) -# return nothing -# end - -# function ifelse_broadcast_simple(a, b, c) -# FT = eltype(a) -# @. a = ifelse(c < b * FT(1), FT(0), c) -# return nothing -# end - -# @testset "Broadcasting ifelse" begin -# FT = Float32 -# device = ClimaComms.CPUSingleThreaded() # broken on gpu -# context = ClimaComms.SingletonCommsContext(device) -# for space in ( -# TU.CenterExtrudedFiniteDifferenceSpace(FT; context), -# TU.ColumnCenterFiniteDifferenceSpace(FT; context), -# ) -# a = Fields.level(fill(FT(0), space), 1) -# b = Fields.level(fill(FT(2), space), 1) -# c = Fields.level(fill(FT(3), space), 1) - -# ifelse_broadcast_allocating(a, b, c) -# p_allocated = @allocated ifelse_broadcast_allocating(a, b, c) -# @test_broken p_allocated == 0 - -# ifelse_broadcast_or(a, b, c) -# p_allocated = @allocated ifelse_broadcast_or(a, b, c) -# @test p_allocated == 0 - -# ifelse_broadcast_simple(a, b, c) -# p_allocated = @allocated ifelse_broadcast_simple(a, b, c) -# @test p_allocated == 0 -# end -# end - -# # Requires `--check-bounds=yes` -# @testset "Constructing & broadcasting over empty fields" begin -# FT = Float32 -# for space in TU.all_spaces(FT) -# f = fill((;), space) -# @. f += f -# end - -# function test_broken_throws(f) -# try -# @. f += 1 -# # we want to throw exception, test is broken -# @test_broken false -# catch -# # we want to throw exception, unexpected pass -# @test_broken true -# end -# end -# empty_field(space) = fill((;), space) - -# # Broadcasting over the wrong size should error -# test_broken_throws(empty_field(TU.PointSpace(FT))) -# test_broken_throws(empty_field(TU.SpectralElementSpace1D(FT))) -# test_broken_throws(empty_field(TU.SpectralElementSpace2D(FT))) -# test_broken_throws(empty_field(TU.ColumnCenterFiniteDifferenceSpace(FT))) -# test_broken_throws(empty_field(TU.ColumnFaceFiniteDifferenceSpace(FT))) -# test_broken_throws(empty_field(TU.SphereSpectralElementSpace(FT))) -# test_broken_throws(empty_field(TU.CenterExtrudedFiniteDifferenceSpace(FT))) -# test_broken_throws(empty_field(TU.FaceExtrudedFiniteDifferenceSpace(FT))) - -# # TODO: performance optimization: shouldn't we do -# # nothing when broadcasting over empty fields? -# # This is otherwise a performance penalty if -# # users regularly rely on empty fields. In particular: -# # - does iterating over empty fields load data? -# # - what is the overhead in iterating over empty fields? -# # - what is the use case of anything useful that can be -# # done by iterating over empty fields? -# end - -# @testset "Broadcasting interception for tuple-valued fields" begin -# n1 = n2 = 1 -# Nij = 4 -# space = spectral_space_2D(n1 = n1, n2 = n2, Nij = Nij) - -# nt_field = Fields.Field( -# IJFH{NamedTuple{(:a, :b), Tuple{Float64, Float64}}, Nij}( -# ones(Nij, Nij, 2, n1 * n2), -# ), -# space, -# ) -# nt_sum = sum(nt_field) -# @test nt_sum isa NamedTuple{(:a, :b), Tuple{Float64, Float64}} -# @test nt_sum.a ≈ 8.0 * 10.0 rtol = 10eps() -# @test nt_sum.b ≈ 8.0 * 10.0 rtol = 10eps() -# @test norm(nt_field) ≈ sqrt(2.0) rtol = 10eps() - -# # test scalar asignment -# nt_field.a .= 0.0 -# @test sum(nt_field.a) == 0.0 -# end - -# @testset "Special case handling for broadcased norm to pass through space local geometry" begin -# space = spectral_space_2D() -# u = Geometry.Covariant12Vector.(ones(space), ones(space)) -# @test norm.(u) ≈ hypot(4 / 8 / 2, 4 / 10 / 2) .* ones(space) -# end - -# @testset "FieldVector" begin -# space = spectral_space_2D() -# u = Geometry.Covariant12Vector.(ones(space), ones(space)) -# x = Fields.coordinate_field(space) -# y = [1.0, 2.0, 3.0] -# z = 1.0 -# Y = Fields.FieldVector(u = u, k = (x = x, y = y, z = z)) - -# @test propertynames(Y) == (:u, :k) -# @test propertynames(Y.k) == (:x, :y, :z) -# @test Y.u === u -# @test Y.k.x === x -# @test Y.k.y === y -# @test Y.k.z === z - -# @test deepcopy(Y).u !== u -# @test deepcopy(Y).k.x !== x -# @test deepcopy(Y).k.y !== y - -# @test getfield(deepcopy(Y).u, :space) === space - -# Y1 = 2 .* Y -# @test parent(Y1.u) == 2 .* parent(u) -# @test parent(Y1.k.x) == 2 .* parent(x) -# @test Y1.k.y == 2 .* y -# @test Y1.k.z === 2 * z - -# Y1 .= Y1 .+ 2 .* Y -# @test parent(Y1.u) == 4 .* parent(u) -# @test parent(Y1.k.x) == 4 .* parent(x) -# @test Y1.k.y == 4 .* y -# @test Y1.k.z === 4 * z - -# Y.k.z = 3.0 -# @test Y.k.z === 3.0 -# end - -# @testset "FieldVector array_type" begin -# device = ClimaComms.device() -# context = ClimaComms.SingletonCommsContext(device) -# space = TU.PointSpace(Float32; context) -# xcenters = Fields.coordinate_field(space).x -# y = Fields.FieldVector(x = xcenters) -# @test ClimaComms.array_type(y) == ClimaComms.array_type(device) -# y = Fields.FieldVector(x = xcenters, y = xcenters) -# @test ClimaComms.array_type(y) == ClimaComms.array_type(device) -# end - -# @testset "FieldVector basetype replacement and deepcopy" begin -# device = ClimaComms.CPUSingleThreaded() # constructing space_vijfh is broken -# context = ClimaComms.SingletonCommsContext(device) -# domain_z = Domains.IntervalDomain( -# Geometry.ZPoint(-1.0) .. Geometry.ZPoint(1.0), -# periodic = true, -# ) -# mesh_z = Meshes.IntervalMesh(domain_z; nelems = 10) -# topology_z = Topologies.IntervalTopology(context, mesh_z) - -# domain_x = Domains.IntervalDomain( -# Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0), -# periodic = true, -# ) -# mesh_x = Meshes.IntervalMesh(domain_x; nelems = 10) -# topology_x = Topologies.IntervalTopology(context, mesh_x) - -# domain_xy = Domains.RectangleDomain( -# Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0), -# Geometry.YPoint(-1.0) .. Geometry.YPoint(1.0), -# x1periodic = true, -# x2periodic = true, -# ) -# mesh_xy = Meshes.RectilinearMesh(domain_xy, 10, 10) -# topology_xy = Topologies.Topology2D(context, mesh_xy) - -# quad = Spaces.Quadratures.GLL{4}() - -# space_vf = Spaces.CenterFiniteDifferenceSpace(topology_z) -# space_ifh = Spaces.SpectralElementSpace1D(topology_x, quad) -# space_ijfh = Spaces.SpectralElementSpace2D(topology_xy, quad) -# space_vifh = Spaces.ExtrudedFiniteDifferenceSpace(space_ifh, space_vf) -# space_vijfh = Spaces.ExtrudedFiniteDifferenceSpace(space_ijfh, space_vf) - -# space2field(space) = map( -# coord -> (coord, Geometry.Covariant12Vector(1.0, 2.0)), -# Fields.coordinate_field(space), -# ) - -# Y = Fields.FieldVector( -# field_vf = space2field(space_vf), -# field_if = slab(space2field(space_ifh), 1), -# field_ifh = space2field(space_ifh), -# field_ijf = slab(space2field(space_ijfh), 1, 1), -# field_ijfh = space2field(space_ijfh), -# field_vifh = space2field(space_vifh), -# field_vijfh = space2field(space_vijfh), -# array = [1.0, 2.0, 3.0], -# scalar = 1.0, -# ) - -# Yf = ForwardDiff.Dual{Nothing}.(Y, 1.0) -# Yf .= Yf .^ 2 .+ Y -# @test all(ForwardDiff.value.(Yf) .== Y .^ 2 .+ Y) -# @test all(ForwardDiff.partials.(Yf, 1) .== 2 .* Y) - -# dual_field = Yf.field_vf -# dual_field_original_basetype = similar(Y.field_vf, eltype(dual_field)) -# @test eltype(dual_field_original_basetype) === eltype(dual_field) -# @test eltype(parent(dual_field_original_basetype)) === Float64 -# @test eltype(parent(dual_field)) === ForwardDiff.Dual{Nothing, Float64, 1} - -# object_that_contains_Yf = (; Yf) -# @test axes(deepcopy(Yf).field_vf) === space_vf -# @test axes(deepcopy(object_that_contains_Yf).Yf.field_vf) === space_vf -# end - -# @testset "Scalar field iterator" begin -# space = spectral_space_2D() -# u = Geometry.Covariant12Vector.(ones(space), ones(space)) -# x = Fields.coordinate_field(space) -# y = [1.0, 2.0, 3.0] -# z = 1.0 -# Y = Fields.FieldVector(u = u, k = (x = x, y = y, z = z)) - -# prop_chains = Fields.property_chains(Y) -# @test length(prop_chains) == 6 -# @test prop_chains[1] == (:u, :components, :data, 1) -# @test prop_chains[2] == (:u, :components, :data, 2) -# @test prop_chains[3] == (:k, :x, :x) -# @test prop_chains[4] == (:k, :x, :y) -# @test prop_chains[5] == (:k, :y) -# @test prop_chains[6] == (:k, :z) - -# FT = Float64 -# nt = -# (; x = FT(0), y = FT(0), tup = ntuple(i -> (; a = FT(1), b = FT(1)), 2)) -# Y = fill(nt, space) - -# prop_chains = Fields.property_chains(Y) -# @test prop_chains[1] == (:x,) -# @test prop_chains[2] == (:y,) -# @test prop_chains[3] == (:tup, 1, :a) -# @test prop_chains[4] == (:tup, 1, :b) -# @test prop_chains[5] == (:tup, 2, :a) -# @test prop_chains[6] == (:tup, 2, :b) - -# @test Fields.single_field(Y, prop_chains[1]) === Y.x -# @test Fields.single_field(Y, prop_chains[2]) === Y.y -# @test Fields.single_field(Y, prop_chains[3]) === getproperty(Y.tup, 1).a -# @test Fields.single_field(Y, prop_chains[4]) === getproperty(Y.tup, 1).b -# @test Fields.single_field(Y, prop_chains[5]) === getproperty(Y.tup, 2).a -# @test Fields.single_field(Y, prop_chains[6]) === getproperty(Y.tup, 2).b - -# for (i, (var, prop_chain)) in enumerate(Fields.field_iterator(Y)) -# @test prop_chains[i] == prop_chain -# @test var === Fields.single_field(Y, prop_chain) -# end -# end - -# # Test truncated field type printing: -# ClimaCore.Fields.truncate_printing_field_types() = true -# @testset "Truncated printing" begin -# nt = (; x = Float64(0), y = Float64(0)) -# Y = fill(nt, spectral_space_2D()) -# @test sprint(show, typeof(Y); context = IOContext(stdout)) == -# "Field{(:x, :y)} (trunc disp)" -# end -# ClimaCore.Fields.truncate_printing_field_types() = false - -# @testset "Standard printing" begin -# nt = (; x = Float64(0), y = Float64(0)) -# Y = fill(nt, spectral_space_2D()) -# s = sprint(show, typeof(Y)) # just make sure this doesn't break -# end - -# @testset "Set!" begin -# space = spectral_space_2D() -# FT = Float64 -# nt = (; x = FT(0), y = FT(0)) -# Y = fill(nt, space) -# foo(local_geom) = -# sin(local_geom.coordinates.x * local_geom.coordinates.y) + 3 -# Fields.set!(foo, Y.x) -# @test all((parent(Y.x) .> 1)) -# end - -# @testset "PointField" begin -# device = ClimaComms.CPUSingleThreaded() # a bunch of cuda pieces are broken -# context = ClimaComms.SingletonCommsContext(device) -# FT = Float64 -# coord = Geometry.XPoint(FT(π)) -# space = Spaces.PointSpace(context, coord) -# @test parent(Spaces.local_geometry_data(space)) == -# FT[Geometry.component(coord, 1), 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] -# field = Fields.coordinate_field(space) -# @test field isa Fields.PointField -# @test Fields.field_values(field)[] == coord - -# if ClimaComms.device(context) isa ClimaComms.AbstractCPUDevice -# @test sum(field.x) == FT(π) -# elseif ClimaComms.device(context) isa ClimaComms.CUDADevice -# # Not yet supported -# # @test sum(field.x) == FT(π) -# end - -# field = ones(space) .* π -# sin_field = sin.(field) -# add_field = field .+ field -# @test isapprox(Fields.field_values(sin_field)[], FT(0.0); atol = √eps(FT)) -# @test isapprox(Fields.field_values(add_field)[], FT(2π)) -# end - -# @testset "Level" begin -# FT = Float64 -# for space in TU.all_spaces(FT) -# TU.levelable(space) || continue -# Y = fill((; x = FT(2)), space) -# lg_space = Spaces.level(space, TU.fc_index(1, space)) -# lg_field_space = axes(Fields.level(Y, TU.fc_index(1, space))) -# @test all( -# lg_space.local_geometry.coordinates === -# lg_field_space.local_geometry.coordinates, -# ) -# @test all(Fields.zeros(lg_space) == Fields.zeros(lg_field_space)) -# end -# end - -# @testset "Points from Columns" begin -# FT = Float64 -# for space in TU.all_spaces(FT) -# if space isa Spaces.SpectralElementSpace1D -# Y = fill((; x = FT(1)), space) -# point_space_from_field = axes(Fields.column(Y.x, 1, 1)) -# point_space = Spaces.column(space, 1, 1) -# @test Fields.ones(point_space) == -# Fields.ones(point_space_from_field) -# end -# if space isa Spaces.SpectralElementSpace2D -# Y = fill((; x = FT(1)), space) -# point_space_from_field = axes(Fields.column(Y.x, 1, 1, 1)) -# point_space = Spaces.column(space, 1, 1, 1) -# @test Fields.ones(point_space) == -# Fields.ones(point_space_from_field) -# end - -# end -# end - -# @testset "(Domain/Column)-surface broadcasting" begin -# FT = Float64 -# function domain_surface_bc!(x, ᶜz_surf, ᶜx_surf) -# @. x = x + ᶜz_surf -# # exercises broadcast_shape(PointSpace, PointSpace) -# @. x = x + (ᶜz_surf * ᶜx_surf) -# nothing -# end -# function column_surface_bc!(x, ᶜz_surf, ᶜx_surf) -# Fields.bycolumn(axes(x)) do colidx -# @. x[colidx] = x[colidx] + ᶜz_surf[colidx] -# # exercises broadcast_shape(PointSpace, PointSpace) -# @. x[colidx] = x[colidx] + (ᶜz_surf[colidx] * ᶜx_surf[colidx]) -# end -# nothing -# end -# for space in TU.all_spaces(FT) -# # Filter out spaces without z coordinates: -# TU.has_z_coordinates(space) || continue -# Y = fill((; x = FT(1)), space) -# ᶜz_surf = -# Spaces.level(Fields.coordinate_field(Y).z, TU.fc_index(1, space)) -# ᶜx_surf = copy(Spaces.level(Y.x, TU.fc_index(1, space))) - -# # Still need to define broadcast rules for surface planes with 3D domains -# domain_surface_bc!(Y.x, ᶜz_surf, ᶜx_surf) - -# # Skip spaces incompatible with Fields.bycolumn: -# TU.bycolumnable(space) || continue -# Yc = fill((; x = FT(1)), space) -# column_surface_bc!(Yc.x, ᶜz_surf, ᶜx_surf) -# @test Y.x == Yc.x -# nothing -# end -# nothing -# end - -# @testset "Broadcasting same spaces different instances" begin -# space1 = spectral_space_2D() -# space2 = spectral_space_2D() -# field1 = ones(space1) -# field2 = 2 .* ones(space2) -# @test Fields.is_diagonalized_spaces(typeof(space1), typeof(space2)) -# @test_throws ErrorException( -# "Broacasted spaces are the same ClimaCore.Spaces type but not the same instance", -# ) field1 .= field2 - -# # turn warning on -# Fields.allow_mismatched_diagonalized_spaces() = true -# @test_warn "Broacasted spaces are the same ClimaCore.Spaces type but not the same instance" field1 .= -# field2 -# @test parent(field1) == parent(field2) -# Fields.allow_mismatched_diagonalized_spaces() = false -# end - -# struct InferenceFoo{FT} -# bar::FT -# end -# Base.broadcastable(x::InferenceFoo) = Ref(x) -# @testset "Inference failure message" begin -# function ics_foo(::Type{FT}, lg, foo) where {FT} -# uv = Geometry.UVVector(FT(0), FT(0)) -# z = Geometry.Covariant12Vector(uv, lg) -# y = foo.bingo -# return (; x = FT(0) + y) -# end -# function ics_foo_with_field(::Type{FT}, lg, foo, f) where {FT} -# uv = Geometry.UVVector(FT(0), FT(0)) -# z = Geometry.Covariant12Vector(uv, lg) -# ζ = f.a -# y = foo.baz -# return (; x = FT(0) + y - ζ) -# end -# function FieldFromNamedTupleBroken( -# space, -# ics::Function, -# ::Type{FT}, -# params..., -# ) where {FT} -# lg = Fields.local_geometry_field(space) -# return ics.(FT, lg, params...) -# end -# FT = Float64 -# foo = InferenceFoo(2.0) -# device = ClimaComms.CPUSingleThreaded() # cuda fill is broken -# context = ClimaComms.SingletonCommsContext(device) -# for space in TU.all_spaces(FT; context) -# Y = fill((; a = FT(0), b = FT(1)), space) -# @test_throws ErrorException("type InferenceFoo has no field bingo") FieldFromNamedTupleBroken( -# space, -# ics_foo, -# FT, -# foo, -# ) -# @test_throws ErrorException("type InferenceFoo has no field baz") FieldFromNamedTupleBroken( -# space, -# ics_foo_with_field, -# FT, -# foo, -# Y, -# ) -# end - -# end - -# @testset "Δz_field" begin -# FT = Float64 -# context = ClimaComms.SingletonCommsContext() -# x = FT(1) -# y = FT(2) -# z = FT(3) -# lat, long = FT(4), FT(5) -# x1 = FT(1) -# x2 = FT(2) -# x3 = FT(3) -# coords = [ -# Geometry.ZPoint(z), -# Geometry.XZPoint(x, z), -# Geometry.XYZPoint(x, y, z), -# Geometry.LatLongZPoint(lat, long, z), -# Geometry.Cartesian3Point(x3), -# Geometry.Cartesian13Point(x1, x3), -# Geometry.Cartesian123Point(x1, x2, x3), -# ] -# all_components = [ -# SMatrix{1, 1, FT}(range(1, 1)...), -# SMatrix{2, 2, FT}(range(1, 4)...), -# SMatrix{3, 3, FT}(range(1, 9)...), -# SMatrix{3, 3, FT}(range(1, 9)...), -# SMatrix{1, 1, FT}(range(1, 1)...), -# SMatrix{2, 2, FT}(range(1, 4)...), -# SMatrix{3, 3, FT}(range(1, 9)...), -# ] - -# expected_dzs = [1.0, 4.0, 9.0, 9.0, 1.0, 4.0, 9.0] - -# for (components, coord, expected_dz) in -# zip(all_components, coords, expected_dzs) -# CoordType = typeof(coord) -# AIdx = Geometry.coordinate_axis(CoordType) -# at = Geometry.AxisTensor( -# (Geometry.LocalAxis{AIdx}(), Geometry.CovariantAxis{AIdx}()), -# components, -# ) -# local_geometry = Geometry.LocalGeometry(coord, FT(1.0), FT(1.0), at) -# space = Spaces.PointSpace(context, local_geometry) -# dz_computed = Array(parent(Fields.Δz_field(space))) -# @test length(dz_computed) == 1 -# @test dz_computed[1] == expected_dz -# end -# end - -# @testset "scalar assignment" begin -# device = ClimaComms.CPUSingleThreaded() # constructing space_vijfh is broken -# context = ClimaComms.SingletonCommsContext(device) -# domain_z = Domains.IntervalDomain( -# Geometry.ZPoint(-1.0) .. Geometry.ZPoint(1.0), -# periodic = true, -# ) -# mesh_z = Meshes.IntervalMesh(domain_z; nelems = 10) -# topology_z = Topologies.IntervalTopology(context, mesh_z) - -# domain_x = Domains.IntervalDomain( -# Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0), -# periodic = true, -# ) -# mesh_x = Meshes.IntervalMesh(domain_x; nelems = 10) -# topology_x = Topologies.IntervalTopology(context, mesh_x) - -# domain_xy = Domains.RectangleDomain( -# Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0), -# Geometry.YPoint(-1.0) .. Geometry.YPoint(1.0), -# x1periodic = true, -# x2periodic = true, -# ) -# mesh_xy = Meshes.RectilinearMesh(domain_xy, 10, 10) -# topology_xy = Topologies.Topology2D(context, mesh_xy) - -# quad = Spaces.Quadratures.GLL{4}() - -# space_vf = Spaces.CenterFiniteDifferenceSpace(topology_z) -# space_ifh = Spaces.SpectralElementSpace1D(topology_x, quad) -# space_ijfh = Spaces.SpectralElementSpace2D(topology_xy, quad) -# space_vifh = Spaces.ExtrudedFiniteDifferenceSpace(space_ifh, space_vf) -# space_vijfh = Spaces.ExtrudedFiniteDifferenceSpace(space_ijfh, space_vf) - -# C = map(x -> Geometry.Covariant12Vector(1.0, 1.0), zeros(space_vifh)) -# @test all(==(1.0), parent(C)) -# C .= Ref(zero(eltype(C))) -# @test all(==(0.0), parent(C)) -# end +function spectral_space_2D(; n1 = 1, n2 = 1, Nij = 4) + domain = Domains.RectangleDomain( + Geometry.XPoint(-3.0) .. Geometry.XPoint(5.0), + Geometry.YPoint(-2.0) .. Geometry.YPoint(8.0), + x1periodic = false, + x2periodic = false, + x1boundary = (:east, :west), + x2boundary = (:south, :north), + ) + mesh = Meshes.RectilinearMesh(domain, n1, n2) + device = ClimaComms.CPUSingleThreaded() + grid_topology = + Topologies.Topology2D(ClimaComms.SingletonCommsContext(device), mesh) + + quad = Spaces.Quadratures.GLL{Nij}() + space = Spaces.SpectralElementSpace2D(grid_topology, quad) + return space +end + +@testset "1×1 2D domain space" begin + Nij = 4 + n1 = n2 = 1 + space = spectral_space_2D(n1 = n1, n2 = n2, Nij = Nij) + + field = + Fields.Field(IJFH{ComplexF64, Nij}(ones(Nij, Nij, 2, n1 * n2)), space) + + @test sum(field) ≈ Complex(1.0, 1.0) * 8.0 * 10.0 rtol = 10eps() + @test sum(x -> 3.0, field) ≈ 3 * 8.0 * 10.0 rtol = 10eps() + @test mean(field) ≈ Complex(1.0, 1.0) rtol = 10eps() + @test mean(x -> 3.0, field) ≈ 3 rtol = 10eps() + @test norm(field) ≈ sqrt(2.0) rtol = 10eps() + @test norm(field, 1) ≈ norm(Complex(1.0, 1.0)) rtol = 10eps() + @test norm(field, Inf) ≈ norm(Complex(1.0, 1.0)) rtol = 10eps() + @test norm(field; normalize = false) ≈ sqrt(2.0 * 8.0 * 10.0) rtol = 10eps() + @test norm(field, 1; normalize = false) ≈ + norm(Complex(1.0, 1.0)) * 8.0 * 10.0 rtol = 10eps() + @test norm(field, Inf; normalize = false) ≈ norm(Complex(1.0, 1.0)) rtol = + 10eps() + + @test extrema(real, field) == (1.0, 1.0) + + @test Operators.matrix_interpolate(field, 4) ≈ + [Complex(1.0, 1.0) for i in 1:(4 * n1), j in 1:(4 * n2)] + + + field_sin = map(x -> sin((x.x) / 2), Fields.coordinate_field(space)) + M = Operators.matrix_interpolate(field_sin, 20) + @test size(M) == (20, 20) # 20 x 20 for a 1 element field + + real_field = field.re + + # test broadcasting + res = field .+ 1 + @test parent(Fields.field_values(res)) == Float64[ + f == 1 ? 2 : 1 for i in 1:Nij, j in 1:Nij, f in 1:2, h in 1:(n1 * n2) + ] + + res = field.re .+ 1 + @test parent(Fields.field_values(res)) == + Float64[2 for i in 1:Nij, j in 1:Nij, f in 1:1, h in 1:(n1 * n2)] + + # test field slab broadcasting + f1 = ones(space) + f2 = ones(space) + + for h in 1:(n1 * n2) + f1_slab = Fields.slab(f1, h) + f2_slab = Fields.slab(f2, h) + q = f1_slab .+ f2_slab + f1_slab .= q .+ f2_slab + end + @test all(parent(f1) .== 3) + + point_field = Fields.column(field, 1, 1, 1) + @test axes(point_field) isa Spaces.PointSpace +end + +# https://github.com/CliMA/ClimaCore.jl/issues/1126 +function pow_n(f) + @. f.x = f.x^2 + return nothing +end +@testset "Broadcasting with ^n" begin + FT = Float32 + device = ClimaComms.CPUSingleThreaded() # fill is broken on gpu + context = ClimaComms.SingletonCommsContext(device) + for space in TU.all_spaces(FT; context) + f = fill((; x = FT(1)), space) + pow_n(f) # Compile first + p_allocated = @allocated pow_n(f) + if space isa Spaces.SpectralElementSpace1D + @test p_allocated == 0 + else + @test p_allocated == 0 broken = (device isa ClimaComms.CUDADevice) + end + end +end + +function ifelse_broadcast_allocating(a, b, c) + FT = eltype(a) + @. a = ifelse(true || c < b * FT(1), FT(0), c) + return nothing +end + +function ifelse_broadcast_or(a, b, c) + FT = eltype(a) + val = FT(1) + @. a = ifelse(true || c < b * val, FT(0), c) + return nothing +end + +function ifelse_broadcast_simple(a, b, c) + FT = eltype(a) + @. a = ifelse(c < b * FT(1), FT(0), c) + return nothing +end + +@testset "Broadcasting ifelse" begin + FT = Float32 + device = ClimaComms.CPUSingleThreaded() # broken on gpu + context = ClimaComms.SingletonCommsContext(device) + for space in ( + TU.CenterExtrudedFiniteDifferenceSpace(FT; context), + TU.ColumnCenterFiniteDifferenceSpace(FT; context), + ) + a = Fields.level(fill(FT(0), space), 1) + b = Fields.level(fill(FT(2), space), 1) + c = Fields.level(fill(FT(3), space), 1) + + ifelse_broadcast_allocating(a, b, c) + p_allocated = @allocated ifelse_broadcast_allocating(a, b, c) + @test_broken p_allocated == 0 + + ifelse_broadcast_or(a, b, c) + p_allocated = @allocated ifelse_broadcast_or(a, b, c) + @test p_allocated == 0 + + ifelse_broadcast_simple(a, b, c) + p_allocated = @allocated ifelse_broadcast_simple(a, b, c) + @test p_allocated == 0 + end +end + +# Requires `--check-bounds=yes` +@testset "Constructing & broadcasting over empty fields" begin + FT = Float32 + for space in TU.all_spaces(FT) + f = fill((;), space) + @. f += f + end + + function test_broken_throws(f) + try + @. f += 1 + # we want to throw exception, test is broken + @test_broken false + catch + # we want to throw exception, unexpected pass + @test_broken true + end + end + empty_field(space) = fill((;), space) + + # Broadcasting over the wrong size should error + test_broken_throws(empty_field(TU.PointSpace(FT))) + test_broken_throws(empty_field(TU.SpectralElementSpace1D(FT))) + test_broken_throws(empty_field(TU.SpectralElementSpace2D(FT))) + test_broken_throws(empty_field(TU.ColumnCenterFiniteDifferenceSpace(FT))) + test_broken_throws(empty_field(TU.ColumnFaceFiniteDifferenceSpace(FT))) + test_broken_throws(empty_field(TU.SphereSpectralElementSpace(FT))) + test_broken_throws(empty_field(TU.CenterExtrudedFiniteDifferenceSpace(FT))) + test_broken_throws(empty_field(TU.FaceExtrudedFiniteDifferenceSpace(FT))) + + # TODO: performance optimization: shouldn't we do + # nothing when broadcasting over empty fields? + # This is otherwise a performance penalty if + # users regularly rely on empty fields. In particular: + # - does iterating over empty fields load data? + # - what is the overhead in iterating over empty fields? + # - what is the use case of anything useful that can be + # done by iterating over empty fields? +end + +@testset "Broadcasting interception for tuple-valued fields" begin + n1 = n2 = 1 + Nij = 4 + space = spectral_space_2D(n1 = n1, n2 = n2, Nij = Nij) + + nt_field = Fields.Field( + IJFH{NamedTuple{(:a, :b), Tuple{Float64, Float64}}, Nij}( + ones(Nij, Nij, 2, n1 * n2), + ), + space, + ) + nt_sum = sum(nt_field) + @test nt_sum isa NamedTuple{(:a, :b), Tuple{Float64, Float64}} + @test nt_sum.a ≈ 8.0 * 10.0 rtol = 10eps() + @test nt_sum.b ≈ 8.0 * 10.0 rtol = 10eps() + @test norm(nt_field) ≈ sqrt(2.0) rtol = 10eps() + + # test scalar asignment + nt_field.a .= 0.0 + @test sum(nt_field.a) == 0.0 +end + +@testset "Special case handling for broadcased norm to pass through space local geometry" begin + space = spectral_space_2D() + u = Geometry.Covariant12Vector.(ones(space), ones(space)) + @test norm.(u) ≈ hypot(4 / 8 / 2, 4 / 10 / 2) .* ones(space) +end + +@testset "FieldVector" begin + space = spectral_space_2D() + u = Geometry.Covariant12Vector.(ones(space), ones(space)) + x = Fields.coordinate_field(space) + y = [1.0, 2.0, 3.0] + z = 1.0 + Y = Fields.FieldVector(u = u, k = (x = x, y = y, z = z)) + + @test propertynames(Y) == (:u, :k) + @test propertynames(Y.k) == (:x, :y, :z) + @test Y.u === u + @test Y.k.x === x + @test Y.k.y === y + @test Y.k.z === z + + @test deepcopy(Y).u !== u + @test deepcopy(Y).k.x !== x + @test deepcopy(Y).k.y !== y + + @test getfield(deepcopy(Y).u, :space) === space + + Y1 = 2 .* Y + @test parent(Y1.u) == 2 .* parent(u) + @test parent(Y1.k.x) == 2 .* parent(x) + @test Y1.k.y == 2 .* y + @test Y1.k.z === 2 * z + + Y1 .= Y1 .+ 2 .* Y + @test parent(Y1.u) == 4 .* parent(u) + @test parent(Y1.k.x) == 4 .* parent(x) + @test Y1.k.y == 4 .* y + @test Y1.k.z === 4 * z + + Y.k.z = 3.0 + @test Y.k.z === 3.0 +end + +@testset "FieldVector array_type" begin + device = ClimaComms.device() + context = ClimaComms.SingletonCommsContext(device) + space = TU.PointSpace(Float32; context) + xcenters = Fields.coordinate_field(space).x + y = Fields.FieldVector(x = xcenters) + @test ClimaComms.array_type(y) == ClimaComms.array_type(device) + y = Fields.FieldVector(x = xcenters, y = xcenters) + @test ClimaComms.array_type(y) == ClimaComms.array_type(device) +end + +@testset "FieldVector basetype replacement and deepcopy" begin + device = ClimaComms.CPUSingleThreaded() # constructing space_vijfh is broken + context = ClimaComms.SingletonCommsContext(device) + domain_z = Domains.IntervalDomain( + Geometry.ZPoint(-1.0) .. Geometry.ZPoint(1.0), + periodic = true, + ) + mesh_z = Meshes.IntervalMesh(domain_z; nelems = 10) + topology_z = Topologies.IntervalTopology(context, mesh_z) + + domain_x = Domains.IntervalDomain( + Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0), + periodic = true, + ) + mesh_x = Meshes.IntervalMesh(domain_x; nelems = 10) + topology_x = Topologies.IntervalTopology(context, mesh_x) + + domain_xy = Domains.RectangleDomain( + Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0), + Geometry.YPoint(-1.0) .. Geometry.YPoint(1.0), + x1periodic = true, + x2periodic = true, + ) + mesh_xy = Meshes.RectilinearMesh(domain_xy, 10, 10) + topology_xy = Topologies.Topology2D(context, mesh_xy) + + quad = Spaces.Quadratures.GLL{4}() + + space_vf = Spaces.CenterFiniteDifferenceSpace(topology_z) + space_ifh = Spaces.SpectralElementSpace1D(topology_x, quad) + space_ijfh = Spaces.SpectralElementSpace2D(topology_xy, quad) + space_vifh = Spaces.ExtrudedFiniteDifferenceSpace(space_ifh, space_vf) + space_vijfh = Spaces.ExtrudedFiniteDifferenceSpace(space_ijfh, space_vf) + + space2field(space) = map( + coord -> (coord, Geometry.Covariant12Vector(1.0, 2.0)), + Fields.coordinate_field(space), + ) + + Y = Fields.FieldVector( + field_vf = space2field(space_vf), + field_if = slab(space2field(space_ifh), 1), + field_ifh = space2field(space_ifh), + field_ijf = slab(space2field(space_ijfh), 1, 1), + field_ijfh = space2field(space_ijfh), + field_vifh = space2field(space_vifh), + field_vijfh = space2field(space_vijfh), + array = [1.0, 2.0, 3.0], + scalar = 1.0, + ) + + Yf = ForwardDiff.Dual{Nothing}.(Y, 1.0) + Yf .= Yf .^ 2 .+ Y + @test all(ForwardDiff.value.(Yf) .== Y .^ 2 .+ Y) + @test all(ForwardDiff.partials.(Yf, 1) .== 2 .* Y) + + dual_field = Yf.field_vf + dual_field_original_basetype = similar(Y.field_vf, eltype(dual_field)) + @test eltype(dual_field_original_basetype) === eltype(dual_field) + @test eltype(parent(dual_field_original_basetype)) === Float64 + @test eltype(parent(dual_field)) === ForwardDiff.Dual{Nothing, Float64, 1} + + object_that_contains_Yf = (; Yf) + @test axes(deepcopy(Yf).field_vf) === space_vf + @test axes(deepcopy(object_that_contains_Yf).Yf.field_vf) === space_vf +end + +@testset "Scalar field iterator" begin + space = spectral_space_2D() + u = Geometry.Covariant12Vector.(ones(space), ones(space)) + x = Fields.coordinate_field(space) + y = [1.0, 2.0, 3.0] + z = 1.0 + Y = Fields.FieldVector(u = u, k = (x = x, y = y, z = z)) + + prop_chains = Fields.property_chains(Y) + @test length(prop_chains) == 6 + @test prop_chains[1] == (:u, :components, :data, 1) + @test prop_chains[2] == (:u, :components, :data, 2) + @test prop_chains[3] == (:k, :x, :x) + @test prop_chains[4] == (:k, :x, :y) + @test prop_chains[5] == (:k, :y) + @test prop_chains[6] == (:k, :z) + + FT = Float64 + nt = + (; x = FT(0), y = FT(0), tup = ntuple(i -> (; a = FT(1), b = FT(1)), 2)) + Y = fill(nt, space) + + prop_chains = Fields.property_chains(Y) + @test prop_chains[1] == (:x,) + @test prop_chains[2] == (:y,) + @test prop_chains[3] == (:tup, 1, :a) + @test prop_chains[4] == (:tup, 1, :b) + @test prop_chains[5] == (:tup, 2, :a) + @test prop_chains[6] == (:tup, 2, :b) + + @test Fields.single_field(Y, prop_chains[1]) === Y.x + @test Fields.single_field(Y, prop_chains[2]) === Y.y + @test Fields.single_field(Y, prop_chains[3]) === getproperty(Y.tup, 1).a + @test Fields.single_field(Y, prop_chains[4]) === getproperty(Y.tup, 1).b + @test Fields.single_field(Y, prop_chains[5]) === getproperty(Y.tup, 2).a + @test Fields.single_field(Y, prop_chains[6]) === getproperty(Y.tup, 2).b + + for (i, (var, prop_chain)) in enumerate(Fields.field_iterator(Y)) + @test prop_chains[i] == prop_chain + @test var === Fields.single_field(Y, prop_chain) + end +end + +# Test truncated field type printing: +ClimaCore.Fields.truncate_printing_field_types() = true +@testset "Truncated printing" begin + nt = (; x = Float64(0), y = Float64(0)) + Y = fill(nt, spectral_space_2D()) + @test sprint(show, typeof(Y); context = IOContext(stdout)) == + "Field{(:x, :y)} (trunc disp)" +end +ClimaCore.Fields.truncate_printing_field_types() = false + +@testset "Standard printing" begin + nt = (; x = Float64(0), y = Float64(0)) + Y = fill(nt, spectral_space_2D()) + s = sprint(show, typeof(Y)) # just make sure this doesn't break +end + +@testset "Set!" begin + space = spectral_space_2D() + FT = Float64 + nt = (; x = FT(0), y = FT(0)) + Y = fill(nt, space) + foo(local_geom) = + sin(local_geom.coordinates.x * local_geom.coordinates.y) + 3 + Fields.set!(foo, Y.x) + @test all((parent(Y.x) .> 1)) +end + +@testset "PointField" begin + device = ClimaComms.CPUSingleThreaded() # a bunch of cuda pieces are broken + context = ClimaComms.SingletonCommsContext(device) + FT = Float64 + coord = Geometry.XPoint(FT(π)) + space = Spaces.PointSpace(context, coord) + @test parent(Spaces.local_geometry_data(space)) == + FT[Geometry.component(coord, 1), 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] + field = Fields.coordinate_field(space) + @test field isa Fields.PointField + @test Fields.field_values(field)[] == coord + + if ClimaComms.device(context) isa ClimaComms.AbstractCPUDevice + @test sum(field.x) == FT(π) + elseif ClimaComms.device(context) isa ClimaComms.CUDADevice + # Not yet supported + # @test sum(field.x) == FT(π) + end + + field = ones(space) .* π + sin_field = sin.(field) + add_field = field .+ field + @test isapprox(Fields.field_values(sin_field)[], FT(0.0); atol = √eps(FT)) + @test isapprox(Fields.field_values(add_field)[], FT(2π)) +end + +@testset "Level" begin + FT = Float64 + for space in TU.all_spaces(FT) + TU.levelable(space) || continue + Y = fill((; x = FT(2)), space) + lg_space = Spaces.level(space, TU.fc_index(1, space)) + lg_field_space = axes(Fields.level(Y, TU.fc_index(1, space))) + @test all( + lg_space.local_geometry.coordinates === + lg_field_space.local_geometry.coordinates, + ) + @test all(Fields.zeros(lg_space) == Fields.zeros(lg_field_space)) + end +end + +@testset "Points from Columns" begin + FT = Float64 + for space in TU.all_spaces(FT) + if space isa Spaces.SpectralElementSpace1D + Y = fill((; x = FT(1)), space) + point_space_from_field = axes(Fields.column(Y.x, 1, 1)) + point_space = Spaces.column(space, 1, 1) + @test Fields.ones(point_space) == + Fields.ones(point_space_from_field) + end + if space isa Spaces.SpectralElementSpace2D + Y = fill((; x = FT(1)), space) + point_space_from_field = axes(Fields.column(Y.x, 1, 1, 1)) + point_space = Spaces.column(space, 1, 1, 1) + @test Fields.ones(point_space) == + Fields.ones(point_space_from_field) + end + + end +end + +@testset "(Domain/Column)-surface broadcasting" begin + FT = Float64 + function domain_surface_bc!(x, ᶜz_surf, ᶜx_surf) + @. x = x + ᶜz_surf + # exercises broadcast_shape(PointSpace, PointSpace) + @. x = x + (ᶜz_surf * ᶜx_surf) + nothing + end + function column_surface_bc!(x, ᶜz_surf, ᶜx_surf) + Fields.bycolumn(axes(x)) do colidx + @. x[colidx] = x[colidx] + ᶜz_surf[colidx] + # exercises broadcast_shape(PointSpace, PointSpace) + @. x[colidx] = x[colidx] + (ᶜz_surf[colidx] * ᶜx_surf[colidx]) + end + nothing + end + for space in TU.all_spaces(FT) + # Filter out spaces without z coordinates: + TU.has_z_coordinates(space) || continue + Y = fill((; x = FT(1)), space) + ᶜz_surf = + Spaces.level(Fields.coordinate_field(Y).z, TU.fc_index(1, space)) + ᶜx_surf = copy(Spaces.level(Y.x, TU.fc_index(1, space))) + + # Still need to define broadcast rules for surface planes with 3D domains + domain_surface_bc!(Y.x, ᶜz_surf, ᶜx_surf) + + # Skip spaces incompatible with Fields.bycolumn: + TU.bycolumnable(space) || continue + Yc = fill((; x = FT(1)), space) + column_surface_bc!(Yc.x, ᶜz_surf, ᶜx_surf) + @test Y.x == Yc.x + nothing + end + nothing +end + +@testset "Broadcasting same spaces different instances" begin + space1 = spectral_space_2D() + space2 = spectral_space_2D() + field1 = ones(space1) + field2 = 2 .* ones(space2) + @test Fields.is_diagonalized_spaces(typeof(space1), typeof(space2)) + @test_throws ErrorException( + "Broacasted spaces are the same ClimaCore.Spaces type but not the same instance", + ) field1 .= field2 + + # turn warning on + Fields.allow_mismatched_diagonalized_spaces() = true + @test_warn "Broacasted spaces are the same ClimaCore.Spaces type but not the same instance" field1 .= + field2 + @test parent(field1) == parent(field2) + Fields.allow_mismatched_diagonalized_spaces() = false +end + +struct InferenceFoo{FT} + bar::FT +end +Base.broadcastable(x::InferenceFoo) = Ref(x) +@testset "Inference failure message" begin + function ics_foo(::Type{FT}, lg, foo) where {FT} + uv = Geometry.UVVector(FT(0), FT(0)) + z = Geometry.Covariant12Vector(uv, lg) + y = foo.bingo + return (; x = FT(0) + y) + end + function ics_foo_with_field(::Type{FT}, lg, foo, f) where {FT} + uv = Geometry.UVVector(FT(0), FT(0)) + z = Geometry.Covariant12Vector(uv, lg) + ζ = f.a + y = foo.baz + return (; x = FT(0) + y - ζ) + end + function FieldFromNamedTupleBroken( + space, + ics::Function, + ::Type{FT}, + params..., + ) where {FT} + lg = Fields.local_geometry_field(space) + return ics.(FT, lg, params...) + end + FT = Float64 + foo = InferenceFoo(2.0) + device = ClimaComms.CPUSingleThreaded() # cuda fill is broken + context = ClimaComms.SingletonCommsContext(device) + for space in TU.all_spaces(FT; context) + Y = fill((; a = FT(0), b = FT(1)), space) + @test_throws ErrorException("type InferenceFoo has no field bingo") FieldFromNamedTupleBroken( + space, + ics_foo, + FT, + foo, + ) + @test_throws ErrorException("type InferenceFoo has no field baz") FieldFromNamedTupleBroken( + space, + ics_foo_with_field, + FT, + foo, + Y, + ) + end + +end + +@testset "Δz_field" begin + FT = Float64 + context = ClimaComms.SingletonCommsContext() + x = FT(1) + y = FT(2) + z = FT(3) + lat, long = FT(4), FT(5) + x1 = FT(1) + x2 = FT(2) + x3 = FT(3) + coords = [ + Geometry.ZPoint(z), + Geometry.XZPoint(x, z), + Geometry.XYZPoint(x, y, z), + Geometry.LatLongZPoint(lat, long, z), + Geometry.Cartesian3Point(x3), + Geometry.Cartesian13Point(x1, x3), + Geometry.Cartesian123Point(x1, x2, x3), + ] + all_components = [ + SMatrix{1, 1, FT}(range(1, 1)...), + SMatrix{2, 2, FT}(range(1, 4)...), + SMatrix{3, 3, FT}(range(1, 9)...), + SMatrix{3, 3, FT}(range(1, 9)...), + SMatrix{1, 1, FT}(range(1, 1)...), + SMatrix{2, 2, FT}(range(1, 4)...), + SMatrix{3, 3, FT}(range(1, 9)...), + ] + + expected_dzs = [1.0, 4.0, 9.0, 9.0, 1.0, 4.0, 9.0] + + for (components, coord, expected_dz) in + zip(all_components, coords, expected_dzs) + CoordType = typeof(coord) + AIdx = Geometry.coordinate_axis(CoordType) + at = Geometry.AxisTensor( + (Geometry.LocalAxis{AIdx}(), Geometry.CovariantAxis{AIdx}()), + components, + ) + local_geometry = Geometry.LocalGeometry(coord, FT(1.0), FT(1.0), at) + space = Spaces.PointSpace(context, local_geometry) + dz_computed = Array(parent(Fields.Δz_field(space))) + @test length(dz_computed) == 1 + @test dz_computed[1] == expected_dz + end +end + +@testset "scalar assignment" begin + device = ClimaComms.CPUSingleThreaded() # constructing space_vijfh is broken + context = ClimaComms.SingletonCommsContext(device) + domain_z = Domains.IntervalDomain( + Geometry.ZPoint(-1.0) .. Geometry.ZPoint(1.0), + periodic = true, + ) + mesh_z = Meshes.IntervalMesh(domain_z; nelems = 10) + topology_z = Topologies.IntervalTopology(context, mesh_z) + + domain_x = Domains.IntervalDomain( + Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0), + periodic = true, + ) + mesh_x = Meshes.IntervalMesh(domain_x; nelems = 10) + topology_x = Topologies.IntervalTopology(context, mesh_x) + + domain_xy = Domains.RectangleDomain( + Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0), + Geometry.YPoint(-1.0) .. Geometry.YPoint(1.0), + x1periodic = true, + x2periodic = true, + ) + mesh_xy = Meshes.RectilinearMesh(domain_xy, 10, 10) + topology_xy = Topologies.Topology2D(context, mesh_xy) + + quad = Spaces.Quadratures.GLL{4}() + + space_vf = Spaces.CenterFiniteDifferenceSpace(topology_z) + space_ifh = Spaces.SpectralElementSpace1D(topology_x, quad) + space_ijfh = Spaces.SpectralElementSpace2D(topology_xy, quad) + space_vifh = Spaces.ExtrudedFiniteDifferenceSpace(space_ifh, space_vf) + space_vijfh = Spaces.ExtrudedFiniteDifferenceSpace(space_ijfh, space_vf) + + C = map(x -> Geometry.Covariant12Vector(1.0, 1.0), zeros(space_vifh)) + @test all(==(1.0), parent(C)) + C .= Ref(zero(eltype(C))) + @test all(==(0.0), parent(C)) +end function integrate_bycolumn!(∫y, Y) Fields.bycolumn(axes(Y.y)) do colidx Operators.column_integral_definite!(∫y[colidx], Y.y[colidx]) @@ -698,7 +698,7 @@ convergence_rate(err, Δh) = col_copy = similar(y[Fields.ColumnIndex((1, 1), 1)]) return Fields.Field(Fields.field_values(col_copy), axes(col_copy)) end - device = ClimaComms.device() + device = ClimaComms.CPUSingleThreaded() context = ClimaComms.SingletonCommsContext(device) for zelem in (2^2, 2^3, 2^4, 2^5) for space in TU.all_spaces(FT; zelem, context) @@ -706,7 +706,6 @@ convergence_rate(err, Δh) = TU.has_z_coordinates(space) || continue # Skip spaces incompatible with Fields.bycolumn: TU.bycolumnable(space) || continue - @show nameof(typeof(space)) Y = fill((; y = FT(1)), space) zcf = Fields.coordinate_field(Y.y).z From d8aa6bc45948e8816995f10f1c43d8b634c4075a Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Wed, 19 Jul 2023 15:03:04 -0700 Subject: [PATCH 06/11] Test on gpu --- test/Fields/field.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/Fields/field.jl b/test/Fields/field.jl index cb711f0306..b2501279e4 100644 --- a/test/Fields/field.jl +++ b/test/Fields/field.jl @@ -698,7 +698,7 @@ convergence_rate(err, Δh) = col_copy = similar(y[Fields.ColumnIndex((1, 1), 1)]) return Fields.Field(Fields.field_values(col_copy), axes(col_copy)) end - device = ClimaComms.CPUSingleThreaded() + device = ClimaComms.device() context = ClimaComms.SingletonCommsContext(device) for zelem in (2^2, 2^3, 2^4, 2^5) for space in TU.all_spaces(FT; zelem, context) From 1255b9586328ac3a1d51ecd7500d1292d76dae8b Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Wed, 19 Jul 2023 16:00:27 -0700 Subject: [PATCH 07/11] Dont use PointSpace in array type test --- test/Fields/field.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/Fields/field.jl b/test/Fields/field.jl index b2501279e4..3fdb57d322 100644 --- a/test/Fields/field.jl +++ b/test/Fields/field.jl @@ -269,7 +269,7 @@ end @testset "FieldVector array_type" begin device = ClimaComms.device() context = ClimaComms.SingletonCommsContext(device) - space = TU.PointSpace(Float32; context) + space = TU.SpectralElementSpace1D(Float32; context) xcenters = Fields.coordinate_field(space).x y = Fields.FieldVector(x = xcenters) @test ClimaComms.array_type(y) == ClimaComms.array_type(device) From 32e372ff106115cf2b729b433510619b40c8b87a Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Thu, 20 Jul 2023 09:01:34 -0700 Subject: [PATCH 08/11] define device --- src/Spaces/pointspace.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/Spaces/pointspace.jl b/src/Spaces/pointspace.jl index 82b476c53d..4aee1a8ff4 100644 --- a/src/Spaces/pointspace.jl +++ b/src/Spaces/pointspace.jl @@ -11,7 +11,9 @@ struct PointSpace{LG <: DataLayouts.Data0D} <: AbstractPointSpace local_geometry::LG end -ClimaComms.device(space::PointSpace) = ClimaComms.device(space.context) +# not strictly correct, but it is needed for mapreduce +ClimaComms.device(space::PointSpace) = + ClimaComms.device(ClimaComms.context(space)) ClimaComms.context(space::PointSpace) = ClimaComms.SingletonCommsContext(ClimaComms.CPUSingleThreaded()) From 6af87665eda1ca511fabc283e138fe1efe256606 Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Thu, 20 Jul 2023 10:45:33 -0700 Subject: [PATCH 09/11] allowscalar --- test/Fields/field.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/test/Fields/field.jl b/test/Fields/field.jl index 3fdb57d322..06e8b62c82 100644 --- a/test/Fields/field.jl +++ b/test/Fields/field.jl @@ -11,6 +11,7 @@ import ClimaCore: using LinearAlgebra: norm using Statistics: mean using ForwardDiff +using CUDA include( joinpath(pkgdir(ClimaCore), "test", "TestUtilities", "TestUtilities.jl"), @@ -35,7 +36,7 @@ function spectral_space_2D(; n1 = 1, n2 = 1, Nij = 4) space = Spaces.SpectralElementSpace2D(grid_topology, quad) return space end - +#= @testset "1×1 2D domain space" begin Nij = 4 n1 = n2 = 1 @@ -667,13 +668,13 @@ end C .= Ref(zero(eltype(C))) @test all(==(0.0), parent(C)) end +=# function integrate_bycolumn!(∫y, Y) Fields.bycolumn(axes(Y.y)) do colidx Operators.column_integral_definite!(∫y[colidx], Y.y[colidx]) nothing end end - """ convergence_rate(err, Δh) @@ -711,7 +712,9 @@ convergence_rate(err, Δh) = zcf = Fields.coordinate_field(Y.y).z Δz = Fields.Δz_field(axes(zcf)) Δz_col = Δz[Fields.ColumnIndex((1, 1), 1)] - Δz_1 = parent(Δz_col)[1] + Δz_1 = CUDA.allowscalar() do + parent(Δz_col)[1] + end key = (space, zelem) if !haskey(results, key) results[key] = From dff77b0727f7ebfc5f04a0c03b7d0d50538bff44 Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Thu, 20 Jul 2023 15:31:08 -0700 Subject: [PATCH 10/11] support SubArray on fill --- src/DataLayouts/cuda.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/DataLayouts/cuda.jl b/src/DataLayouts/cuda.jl index 7f67b5de66..57e3dd4dd5 100644 --- a/src/DataLayouts/cuda.jl +++ b/src/DataLayouts/cuda.jl @@ -97,7 +97,11 @@ end function Base.fill!( dest::IJFH{S, Nij, A}, val, -) where {S, Nij, A <: CUDA.CuArray} +) where { + S, + Nij, + A <: Union{CUDA.CuArray, SubArray{<:Any, <:Any, <:CUDA.CuArray}}, +} _, _, _, _, Nh = size(dest) if Nh > 0 CUDA.@cuda threads = (Nij, Nij) blocks = (Nh, 1) knl_fill!(dest, val) From bd87cbc5b484bcbf2be888dc14c51c734f6b86b4 Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Thu, 20 Jul 2023 21:16:00 -0700 Subject: [PATCH 11/11] disable nsys --- .buildkite/pipeline.yml | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 85a14e83fa..fc76fa3b31 100755 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -939,7 +939,8 @@ steps: - label: ":flower_playing_cards: Rising Bubble 3D hybrid invariant (ρe)" key: "gpu_rising_bubble_3d_hybrid_invariant_rhoe" command: - - "nsys profile --trace=nvtx,cuda --output=output/$$BUILDKITE_STEP_KEY/report julia --color=yes --project=examples examples/hybrid/box/bubble_3d_invariant_rhoe.jl" +# - "nsys profile --trace=nvtx,cuda --output=output/$$BUILDKITE_STEP_KEY/report julia --color=yes --project=examples examples/hybrid/box/bubble_3d_invariant_rhoe.jl" + - "julia --color=yes --project=examples examples/hybrid/box/bubble_3d_invariant_rhoe.jl" artifact_paths: - "examples/hybrid/box/output/gpu_bubble_3d_invariant_rhoe/*_low_*" agents: @@ -948,7 +949,8 @@ steps: - label: ":flower_playing_cards: Rising Bubble 3D hybrid invariant (ρe), custom resolution" key: "gpu_rising_bubble_3d_hybrid_invariant_rhoe_custom" command: - - "nsys profile --trace=nvtx,cuda --output=output/$$BUILDKITE_STEP_KEY/report julia --color=yes --project=examples examples/hybrid/box/bubble_3d_invariant_rhoe.jl Float64 custom 1000 1000 4 16 3 0.05 700.0" +# - "nsys profile --trace=nvtx,cuda --output=output/$$BUILDKITE_STEP_KEY/report julia --color=yes --project=examples examples/hybrid/box/bubble_3d_invariant_rhoe.jl Float64 custom 1000 1000 4 16 3 0.05 700.0" + - "julia --color=yes --project=examples examples/hybrid/box/bubble_3d_invariant_rhoe.jl Float64 custom 1000 1000 4 16 3 0.05 700.0" artifact_paths: - "examples/hybrid/box/output/gpu_bubble_3d_invariant_rhoe/*_custom_*" agents: @@ -1097,7 +1099,8 @@ steps: - label: ":computer: MPI Shallow-water 2D sphere barotropic instability alpha0" key: "cpu_mpi_shallowwater_2d_cg_sphere_barotropic_alpha0" command: - - "nsys profile --trace=nvtx,mpi --mpi-impl=openmpi --output=examples/sphere/output/cg_sphere_shallowwater_barotropic_instability_alpha0/report.%q{NPROCS} mpiexec julia --color=yes --project=examples examples/sphere/shallow_water.jl barotropic_instability" +# - "nsys profile --trace=nvtx,mpi --mpi-impl=openmpi --output=examples/sphere/output/cg_sphere_shallowwater_barotropic_instability_alpha0/report.%q{NPROCS} mpiexec julia --color=yes --project=examples examples/sphere/shallow_water.jl barotropic_instability" + - "mpiexec julia --color=yes --project=examples examples/sphere/shallow_water.jl barotropic_instability" artifact_paths: - "examples/sphere/output/cg_sphere_shallowwater_barotropic_instability_alpha0/*" env: @@ -1132,7 +1135,8 @@ steps: key: "cuda_shallowwater_2d_cg_sphere" command: - mkdir -p output/$$BUILDKITE_STEP_KEY - - nsys profile --trace=nvtx,cuda --output=output/$$BUILDKITE_STEP_KEY/report julia --color=yes --project=examples examples/sphere/shallow_water_cuda.jl +# - nsys profile --trace=nvtx,cuda --output=output/$$BUILDKITE_STEP_KEY/report julia --color=yes --project=examples examples/sphere/shallow_water_cuda.jl + - julia --color=yes --project=examples examples/sphere/shallow_water_cuda.jl artifact_paths: - output/cuda_shallowwater_2d_cg_sphere agents: