From 6448041097ddb3058d2586296f86fec08ad7dc87 Mon Sep 17 00:00:00 2001 From: gaelforget Date: Thu, 22 Apr 2021 21:22:30 -0400 Subject: [PATCH 1/4] update dependency compats --- Project.toml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 616e80c0..6aedea83 100644 --- a/Project.toml +++ b/Project.toml @@ -22,9 +22,9 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" CFTime = "0.1" CSV = "0.6, 0.7, 0.8" CyclicArrays = "0.2, 0.3, 0.4, 0.5" -DataFrames = "0.21, 0.22" -MITgcmTools = "0.1.6" -MeshArrays = "0.2.14" +DataFrames = "0.21, 0.22, 1.0" +MITgcmTools = "0.1.20" +MeshArrays = "0.2.16" NetCDF = "0.10, 0.11" OceanStateEstimation = "0.1.5" OrdinaryDiffEq = "5" From 9417d3e2258d97138cc88e522294fe7aa4923ea7 Mon Sep 17 00:00:00 2001 From: gaelforget Date: Thu, 20 May 2021 08:43:57 -0400 Subject: [PATCH 2/4] update to MeshArrays v0.2.18 (Dict -> NamedTuple) --- Project.toml | 2 +- examples/basics/detailed_look.jl | 26 +++++------ examples/example123.jl | 2 + examples/flow_fields.jl | 45 +++++++++---------- examples/helper_functions.jl | 35 +++++++-------- examples/recipes_plots.jl | 2 +- examples/worldwide/three_dimensional_ocean.jl | 6 +-- src/data_wrangling.jl | 8 ++-- 8 files changed, 63 insertions(+), 63 deletions(-) diff --git a/Project.toml b/Project.toml index 6aedea83..70806240 100644 --- a/Project.toml +++ b/Project.toml @@ -24,7 +24,7 @@ CSV = "0.6, 0.7, 0.8" CyclicArrays = "0.2, 0.3, 0.4, 0.5" DataFrames = "0.21, 0.22, 1.0" MITgcmTools = "0.1.20" -MeshArrays = "0.2.16" +MeshArrays = "0.2.18" NetCDF = "0.10, 0.11" OceanStateEstimation = "0.1.5" OrdinaryDiffEq = "5" diff --git a/examples/basics/detailed_look.jl b/examples/basics/detailed_look.jl index 685109ed..ee6f3a0a 100644 --- a/examples/basics/detailed_look.jl +++ b/examples/basics/detailed_look.jl @@ -46,8 +46,8 @@ df=read_flt(dirIn,prec); # ## 4. Visualize Velocity Fields # # ``` -# plt=heatmap(Γ["mskW"][1,1].*𝑃.u0,title="U at the start") -# plt=heatmap(Γ["mskW"][1,1].*𝑃.u1-𝑃.u0,title="U end - U start") +# plt=heatmap(Γ.mskW[1,1].*𝑃.u0,title="U at the start") +# plt=heatmap(Γ.mskW[1,1].*𝑃.u1-𝑃.u0,title="U end - U start") # ``` # ## 5. Visualize Trajectories @@ -59,25 +59,25 @@ tmp[1:4,:] # Super-impose trajectory over velocity field (first for u ...) -x=Γ["XG"].f[1][:,1] -y=Γ["YC"].f[1][1,:] -z=transpose(Γ["mskW"][1].*𝑃.u0); +x=Γ.XG.f[1][:,1] +y=Γ.YC.f[1][1,:] +z=transpose(Γ.mskW[1].*𝑃.u0); # plt=contourf(x,y,z,c=:delta) # plot!(tmp[:,:lon],tmp[:,:lat],c=:red,w=4,leg=false) # Super-impose trajectory over velocity field (... then for v) -x=Γ["XC"].f[1][:,1] -y=Γ["YG"].f[1][1,:] -z=transpose(Γ["mskW"][1].*𝑃.v0); +x=Γ.XC.f[1][:,1] +y=Γ.YG.f[1][1,:] +z=transpose(Γ.mskW[1].*𝑃.v0); # plt=contourf(x,y,z,c=:delta) # plot!(tmp[:,:lon],tmp[:,:lat],c=:red,w=4,leg=false) # ## 6. Interpolate Velocities -dx=Γ["dx"] +dx=Γ.dx uInit=[tmp[1,:lon];tmp[1,:lat]]./dx nSteps=Int32(tmp[end,:time]/3600)-2 du=fill(0.0,2); @@ -95,9 +95,9 @@ for i=1:100 end #md plt=plot(tmpx,tmpu,label="u (interp)") -#md plot!(Γ["XG"].f[1][1:10,1]./dx,𝑃.u0[1:10,1],marker=:o,label="u (C-grid)") +#md plot!(Γ.XG.f[1][1:10,1]./dx,𝑃.u0[1:10,1],marker=:o,label="u (C-grid)") #md plot!(tmpx,tmpv,label="v (interp)") -#md plot!(Γ["XG"].f[1][1:10,1]./dx,𝑃.v0[1:10,1],marker=:o,label="v (C-grid)") +#md plot!(Γ.XG.f[1][1:10,1]./dx,𝑃.v0[1:10,1],marker=:o,label="v (C-grid)") # And similarly in the other direction @@ -112,9 +112,9 @@ for i=1:100 end # plt=plot(tmpx,tmpu,label="u (interp)") -# plot!(Γ["YG"].f[1][1,1:10]./dx,𝑃.u0[1,1:10],marker=:o,label="u (C-grid)") +# plot!(Γ.YG.f[1][1,1:10]./dx,𝑃.u0[1,1:10],marker=:o,label="u (C-grid)") # plot!(tmpx,tmpv,label="v (interp)") -# plot!(Γ["YG"].f[1][1,1:10]./dx,𝑃.v0[1,1:10],marker=:o,label="v (C-grid)") +# plot!(Γ.YG.f[1][1,1:10]./dx,𝑃.v0[1,1:10],marker=:o,label="v (C-grid)") # Compare recomputed velocities with those from `pkg/flt` diff --git a/examples/example123.jl b/examples/example123.jl index 492c371c..8f2ef8f8 100644 --- a/examples/example123.jl +++ b/examples/example123.jl @@ -131,6 +131,8 @@ function example2_setup() mskS=read_mds(γ.path*"hFacS",MeshArray(γ,Float32,nr)) mskS=1.0 .+ 0.0 * mask(mskS[:,kk],NaN,0.0) Γ=merge(Γ,Dict("mskW" => mskW, "mskS" => mskS)) + #convert to NamedTuple + Γ=(; zip(Symbol.(keys(Γ)), values(Γ))...) 𝑃=FlowFields(u0[1], u1[1], v0[1], v1[1], [t0,t1]) return 𝑃,Γ diff --git a/examples/flow_fields.jl b/examples/flow_fields.jl index d7b609f2..f26100b9 100644 --- a/examples/flow_fields.jl +++ b/examples/flow_fields.jl @@ -46,17 +46,17 @@ u,v,w=solid_body_rotation(12,4) """ function solid_body_rotation(np,nz) Γ=simple_periodic_domain(np); - γ=Γ["XC"].grid; + γ=Γ.XC.grid; #Solid-body rotation around central location ... i=Int(np/2+1) - u=-(Γ["YG"].-Γ["YG"][1][i,i]) - v=(Γ["XG"].-Γ["XG"][1][i,i]) + u=-(Γ.YG.-Γ.YG[1][i,i]) + v=(Γ.XG.-Γ.XG[1][i,i]) #... plus a convergent term to / from central location d=-0.01 - u=u+d*(Γ["XG"].-Γ["XG"][1][i,i]) - v=v+d*(Γ["YG"].-Γ["YG"][1][i,i]) + u=u+d*(Γ.XG.-Γ.XG[1][i,i]) + v=v+d*(Γ.YG.-Γ.YG[1][i,i]) #Replicate u,v in vertical dimension uu=MeshArray(γ,γ.ioPrec,nz) @@ -89,10 +89,10 @@ function global_ocean_circulation(;k=1,ny=2) p=dirname(pathof(IndividualDisplacements)) γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90) Γ=GridLoad(γ) - Γ=merge(Γ,IndividualDisplacements.NeighborTileIndices_cs(Γ)) + Γ=merge(Γ,NeighborTileIndices_cs(Γ)) - func=(u -> IndividualDisplacements.update_location_llc!(u,𝐷)) - Γ=merge(Γ,Dict("update_location!" => func)) + func=(u -> update_location_llc!(u,𝐷)) + Γ=merge(Γ,(; update_location! = func)) #initialize u0,u1 etc 𝑃,𝐷=set_up_FlowFields(k,Γ,ECCOclim_path); @@ -108,19 +108,19 @@ end """ OCCA_FlowFields(;backward_in_time::Bool=false) -Define gridded variables and return result as Dictionary (`uvetc`). +Define gridded variables and return result as NamedTuple """ function OCCA_FlowFields(;backward_in_time::Bool=false) γ=GridSpec("PeriodicChannel",MeshArrays.GRID_LL360) Γ=GridLoad(γ) - n=length(Γ["RC"]) + n=length(Γ.RC) n=5 - g=Γ["XC"].grid + g=Γ.XC.grid func=(u -> IndividualDisplacements.update_location_dpdo!(u,g)) - delete!.(Ref(Γ), ["hFacC", "hFacW", "hFacS","DXG","DYG","RAC","RAZ","RAS"]); + #delete!.(Ref(Γ), ["hFacC", "hFacW", "hFacS","DXG","DYG","RAC","RAZ","RAS"]); backward_in_time ? s=-1.0 : s=1.0 s=Float32(s) @@ -153,8 +153,8 @@ function OCCA_FlowFields(;backward_in_time::Bool=false) # 𝑆=read(rd(fileIn,"salt",n),MeshArray(γ,Float64,n)) for i in eachindex(u) - u[i]=u[i]./Γ["DXC"][1] - v[i]=v[i]./Γ["DYC"][1] + u[i]=u[i]./Γ.DXC[1] + v[i]=v[i]./Γ.DYC[1] end for i in eachindex(u) @@ -165,18 +165,18 @@ function OCCA_FlowFields(;backward_in_time::Bool=false) end for i in eachindex(w) - w[i]=w[i]./Γ["DRC"][min(i[2]+1,n)] + w[i]=w[i]./Γ.DRC[min(i[2]+1,n)] w[i]=circshift(w[i],[-180 0]) end - tmpx=circshift(Γ["XC"][1],[-180 0]) + tmpx=circshift(Γ.XC[1],[-180 0]) tmpx[1:180,:]=tmpx[1:180,:] .- 360.0 - Γ["XC"][1]=tmpx + Γ.XC[1]=tmpx - tmpx=circshift(Γ["XG"][1],[-180 0]) + tmpx=circshift(Γ.XG[1],[-180 0]) tmpx[1:180,:]=tmpx[1:180,:] .- 360.0 - Γ["XG"][1]=tmpx - Γ["Depth"][1]=circshift(Γ["Depth"][1],[-180 0]) + Γ.XG[1]=tmpx + Γ.Depth[1]=circshift(Γ.Depth[1],[-180 0]) t0=0.0; t1=86400*366*2.0; @@ -192,8 +192,8 @@ function OCCA_FlowFields(;backward_in_time::Bool=false) 𝑃=FlowFields(u,u,v,v,w,w,[t0,t1],func) - 𝐷 = (θ0=θ, θ1=θ, XC=exchange(Γ["XC"]), YC=exchange(Γ["YC"]), - RF=Γ["RF"], RC=Γ["RC"],ioSize=(360,160,n)) + 𝐷 = (θ0=θ, θ1=θ, XC=exchange(Γ.XC), YC=exchange(Γ.YC), + RF=Γ.RF, RC=Γ.RC,ioSize=(360,160,n)) return 𝑃,𝐷,Γ @@ -254,7 +254,6 @@ end function test2_periodic_domain(np = 12, nq = 12) #domain and time parameters Γ = simple_periodic_domain(np, nq) - Γ = IndividualDisplacements.dict_to_nt(Γ) u = 0.1 ./ Γ.DXC v = 0.3 ./ Γ.DYC diff --git a/examples/helper_functions.jl b/examples/helper_functions.jl index ec27ec01..09bdd3ac 100644 --- a/examples/helper_functions.jl +++ b/examples/helper_functions.jl @@ -39,7 +39,7 @@ end isosurface(θ,T,z) ``` -isosurface(𝐼.𝑃.θ0,15,Γ["RC"]) +isosurface(𝐼.𝑃.θ0,15,Γ.RC) ``` """ function isosurface(θ,T,z) @@ -49,17 +49,17 @@ function isosurface(θ,T,z) 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]) - d[j][i]=(1 .- a).*Γ["RC"][k] + a.*Γ["RC"][k+1] + d[j][i]=(1 .- a).*Γ.RC[k] + a.*Γ.RC[k+1] i=findall(isnan.(d[j]).&(θ[j,k].<=T).&(θ[j,k+1].>T)) a=(θ[j,k+1][i] .- T)./(θ[j,k+1][i] .- θ[j,k][i]) - d[j][i]=(1 .- a).*Γ["RC"][k+1] + a.*Γ["RC"][k] + d[j][i]=(1 .- a).*Γ.RC[k+1] + a.*Γ.RC[k] end end return d end """ - set_up_𝑃(k::Int,t::Float64,Γ::Dict,pth::String) + set_up_𝑃(k::Int,t::Float64,Γ::NamedTuple,pth::String) Define `FlowFields` data structure (𝑃) along with ancillary variables (𝐷) for the specified grid (`Γ` dictionnary), vertical level (`k`), and @@ -68,24 +68,24 @@ file location (`pth`). _Note: the initial implementation approximates month durations to 365 days / 12 months for simplicity and sets 𝑃.𝑇 to [-mon/2,mon/2]_ """ -function set_up_FlowFields(k::Int,Γ::Dict,pth::String) - XC=exchange(Γ["XC"]) #add 1 lon point at each edge - YC=exchange(Γ["YC"]) #add 1 lat point at each edge - iDXC=1. ./Γ["DXC"] - iDYC=1. ./Γ["DYC"] - γ=Γ["XC"].grid +function set_up_FlowFields(k::Int,Γ::NamedTuple,pth::String) + XC=exchange(Γ.XC) #add 1 lon point at each edge + YC=exchange(Γ.YC) #add 1 lat point at each edge + iDXC=1. ./Γ.DXC + iDYC=1. ./Γ.DYC + γ=Γ.XC.grid mon=86400.0*365.0/12.0 - func=Γ["update_location!"] + func=Γ.update_location! if k==0 - msk=Γ["hFacC"] + msk=Γ.hFacC (_,nr)=size(msk) 𝑃=FlowFields(MeshArray(γ,Float64,nr),MeshArray(γ,Float64,nr), MeshArray(γ,Float64,nr),MeshArray(γ,Float64,nr), MeshArray(γ,Float64,nr+1),MeshArray(γ,Float64,nr+1), [-mon/2,mon/2],func) else - msk=Γ["hFacC"][:, k] + msk=Γ.hFacC[:, k] 𝑃=FlowFields(MeshArray(γ,Float64),MeshArray(γ,Float64), MeshArray(γ,Float64),MeshArray(γ,Float64),[-mon/2,mon/2],func) end @@ -94,8 +94,7 @@ function set_up_FlowFields(k::Int,Γ::Dict,pth::String) XC=XC, YC=YC, iDXC=iDXC, iDYC=iDYC, k=k, msk=msk, θ0=similar(msk), θ1=similar(msk)) - tmp = IndividualDisplacements.dict_to_nt(IndividualDisplacements.NeighborTileIndices_cs(Γ)) - 𝐷 = merge(𝐷 , tmp) + 𝐷 = merge(𝐷 , NeighborTileIndices_cs(Γ)) return 𝑃,𝐷 end @@ -166,7 +165,7 @@ function update_FlowFields!(𝑃::𝐹_MeshArray3D,𝐷::NamedTuple,t::Float64) m1=mod(m1,12) m1==0 ? m1=12 : nothing - (_,nr)=size(𝐷.Γ["hFacC"]) + (_,nr)=size(𝐷.Γ.hFacC) (U,V)=read_velocities(𝑃.u0.grid,m0,𝐷.pth) u0=U; v0=V @@ -198,9 +197,9 @@ function update_FlowFields!(𝑃::𝐹_MeshArray3D,𝐷::NamedTuple,t::Float64) 𝑃.v1[:,:]=v1[:,:] for k=1:nr tmpw=exchange(-w0[:,k],1) - 𝑃.w0[:,k]=tmpw./𝐷.Γ["DRC"][k] + 𝑃.w0[:,k]=tmpw./𝐷.Γ.DRC[k] tmpw=exchange(-w1[:,k],1) - 𝑃.w1[:,k]=tmpw./𝐷.Γ["DRC"][k] + 𝑃.w1[:,k]=tmpw./𝐷.Γ.DRC[k] end 𝑃.w0[:,1]=0*exchange(-w0[:,1],1) 𝑃.w1[:,1]=0*exchange(-w1[:,1],1) diff --git a/examples/recipes_plots.jl b/examples/recipes_plots.jl index fd513096..40a95274 100644 --- a/examples/recipes_plots.jl +++ b/examples/recipes_plots.jl @@ -108,7 +108,7 @@ Compute Ocean depth logarithm on regular grid. function OceanDepthLog(Γ) lon=[i for i=20.:2.0:380., j=-79.:2.0:89.] lat=[j for i=20.:2.0:380., j=-79.:2.0:89.] - DL=interp_to_lonlat(Γ["Depth"],Γ,lon,lat) + DL=interp_to_lonlat(Γ.Depth,Γ,lon,lat) DL[findall(DL.<0)].=0 DL=transpose(log10.(DL)) DL[findall((!isfinite).(DL))].=NaN diff --git a/examples/worldwide/three_dimensional_ocean.jl b/examples/worldwide/three_dimensional_ocean.jl index d4b02653..5e2c7680 100644 --- a/examples/worldwide/three_dimensional_ocean.jl +++ b/examples/worldwide/three_dimensional_ocean.jl @@ -42,11 +42,11 @@ To convert from longitude,latitude here we take advantage of the regularity of the 1 degree grid being used -- for a more general alternative, see the global ocean example. """ -function initial_positions(Γ::Dict, nf=10000, lon_rng=(-160.0,-159.0), lat_rng=(30.0,31.0)) +function initial_positions(Γ::NamedTuple, nf=10000, lon_rng=(-160.0,-159.0), lat_rng=(30.0,31.0)) lon=lon_rng[1] .+(lon_rng[2]-lon_rng[1]).*rand(nf) lat=lat_rng[1] .+(lat_rng[2]-lat_rng[1]).*rand(nf) - x=lon .+ (21. - Γ["XC"][1][21,1]) - y=lat .+ (111. - Γ["YC"][1][1,111]) + x=lon .+ (21. - Γ.XC[1][21,1]) + y=lat .+ (111. - Γ.YC[1][1,111]) return x,y end diff --git a/src/data_wrangling.jl b/src/data_wrangling.jl index 9379264d..5b864069 100644 --- a/src/data_wrangling.jl +++ b/src/data_wrangling.jl @@ -14,7 +14,7 @@ function convert_to_FlowFields(U::Array{T,2},V::Array{T,2},t1::T) where T np,nq=size(U) Γ=simple_periodic_domain(np,nq) - g=Γ["XC"].grid + g=Γ.XC.grid u=MeshArray(g,[U]) v=MeshArray(g,[V]) (u,v)=exchange(u,v,1) @@ -187,8 +187,8 @@ lat=[j for i=20.:20.0:380., j=-70.:10.0:70.] (f,i,j,w,_,_,_)=InterpolationFactors(𝐷.Γ,vec(lon),vec(lat)) IntFac=(lon=lon,lat=lat,f=f,i=i,j=j,w=w) -tmp1=interp_to_lonlat(𝐷.Γ["Depth"],𝐷.Γ,lon,lat) -tmp2=interp_to_lonlat(𝐷.Γ["Depth"],IntFac) +tmp1=interp_to_lonlat(𝐷.Γ.Depth,𝐷.Γ,lon,lat) +tmp2=interp_to_lonlat(𝐷.Γ.Depth,IntFac) ref=[5896. 5896.] prod(isapprox.([maximum(tmp1) maximum(tmp2)],ref,atol=1.0)) @@ -198,7 +198,7 @@ prod(isapprox.([maximum(tmp1) maximum(tmp2)],ref,atol=1.0)) true ``` """ -function interp_to_lonlat(X::MeshArray,Γ::Dict,lon,lat) +function interp_to_lonlat(X::MeshArray,Γ::NamedTuple,lon,lat) (f,i,j,w,_,_,_)=InterpolationFactors(Γ,vec(lon),vec(lat)) return reshape(Interpolate(X,f,i,j,w),size(lon)) end From e06ff782ffe5080a0233baafce39a41cef786435 Mon Sep 17 00:00:00 2001 From: gaelforget Date: Thu, 20 May 2021 09:23:05 -0400 Subject: [PATCH 3/4] rm src/update_locations.jl (now in MeshArrays.jl) --- Project.toml | 2 +- examples/flow_fields.jl | 5 +- src/IndividualDisplacements.jl | 1 - src/update_locations.jl | 221 --------------------------------- 4 files changed, 5 insertions(+), 224 deletions(-) delete mode 100644 src/update_locations.jl diff --git a/Project.toml b/Project.toml index 70806240..6d5b6f02 100644 --- a/Project.toml +++ b/Project.toml @@ -22,7 +22,7 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" CFTime = "0.1" CSV = "0.6, 0.7, 0.8" CyclicArrays = "0.2, 0.3, 0.4, 0.5" -DataFrames = "0.21, 0.22, 1.0" +DataFrames = "0.21, 0.22, 1.0, 1.1" MITgcmTools = "0.1.20" MeshArrays = "0.2.18" NetCDF = "0.10, 0.11" diff --git a/examples/flow_fields.jl b/examples/flow_fields.jl index f26100b9..5445e495 100644 --- a/examples/flow_fields.jl +++ b/examples/flow_fields.jl @@ -120,7 +120,10 @@ function OCCA_FlowFields(;backward_in_time::Bool=false) g=Γ.XC.grid func=(u -> IndividualDisplacements.update_location_dpdo!(u,g)) - #delete!.(Ref(Γ), ["hFacC", "hFacW", "hFacS","DXG","DYG","RAC","RAZ","RAS"]); + jj=[:hFacC, :hFacW, :hFacS, :DXG, :DYG, :RAC, :RAZ, :RAS] + ii=findall([!in(i,jj) for i in keys(Γ)]) + Γ=(; zip(Symbol.(keys(Γ)[ii]), values(Γ)[ii])...) + backward_in_time ? s=-1.0 : s=1.0 s=Float32(s) diff --git a/src/IndividualDisplacements.jl b/src/IndividualDisplacements.jl index 03b50214..0d3c19e6 100644 --- a/src/IndividualDisplacements.jl +++ b/src/IndividualDisplacements.jl @@ -13,7 +13,6 @@ include("API.jl") include("compute.jl") include("data_wrangling.jl") include("read.jl") -include("update_locations.jl") export Individuals, ∫! export FlowFields, convert_to_FlowFields diff --git a/src/update_locations.jl b/src/update_locations.jl deleted file mode 100644 index cd8e6689..00000000 --- a/src/update_locations.jl +++ /dev/null @@ -1,221 +0,0 @@ -""" - location_is_out(u::AbstractArray{T,1},grid::gcmgrid) - -Test whether location (x,y,fIndex) is out of domain. If true then -one typically needs to call update_location_cs! or update_location_dpdo! -""" - -function location_is_out(u::AbstractArray{T,1},grid::gcmgrid) where T - u[1]<0|| u[1]> grid.fSize[Int(u[end])][1]|| - u[2]<0|| u[2]> grid.fSize[Int(u[end])][2] -end - -""" - NeighborTileIndices_dpdo(ni::Int,nj::Int) - -List of W, E, S, N neighbor tile IDs in the case of a doubly -periodic domain with ni x nj tiles. -""" -function NeighborTileIndices_dpdo(ni::Int,nj::Int) - tmp=fill(0,ni*nj,4) - for i=1:ni - for j=1:nj - k=i+ni*(j-1) - kS=j-1; kS==0 ? kS=nj : nothing; kS=i+ni*(kS-1) - kN=j+1; kN==nj+1 ? kN=1 : nothing; kN=i+ni*(kN-1) - kW=i-1; kW==0 ? kW=ni : nothing; kW=kW+ni*(j-1) - kE=i+1; kE==ni+1 ? kE=1 : nothing; kE=kE+ni*(j-1) - tmp[k,1]=kW - tmp[k,2]=kE - tmp[k,3]=kS - tmp[k,4]=kN - end - end - return tmp -end - -""" - update_location_cs! - -Update location (x,y,fIndex) when out of domain for cube sphere (cs) grid -as implemented by `MeshArrays.jl` (and MITgcm) -""" -function update_location_cs!(u::Array{Float64,1},𝑃::NamedTuple) - x,y = u[1:2] - fIndex = Int(u[end]) - nx,ny=𝑃.XC.fSize[fIndex] - if x<0||x>nx||y<0||y>ny - j = 0 - x<0 ? j=𝑃.aW[fIndex] : nothing - x>nx ? j=𝑃.aE[fIndex] : nothing - y<0 ? j=𝑃.aS[fIndex] : nothing - y>ny ? j=𝑃.aN[fIndex] : nothing - (x,y)=𝑃.RelocFunctions[j,fIndex](x,y) - u[1]=x - u[2]=y - u[end]=j - end - # - return u -end - -""" - update_location_llc! - -Update location (x,y,fIndex) when out of domain for lat-lon-cap (llc) grid -as implemented by `MeshArrays.jl` (and MITgcm) -""" -function update_location_llc!(u::Array{Float64,1},𝑃::NamedTuple) - x,y = u[1:2] - fIndex = Int(u[end]) - nx,ny=𝑃.XC.fSize[fIndex] - if y<0&&(fIndex==1||fIndex==2) - u[2]=eps(y) - elseif x>nx&&(fIndex==4||fIndex==5) - u[1]=nx-eps(x) - elseif x<0||x>nx||y<0||y>ny - j = 0 - x<0 ? j=𝑃.aW[fIndex] : nothing - x>nx ? j=𝑃.aE[fIndex] : nothing - y<0 ? j=𝑃.aS[fIndex] : nothing - y>ny ? j=𝑃.aN[fIndex] : nothing - (x,y)=𝑃.RelocFunctions[j,fIndex](x,y) - u[1]=x - u[2]=y - u[end]=j - end - # - return u -end - - -""" - update_location_dpdo! - -Update location (x,y,fIndex) when out of domain. Note: initially, this -only works for the `dpdo` grid type provided by `MeshArrays.jl`. -""" -function update_location_dpdo!(u::AbstractArray{T,1},grid::gcmgrid) where T - x,y = u[1:2] - fIndex = Int(u[3]) - # - nx,ny=grid.fSize[fIndex] - ni,nj=Int.(transpose(grid.ioSize)./grid.fSize[1]) - WESN=NeighborTileIndices_dpdo(ni,nj) - # - if x<0 - x=x+nx - u[1]=x - fIndex=WESN[fIndex,1] - u[3]=fIndex - elseif x>=nx - x=x-nx - u[1]=x - fIndex=WESN[fIndex,2] - u[3]=fIndex - end - # - if y<0 - y=y+ny - u[2]=y - fIndex=WESN[fIndex,3] - u[3]=fIndex - elseif y>=ny - y=y-ny - u[2]=y - fIndex=WESN[fIndex,4] - u[3]=fIndex - end - # - return u -end - - -""" - NeighborTileIndices_cs(grid::Dict) - -Derive list of neighboring tile indices for a cs or llc grid + functions that -convert indices from one tile to another. Returns a Dict to merge later. -""" -function NeighborTileIndices_cs(grid::Dict) - s = grid["XC"].fSize - nFaces = length(s) - nFaces == 5 ? s = vcat(s, s[3]) : nothing - aW=Array{Int,1}(undef,nFaces) - aE=similar(aW); aS=similar(aW); aN=similar(aW); - for i = 1:nFaces - (aW[i], aE[i], aS[i], aN[i], _, _, _, _) = MeshArrays.exch_cs_sources(i, s, 1) - end - RelocFunctions=RelocationFunctions_cs(grid["XC"]) - return Dict("aW" => aW, "aE" => aE, "aS" => aS, "aN" => aN, "RelocFunctions" => RelocFunctions) -end - -""" - RelocationFunctions_cs(xmpl) - -Define matrix of functions to convert indices across neighboring tiles -""" -function RelocationFunctions_cs(xmpl::MeshArray) - -# f1 : 0-n,0-n => -n-0,0-n for 1->2, 3->4, 5->6 -# f2 : 0-n,0-n => n-0,-n-0 for 2->4, 4->6, 6->2 -# f3 : 0-n,0-n => 0-n,-n-0 for 2->3, 4->5, 6->1 -# f4 : 0-n,0-n => -n-0,n-0 for 1->3, 3->5, 5->1 -# g1, g2, g3, g4 : the reverse connections - - f1(x, y, nx, ny) = (x .- Float64(nx), y) - f2(x, y, nx, ny) = (Float64(ny) .- y , x .- Float64(nx)) - f3(x, y, nx, ny) = (x, y .- Float64(ny)) - f4(x, y, nx, ny) = (y .- Float64(ny), Float64(nx) .- x ) - - g1(x, y, nx, ny) = (x .+ Float64(nx), y) - g2(x, y, nx, ny) = (y .+ Float64(ny), Float64(nx) .- x ) - g3(x, y, nx, ny) = (x, y .+ Float64(ny)) - g4(x, y, nx, ny) = (Float64(ny) .- y , x .+ Float64(nx)) - -# - - s = size.(xmpl.f) - nFaces = length(s) - tmp = Array{Function,2}(undef, 6, 6) - -# f1, f2, f3, f4 : always get nx & ny from the source tile - - tmp[2, 1] = (x, y) -> f1(x, y, s[1][1], s[1][2]) - tmp[4, 3] = (x, y) -> f1(x, y, s[3][1], s[3][2]) - tmp[6, 5] = (x, y) -> f1(x, y, s[5][1], s[5][2]) - - tmp[4, 2] = (x, y) -> f2(x, y, s[2][1], s[2][2]) - tmp[6, 4] = (x, y) -> f2(x, y, s[4][1], s[4][2]) - tmp[2, 6] = (x, y) -> f2(x, y, s[6][1], s[6][2]) - - tmp[3, 2] = (x, y) -> f3(x, y, s[2][1], s[2][2]) - tmp[5, 4] = (x, y) -> f3(x, y, s[4][1], s[4][2]) - tmp[1, 6] = (x, y) -> f3(x, y, s[6][1], s[6][2]) - - tmp[3, 1] = (x, y) -> f4(x, y, s[1][1], s[1][2]) - tmp[5, 3] = (x, y) -> f4(x, y, s[3][1], s[3][2]) - tmp[1, 5] = (x, y) -> f4(x, y, s[5][1], s[5][2]) - -# g1, g2, g3, g4 : nx or ny can come from source or target + notice nx/ny flips - - tmp[1, 2] = (x, y) -> g1(x, y, s[1][1], s[2][2]) - tmp[3, 4] = (x, y) -> g1(x, y, s[3][1], s[4][2]) - tmp[5, 6] = (x, y) -> g1(x, y, s[5][1], s[6][2]) - - tmp[2, 4] = (x, y) -> g2(x, y, s[4][1], s[2][1]) - tmp[4, 6] = (x, y) -> g2(x, y, s[6][1], s[4][1]) - tmp[6, 2] = (x, y) -> g2(x, y, s[2][1], s[6][1]) - - tmp[2, 3] = (x, y) -> g3(x, y, s[3][1], s[2][2]) - tmp[4, 5] = (x, y) -> g3(x, y, s[5][1], s[4][2]) - tmp[6, 1] = (x, y) -> g3(x, y, s[1][1], s[6][2]) - - tmp[1, 3] = (x, y) -> g4(x, y, s[1][2], s[3][2]) - tmp[3, 5] = (x, y) -> g4(x, y, s[3][2], s[5][2]) - tmp[5, 1] = (x, y) -> g4(x, y, s[5][2], s[1][2]) - - return tmp - -end - From 84f918344bbdd10588298440a3879ad74d5dbd9e Mon Sep 17 00:00:00 2001 From: gaelforget Date: Thu, 20 May 2021 09:30:42 -0400 Subject: [PATCH 4/4] bump version --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 6d5b6f02..59674942 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "IndividualDisplacements" uuid = "b92f0c32-5b7e-11e9-1d7b-238b2da8b0e6" authors = ["gaelforget "] -version = "0.3.0" +version = "0.3.1" [deps] CFTime = "179af706-886a-5703-950a-314cd64e0468" @@ -24,7 +24,7 @@ CSV = "0.6, 0.7, 0.8" CyclicArrays = "0.2, 0.3, 0.4, 0.5" DataFrames = "0.21, 0.22, 1.0, 1.1" MITgcmTools = "0.1.20" -MeshArrays = "0.2.18" +MeshArrays = "0.2.19" NetCDF = "0.10, 0.11" OceanStateEstimation = "0.1.5" OrdinaryDiffEq = "5"