From 597c0631205c2350ba05d600cf768c1c0b663e50 Mon Sep 17 00:00:00 2001 From: gaelforget Date: Sat, 18 Jul 2026 01:19:19 +0200 Subject: [PATCH 01/10] hton fix --- src/ReadWrite.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/ReadWrite.jl b/src/ReadWrite.jl index c279b10..34583e8 100644 --- a/src/ReadWrite.jl +++ b/src/ReadWrite.jl @@ -19,7 +19,8 @@ function _load_binary(fil::String, y::AbstractMeshArray) fid = open(fil) xx = Array{eltype(y),1}(undef, n1*n2*n3*n4) read!(fid, xx) - xx = reshape(hton.(xx), (n1,n2,n3,n4)) + map!(hton, xx, xx) + xx = reshape(xx, (n1,n2,n3,n4)) close(fid) xx end @@ -236,7 +237,8 @@ function read_tiles(fil::String,x::AbstractMeshArray) fid = open(fil) xx = Array{eltype(x),2}(undef,(n1*n2,n3)) read!(fid,xx) - xx = reshape(hton.(xx),(n1,n2,n3)) + map!(hton, xx, xx) + xx = reshape(xx,(n1,n2,n3)) close(fid) return x.grid.read(xx,x) From 82f6c4c7fdde4ae62dce1ce8d019f1d537ef7eb6 Mon Sep 17 00:00:00 2001 From: gaelforget Date: Sat, 18 Jul 2026 17:51:48 +0200 Subject: [PATCH 02/10] streamline read operations --- src/ReadWrite.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ReadWrite.jl b/src/ReadWrite.jl index 34583e8..e090218 100644 --- a/src/ReadWrite.jl +++ b/src/ReadWrite.jl @@ -241,7 +241,7 @@ function read_tiles(fil::String,x::AbstractMeshArray) xx = reshape(xx,(n1,n2,n3)) close(fid) - return x.grid.read(xx,x) + return read_tiles(xx,x) end From 65cd5a6f2048ce6c253f4cfe2563d7f209539a55 Mon Sep 17 00:00:00 2001 From: gaelforget Date: Sat, 25 Jul 2026 13:44:40 +0200 Subject: [PATCH 03/10] fix : add an explicit Float64 type prefix, which routes through typed_hcat instead of generic hcat and avoids the SparseArrays dispatch --- src/Polygons.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Polygons.jl b/src/Polygons.jl index e0c7a2b..f8e56a2 100644 --- a/src/Polygons.jl +++ b/src/Polygons.jl @@ -48,9 +48,9 @@ function to_Polygons(XG,YG,ff=1) IJ1=1:n1; IJ2=1:n2 for i in IJ1 for j in IJ2 - x=[XG[ff][i,j] XG[ff][i+1,j] XG[ff][i+1,j+1] XG[ff][i,j+1] XG[ff][i,j]] + x=Float64[XG[ff][i,j] XG[ff][i+1,j] XG[ff][i+1,j+1] XG[ff][i,j+1] XG[ff][i,j]] treat_180lon!(x) - y=[YG[ff][i,j] YG[ff][i+1,j] YG[ff][i+1,j+1] YG[ff][i,j+1] YG[ff][i,j]] + y=Float64[YG[ff][i,j] YG[ff][i+1,j] YG[ff][i+1,j+1] YG[ff][i,j+1] YG[ff][i,j]] arr2[i,j]=GI.Polygon([ GI.LinearRing(GI.Point.(zip(vec(x),vec(y)))) ]) # arr2[i,j]=GI.LineString(GI.Point.(zip(vec(x),vec(y)))) end @@ -64,9 +64,9 @@ function to_LineStrings2D(XG,YG,ff=1;do_sphere=true) IJ1=1:n1; IJ2=1:n2 for i in IJ1 for j in IJ2 - x=[XG[ff][i,j] XG[ff][i+1,j] XG[ff][i+1,j+1] XG[ff][i,j+1] XG[ff][i,j]] + x=Float64[XG[ff][i,j] XG[ff][i+1,j] XG[ff][i+1,j+1] XG[ff][i,j+1] XG[ff][i,j]] treat_180lon!(x) - y=[YG[ff][i,j] YG[ff][i+1,j] YG[ff][i+1,j+1] YG[ff][i,j+1] YG[ff][i,j]] + y=Float64[YG[ff][i,j] YG[ff][i+1,j] YG[ff][i+1,j+1] YG[ff][i,j+1] YG[ff][i,j]] arr2[i,j]=GI.LineString(GI.Point.(zip(vec(x),vec(y)))) end end From 10de5acc68d03dee2a0b1409f7c3a93d5a75890f Mon Sep 17 00:00:00 2001 From: gaelforget Date: Sat, 25 Jul 2026 16:05:35 +0200 Subject: [PATCH 04/10] Fix CatView deep nesting in view(gcmarray,gcmvector) and bypass it in define_regions \n Replace element-by-element broadcast copyto! with face-by-face broadcast! for gcmarray/gcmvector --- src/Integration.jl | 4 ++-- src/types/gcmarray.jl | 30 +++++++++++++++++++----------- src/types/gcmvector.jl | 35 +++++++++++++++++------------------ 3 files changed, 38 insertions(+), 31 deletions(-) diff --git a/src/Integration.jl b/src/Integration.jl index 363b733..4959a1e 100644 --- a/src/Integration.jl +++ b/src/Integration.jl @@ -48,7 +48,7 @@ function define_regions(;option=:global,grid::NamedTuple) lats=[-90 ; -75:10:75 ; 90] nl=length(lats)-1 name=[Symbol("lat_$(lats[l])_to_$(lats[l+1])") for l in 1:nl] - [mask[findall((mask.>0)*(la.>=lats[l])*(la.0)*(la.>=lats[l])*(la.0)*(la.>=lats[i_a])*(la.=lons[i_o])*(lo.=lons[i_o])*(lo. _bc_face(a, I), bc.args) + Broadcast.Broadcasted(bc.f, new_args) # axes=nothing → recomputed from face arrays +end +_bc_face(a::AbstractMeshArray, I) = ndims(a.f) == 1 ? a.f[I[1]] : a.f[I] +_bc_face(a, I) = a # scalar, plain array, etc. + @inline function copyto!(dest::AbstractMeshArray, bc::Broadcast.Broadcasted{Nothing}) axes(dest) == axes(bc) || throwdm(axes(dest), axes(bc)) - # Performance optimization: broadcast!(identity, dest, A) is equivalent to copyto!(dest, A) if indices match - if bc.f === identity && bc.args isa Tuple{AbstractArray} # only a single input argument to broadcast! - A = bc.args[1] - if axes(dest) == axes(A) - return copyto!(dest, A) + for I in CartesianIndices(dest.f) + face_args = map(a -> _bc_face(a, I), bc.args) + array_args = filter(a -> isa(a, AbstractArray), face_args) + if !isempty(array_args) + result_axes = Base.Broadcast.combine_axes(array_args...) + result_size = map(length, result_axes) + if size(dest.f[I]) != result_size + dest.f[I] = similar(dest.f[I], result_size) + end end - end - bc′ = Broadcast.preprocess(dest, bc) - @simd for I in eachindex(bc′) - #@inbounds dest[I] = bc′[I] - @inbounds dest[I] = gcmarray_getindex_evalf(bc′,I) + broadcast!(bc.f, dest.f[I], face_args...) end return dest end diff --git a/src/types/gcmvector.jl b/src/types/gcmvector.jl index b4b6c96..06068bd 100644 --- a/src/types/gcmvector.jl +++ b/src/types/gcmvector.jl @@ -62,10 +62,12 @@ end Base.BroadcastStyle(::Type{<:gcmvector}) = Broadcast.ArrayStyle{gcmvector}() function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{gcmvector}}, ::Type{ElType}) where ElType - # Scan the inputs for the gcmarray: A = find_gcmvector(bc) - # Create the gcmvector output: - return gcmvector{ElType,ndims(A)}(A.grid,similar(A.f),copy(A.fSize),copy(A.fIndex)) + f = similar(A.f, Array{ElType,1}) + for I in eachindex(A.f) + f[I] = similar(A.f[I], ElType) + end + return gcmvector{ElType,ndims(A)}(A.grid,f,copy(A.fSize),copy(A.fIndex)) end find_gcmvector(bc::Base.Broadcast.Broadcasted) = find_gcmvector(bc.args) @@ -77,17 +79,17 @@ find_gcmvector(::Any, rest) = find_gcmvector(rest) # Specialize this method if all you want to do is specialize on typeof(dest) @inline function copyto!(dest::gcmvector, bc::Broadcast.Broadcasted{Nothing}) axes(dest) == axes(bc) || throwdm(axes(dest), axes(bc)) - # Performance optimization: broadcast!(identity, dest, A) is equivalent to copyto!(dest, A) if indices match - if bc.f === identity && bc.args isa Tuple{AbstractArray} # only a single input argument to broadcast! - A = bc.args[1] - if axes(dest) == axes(A) - return copyto!(dest, A) + for I in CartesianIndices(dest.f) + face_args = map(a -> _bc_face(a, I), bc.args) + array_args = filter(a -> isa(a, AbstractArray), face_args) + if !isempty(array_args) + result_axes = Base.Broadcast.combine_axes(array_args...) + result_size = map(length, result_axes) + if size(dest.f[I]) != result_size + dest.f[I] = similar(dest.f[I], result_size) + end end - end - bc′ = Broadcast.preprocess(dest, bc) - @simd for I in eachindex(bc′) - #@inbounds dest[I] = bc′[I] - @inbounds dest[I] = gcmarray_getindex_evalf(bc′,I) + broadcast!(bc.f, dest.f[I], face_args...) end return dest end @@ -126,11 +128,8 @@ function Base.getindex(A::gcmarray{T,N}, B::gcmvector{CartesianIndex{2},N}) wher end function Base.view(A::gcmarray{T,N}, B::gcmvector{CartesianIndex{2},N}) where {T,N} - tmpOut=missing - for a in eachindex(A) - ismissing(tmpOut) ? tmpOut=view(A[a],B[a]) : tmpOut=CatView(tmpOut,view(A[a],B[a])) - end - return tmpOut + views = [view(A[a], B[a]) for a in eachindex(A)] + return CatView(views...) end function Base.setindex!(A::gcmarray{T,N}, B::gcmvector{T,N}, C::gcmvector{CartesianIndex{2},N}) where {T,N} From c8969d6265027801472263761db964cd296e6c06 Mon Sep 17 00:00:00 2001 From: gaelforget Date: Sat, 25 Jul 2026 19:06:47 +0200 Subject: [PATCH 05/10] =?UTF-8?q?=20func=20no=20longer=20allocates=20a=20V?= =?UTF-8?q?ector=20per=20call=20=E2=80=94=20accumulates=20directly=20with?= =?UTF-8?q?=20a=20for=20loop,=20reusing=20xymsk(b)=20mask.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/Integration.jl | 10 ++++++++-- test/testsets/transport.jl | 2 +- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/src/Integration.jl b/src/Integration.jl index 4959a1e..9c5dee9 100644 --- a/src/Integration.jl +++ b/src/Integration.jl @@ -108,8 +108,14 @@ function define_sums(;option=:loops, grid::NamedTuple, regions=:global, depths=[ tmp2d=MeshArray(grid.XC.grid,Float32) zmsk(d0,d1,k) = layer_mask(grid.RF,d0,d1)[k] - func(X,b,d0,d1)=sum([sum(xymsk(b)*zmsk(d0,d1,k)*X[:,k]* - grid.DRF[k]*grid.hFacC[:,k]*grid.RAC) for k in 1:nr]) + func(X,b,d0,d1)=begin + s=0.0 + msk=xymsk(b) + for k in 1:nr + s+=sum(msk*zmsk(d0,d1,k)*X[:,k]*grid.DRF[k]*grid.hFacC[:,k]*grid.RAC) + end + s + end tmp3d=MeshArray(grid.XC.grid,Float32,nr) function func_v(X,d0,d1) diff --git a/test/testsets/transport.jl b/test/testsets/transport.jl index d348fb6..eeded10 100644 --- a/test/testsets/transport.jl +++ b/test/testsets/transport.jl @@ -36,7 +36,7 @@ uv=Dict("U"=>Tx,"V"=>Ty,"dimensions"=>["x","y"]) L=-85.0:5.0:85.0; LC=LatitudeCircles(L,Γ,format=:gridpath) T=Array{Float64,1}(undef,length(LC)) - [T[i]=1e-6*ThroughFlow(uv,LC[i],Γ) for i=1:length(LC)] + for i=1:length(LC); T[i]=1e-6*ThroughFlow(uv,LC[i],Γ); end plot(LC) plot(LC[1]) From b16ef62dd06dd5a7942d256ba7581fb1729879d2 Mon Sep 17 00:00:00 2001 From: gaelforget Date: Sat, 25 Jul 2026 19:39:01 +0200 Subject: [PATCH 06/10] Operations.land_mask: use setindex! (=) not broadcast (.=) on gcmvector results to bypass CatView materialize! and its ~40-level recursion --- src/Operations.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Operations.jl b/src/Operations.jl index fb10b6d..6cd83e1 100644 --- a/src/Operations.jl +++ b/src/Operations.jl @@ -277,8 +277,8 @@ function land_mask(m::AbstractMeshArray) for I in eachindex(m.f) μ.f[I]=copy(m.f[I]) end - μ[findall(μ.>0.0)].=1.0 - μ[findall(μ.==0.0)].=NaN + μ[findall(μ.>0.0)]=1.0 + μ[findall(μ.==0.0)]=NaN μ end From cbcbc58d6f01eddb086f3ee0516f5aaa07bef72a Mon Sep 17 00:00:00 2001 From: gaelforget Date: Sun, 26 Jul 2026 00:10:42 +0200 Subject: [PATCH 07/10] 1) pre-compute basin mask in define_sums; replace func/func_h with face-level loops to eliminate per-call gcmarray allocations ; 2) replace gcmarray broadcast+findall pattern in define_regions with face-level loop, eliminating CatView nesting --- src/Integration.jl | 156 ++++++++++++++++++++++++++++++--------------- 1 file changed, 106 insertions(+), 50 deletions(-) diff --git a/src/Integration.jl b/src/Integration.jl index 9c5dee9..07099b3 100644 --- a/src/Integration.jl +++ b/src/Integration.jl @@ -48,7 +48,19 @@ function define_regions(;option=:global,grid::NamedTuple) lats=[-90 ; -75:10:75 ; 90] nl=length(lats)-1 name=[Symbol("lat_$(lats[l])_to_$(lats[l+1])") for l in 1:nl] - [mask[findall((mask.>0)*(la.>=lats[l])*(la. 0 && laf[i,j] >= lats[l] && laf[i,j] < lats[l+1] + mf[i,j] = Float64(l) + end + end + end + end + (mask=mask,name=name) elseif isa(option,Tuple) dlo=option[1]; dla=option[2] @@ -90,81 +102,125 @@ layer_mask(dF,d0,d1)=begin md end +""" + func_h_inner(X, msk) + +Zero-allocation horizontal sum: replaces sum(xymsk(b)*X). +msk is pre-computed once per basin in define_sums. +""" +function func_h_inner(X::AbstractMeshArray, msk::AbstractMeshArray) + s = 0.0 + nf = length(X.fIndex) + for f in 1:nf + Xf = X.f[f] + mf = msk.f[f] + @inbounds for j in axes(Xf,2), i in axes(Xf,1) + s += mf[i,j] * Xf[i,j] + end + end + s +end + +""" + func_inner(X, msk, lmsk_d, grid, nr) + +Zero-allocation depth-integrated sum. +- msk: pre-computed horizontal mask, once per basin +- lmsk_d: pre-computed layer_mask vector for one depth range (length nr) +Replaces the loop over k with X[:,k] and hFacC[:,k] gcmarray slicing. +""" +function func_inner(X::AbstractMeshArray, msk::AbstractMeshArray, + lmsk_d, grid::NamedTuple, nr::Int) + s = 0.0 + nf = length(X.fIndex) + for k in 1:nr + zf = lmsk_d[k] * grid.DRF[k] + iszero(zf) && continue + for f in 1:nf + Xf = X.f[f, k] + hf = grid.hFacC.f[f, k] + mf = msk.f[f] + Rf = grid.RAC.f[f] + @inbounds for j in axes(Xf,2), i in axes(Xf,1) + s += zf * mf[i,j] * Xf[i,j] * hf[i,j] * Rf[i,j] + end + end + end + s +end + """ define_sums(;option=:loops, grid::NamedTuple, regions=:global, depths=[(0,7000)]) Define regional integration function for each basin and depth range. """ function define_sums(;option=:loops, grid::NamedTuple, regions=:global, depths=[(0,7000)]) - dep=(isa(depths,Tuple) ? [depths] : depths) - nd=length(dep) - rgns=define_regions(option=regions,grid=grid) - nb=length(rgns.name) - allones=1.0 .+0*grid.hFacC - nr=length(grid.RC) + dep = (isa(depths,Tuple) ? [depths] : depths) + nd = length(dep) + rgns = define_regions(option=regions, grid=grid) + nb = length(rgns.name) + allones = 1.0 .+ 0*grid.hFacC + nr = length(grid.RC) xymsk(b) = 1.0*(rgns.mask.==b) - func_h(X,b)=sum(xymsk(b)*X) - tmp2d=MeshArray(grid.XC.grid,Float32) - zmsk(d0,d1,k) = layer_mask(grid.RF,d0,d1)[k] - func(X,b,d0,d1)=begin - s=0.0 - msk=xymsk(b) - for k in 1:nr - s+=sum(msk*zmsk(d0,d1,k)*X[:,k]*grid.DRF[k]*grid.hFacC[:,k]*grid.RAC) - end - s - end - tmp3d=MeshArray(grid.XC.grid,Float32,nr) + # pre-compute layer masks once per depth range — avoids re-allocation in hot loop + lmsk_cache = [layer_mask(grid.RF, d0, d1) for (d0,d1) in dep] - function func_v(X,d0,d1) - tmp2d.=0.0 + zmsk(d0,d1,k) = layer_mask(grid.RF,d0,d1)[k] # kept for func_v only + + tmp2d = MeshArray(grid.XC.grid, Float32) + + function func_v(X, d0, d1) # not yet optimized + tmp2d .= 0.0 for k in 1:nr - tmp2d.+=zmsk(d0,d1,k)*X[:,k]*grid.DRF[k]*grid.hFacC[:,k]*grid.RAC + tmp2d .+= zmsk(d0,d1,k)*X[:,k]*grid.DRF[k]*grid.hFacC[:,k]*grid.RAC end tmp2d end + tmp3d = MeshArray(grid.XC.grid, Float32, nr) + if option==:streamlined_loop - #ocn_surf=[sum(xymsk(b)*(grid.hFacC[:,1].>0)*grid.RAC) for b in 1:nb] - BX=(name=String[],volsum=Function[],volume=Float64[], - ocn_surf=Float64[],tmp2d=tmp2d,tmp3d=tmp3d) - for b in 1:nb - for d in 1:nd - (d0,d1)=dep[d] - n=string(rgns.name[b])*"_dep_$(d0)_to_$(d1)" - @inline f=X->func(X,b,d0,d1) - v=f(allones) - push!(BX.name,n) - push!(BX.volsum,f) - push!(BX.volume,v) - #push!(BX.ocn_surf,ocn_surf[b]) - end - end + BX = (name=String[], volsum=Function[], volume=Float64[], + ocn_surf=Float64[], tmp2d=tmp2d, tmp3d=tmp3d) + for b in 1:nb + msk_b = xymsk(b) # once per basin + for d in 1:nd + (d0,d1) = dep[d] + lmsk_d = lmsk_cache[d] + n = string(rgns.name[b])*"_dep_$(d0)_to_$(d1)" + @inline f = X -> func_inner(X, msk_b, lmsk_d, grid, nr) + v = f(allones) + push!(BX.name, n) + push!(BX.volsum, f) + push!(BX.volume, v) + end + end end - BXh=(name=String[],hsum=Function[],tmp2d=tmp2d,tmp3d=tmp3d) + BXh = (name=String[], hsum=Function[], tmp2d=tmp2d, tmp3d=tmp3d) for b in 1:nb - n=string(rgns.name[b]) - @inline f=X->func_h(X,b) - push!(BXh.name,n) - push!(BXh.hsum,f) + msk_b = xymsk(b) # once per basin + n = string(rgns.name[b]) + @inline f = X -> func_h_inner(X, msk_b) + push!(BXh.name, n) + push!(BXh.hsum, f) end - BXv=(name=String[],vint=Function[],tmp2d=tmp2d,tmp3d=tmp3d) + BXv = (name=String[], vint=Function[], tmp2d=tmp2d, tmp3d=tmp3d) for d in 1:nd - (d0,d1)=dep[d] - n="$(d0)-$(d1)m" - @inline f=X->func_v(X,d0,d1) - push!(BXv.name,n) - push!(BXv.vint,f) + (d0,d1) = dep[d] + n = "$(d0)-$(d1)m" + @inline f = X -> func_v(X, d0, d1) + push!(BXv.name, n) + push!(BXv.vint, f) end if option==:loops - gridmask(rgns.mask,BXh.name,depths,BXh.hsum,BXv.vint,tmp2d,tmp3d) + gridmask(rgns.mask, BXh.name, depths, BXh.hsum, BXv.vint, tmp2d, tmp3d) elseif option==:streamlined_loop - gridmask(rgns.mask,BX.name,depths,BX.volsum,[],tmp2d,tmp3d) + gridmask(rgns.mask, BX.name, depths, BX.volsum, [], tmp2d, tmp3d) else error("unknown option") end From f24109c9f8df9cc2c1357b7a07264429ea268382 Mon Sep 17 00:00:00 2001 From: gaelforget Date: Mon, 27 Jul 2026 08:25:50 +0200 Subject: [PATCH 08/10] func_v_inner instead of func_v --- src/Integration.jl | 39 +++++++++++++++++++++++++++------------ 1 file changed, 27 insertions(+), 12 deletions(-) diff --git a/src/Integration.jl b/src/Integration.jl index 07099b3..bec0864 100644 --- a/src/Integration.jl +++ b/src/Integration.jl @@ -167,17 +167,31 @@ function define_sums(;option=:loops, grid::NamedTuple, regions=:global, depths=[ # pre-compute layer masks once per depth range — avoids re-allocation in hot loop lmsk_cache = [layer_mask(grid.RF, d0, d1) for (d0,d1) in dep] - zmsk(d0,d1,k) = layer_mask(grid.RF,d0,d1)[k] # kept for func_v only - tmp2d = MeshArray(grid.XC.grid, Float32) - function func_v(X, d0, d1) # not yet optimized - tmp2d .= 0.0 + function func_v_inner!(out::AbstractMeshArray, X::AbstractMeshArray, + lmsk_d, grid::NamedTuple, nr::Int) + for f in 1:length(out.fIndex) + of = out.f[f] + @inbounds for j in axes(of,2), i in axes(of,1) + of[i,j] = 0.0 + end + end for k in 1:nr - tmp2d .+= zmsk(d0,d1,k)*X[:,k]*grid.DRF[k]*grid.hFacC[:,k]*grid.RAC + zf = lmsk_d[k] * grid.DRF[k] + iszero(zf) && continue + for f in 1:length(out.fIndex) + of = out.f[f] + Xf = X.f[f,k] + hf = grid.hFacC.f[f,k] + Rf = grid.RAC.f[f] + @inbounds for j in axes(Xf,2), i in axes(Xf,1) + of[i,j] += zf * Xf[i,j] * hf[i,j] * Rf[i,j] + end + end end - tmp2d - end + out +end tmp3d = MeshArray(grid.XC.grid, Float32, nr) @@ -210,11 +224,12 @@ function define_sums(;option=:loops, grid::NamedTuple, regions=:global, depths=[ BXv = (name=String[], vint=Function[], tmp2d=tmp2d, tmp3d=tmp3d) for d in 1:nd - (d0,d1) = dep[d] - n = "$(d0)-$(d1)m" - @inline f = X -> func_v(X, d0, d1) - push!(BXv.name, n) - push!(BXv.vint, f) + (d0,d1) = dep[d] + lmsk_d = lmsk_cache[d] # already computed above + n = "$(d0)-$(d1)m" + @inline f = X -> func_v_inner!(tmp2d, X, lmsk_d, grid, nr) + push!(BXv.name, n) + push!(BXv.vint, f) end if option==:loops From 119466737d85e06449b121b6352dd52cc39d4a78 Mon Sep 17 00:00:00 2001 From: gaelforget Date: Mon, 27 Jul 2026 22:44:46 +0200 Subject: [PATCH 09/10] reduce allocs --- src/Operations.jl | 99 ++++++++++++++++++++++++++----------------- src/types/gcmarray.jl | 22 +++++----- 2 files changed, 70 insertions(+), 51 deletions(-) diff --git a/src/Operations.jl b/src/Operations.jl index 6cd83e1..5cadbbc 100644 --- a/src/Operations.jl +++ b/src/Operations.jl @@ -488,8 +488,16 @@ Returns a `gridpath` by default, or a `NamedTuple` with fields `lat`, `name`, """ function LatitudeCircle(lat,Γ::NamedTuple; format=:gridpath, range=(0.0,360.0)) - mskCint=1*(Γ.YC .>= lat) - mskC,mskW,mskS=edge_mask(mskCint) + + mskCint = similar(Γ.YC) + 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) + mf[i,j] = Float64(yf[i,j] >= lat) + end + end + + mskC,mskW,mskS=edge_mask(mskCint) restrict_longitudes!(mskC,Γ.XC,range=range) restrict_longitudes!(mskS,Γ.XS,range=range) restrict_longitudes!(mskW,Γ.XW,range=range) @@ -505,22 +513,40 @@ end is_in_lon_range(x,range)=(range[2].-range[1]>=360)|| (mod(x-range[1],360).0).&(tmp1[2:end-1,2:end-1].==0)) - # velocity masks: - mskW[i]=tmp1[2:end-1,2:end-1] - tmp1[1:end-2,2:end-1] - mskS[i]=tmp1[2:end-1,2:end-1] - tmp1[2:end-1,1:end-2] + mskC = similar(mskCint) + mskW = similar(mskCint) + mskS = similar(mskCint) + + exFLD = exchange(mskCint).MA # still needed for halo + + for i in eachindex(mskC.fIndex) + tmp1 = exFLD.f[i] # halo array, size (s1+2, s2+2) + mCf = mskC.f[i] + mWf = mskW.f[i] + mSf = mskS.f[i] + s1, s2 = size(mCf) + @inbounds for j in 1:s2, ii in 1:s1 + c = tmp1[ii+1, j+1] + neigh = tmp1[ii+1, j ] + tmp1[ii+1, j+2] + + tmp1[ii, j+1] + tmp1[ii+2, j+1] + mCf[ii,j] = Float64((neigh > 0) & (c == 0)) + mWf[ii,j] = c - tmp1[ii, j+1] + mSf[ii,j] = c - tmp1[ii+1, j ] + end end - #treat the case of blank tiles: - #mskC[findall(isnan.(mskC))].=0.0 - #mskW[findall(isnan.(mskW))].=0.0 - #mskS[findall(isnan.(mskS))].=0.0 - - return mskC,mskW,mskS + return mskC, mskW, mskS end ## diff --git a/src/types/gcmarray.jl b/src/types/gcmarray.jl index 4aad11d..9cdb7a3 100644 --- a/src/types/gcmarray.jl +++ b/src/types/gcmarray.jl @@ -241,11 +241,11 @@ end import Base: display; display(X::AbstractMeshArray)=show(X) -function Base.similar(A::gcmarray;m::varmeta=defaultmeta) +function Base.similar(A::gcmarray; m::varmeta=defaultmeta) if ndims(A)==1 - B=gcmarray(similar(A.grid),eltype(A),copy(A.fSize),copy(A.fIndex); meta=m) + B = gcmarray(A.grid, eltype(A), A.fSize, A.fIndex; meta=m) else - B=gcmarray(similar(A.grid),eltype(A),copy(A.fSize),copy(A.fIndex),size(A)[2:end]...; meta=m) + B = gcmarray(A.grid, eltype(A), A.fSize, A.fIndex, size(A)[2:end]...; meta=m) end return B end @@ -255,15 +255,13 @@ end Base.BroadcastStyle(::Type{<:AbstractMeshArray}) = Broadcast.ArrayStyle{AbstractMeshArray}() function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{AbstractMeshArray}}, ::Type{ElType}) where ElType - # Scan the inputs for the gcmarray: - A = find_gcmarray(bc) - # Create the gcmarray output: - if ndims(A)==1 - B=gcmarray(similar(A.grid),ElType,copy(A.fSize),copy(A.fIndex)) - else - B=gcmarray(similar(A.grid),ElType,copy(A.fSize),copy(A.fIndex),size(A)[2:end]...) - end - return B + 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]...) + end + return B end find_gcmarray(bc::Base.Broadcast.Broadcasted) = find_gcmarray(bc.args) From 124cc5a9d51917fc3289e51b74aaf3e7ab889e8a Mon Sep 17 00:00:00 2001 From: gaelforget Date: Mon, 27 Jul 2026 22:45:00 +0200 Subject: [PATCH 10/10] bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 18f588d..31a7e62 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "MeshArrays" uuid = "cb8c808f-1acf-59a3-9d2b-6e38d009f683" -version = "0.5.8" +version = "0.5.9" authors = ["gaelforget "] [deps]