gcmarray internals revamp, skip allocs by default in read and similar

Project.toml

@@ -1,6 +1,6 @@
 name = "MeshArrays"
 uuid = "cb8c808f-1acf-59a3-9d2b-6e38d009f683"
-version = "0.5.9"
+version = "0.5.10"
 authors = ["gaelforget <gforget@mit.edu>"]
 
 [deps]

src/Operations.jl

@@ -489,7 +489,7 @@ Returns a `gridpath` by default, or a `NamedTuple` with fields `lat`, `name`,
 function LatitudeCircle(lat,Γ::NamedTuple;
   format=:gridpath, range=(0.0,360.0))
 
-  mskCint = similar(Γ.YC)
+  mskCint = similar(Γ.YC,allocate=true)
   for f in 1:mskCint.grid.nFaces
       yf = Γ.YC.f[f]; mf = mskCint.f[f]
       @inbounds for j in axes(mf,2), i in axes(mf,1)
@@ -558,9 +558,11 @@ Compute edge mask (mskC,mskW,mskS) from domain interior mask (mskCint).
 This is used in `LatitudeCircles` and `Transect`.
 """
 function edge_mask(mskCint::AbstractMeshArray)
-  mskC = similar(mskCint)
-  mskW = similar(mskCint)
-  mskS = similar(mskCint)
+  mskC = similar(mskCint,allocate=true)
+  mskW = similar(mskCint,allocate=true)
+  mskS = similar(mskCint,allocate=true)
+
+  #need to allocate here
 
   exFLD = exchange(mskCint).MA   # still needed for halo
 

src/ReadWrite.jl

@@ -26,15 +26,15 @@ function _load_binary(fil::String, y::AbstractMeshArray)
 end
 
 """
-    read(fil::String,x::AbstractMeshArray)
+    read(fil::String, x::AbstractMeshArray)
 
-Read array from file and return as a MeshArray.
+Read binary file into a new MeshArray with the same structure as `x`.
 
-_The second argument (MeshArray or gcmgrid) provides the grid specifications (x.grid.ioSize)._
-```
+Creates a new MeshArray (allocating all face arrays) and populates it from `fil`.
+The template `x` provides grid specifications via `x.grid.ioSize`.
 """
 function read(fil::String, x::AbstractMeshArray)
-  y = similar(x; m=x.meta)
+  y = similar(x;m=x.meta,allocate=true)
   read!(fil, y)
   return y
 end
@@ -85,12 +85,14 @@ function read(xx::Array,γ::gcmgrid; verbose=false)
 end
 
 """
-    read(xx::Array,x::AbstractMeshArray)
+    read(xx::Array, x::AbstractMeshArray)
+
+Reformat Array data into a new MeshArray with the same structure as `x`.
 
-Reformat Array data into a MeshArray similar to `x`. 
+Creates a new MeshArray (allocating all face arrays) and populates it from the Array `xx`.
 """
 function read(xx::Array,x::AbstractMeshArray)
-  y=similar(x; m=x.meta)
+  y=similar(x;m=x.meta,allocate=true)
   read!(xx,y)
   return y
 end
@@ -115,6 +117,7 @@ function read!(xx::Array,x::AbstractMeshArray)
         elseif n3>1
           x[f,i3].=tmp[f]
         else
+          !isempty(x[f]) ? nothing : (x[f]=zeros(facesSize[f][1],facesSize[f][2]))
           x[f].=tmp[f]
         end
       end
@@ -251,7 +254,7 @@ function read_tiles(xx::Array,γ::gcmgrid)
 end
 
 function read_tiles(xx::Array,x::AbstractMeshArray)
-	tmp=similar(x)
+	tmp=similar(x,allocate=true)
 	s=x.grid.ioSize
 	(n1,n2)=x.grid.fSize[1]
 	ni=Int(s[1]/n1)

src/VerticalDimension.jl

@@ -7,6 +7,7 @@ function isosurface(θ,T,Γ)
     d=NaN*similar(θ[:,1])
     nr=size(θ,2)
     for j=1:size(d,1)
+        !isempty(d[j]) ? nothing : (d[j]=NaN*zeros(d.fSize[j][1],d.fSize[j][2]))
         for k=1:nr-1
             i=findall(isnan.(d[j]).&(θ[j,k].>T).&(θ[j,k+1].<=T))
             a=(θ[j,k][i] .- T)./(θ[j,k][i] .- θ[j,k+1][i])

src/grids/tiles.jl

@@ -52,7 +52,7 @@ read!(f,d)
 τ=Tiles(γ,30,30)
 td=Tiles(τ,d)
 
-D=similar(d)
+D=similar(d,allocate=true)
 Tiles!(τ,td,D)
 
 isa(td[1],Array)

src/types/gcmarray.jl

@@ -153,22 +153,40 @@ end
 Base.dataids(A::gcmarray) = (Base.dataids(A.f)..., Base.dataids(A.fSize)..., Base.dataids(A.fIndex)...)
 
 # +
-function Base.getindex(A::AbstractMeshArray{T, N}, I::Vararg{Union{Int,Array{Int},AbstractUnitRange,Colon}, N}) where {T,N}
-  J=1:length(A.fIndex)
-  !isa(I[1],Colon) ? J=J[I[1]] : nothing
-  nFaces=length(J)
-
-  tmpf=A.f[I...]
-  if isa(tmpf,Array{eltype(A),2})
-    tmp=tmpf
-  else
-    n3=Int(length(tmpf)/nFaces)
-    K=(A.grid,eltype(A),A.fSize[J],A.fIndex[J])
-    n3>1 ? tmp=gcmarray(K...,n3) : tmp=gcmarray(K...)
-    for I in eachindex(tmpf); tmp.f[I] = view(tmpf[I],:,:); end
-  end
 
-  return tmp
+function Base.getindex(A::AbstractMeshArray{T, N}, I::Vararg{Union{Int,Array{Int},AbstractUnitRange,Colon}, N}) where {T,N}
+    J = 1:length(A.fIndex)
+    !isa(I[1], Colon) ? J = J[I[1]] : nothing
+    nFaces = length(J)
+
+    tmpf = A.f[I...]
+
+    if isa(tmpf, Array{eltype(A), 2})
+        return tmpf  # Scalar index → single face
+    else
+        fSize_sub = A.fSize[J]
+        fIndex_sub = A.fIndex[J]
+        n3 = Int(length(tmpf) / nFaces)
+
+        if n3 == 1
+            # 1D result
+            f_new = OuterArray{InnerArray{T,2},1}(undef, nFaces)
+            for I_iter in eachindex(tmpf)
+                f_new[I_iter] = view(tmpf[I_iter], :, :)
+            end
+            B = gcmarray{T, 1, InnerArray{T,2}}(
+                A.grid, A.meta, f_new, fSize_sub, fIndex_sub, thisversion)
+        else
+            # 2D result
+            f_new = OuterArray{InnerArray{T,2}, 2}(undef, nFaces, n3)
+            for I_iter in eachindex(tmpf)
+                f_new[I_iter] = view(tmpf[I_iter], :, :)
+            end
+            B = gcmarray{T, 2, InnerArray{T,2}}(
+                A.grid, A.meta, f_new, fSize_sub, fIndex_sub, thisversion)
+        end
+        return B
+    end
 end
 
 """
@@ -189,18 +207,34 @@ function Base.setindex!(A::AbstractMeshArray{T, N}, v, I::Vararg{Int, N}) where
 end
 
 function Base.view(A::AbstractMeshArray{T, N}, I::Vararg{Union{Int,AbstractUnitRange,Colon}, N}) where {T,N}
-  J=1:length(A.fIndex)
-  !isa(I[1],Colon) ? J=J[I[1]] : nothing
-  nFaces=length(J)
-
-  tmpf=view(A.f,I...)
-  n3=Int(length(tmpf)/nFaces) #length(tmpf)>nFaces ? n3=Int(length(tmpf)/nFaces) : n3=1
-
-  K=(A.grid,eltype(A),A.fSize[J],A.fIndex[J])
-  n3>1 ? tmp=gcmarray(K...,n3) : tmp=gcmarray(K...)
-  for I in eachindex(tmpf); tmp.f[I] = view(tmpf[I],:,:); end
-
-  return tmp
+    J = 1:length(A.fIndex)
+    !isa(I[1], Colon) ? J = J[I[1]] : nothing
+    nFaces = length(J)
+
+    tmpf = view(A.f, I...)
+    n3 = Int(length(tmpf) / nFaces)
+
+    fSize_sub = A.fSize[J]
+    fIndex_sub = A.fIndex[J]
+
+    if n3 == 1
+        # 1D result
+        f_new = OuterArray{InnerArray{T,2},1}(undef, nFaces)
+        for I_iter in eachindex(tmpf)
+            f_new[I_iter] = view(tmpf[I_iter], :, :)
+        end
+        B = gcmarray{T, 1, InnerArray{T,2}}(
+            A.grid, A.meta, f_new, fSize_sub, fIndex_sub, thisversion)
+    else
+        # 2D result
+        f_new = OuterArray{InnerArray{T,2}, 2}(undef, nFaces, n3)
+        for I_iter in eachindex(tmpf)
+            f_new[I_iter] = view(tmpf[I_iter], :, :)
+        end
+        B = gcmarray{T, 2, InnerArray{T,2}}(
+            A.grid, A.meta, f_new, fSize_sub, fIndex_sub, thisversion)
+    end
+    return B
 end
 
 # ### Custom pretty-printing, similar, and broadcast
@@ -241,11 +275,65 @@ end
 
 import Base: display; display(X::AbstractMeshArray)=show(X)
 
-function Base.similar(A::gcmarray; m::varmeta=defaultmeta)
-    if ndims(A)==1
-        B = gcmarray(A.grid, eltype(A), A.fSize, A.fIndex; meta=m)
-    else
-        B = gcmarray(A.grid, eltype(A), A.fSize, A.fIndex, size(A)[2:end]...; meta=m)
+"""
+    similar(A::gcmarray; m::varmeta=defaultmeta, allocate=false)
+
+Create a gcmarray with the same structure and type as `A`.
+
+If `allocate=true`, eagerly allocates all face arrays sized according to `A.fSize`.
+If `allocate=false` (default), creates lazy empty placeholders; use when the array
+will be filled in-place (e.g., `read!`, `readtiles`).
+"""
+function Base.similar(A::gcmarray; m::varmeta=defaultmeta,allocate=false)
+    ElType = eltype(A)
+    nFaces = length(A.fIndex)
+
+    if ndims(A) == 1
+      if allocate
+        # 1D case: (nFaces,)
+        f = OuterArray{InnerArray{ElType,2},1}(undef, nFaces)
+        for a in 1:nFaces
+            f[a] = InnerArray{ElType}(undef, A.fSize[a]...)
+        end
+      else
+        f = fill(InnerArray{ElType}(undef, 0,0),nFaces)
+      end
+      B = gcmarray{ElType, 1, InnerArray{ElType,2}}(
+          A.grid, m, f, A.fSize, A.fIndex, thisversion)
+    elseif ndims(A) == 2
+        # 2D case: (nFaces, n3)
+        n3 = size(A, 2)
+        if allocate
+        f = OuterArray{InnerArray{ElType,2}, 2}(undef, nFaces, n3)
+        for a in 1:nFaces
+            for i3 in 1:n3
+                f[a, i3] = InnerArray{ElType}(undef, A.fSize[a]...)
+            end
+        end
+        else
+        f = fill(InnerArray{ElType}(undef, 0,0),nFaces,n3)
+        end
+        B = gcmarray{ElType, 2, InnerArray{ElType,2}}(
+            A.grid, m, f, A.fSize, A.fIndex, thisversion)
+    else  # ndims(A) == 3
+        # 3D case: (nFaces, n3, n4)
+        n3 = size(A, 2)
+        n4 = size(A, 3)
+        if allocate
+        f = OuterArray{InnerArray{ElType,2}, 3}(undef, nFaces, n3, n4)
+        for a in 1:nFaces
+            for i4 in 1:n4
+                for i3 in 1:n3
+                    f[a, i3, i4] = InnerArray{ElType}(undef, A.fSize[a]...)
+                end
+            end
+        end
+        else
+        f = fill(InnerArray{ElType}(undef, 0,0),nFaces,n3,n4)
+        end
+        B = gcmarray{ElType, 3, InnerArray{ElType,2}}(
+            A.grid, m, f, A.fSize, A.fIndex, thisversion)
+
     end
     return B
 end
@@ -256,10 +344,26 @@ Base.BroadcastStyle(::Type{<:AbstractMeshArray}) = Broadcast.ArrayStyle{Abstract
 
 function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{AbstractMeshArray}}, ::Type{ElType}) where ElType
     A = find_gcmarray(bc)
-    if ndims(A)==1
-        B = gcmarray(A.grid, ElType, A.fSize, A.fIndex)
-    else
-        B = gcmarray(A.grid, ElType, A.fSize, A.fIndex, size(A)[2:end]...)
+    nFaces = length(A.fIndex)
+
+    if ndims(A) == 1
+        # 1D case: lazy init — copyto! will allocate on first write
+        f = fill(InnerArray{ElType}(undef, 0,0), nFaces)
+        B = gcmarray{ElType, 1, InnerArray{ElType,2}}(
+            A.grid, defaultmeta, f, A.fSize, A.fIndex, thisversion)
+    elseif ndims(A) == 2
+        # 2D case: lazy init — copyto! will allocate on first write
+        n3 = size(A, 2)
+        f = fill(InnerArray{ElType}(undef, 0,0), nFaces, n3)
+        B = gcmarray{ElType, 2, InnerArray{ElType,2}}(
+            A.grid, defaultmeta, f, A.fSize, A.fIndex, thisversion)
+    else  # ndims(A) == 3
+        # 3D case: lazy init — copyto! will allocate on first write
+        n3 = size(A, 2)
+        n4 = size(A, 3)
+        f = fill(InnerArray{ElType}(undef, 0,0), nFaces, n3, n4)
+        B = gcmarray{ElType, 3, InnerArray{ElType,2}}(
+            A.grid, defaultmeta, f, A.fSize, A.fIndex, thisversion)
     end
     return B
 end

test/testsets/unitgrid.jl

@@ -27,4 +27,50 @@
     gr=Grids_simple.grid_add_z(gr,dep,msk)
 
     @test haskey(gr,:hFacC)
+
+    # Test read(xx::Array, x::AbstractMeshArray) — untested pathway
+    # Create a 4D array with shape (ioSize..., n3, n4) for PeriodicDomain grid
+    xx_data = randn(10, 10, 1, 1)  # 10×10 ioSize, no extra dims
+    A_template = MeshArray(γ, Float64)
+    B = read(xx_data, A_template)
+    @test isa(B, MeshArray)
+    @test eltype(B) == Float64
+    @test size(B) == size(A_template)
+
+    # Test read(xx::Array, x::AbstractMeshArray) with n3 > 1
+    xx_3d = randn(10, 10, 2, 1)
+    C_template = MeshArray(γ, Float64, 2)
+    C = read(xx_3d, C_template)
+    @test isa(C, MeshArray)
+    @test size(C) == size(C_template)
+    @test size(C, 2) == 2
+
+    # Test read(fil::String, x::AbstractMeshArray) with allocate=true
+    tmp_file = tempname()
+    write(tmp_file, Γ.XC)
+    A_template2 = MeshArray(γ, Float64)
+    D = read(tmp_file, A_template2)
+    @test isa(D, MeshArray)
+    @test size(D) == size(A_template2)
+    @test eltype(D) == Float64
+
+    # Test similar(x, allocate=true) pathway
+    A_1d = MeshArray(γ)
+    A_1d_full = similar(A_1d; allocate=true)
+    @test isa(A_1d_full, MeshArray)
+    @test size(A_1d_full) == size(A_1d)
+    # Verify faces are allocated (not 0×0)
+    @test size(A_1d_full.f[1]) == size(A_1d.f[1])
+
+    A_2d = MeshArray(γ, Float64, 3)
+    A_2d_full = similar(A_2d; allocate=true)
+    @test isa(A_2d_full, MeshArray)
+    @test size(A_2d_full) == size(A_2d)
+    @test size(A_2d_full.f[1]) == size(A_2d.f[1])
+
+    A_3d = MeshArray(γ, Float64, 3, 4)
+    A_3d_full = similar(A_3d; allocate=true)
+    @test isa(A_3d_full, MeshArray)
+    @test size(A_3d_full) == size(A_3d)
+    @test size(A_3d_full.f[1]) == size(A_3d.f[1])
 end