Add support for HDF5 read and write with PointSpaces (#2162)

* Add support for InputOuput with PointSpaces

`write!` on Fields is now dispatched off the input
field's space. When it is a point space, the local
geometry data is written. When read_field is called,
and no grid is present, then the field is assumed
to be on a point space.

Add function comments

Simplify unit_point read/write test

Make suggested changes to unit_point and docstrings

clean up point space test

* Add `read_type` to read field eltypes

When the reader recreates a field from a HDF5, it called
eval on a parsed input string. This replaces that with a function
that ensures there are no function calls inside the parsed string
before eval'ing it. This is based on old HDF5.jl code.

make suggested read_type changes

Revert field_eltype to value_type

* Change value_type attribute to field_eltype

Change the value_type attribute, used by the HDF5 reader
and writer, to field_eltype. Backwards support is maintained
and tested.

Remove HDF5 import from test

import climacomms backend

Author: imreddyTeja

.buildkite/pipeline.yml

@@ -1149,6 +1149,15 @@ steps:
   - group: "Unit: InputOutput"
     steps:
 
+      - label: "Unit: read_type"
+        key: unit_read_type
+        command: "julia --color=yes --check-bounds=yes --project=.buildkite test/InputOutput/unit_read_type.jl"
+        env:
+          CLIMACOMMS_DEVICE: "CUDA"
+        agents:
+          slurm_gpus: 1
+          slurm_mem: 3GB
+
       - label: "Unit: spectralelement2d"
         key: unit_spectralelement2d
         command: "julia --color=yes --check-bounds=yes --project=.buildkite test/InputOutput/unit_spectralelement2d.jl"

NEWS.md

@@ -4,6 +4,13 @@ ClimaCore.jl Release Notes
 main
 -------
 
+- Added support for InputOutput with PointSpaces
+  PR [2162](https://github.com/CliMA/ClimaCore.jl/pull/2162).
+
+- Reading a `Field` with `HDF5Reader` will not call `eval ∘ Meta.parse` on the raw
+`value_type` attribute. Instead, the new `read_type` function is used, which prevents
+the execution of arbitrary code. PR [2162](https://github.com/CliMA/ClimaCore.jl/pull/2162).
+
  - Improved `show` for spaces, and added `show` for grids. PR [2202](https://github.com/CliMA/ClimaCore.jl/pull/2202).
 
 v0.14.26

src/InputOutput/readers.jl

@@ -41,6 +41,42 @@ function get_key(file, key)
     end
 end
 
+"""
+    read_type(ts::AbstractString)
+
+Parse a string `ts` into an expression, and then evaluate it into a type if it is a valid type expression.
+"""
+function read_type(ts::AbstractString)
+    type_expr = Meta.parse(ts)
+    if is_type_expr(type_expr)
+        return eval(type_expr)
+    end
+    error("$ts cannot be parsed into a valid type")
+end
+
+"""
+    is_type_expr(expr)
+
+Check if an expression is a type expression, with no function calls or other non-type
+expressions, with the expection of @NamedTuple.
+This function is based on the JLD.jl `is_valid_type_exp`.
+See https://github.com/JuliaIO/JLD.jl/blob/80ac89643e3ad87545e48f4d361a00a29cdf4e2f/src/JLD.jl#L922
+"""
+is_type_expr(s::Symbol) = true
+is_type_expr(q::QuoteNode) = is_type_expr(q.value)
+function is_type_expr(expr::Expr)
+    if expr.head == :macrocall && expr.args[1] == Symbol("@NamedTuple")
+        # skip the LineNumberNode, as eval does nothing for it
+        if length(expr.args) > 2
+            return all(map(is_type_expr, expr.args[3:end]))
+        end
+        return true
+    end
+    return expr.head in (:curly, :., :tuple, :braces, Symbol("::")) &&
+           all(map(is_type_expr, expr.args))
+end
+is_type_expr(t) = isbits(t)
+
 """
     HDF5Reader(filename::AbstractString[, context::ClimaComms.AbstractCommsContext])
     HDF5Reader(::Function, filename::AbstractString[, context::ClimaComms.AbstractCommsContext])
@@ -182,6 +218,7 @@ function _scan_data_layout(layoutstring::AbstractString)
         "VIJHF",
         "VIFH",
         "VIHF",
+        "DataF",
     ) "datalayout is $layoutstring"
     layoutstring == "IJFH" && return DataLayouts.IJFH
     layoutstring == "IJHF" && return DataLayouts.IJHF
@@ -192,6 +229,7 @@ function _scan_data_layout(layoutstring::AbstractString)
     layoutstring == "VF" && return DataLayouts.VF
     layoutstring == "VIJFH" && return DataLayouts.VIJFH
     layoutstring == "VIJHF" && return DataLayouts.VIJHF
+    layoutstring == "DataF" && return DataLayouts.DataF
     return DataLayouts.VIFH
 end
 
@@ -513,11 +551,23 @@ function read_field(reader::HDF5Reader, name::AbstractString)
                 staggering = Grids.CellFace()
             end
             space = Spaces.space(grid, staggering)
-        else
+        elseif haskey(attrs(obj), "space")
             space = read_space(reader, attrs(obj)["space"])
+        else
+            # if the there is no grid, then the field is on a PointSpace
+            lg_obj = reader.file["local_geometry_data/$name"]
+            ArrayType = ClimaComms.array_type(ClimaComms.device(reader.context))
+            lg_data = ArrayType(read(lg_obj))
+            # because it is a point space, the data layout of local_geometry_data is always DataF
+            lg_type = read_type(attrs(lg_obj)["local_geometry_type"])
+            local_geometry_data = DataLayouts.DataF{lg_type}(lg_data)
+            space = Spaces.PointSpace(local_geometry_data)
+            topology = nothing
+        end
+        if !(space isa Spaces.AbstractPointSpace)
+            topology = Spaces.topology(space)
+            ArrayType = ClimaComms.array_type(topology)
         end
-        topology = Spaces.topology(space)
-        ArrayType = ClimaComms.array_type(topology)
         data_layout = attrs(obj)["data_layout"]
         has_horizontal = occursin('I', data_layout)
         DataLayout = _scan_data_layout(data_layout)
@@ -535,14 +585,20 @@ function read_field(reader::HDF5Reader, name::AbstractString)
         # For when `Nh` is added back to the type space
         #     Nhd = Nh_dim(data_layout)
         #     Nht = Nhd == -1 ? () : (size(data, Nhd),)
-        ElType = eval(Meta.parse(attrs(obj)["value_type"]))
+        # The `value_type` attribute is deprecated. here we mantain backwards compatibility
+        ElType = read_type(
+            haskey(attrs(obj), "field_eltype") ?
+            attrs(obj)["field_eltype"] : attrs(obj)["value_type"],
+        )
         if data_layout in ("VIJFH", "VIFH")
             Nv = size(data, 1)
             # values = DataLayout{ElType, Nv, Nij, Nht...}(data) # when Nh is in type-domain
             values = DataLayout{ElType, Nv, Nij}(data)
         elseif data_layout in ("VF",)
             Nv = size(data, 1)
             values = DataLayout{ElType, Nv}(data)
+        elseif data_layout in ("DataF",)
+            values = DataLayout{ElType}(data)
         else
             # values = DataLayout{ElType, Nij, Nht...}(data) # when Nh is in type-domain
             values = DataLayout{ElType, Nij}(data)

src/InputOutput/writers.jl

@@ -451,16 +451,91 @@ function write_new!(
 end
 
 # write fields
+"""
+    write!(writer::HDF5Writer, field::Fields.Field, name::AbstractString)
+
+Write the `field` to the HDF5 in `writer` and assign it the given `name`.
+"""
 function write!(writer::HDF5Writer, field::Fields.Field, name::AbstractString)
+    write!(writer, field, name, axes(field))
+end
+
+"""
+    write!(
+        writer::HDF5Writer,
+        field::Fields.Field,
+        name::AbstractString,
+        space::Spaces.AbstractPointSpace,
+    )
+
+Write a `Field`, with `axes` of type `PointSpace`,  to the HDF5 file. The field
+is written to the `fields` group in the file, with the name `name`. The local
+geometry data of the `PointSpace` is written to the `local_geometry_data` group
+with name `name`.
+"""
+function write!(
+    writer::HDF5Writer,
+    field::Fields.Field,
+    name::AbstractString,
+    space::Spaces.AbstractPointSpace,
+)
+    array = parent(field)
+    lg_data = Grids.local_geometry_data(space)
+    lg_type = Grids.local_geometry_type(typeof(space))
+    lg_array = parent(lg_data)
+    dataset = create_dataset(
+        writer.file,
+        "fields/$name",
+        datatype(eltype(array)),
+        dataspace(size(array)),
+    )
+    dataset[:] = array
+    write_attribute(dataset, "type", "Field")
+    write_attribute(
+        dataset,
+        "data_layout",
+        string(nameof(typeof(Fields.field_values(field)))),
+    )
+    write_attribute(dataset, "field_eltype", string(eltype(field)))
+    local_geometry_dataset = create_dataset(
+        writer.file,
+        "local_geometry_data/$name",
+        datatype(eltype(array)),
+        dataspace(size(lg_array)),
+    )
+    local_geometry_dataset[:] = lg_array
+    write_attribute(
+        local_geometry_dataset,
+        "local_geometry_type",
+        string(lg_type),
+    )
+end
+
+"""
+    write!(
+        writer::HDF5Writer,
+        field::Fields.Field,
+        name::AbstractString,
+        space::Spaces.AbstractSpace,
+    )
+
+Write an object of type 'Field' and name 'name' to the HDF5 file.
+"""
+function write!(
+    writer::HDF5Writer,
+    field::Fields.Field,
+    name::AbstractString,
+    space::Spaces.AbstractSpace,
+)
     values = Fields.field_values(field)
-    space = axes(field)
+    array = parent(field)
+    nd = ndims(array)
+
     staggering = Spaces.staggering(space)
+    topology = Spaces.topology(space)
     grid = Spaces.grid(space)
     grid_name = write!(writer, grid)
 
-    array = parent(field)
-    topology = Spaces.topology(space)
-    nd = ndims(array)
     if topology isa Topologies.Topology2D &&
        !(writer.context isa ClimaComms.SingletonCommsContext)
         nelems = Topologies.nelems(topology)
@@ -491,7 +566,7 @@ function write!(writer::HDF5Writer, field::Fields.Field, name::AbstractString)
         "data_layout",
         string(nameof(typeof(Fields.field_values(field)))),
     )
-    write_attribute(dataset, "value_type", string(eltype(field)))
+    write_attribute(dataset, "field_eltype", string(eltype(field)))
     write_attribute(dataset, "grid", grid_name)
     if !isnothing(staggering)
         write_attribute(

test/InputOutput/unit_point.jl

@@ -0,0 +1,34 @@
+using Test
+import ClimaCore: Spaces, Fields, InputOutput, Geometry
+using ClimaComms
+const comms_ctx = ClimaComms.context()
+
+
+@testset "HDF5 write/read test for 0d PointSpace" begin
+
+    FT = Float32
+
+    space = Spaces.PointSpace(comms_ctx, Geometry.ZPoint(FT(1)))
+    field_0d = Fields.local_geometry_field(space)
+    Y = Fields.FieldVector(; p = field_0d)
+
+    filename = tempname()
+
+    InputOutput.HDF5Writer(filename, comms_ctx) do writer
+        InputOutput.write!(writer, "Y" => Y) # write field vector from hdf5 file
+    end
+
+    InputOutput.HDF5Reader(filename, comms_ctx) do reader
+        restart_Y = InputOutput.read_field(reader, "Y") # read fieldvector from hdf5 file
+        @test restart_Y == Y # test if restart is exact
+        # test if space is the same by comparing local geometry data
+        # note that the read spaces are not "===" to the written spaces because there is no grid to cache
+        lg_point = Spaces.local_geometry_data(axes(Y.p))[]
+        lg_restart = Spaces.local_geometry_data(axes(restart_Y.p))[]
+        @test typeof(lg_point) == typeof(lg_restart)
+        @test all(
+            getproperty(lg_point, p) == getproperty(lg_restart, p) for
+            p in propertynames(lg_point)
+        )
+    end
+end

test/InputOutput/unit_read_type.jl

@@ -0,0 +1,48 @@
+using Test
+import ClimaCore
+import ClimaCore: Fields, InputOutput
+using ClimaComms
+ClimaComms.@import_required_backends
+include(
+    joinpath(pkgdir(ClimaCore), "test", "TestUtilities", "TestUtilities.jl"),
+)
+import .TestUtilities as TU
+
+compare_read_type(x) = InputOutput.read_type(string(eltype(x))) == eltype(x)
+@testset "Read field element types" begin
+    context = ClimaComms.context()
+    FT = Float64
+
+    for space in TU.all_spaces(FT; context)
+        lg_field = Fields.local_geometry_field(space)
+        @test compare_read_type(lg_field)
+        coord_field = Fields.coordinate_field(space)
+        @test compare_read_type(coord_field)
+    end
+
+    space = TU.ColumnCenterFiniteDifferenceSpace(FT; context)
+    @test compare_read_type(fill((1.0, 2.0, (3.0, 4.0, (5.0,))), space))
+    @test compare_read_type(fill((; a = FT(0), b = (; c = FT(1))), space))
+    # test that attempting to read a type from a string with an expression that executes code throws an error
+    @test_throws ErrorException InputOutput.read_type("1+2")
+    @test_throws ErrorException InputOutput.read_type("Base.rm(\"foo.jl\")")
+end
+
+@testset "backwards compatibility with `value_type` attribute" begin
+    context = ClimaComms.context()
+    FT = Float64
+
+    space = TU.ColumnCenterFiniteDifferenceSpace(FT; context)
+    Y = fill((1.0, 2.0, (3.0, 4.0, (5.0,))), space)
+    filename = tempname()
+    InputOutput.HDF5Writer(filename, context) do writer
+        InputOutput.write!(writer, "Y" => Y) # write field vector from hdf5 file
+        obj = writer.file["fields/Y"]
+        InputOutput.HDF5.delete_attribute(obj, "field_eltype") # remove the field_eltype attribute
+        InputOutput.HDF5.write_attribute(obj, "value_type", string(eltype(Y))) # add the deprecated value_type attribute
+    end
+    InputOutput.HDF5Reader(filename, context) do reader
+        restart_Y = InputOutput.read_field(reader, "Y")
+        @test eltype(restart_Y) == eltype(Y)
+    end
+end

test/runtests.jl

@@ -96,6 +96,7 @@ UnitTest("Hypsography - 3d sphere"                 ,"Hypsography/3dsphere.jl"),
 UnitTest("Remapping"                               ,"Operators/remapping.jl"),
 UnitTest("Limiter"                                 ,"Limiters/limiter.jl"),
 UnitTest("InputOutput - hdf5"                      ,"InputOutput/unit_hdf5.jl"),
+UnitTest("InputOutput - parse_type"                ,"InputOutput/unit_read_type.jl"),
 UnitTest("InputOutput - spectralelement2d"         ,"InputOutput/unit_spectralelement2d.jl"),
 UnitTest("InputOutput - hybrid2dbox"               ,"InputOutput/unit_hybrid2dbox.jl"),
 UnitTest("InputOutput - hybrid2dbox_topography"    ,"InputOutput/unit_hybrid2dbox_topography.jl"),
@@ -104,6 +105,7 @@ UnitTest("InputOutput - hybrid3dbox"               ,"InputOutput/unit_hybrid3dbo
 UnitTest("InputOutput - hybrid3dcubedsphere"       ,"InputOutput/unit_hybrid3dcubedsphere.jl"),
 UnitTest("InputOutput - hybrid3dcubedsphere_topo"  ,"InputOutput/unit_hybrid3dcubedsphere_topography.jl"),
 UnitTest("InputOutput - finitedifferences"         ,"InputOutput/unit_finitedifference.jl"),
+UnitTest("InputOutput - pointspaces"               ,"InputOutput/unit_point.jl"),
 UnitTest("Array interpolation"                     ,"Remapping/interpolate_array.jl"),
 UnitTest("Array interpolation"                     ,"Remapping/distributed_remapping.jl"),
 UnitTest("Aqua"                                    ,"aqua.jl"),