Support DataLayouts.VF in InputOutput.read_field()

[valentinaschueller]

The read_field() function does not support reading in VF fields at the moment. This is the type of data layout I get when storing values from 1D FiniteDifferenceSpace data, so I would like to be able to load it again.

The function should be adapted to support this type of DataLayout. Things that need to be fixed for sure to make this work:

• add a VF case to _scan_data_layout()
  
• extend h_dim() for VFSingleton (needed here
• after that, this line fails for me with the following error, this needs to be fixed:

ERROR: MethodError: no method matching getindex(::Nothing, ::Int64)
The function `getindex` exists, but no method is defined for this combination of argument types.

Here I do not know what should be fixed.

For the h_dim() function, I was not sure what value should be put. 1? 0?

[Sbozzolo]

Here's a broken test for future use:

using Test
import ClimaCore

using ClimaComms
const comms_ctx = ClimaComms.context(ClimaComms.CPUSingleThreaded())
pid, nprocs = ClimaComms.init(comms_ctx)
filename = ClimaComms.bcast(comms_ctx, tempname(pwd()))
if ClimaComms.iamroot(comms_ctx)
    @info "Comms context" comms_ctx nprocs filename
end

@testset "HDF5 restart test for 1d finite difference space" begin
    FT = Float32

    z_min = FT(0)
    z_max = FT(30e3)
    z_elem = 10
    center_staggering = ClimaCore.Grids.CellCenter()
    face_staggering = ClimaCore.Grids.CellFace()

    center_space = ClimaCore.CommonSpaces.ColumnSpace(;
                                                      z_min,
                                                      z_max,
                                                      z_elem,
                                                      staggering = center_staggering)

    face_space = ClimaCore.CommonSpaces.ColumnSpace(;
                                                      z_min,
                                                      z_max,
                                                      z_elem,
                                                      staggering = face_staggering)

    center_field = ClimaCore.Fields.coordinate_field(center_space)
    face_field = ClimaCore.Fields.coordinate_field(face_space)

    Y = ClimaCore.Fields.FieldVector(; c = center_field, f = face_field)

    # write field vector to hdf5 file
    writer = ClimaCore.InputOutput.HDF5Writer(filename, comms_ctx)
    ClimaCore.InputOutput.write!(writer, Y, "Y")
    close(writer)

    reader = ClimaCore.InputOutput.HDF5Reader(filename, comms_ctx)
    restart_Y = ClimaCore.InputOutput.read_field(reader, "Y") # read fieldvector from hdf5 file
    close(reader)
    @test restart_Y == Y # test if restart is exact
end

I sort of got it to work with, but a little bit more work is needed

        # h_dim = DataLayouts.h_dim(DataLayouts.singleton(DataLayout))
        if topology isa Topologies.Topology2D
            nd = ndims(obj)
            localidx =
                ntuple(d -> d == h_dim ? topology.local_elem_gidx : (:), nd)
            data = ArrayType(obj[localidx...])
        else
            data = ArrayType(read(obj))
        end
        # Nij = size(data, findfirst("I", data_layout)[1])
        # 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"]))
        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)
        else
            # values = DataLayout{ElType, Nij, Nht...}(data) # when Nh is in type-domain
            values = DataLayout{ElType, Nij}(data)
        end