diff --git a/Project.toml b/Project.toml index a176e2c..f0982ff 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "MeshArrays" uuid = "cb8c808f-1acf-59a3-9d2b-6e38d009f683" -version = "0.5.6" +version = "0.5.7" authors = ["gaelforget "] [deps] diff --git a/src/GridPaths.jl b/src/GridPaths.jl index a84d00f..7a823ff 100644 --- a/src/GridPaths.jl +++ b/src/GridPaths.jl @@ -64,9 +64,9 @@ function shorter_paths!(xyz,xyz0,msk_in) theta[2]=atan(y0[2],x0[2]) γ=msk_in[1].grid - tmpx=γ.write(x); tmpy=γ.write(y); tmpz=γ.write(z); + tmpx=γ.write(x); tmpy=γ.write(y); tmpz=γ.write(z) tmptheta=atan.(tmpy,tmpx) - if theta[2]<0; + if theta[2]<0 tmp00=findall(tmptheta.<=theta[2]) tmptheta[tmp00].=tmptheta[tmp00].+2*pi theta[2]=theta[2]+2*pi @@ -83,7 +83,7 @@ function shorter_paths!(xyz,xyz0,msk_in) else tmpm[findall( (tmptheta.<=theta[2]).&(tmptheta.>=theta[1]) )].=0.0 end - mm=γ.read(tmpm,mm); + mm=γ.read(tmpm,mm) #store result: push!(msk_out,mm) end diff --git a/src/Integration.jl b/src/Integration.jl index ae58ab0..363b733 100644 --- a/src/Integration.jl +++ b/src/Integration.jl @@ -175,7 +175,7 @@ begin @everywhere rd(F,var,tim,tmp)=read(read_mdsio(F,var),tmp) @everywhere G,M,files=Integration.example() #,regions=(30,10),depths=Integration.DEPTHS) -end; +end H=Integration.loops(M,files=files,rd=rd) # Hbis=Integration.streamlined_loop(M,files=files,rd=rd) diff --git a/src/Interpolation.jl b/src/Interpolation.jl index 3f68789..4cb8328 100644 --- a/src/Interpolation.jl +++ b/src/Interpolation.jl @@ -8,7 +8,7 @@ import NearestNeighbors: knn Find k nearest neighbors to each point in x,y on xgrid,ygrid ``` -lon=collect(0.1:0.5:2.1); lat=collect(0.1:0.5:2.1); +lon=collect(0.1:0.5:2.1); lat=collect(0.1:0.5:2.1) (f,i,j,c)=knn(Γ.XC,Γ.YC,lon,lat) ``` """ @@ -32,12 +32,12 @@ function knn(xgrid::AbstractMeshArray,ygrid::AbstractMeshArray, kk=findall(isfinite.(a_x)) x=sin.(pi/2 .-a_y[kk]*pi/180).*cos.(a_x[kk]*pi/180) y=sin.(pi/2 .-a_y[kk]*pi/180).*sin.(a_x[kk]*pi/180) - z=cos.(pi/2 .-a_y[kk]*pi/180); + z=cos.(pi/2 .-a_y[kk]*pi/180) #vector of target points in Cartesian, 3D, coordinates - xx=sin.(pi/2 .-yvec*pi/180).*cos.(xvec*pi/180); - yy=sin.(pi/2 .-yvec*pi/180).*sin.(xvec*pi/180); - zz=cos.(pi/2 .-yvec*pi/180); + xx=sin.(pi/2 .-yvec*pi/180).*cos.(xvec*pi/180) + yy=sin.(pi/2 .-yvec*pi/180).*sin.(xvec*pi/180) + zz=cos.(pi/2 .-yvec*pi/180) #define tree kdtree = KDTree([x y z]') @@ -152,7 +152,7 @@ function InterpolationFactors(Γ,lon::Array{T,1},lat::Array{T,1}) where {T} s=fill(0,2*length(fs)) [s[collect(1:2) .+ (i-1)*2]=collect(fs[i]) for i in 1:length(fs)] ni=gcd(s); nj=gcd(s); γ=Γ.XC.grid - τ=Tiles(γ,ni,nj); tiles=MeshArray(γ,Int); + τ=Tiles(γ,ni,nj); tiles=MeshArray(γ,Int) [tiles[τ[ii].face][τ[ii].i,τ[ii].j].=ii for ii in 1:length(τ)] #2. t_XC, t_XC, t_f, t_i, t_j @@ -239,7 +239,7 @@ end interpolation_setup(γ::gcmgrid) Download/read precomputed interpolation coefficients for gcmgrid `γ` if they exist, - currently for LLC90 and LLC270, otherwise compute interpolation coefficients + currently for LLC90 and LLC270, otherwise compute interpolation coefficients defaulting to 2° grid. """ function interpolation_setup(γ::gcmgrid) @@ -255,7 +255,7 @@ function interpolation_setup(γ::gcmgrid) # Default to 2 degree interpolation weights for other grids Γ=GridLoad(γ; option="full") interpolation_setup(; Γ) - end + end end """ @@ -301,32 +301,32 @@ x,y=StereographicProjection(45.,60.,lon,lat) """ function StereographicProjection(XC0::Number,YC0::Number,XC,YC) #compute spherical coordinates: - phi=XC; theta=90 .-YC; - phi0=XC0; theta0=90-YC0; + phi=XC; theta=90 .-YC + phi0=XC0; theta0=90-YC0 #compute cartesian coordinates: - X=sind.(theta).*cosd.(phi); - Y=sind.(theta).*sind.(phi); - Z=cosd.(theta); + X=sind.(theta).*cosd.(phi) + Y=sind.(theta).*sind.(phi) + Z=cosd.(theta) - x=X; y=Y; z=Z; + x=X; y=Y; z=Z #bring chosen point to the north pole: - xx=x; yy=y; zz=z; - x=cosd(phi0).*xx+sind(phi0).*yy; - y=-sind(phi0).*xx+cosd(phi0).*yy; - z=zz; + xx=x; yy=y; zz=z + x=cosd(phi0).*xx+sind(phi0).*yy + y=-sind(phi0).*xx+cosd(phi0).*yy + z=zz - xx=x; yy=y; zz=z; - x=cosd(theta0)*xx-sind(theta0)*zz; - y=yy; - z=sind(theta0)*xx+cosd(theta0)*zz; + xx=x; yy=y; zz=z + x=cosd(theta0)*xx-sind(theta0)*zz + y=yy + z=sind(theta0)*xx+cosd(theta0)*zz #stereographic projection from the south pole: - xx=x./(1 .+z); - yy=y./(1 .+z); + xx=x./(1 .+z) + yy=y./(1 .+z) - #nrm=sqrt(xx.^2+yy.^2); + #nrm=sqrt(xx.^2+yy.^2) #msk=1+0*nrm; msk(nrm>tan(pi/4/2))=NaN;%mask points outside of pi/4 cone return xx,yy @@ -367,8 +367,8 @@ specifies one polygon. (optional) `x,y` are position vectors. ```jldoctest; output = false using MeshArrays -px=[0. 0. 1. 1.]; py=[0. 1. 1. 0.]; -x=collect(-1.0:0.25:2.0); y=x; +px=[0. 0. 1. 1.]; py=[0. 1. 1. 0.] +x=collect(-1.0:0.25:2.0); y=x angsum=fill(0.0,1,length(x)) MeshArrays.PolygonAngle(px,py,x,y,angsum) @@ -439,7 +439,7 @@ function QuadArrays(x_grid::Array{T,2},y_grid::Array{T,2}) where {T} tmp=collect(0+di:ni+di)*ones(1,nj+1) i_quad[:,pp]=vec(tmp) - tmp=ones(ni+1,1)*transpose(collect(0+dj:nj+dj)); + tmp=ones(ni+1,1)*transpose(collect(0+dj:nj+dj)) j_quad[:,pp]=vec(tmp) end @@ -469,13 +469,13 @@ function QuadCoeffs(px,py,ox=[],oy=[],ow=[]) # ParaCoeffs([0., 2., 3., 1.]',[0., 0., 1., 1.]',0.1,0.1) #1. solve linear problem (`a,b` vectors from `px,py`) - # A=[1 0 0 0;1 1 0 0;1 1 1 1;1 0 1 0]; AI = inv(A); - # AI=[1 0 0 0;-1 1 0 0;-1 0 0 1; 1 -1 1 -1]; - # a = AI*px'; - # b = AI*py'; + # A=[1 0 0 0;1 1 0 0;1 1 1 1;1 0 1 0]; AI = inv(A) + # AI=[1 0 0 0;-1 1 0 0;-1 0 0 1; 1 -1 1 -1] + # a = AI*px' + # b = AI*py' #This defines the mapping from logical `l,m` to physical `x,y` as - # x=a(1)+a(2)*l+a(3)*m+a(2)*l*m; - # y=b(1)+b(2)*l+b(3)*m+b(2)*l*m; + # x=a(1)+a(2)*l+a(3)*m+a(2)*l*m + # y=b(1)+b(2)*l+b(3)*m+b(2)*l*m a=[px[1] -px[1]+px[2] -px[1]+px[4] px[1]-px[2]+px[3]-px[4]] a[findall(abs.(a).<1e-8)].=0.0 @@ -496,18 +496,18 @@ function QuadCoeffs(px,py,ox=[],oy=[],ow=[]) #3. solve non-linear problem for `pl,pm` from `px,py` & `a,b` # This defines the mapping from physical `x,y` to logical `l,m` - a=reshape(a,(size(a,1),1,size(a,2))); - b=reshape(b,(size(b,1),1,size(b,2))); + a=reshape(a,(size(a,1),1,size(a,2))) + b=reshape(b,(size(b,1),1,size(b,2))) # quadratic equation coeffs, `aa*mm^2+bb*m+cc=0` - if ~isempty(ox); - x=ox; y=oy; - else; - x=px; y=py; - a=repeat(a,1,size(x,2),1); - b=repeat(b,1,size(x,2),1); - sgn=repeat(sgn,1,size(x,2)); - end; + if ~isempty(ox) + x=ox; y=oy + else + x=px; y=py + a=repeat(a,1,size(x,2),1) + b=repeat(b,1,size(x,2),1) + sgn=repeat(sgn,1,size(x,2)) + end # det=fill(0.0,size(x)) pm=fill(0.0,size(x)) @@ -523,7 +523,7 @@ function QuadCoeffs(px,py,ox=[],oy=[],ow=[]) pl[ii,jj] = (x[ii,jj]-a[ii,jj,1]-a[ii,jj,3]*pm[ii,jj])/(a[ii,jj,2]+a[ii,jj,4]*pm[ii,jj]) end - if ~isempty(ox); + if ~isempty(ox) tmp1=(1 .-pl).*(1 .-pm) tmp2=pl.*(1 .-pm) tmp3=pl.*pm @@ -556,21 +556,21 @@ println(vec(QuadCoeffs([0., 2.01, 3., 1.]',[0., 0., 1., 1.]',x,y))) """ function ParaCoeffs(px,py,ox=[],oy=[]) - tmp1=px[:,1]; - tmp2=-px[:,1]+px[:,2]; - tmp3=-px[:,2]+px[:,3]; - a=[tmp1 tmp2 tmp3]; + tmp1=px[:,1] + tmp2=-px[:,1]+px[:,2] + tmp3=-px[:,2]+px[:,3] + a=[tmp1 tmp2 tmp3] - tmp1=py[:,1]; - tmp2=-py[:,1]+py[:,2]; - tmp3=-py[:,2]+py[:,3]; - b=[tmp1 tmp2 tmp3]; + tmp1=py[:,1] + tmp2=-py[:,1]+py[:,2] + tmp3=-py[:,2]+py[:,3] + b=[tmp1 tmp2 tmp3] # (m,l)=inv([a[1,2] a[1,3];b[1,2] b[1,3]])*[ox[1]-a[1,1]; oy[1]-b[1,1]] m=( b[:,3].*(ox-a[:,1])-a[:,3].*(oy-b[:,1]) ) ./(a[:,2].*b[:,3]-a[:,3].*b[:,2]) l=( -b[:,2].*(ox-a[:,1])+a[:,2].*(oy-b[:,1]) ) ./(a[:,2].*b[:,3]-a[:,3].*b[:,2]) - ow=[]; + ow=[] tmp1=(1 .-l).*(1 .-m) tmp4=l.*(1 .-m) tmp3=l.*m diff --git a/src/Operations.jl b/src/Operations.jl index d77faf1..2f40753 100644 --- a/src/Operations.jl +++ b/src/Operations.jl @@ -167,8 +167,6 @@ end -export f, β - Ω = 7.2921 * 10^(-5) #rad/s R = 6.3781 * 10^(6) #m #Coriolis Parameter @@ -191,8 +189,8 @@ function EkmanTrsp(u::AbstractMeshArray,v::AbstractMeshArray,Γ::NamedTuple) fcur=f_Coriolis.(Γ.YS[i]) fcur[abs.(Γ.YS[i]).<10].=NaN EkY[i]=-ucur./fcur./1029.0 - ucur=1/4* (V.MA[i][1:end-2,2:end-1]+V.MA[i][1:end-2,2:end-1] - +V.MA[i][2:end-1,3:end]+V.MA[i][2:end-1,3:end]) + ucur=1/4* (V.MA[i][1:end-2,2:end-1]+V.MA[i][1:end-2,3:end] + +V.MA[i][2:end-1,2:end-1]+V.MA[i][2:end-1,3:end]) fcur=f_Coriolis.(Γ.YW[i]) fcur[abs.(Γ.YW[i]).<10].=NaN EkX[i]=ucur./fcur./1029.0 @@ -268,7 +266,10 @@ end Define land mask from `m` (1 if m>0; NaN if otherwise). """ function land_mask(m::AbstractMeshArray) - μ=m + μ=similar(m) + for I in eachindex(m.f) + μ.f[I]=copy(m.f[I]) + end μ[findall(μ.>0.0)].=1.0 μ[findall(μ.==0.0)].=NaN μ @@ -286,9 +287,9 @@ Compute convergence of a vector field function convergence(uFLD::AbstractMeshArray,vFLD::AbstractMeshArray) #important note: -# Normally uFLD, vFLD should not contain any NaN; +# Normally uFLD, vFLD should not contain any NaN # if otherwise then something this may be needed: -# uFLD=mask(uFLD,0.0); vFLD=mask(vFLD,0.0); +# uFLD=mask(uFLD,0.0); vFLD=mask(vFLD,0.0) CONV=similar(uFLD) @@ -338,15 +339,15 @@ for a=1:FLD.grid.nFaces end #Before scaling the diffusive operator ... -tmp0=DXCsm*iDXC*mskW; -tmp00=maximum(tmp0); -tmp0=DYCsm*iDYC*mskS; -tmp00=max(tmp00,maximum(tmp0)); +tmp0=DXCsm*iDXC*mskW +tmp00=maximum(tmp0) +tmp0=DYCsm*iDYC*mskS +tmp00=max(tmp00,maximum(tmp0)) #... determine a suitable time period: -nbt=ceil(1.1*2*tmp00^2); -dt=1.; -T=nbt*dt; +nbt=ceil(1.1*2*tmp00^2) +dt=1. +T=nbt*dt #println("nbt="*"$nbt") #diffusion operator times DYG / DXG @@ -358,16 +359,16 @@ dtFac=dt*mskC/Γ.RAC #loop: for it=1:nbt - (dTdxAtU,dTdyAtV)=gradient(FLD,iDXC,iDYC); + (dTdxAtU,dTdyAtV)=gradient(FLD,iDXC,iDYC) tmpU=similar(FLD) tmpV=similar(FLD) for a=1:FLD.grid.nFaces - tmpU.f[a]=dTdxAtU.f[a].*KuxFac.f[a]; - tmpV.f[a]=dTdyAtV.f[a].*KvyFac.f[a]; + tmpU.f[a]=dTdxAtU.f[a].*KuxFac.f[a] + tmpV.f[a]=dTdyAtV.f[a].*KvyFac.f[a] end - tmpC=convergence(tmpU,tmpV); + tmpC=convergence(tmpU,tmpV) for a=1:FLD.grid.nFaces - FLD.f[a]=FLD.f[a]-dtFac.f[a].*tmpC.f[a]; + FLD.f[a]=FLD.f[a]-dtFac.f[a].*tmpC.f[a] end end @@ -421,8 +422,8 @@ function ThroughFlow(VectorField,IntegralPath,Γ::NamedTuple,msk=[]) #do_dz==1 ? mskS=Γ.DRF[i3]*mskS : nothing # #method 2: less slow - tabW=(isa(IntegralPath,NamedTuple) ? IntegralPath.tabW : IntegralPath.W) - tabS=(isa(IntegralPath,NamedTuple) ? IntegralPath.tabS : IntegralPath.S) + tabW=IntegralPath.W + tabS=IntegralPath.S for i4=1:n[4] #method 1: quite slow #trsp[1,i3,i4]=sum(mskW*U[:,:,i3,i4])+sum(mskS*V[:,:,i3,i4]) @@ -460,18 +461,18 @@ end Compute integration paths that follow latitude circles, within the specified longitude `range`. """ -function LatitudeCircles(LatValues,Γ::NamedTuple; +function LatitudeCircles(LatValues,Γ::NamedTuple; format=:gridpath, range=(0.0,360.0)) T=(format==:NamedTuple ? NamedTuple : gridpath) LatitudeCircles=Array{T}(undef,length(LatValues)) for j=1:length(LatValues) - LatitudeCircles[j]=LatitudeCircle(LatValues[j],Γ; + LatitudeCircles[j]=LatitudeCircle(LatValues[j],Γ; format=format,range=range) end (length(LatValues)==1 ? LatitudeCircles[1] : LatitudeCircles) end -function LatitudeCircle(lat,Γ::NamedTuple; +function LatitudeCircle(lat,Γ::NamedTuple; format=:gridpath, range=(0.0,360.0)) mskCint=1*(Γ.YC .>= lat) mskC,mskW,mskS=edge_mask(mskCint) @@ -479,13 +480,12 @@ function LatitudeCircle(lat,Γ::NamedTuple; restrict_longitudes!(mskS,Γ.XS,range=range) restrict_longitudes!(mskW,Γ.XW,range=range) LC=if format==:NamedTuple - (lat=LatValues[j],tabC=MskToTab(mskC), - tabW=MskToTab(mskW),tabS=MskToTab(mskS)) + (lat=lat,name="Parallel $lat", grid=Γ, + C=MskToTab(mskC),W=MskToTab(mskW),S=MskToTab(mskS)) else gridpath(name="Parallel $lat", grid=Γ, C=MskToTab(mskC),W=MskToTab(mskW),S=MskToTab(mskS)) end - end is_in_lon_range(x,range)=(range[2].-range[1]>=360)|| @@ -577,7 +577,7 @@ function Transect(name,lons,lats,Γ; segment=:short, format=:gridpath) tabS=MskToTab(mskSedge) if format==:NamedTuple - (name=name,tabC=tabC,tabW=tabW,tabS=tabS) + (name=name,C=tabC,W=tabW,S=tabS) else gridpath(name=name,grid=Γ,C=tabC,W=tabW,S=tabS) end @@ -678,29 +678,29 @@ compute bolus velocty field (bolusU,bolusV,bolusW) from gm streamfunction (GM_PsiX,GM_PsiY) """ function calc_bolus(GM_PsiX,GM_PsiY, Γ) - nr=length(Γ.RC); - mskW = 0 .*Γ.hFacW; mskS = 0 .*Γ.hFacS; - mskC = 0 .*Γ.hFacC; + nr=length(Γ.RC) + mskW = 0 .*Γ.hFacW; mskS = 0 .*Γ.hFacS + mskC = 0 .*Γ.hFacC for ff in eachindex(mskC) mskW.f[ff][Γ.hFacW.f[ff] .> 0] .= 1 - mskS.f[ff][Γ.hFacS.f[ff] .> 0] .= 1; - mskC.f[ff][Γ.hFacC.f[ff] .> 0] .= 1; + mskS.f[ff][Γ.hFacS.f[ff] .> 0] .= 1 + mskC.f[ff][Γ.hFacC.f[ff] .> 0] .= 1 end - GM_PsiX[findall((!isfinite).(GM_PsiX))]=0; - GM_PsiY[findall((!isfinite).(GM_PsiY))]=0; + GM_PsiX[findall((!isfinite).(GM_PsiX))]=0 + GM_PsiY[findall((!isfinite).(GM_PsiY))]=0 - bolusU=0*Γ.hFacW; - bolusV=0*Γ.hFacS; - for k=1:nr-1; - bolusU.f[:,k].=(GM_PsiX.f[:,k+1].-GM_PsiX.f[:,k])/Γ.DRF[k]; - bolusV.f[:,k].=(GM_PsiY.f[:,k+1].-GM_PsiY.f[:,k])/Γ.DRF[k]; - end; - bolusU.f[:, nr] .= 0*GM_PsiX.f[:,nr] .-GM_PsiX.f[:,nr]./Γ.DRF[nr]; - bolusV.f[:, nr] .= 0*GM_PsiY.f[:,nr] .-GM_PsiY.f[:,nr]./Γ.DRF[nr]; - - bolusU=bolusU.*mskW; - bolusV=bolusV.*mskS; + bolusU=0*Γ.hFacW + bolusV=0*Γ.hFacS + for k=1:nr-1 + bolusU.f[:,k].=(GM_PsiX.f[:,k+1].-GM_PsiX.f[:,k])/Γ.DRF[k] + bolusV.f[:,k].=(GM_PsiY.f[:,k+1].-GM_PsiY.f[:,k])/Γ.DRF[k] + end + bolusU.f[:, nr] .= 0*GM_PsiX.f[:,nr] .-GM_PsiX.f[:,nr]./Γ.DRF[nr] + bolusV.f[:, nr] .= 0*GM_PsiY.f[:,nr] .-GM_PsiY.f[:,nr]./Γ.DRF[nr] + + bolusU=bolusU.*mskW + bolusV=bolusV.*mskS #and its vertical part # (seems correct, leading to 0 divergence) @@ -721,7 +721,7 @@ function calc_bolus(GM_PsiX,GM_PsiY, Γ) end end - bolusW=tmp_w.*mskC; + bolusW=tmp_w.*mskC return bolusU, bolusV, bolusW end diff --git a/src/Polygons.jl b/src/Polygons.jl index a433771..3101406 100644 --- a/src/Polygons.jl +++ b/src/Polygons.jl @@ -7,9 +7,9 @@ using Glob ## functions to derive polygons from XG,YG read from native grid function to_sphere(XC,YC) - x=sin.(pi/2 .-YC*pi/180).*cos.(XC*pi/180); - y=sin.(pi/2 .-YC*pi/180).*sin.(XC*pi/180); - z=cos.(pi/2 .-YC*pi/180); + x=sin.(pi/2 .-YC*pi/180).*cos.(XC*pi/180) + y=sin.(pi/2 .-YC*pi/180).*sin.(XC*pi/180) + z=cos.(pi/2 .-YC*pi/180) (x,y,z) end diff --git a/src/ReIndexing.jl b/src/ReIndexing.jl index 04e9a0c..5da3306 100644 --- a/src/ReIndexing.jl +++ b/src/ReIndexing.jl @@ -264,7 +264,7 @@ function NeighborTileIndices_cs(grid::NamedTuple) 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); + 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 diff --git a/src/ReadWrite.jl b/src/ReadWrite.jl index 1779fb4..f53ead5 100644 --- a/src/ReadWrite.jl +++ b/src/ReadWrite.jl @@ -121,10 +121,10 @@ function read_one!(xx::Array,x::AbstractMeshArray; verbose=false) facesSize=x.grid.fSize (n1,n2)=x.grid.ioSize (nFaces,n3,n4)=nFacesEtc(x) - i0=1; i1=0; - j0=1; j1=0; + i0=1; i1=0 + j0=1; j1=0 for iFace=1:nFaces - nn=facesSize[iFace][1]; mm=facesSize[iFace][2]; + nn=facesSize[iFace][1]; mm=facesSize[iFace][2] if format==:compact i0=i1+1 i1=i1+nn*mm @@ -170,11 +170,11 @@ function write(x::AbstractMeshArray; verbose=false) y = Array{eltype(x),4}(undef,(n1,n2,n3,n4)) end - i0=1; i1=0; - j0=1; j1=0; - for iFace=1:nFaces; - nn=facesSize[iFace][1]; - mm=facesSize[iFace][2]; + i0=1; i1=0 + j0=1; j1=0 + for iFace=1:nFaces + nn=facesSize[iFace][1] + mm=facesSize[iFace][2] if format==:compact i0=i1+1 i1=i1+nn*mm diff --git a/src/Solvers.jl b/src/Solvers.jl index 8d9a33e..678ae12 100644 --- a/src/Solvers.jl +++ b/src/Solvers.jl @@ -20,14 +20,14 @@ function MaskWetPoints(TrspCon) tmp1=fill(1.0,mskWet) tmp2=exchange(tmp1).MA for I=1:size(tmp1.f,1) - tmp3=mskWet[I]; tmp4=tmp2[I]; - tmp4=tmp4[2:end-1,1:end-2]+tmp4[2:end-1,3:end]+tmp4[1:end-2,2:end-1]+tmp4[3:end,2:end-1]; + tmp3=mskWet[I]; tmp4=tmp2[I] + tmp4=tmp4[2:end-1,1:end-2]+tmp4[2:end-1,3:end]+tmp4[1:end-2,2:end-1]+tmp4[3:end,2:end-1] !isempty(findall(isnan.(tmp4) .& (!isnan).(tmp3))) ? println("warning: modified mask") : nothing tmp3[findall(isnan.(tmp4))] .= NaN mskWet[I]=tmp3 end # - TrspCon=mask(TrspCon,0.0)*mskWet; + TrspCon=mask(TrspCon,0.0)*mskWet return TrspCon, mskWet, mskDry end @@ -42,7 +42,7 @@ Mapping from global array to global ocean vector. function MapWetPoints(mskWet) tmp1=write(mskWet)[:] kk=findall((!isnan).(tmp1)) - nn=length(kk); s0=size(tmp1); s1=mskWet.grid.ioSize; + nn=length(kk); s0=size(tmp1); s1=mskWet.grid.ioSize Kvec=fill(0.0,s0...); Kvec[kk]=kk; Kmap=read(reshape(Kvec,s1...),mskWet) #global array indices Lvec=fill(0.0,s0...); Lvec[kk]=1:nn; Lmap=read(reshape(Lvec,s1...),mskWet) #global vector indices return Kvec,Lvec,Kmap,Lmap @@ -100,28 +100,28 @@ function MatrixForPoisson(TrspCon,mskWet,mskDry,Kvec,Lvec,Kmap,Lmap) J=Array{Int}(undef,0) V=Array{Float64}(undef,0) - for aa=1:TrspCon.grid.nFaces; - for ii=1:3; for jj=1:3; + for aa=1:TrspCon.grid.nFaces + for ii=1:3; for jj=1:3 #1) compute effect of each point on neighboring target point: - (FLDones,FLDkkFROM)=SeedWetPoints(TrspCon,Kmap,Lmap,aa,ii,jj); + (FLDones,FLDkkFROM)=SeedWetPoints(TrspCon,Kmap,Lmap,aa,ii,jj) (tmpU,tmpV)=gradient(FLDones,(;),false) - dFLDdt=convergence(tmpU,tmpV); + dFLDdt=convergence(tmpU,tmpV) #2) mask `dFLDdt` since we use a **Neumann** boundary condition. #Extrapolation uses a **Dirichlet** boundary condition, so mskFreeze should not be applied then. - isa(FLDkkFROM,MeshArray) ? FLDkkFROM=write(FLDkkFROM)[:] : nothing; + isa(FLDkkFROM,MeshArray) ? FLDkkFROM=write(FLDkkFROM)[:] : nothing #3.1) For wet points -- add contributions in matrix: - dFLDdtWet=write(dFLDdt.*mskWet)[:]; - tmp1=findall( (dFLDdtWet .!= 0.0) .* (!isnan).(dFLDdtWet)); - tmpV=dFLDdtWet[tmp1]; tmpJ=FLDkkFROM[tmp1]; tmpI=Kvec[tmp1]; - I=[I;Lvec[Int.(tmpI)]]; J=[J;Lvec[Int.(tmpJ)]]; V=[V;tmpV]; + dFLDdtWet=write(dFLDdt.*mskWet)[:] + tmp1=findall( (dFLDdtWet .!= 0.0) .* (!isnan).(dFLDdtWet)) + tmpV=dFLDdtWet[tmp1]; tmpJ=FLDkkFROM[tmp1]; tmpI=Kvec[tmp1] + I=[I;Lvec[Int.(tmpI)]]; J=[J;Lvec[Int.(tmpJ)]]; V=[V;tmpV] size(tmpV) #3.2) For dry points -- This part reflects the `Neumann` boundary condition: - dFLDdtDry=write(dFLDdt.*mskDry)[:]; - tmp1=findall( (dFLDdtDry .!= 0.0) .* (!isnan).(dFLDdtDry)); - tmpV=dFLDdtDry[tmp1]; tmpIJ=FLDkkFROM[tmp1]; - I=[I;Lvec[Int.(tmpIJ)]]; J=[J;Lvec[Int.(tmpIJ)]]; V=[V;tmpV]; - end; end; - end; + dFLDdtDry=write(dFLDdt.*mskDry)[:] + tmp1=findall( (dFLDdtDry .!= 0.0) .* (!isnan).(dFLDdtDry)) + tmpV=dFLDdtDry[tmp1]; tmpIJ=FLDkkFROM[tmp1] + I=[I;Lvec[Int.(tmpIJ)]]; J=[J;Lvec[Int.(tmpIJ)]]; V=[V;tmpV] + end; end + end nn=sum((!isnan).(mskWet)) A=sparse(I,J,V,nn,nn) @@ -165,8 +165,8 @@ function VectorPotential(TrspX::AbstractMeshArray,TrspY::AbstractMeshArray,Γ::N # 1) streamfunction face by face: - (fldU,fldV)=exch_UV(TrspX,TrspY); - fldU=mask(fldU,0.0); fldV=mask(fldV,0.0); + (fldU,fldV)=exch_UV(TrspX,TrspY) + fldU=mask(fldU,0.0); fldV=mask(fldV,0.0) psi=similar(fldV) for I in eachindex(fldV) @@ -187,7 +187,7 @@ function VectorPotential(TrspX::AbstractMeshArray,TrspY::AbstractMeshArray,Γ::N tmp3b=[tmp2;ones(1,size(tmp2,2))] tmpB=tmp3a.*tmp3b # - tmpA=psi[iF]; + tmpA=psi[iF] I=findall(tmpB.==0) ii=I[1][1]; jj=I[1][2] tmpA[:,jj]=tmpA[:,jj] .- tmpB[ii,jj] @@ -204,7 +204,7 @@ function VectorPotential(TrspX::AbstractMeshArray,TrspY::AbstractMeshArray,Γ::N tmpE+=fldU[iF][:,kk-1] tmpE=median(tmpE[findall((!isnan).(tmpE))]) tmpA[:,kk]=tmpA[:,kk] .- tmpE - end; + end # psi[iF]=tmpA end @@ -264,7 +264,7 @@ function VectorPotential(TrspX::AbstractMeshArray,TrspY::AbstractMeshArray,Γ::N tmp2=tmp2[findall(tmp_dis .== minimum(tmp_dis))] #set that point to zero - psi=psi .- tmp1[tmp2[1]]; + psi=psi .- tmp1[tmp2[1]] return psi end diff --git a/src/demo.jl b/src/demo.jl index d7ba311..9ad2cc0 100644 --- a/src/demo.jl +++ b/src/demo.jl @@ -20,29 +20,29 @@ module demo latPairs=[] namPairs=[] - push!(lonPairs,[-173 -164]); push!(latPairs,[65.5 65.5]); push!(namPairs,"Bering Strait"); - push!(lonPairs,[-5 -5]); push!(latPairs,[34 40]); push!(namPairs,"Gibraltar"); - push!(lonPairs,[-81 -77]); push!(latPairs,[28 26]); push!(namPairs,"Florida Strait"); - push!(lonPairs,[-81 -79]); push!(latPairs,[28 22]); push!(namPairs,"Florida Strait W1"); - push!(lonPairs,[-76 -76]); push!(latPairs,[21 8]); push!(namPairs,"Florida Strait S1"); - push!(lonPairs,[-77 -77]); push!(latPairs,[26 24]); push!(namPairs,"Florida Strait E1"); - push!(lonPairs,[-77 -77]); push!(latPairs,[24 22]); push!(namPairs,"Florida Strait E2"); - push!(lonPairs,[-65 -50]); push!(latPairs,[66 66]); push!(namPairs,"Davis Strait"); - push!(lonPairs,[-35 -20]); push!(latPairs,[67 65]); push!(namPairs,"Denmark Strait"); - push!(lonPairs,[-16 -7]); push!(latPairs,[65 62.5]); push!(namPairs,"Iceland Faroe"); - push!(lonPairs,[-6.5 -4]); push!(latPairs,[62.5 57]); push!(namPairs,"Faroe Scotland"); - push!(lonPairs,[-4 8]); push!(latPairs,[57 62]); push!(namPairs,"Scotland Norway"); - push!(lonPairs,[-68 -63]); push!(latPairs,[-54 -66]); push!(namPairs,"Drake Passage"); - push!(lonPairs,[103 103]); push!(latPairs,[4 -1]); push!(namPairs,"Indonesia W1"); - push!(lonPairs,[104 109]); push!(latPairs,[-3 -8]); push!(namPairs,"Indonesia W2"); - push!(lonPairs,[113 118]); push!(latPairs,[-8.5 -8.5]); push!(namPairs,"Indonesia W3"); - push!(lonPairs,[118 127 ]); push!(latPairs,[-8.5 -15]); push!(namPairs,"Indonesia W4"); - push!(lonPairs,[127 127]); push!(latPairs,[-25 -68]); push!(namPairs,"Australia Antarctica"); - push!(lonPairs,[38 46]); push!(latPairs,[-10 -22]); push!(namPairs,"Madagascar Channel"); - push!(lonPairs,[46 46]); push!(latPairs,[-22 -69]); push!(namPairs,"Madagascar Antarctica"); - push!(lonPairs,[20 20]); push!(latPairs,[-30 -69.5]); push!(namPairs,"South Africa Antarctica"); - push!(lonPairs,[-76 -72]); push!(latPairs,[21 18.5]); push!(namPairs,"Florida Strait E3"); - push!(lonPairs,[-72 -72]); push!(latPairs,[18.5 10]); push!(namPairs,"Florida Strait E4"); + push!(lonPairs,[-173 -164]); push!(latPairs,[65.5 65.5]); push!(namPairs,"Bering Strait") + push!(lonPairs,[-5 -5]); push!(latPairs,[34 40]); push!(namPairs,"Gibraltar") + push!(lonPairs,[-81 -77]); push!(latPairs,[28 26]); push!(namPairs,"Florida Strait") + push!(lonPairs,[-81 -79]); push!(latPairs,[28 22]); push!(namPairs,"Florida Strait W1") + push!(lonPairs,[-76 -76]); push!(latPairs,[21 8]); push!(namPairs,"Florida Strait S1") + push!(lonPairs,[-77 -77]); push!(latPairs,[26 24]); push!(namPairs,"Florida Strait E1") + push!(lonPairs,[-77 -77]); push!(latPairs,[24 22]); push!(namPairs,"Florida Strait E2") + push!(lonPairs,[-65 -50]); push!(latPairs,[66 66]); push!(namPairs,"Davis Strait") + push!(lonPairs,[-35 -20]); push!(latPairs,[67 65]); push!(namPairs,"Denmark Strait") + push!(lonPairs,[-16 -7]); push!(latPairs,[65 62.5]); push!(namPairs,"Iceland Faroe") + push!(lonPairs,[-6.5 -4]); push!(latPairs,[62.5 57]); push!(namPairs,"Faroe Scotland") + push!(lonPairs,[-4 8]); push!(latPairs,[57 62]); push!(namPairs,"Scotland Norway") + push!(lonPairs,[-68 -63]); push!(latPairs,[-54 -66]); push!(namPairs,"Drake Passage") + push!(lonPairs,[103 103]); push!(latPairs,[4 -1]); push!(namPairs,"Indonesia W1") + push!(lonPairs,[104 109]); push!(latPairs,[-3 -8]); push!(namPairs,"Indonesia W2") + push!(lonPairs,[113 118]); push!(latPairs,[-8.5 -8.5]); push!(namPairs,"Indonesia W3") + push!(lonPairs,[118 127 ]); push!(latPairs,[-8.5 -15]); push!(namPairs,"Indonesia W4") + push!(lonPairs,[127 127]); push!(latPairs,[-25 -68]); push!(namPairs,"Australia Antarctica") + push!(lonPairs,[38 46]); push!(latPairs,[-10 -22]); push!(namPairs,"Madagascar Channel") + push!(lonPairs,[46 46]); push!(latPairs,[-22 -69]); push!(namPairs,"Madagascar Antarctica") + push!(lonPairs,[20 20]); push!(latPairs,[-30 -69.5]); push!(namPairs,"South Africa Antarctica") + push!(lonPairs,[-76 -72]); push!(latPairs,[21 18.5]); push!(namPairs,"Florida Strait E3") + push!(lonPairs,[-72 -72]); push!(latPairs,[18.5 10]); push!(namPairs,"Florida Strait E4") (name=namPairs,lon=lonPairs,lat=latPairs) end @@ -194,8 +194,8 @@ module demo """ function get_basemap() dx=0.1 - lat=[j for i=-0.05:dx:359.95, j=-89.95:dx:89.95]; - lon=[i for i=-0.05:dx:359.95, j=-89.95:dx:89.95]; + lat=[j for i=-0.05:dx:359.95, j=-89.95:dx:89.95] + lon=[i for i=-0.05:dx:359.95, j=-89.95:dx:89.95] earth_jpg=MeshArrays.Dataset("basemap_jpg1") earth_img=read_jld2(earth_jpg) diff --git a/src/exchanges/CubeSphere.jl b/src/exchanges/CubeSphere.jl index 59f4e02..f366ee7 100644 --- a/src/exchanges/CubeSphere.jl +++ b/src/exchanges/CubeSphere.jl @@ -15,15 +15,15 @@ nf==5 ? s=vcat(s,s[3]) : nothing tp=fld.grid.class FLD=similar(fld;m=fld.meta) -for i=1:nf; FLD.f[i]=fill(fillval,s[i].+2N); end; +for i=1:nf; FLD.f[i]=fill(fillval,s[i].+2N); end #code below yields strange, seemingly incorrect results: -#for i=1:nf; FLD.f[i]=Array{eltype(fld.f[i])}(undef,s[i].+2N); end; +#for i=1:nf; FLD.f[i]=Array{eltype(fld.f[i])}(undef,s[i].+2N); end #all versions below yield same @time and memory (despite diff in allocs) -for i=1:nf; -# FLD.f[i][N+1:end-N,N+1:end-N]=fld.f[i]; - @views FLD.f[i][N+1:N+s[i][1],N+1:N+s[i][2]]=fld.f[i]; -end; +for i=1:nf +# FLD.f[i][N+1:end-N,N+1:end-N]=fld.f[i] + @views FLD.f[i][N+1:N+s[i][1],N+1:N+s[i][2]]=fld.f[i] +end #step 2 @@ -62,12 +62,12 @@ tp=fldU.grid.class FLDU=similar(fldU;m=fldU.meta) FLDV=similar(fldV;m=fldV.meta) -for i=1:nf; - FLDU.f[i]=fill(fillval,s[i].+2N); - FLDV.f[i]=fill(fillval,s[i].+2N); - @views FLDU.f[i][N+1:N+s[i][1],N+1:N+s[i][2]]=fldU.f[i]; - @views FLDV.f[i][N+1:N+s[i][1],N+1:N+s[i][2]]=fldV.f[i]; -end; +for i=1:nf + FLDU.f[i]=fill(fillval,s[i].+2N) + FLDV.f[i]=fill(fillval,s[i].+2N) + @views FLDU.f[i][N+1:N+s[i][1],N+1:N+s[i][2]]=fldU.f[i] + @views FLDV.f[i][N+1:N+s[i][1],N+1:N+s[i][2]]=fldV.f[i] +end #step 2 @@ -115,10 +115,10 @@ FLDU=similar(fldU;m=fldU.meta) FLDV=similar(fldV;m=fldV.meta) for i=1:nf - FLDU.f[i]=fill(fillval,s[i][1]+1,s[i][2]); - FLDV.f[i]=fill(fillval,s[i][1],s[i][2]+1); - @views FLDU.f[i][1:s[i][1],1:s[i][2]]=fldU.f[i]; - @views FLDV.f[i][1:s[i][1],1:s[i][2]]=fldV.f[i]; + FLDU.f[i]=fill(fillval,s[i][1]+1,s[i][2]) + FLDV.f[i]=fill(fillval,s[i][1],s[i][2]+1) + @views FLDU.f[i][1:s[i][1],1:s[i][2]]=fldU.f[i] + @views FLDV.f[i][1:s[i][1],1:s[i][2]]=fldV.f[i] end #step 2 @@ -144,10 +144,10 @@ end function exch_cs_target(sa::Tuple{Int64,Int64},N::Integer) #target array indices - jW=(1:N,N+1:N+sa[2]); - jE=(N+1+sa[1]:2N+sa[1],N+1:N+sa[2]); - jS=(N+1:N+sa[1],1:N); - jN=(N+1:N+sa[1],N+1+sa[2]:2N+sa[2]); + jW=(1:N,N+1:N+sa[2]) + jE=(N+1+sa[1]:2N+sa[1],N+1:N+sa[2]) + jS=(N+1:N+sa[1],1:N) + jN=(N+1:N+sa[1],N+1+sa[2]:2N+sa[2]) return jW, jE, jS, jN @@ -156,28 +156,28 @@ end function exch_cs_sources(a::Integer,s::Array{Tuple{Int64,Int64},1},N::Integer) #source array IDs -aW=0; aE=0; aS=0; aN=0; -if a==1; aW=5; aE=2; aS=6; aN=3; -elseif a==2; aW=1; aE=4; aS=6; aN=3; -elseif a==3; aW=1; aE=4; aS=2; aN=5; -elseif a==4; aW=3; aE=6; aS=2; aN=5; -elseif a==5; aW=3; aE=6; aS=4; aN=1; -elseif a==6; aW=5; aE=2; aS=4; aN=1; -else; error("Array index is out of bounds."); -end; +aW=0; aE=0; aS=0; aN=0 +if a==1; aW=5; aE=2; aS=6; aN=3 +elseif a==2; aW=1; aE=4; aS=6; aN=3 +elseif a==3; aW=1; aE=4; aS=2; aN=5 +elseif a==4; aW=3; aE=6; aS=2; aN=5 +elseif a==5; aW=3; aE=6; aS=4; aN=1 +elseif a==6; aW=5; aE=2; aS=4; aN=1 +else; error("Array index is out of bounds.") +end if !iseven(a) #source array indices - iW=(1:s[aW][1],s[aW][2]-N+1:s[aW][2]); - iE=(1:N,1:s[aE][2]); - iS=(1:s[aS][1],s[aS][2]-N+1:s[aS][2]); - iN=(1:N,1:s[aN][2]); + iW=(1:s[aW][1],s[aW][2]-N+1:s[aW][2]) + iE=(1:N,1:s[aE][2]) + iS=(1:s[aS][1],s[aS][2]-N+1:s[aS][2]) + iN=(1:N,1:s[aN][2]) else #source array indices - iW=(s[aW][1]-N+1:s[aW][1],1:s[aW][2]); - iE=(1:s[aE][1],1:N); - iS=(s[aS][1]-N+1:s[aS][1],1:s[aS][2]); - iN=(1:s[aN][1],1:N); + iW=(s[aW][1]-N+1:s[aW][1],1:s[aW][2]) + iE=(1:s[aE][1],1:N) + iS=(s[aS][1]-N+1:s[aS][1],1:s[aS][2]) + iN=(1:s[aN][1],1:N) end return aW,aE,aS,aN,iW,iE,iS,iN diff --git a/src/exchanges/PeriodicChannel.jl b/src/exchanges/PeriodicChannel.jl index 69bae65..5840498 100644 --- a/src/exchanges/PeriodicChannel.jl +++ b/src/exchanges/PeriodicChannel.jl @@ -6,17 +6,17 @@ fillval=0.0 #step 1 -s=size.(fld.f); -FLD=similar(fld;m=fld.meta); -FLD.f[1]=fill(fillval,s[1].+2N); -@views FLD.f[1][N+1:N+s[1][1],N+1:N+s[1][2]]=fld.f[1]; +s=size.(fld.f) +FLD=similar(fld;m=fld.meta) +FLD.f[1]=fill(fillval,s[1].+2N) +@views FLD.f[1][N+1:N+s[1][1],N+1:N+s[1][2]]=fld.f[1] #step 2 -iW=(s[1][1]-N+1:s[1][1],1:s[1][2]); -iE=(1:N,1:s[1][2]); -jW=(1:N,N+1:N+s[1][2]); -jE=(N+1+s[1][1]:2N+s[1][1],N+1:N+s[1][2]); +iW=(s[1][1]-N+1:s[1][1],1:s[1][2]) +iE=(1:N,1:s[1][2]) +jW=(1:N,N+1:N+s[1][2]) +jE=(N+1+s[1][1]:2N+s[1][1],N+1:N+s[1][2]) FLD.f[1][jW[1],jW[2]]=view(fld.f[1],iW[1],iW[2]) FLD.f[1][jE[1],jE[2]]=view(fld.f[1],iE[1],iE[2]) @@ -26,10 +26,10 @@ end ## -function exch_UV_N_PeriodicChannel(fldU,fldV,N); +function exch_UV_N_PeriodicChannel(fldU,fldV,N) -FLDU=exch_T_N_PeriodicChannel(fldU,N); -FLDV=exch_T_N_PeriodicChannel(fldV,N); +FLDU=exch_T_N_PeriodicChannel(fldU,N) +FLDV=exch_T_N_PeriodicChannel(fldV,N) return FLDU,FLDV @@ -37,20 +37,20 @@ end ## -function exch_UV_PeriodicChannel(fldU,fldV); +function exch_UV_PeriodicChannel(fldU,fldV) fillval=0.0 #step 1 -s=size.(fldU.f); -FLDU=similar(fldU;m=fldU.meta); -FLDV=similar(fldV;m=fldV.meta); +s=size.(fldU.f) +FLDU=similar(fldU;m=fldU.meta) +FLDV=similar(fldV;m=fldV.meta) -FLDU.f[1]=fill(fillval,s[1][1]+1,s[1][2]); -FLDV.f[1]=fill(fillval,s[1][1],s[1][2]+1); -@views FLDU.f[1][1:s[1][1],1:s[1][2]]=fldU.f[1]; -@views FLDV.f[1][1:s[1][1],1:s[1][2]]=fldV.f[1]; +FLDU.f[1]=fill(fillval,s[1][1]+1,s[1][2]) +FLDV.f[1]=fill(fillval,s[1][1],s[1][2]+1) +@views FLDU.f[1][1:s[1][1],1:s[1][2]]=fldU.f[1] +@views FLDV.f[1][1:s[1][1],1:s[1][2]]=fldV.f[1] #step 2 diff --git a/src/exchanges/PeriodicDomain.jl b/src/exchanges/PeriodicDomain.jl index 3779a2f..14bc8d1 100644 --- a/src/exchanges/PeriodicDomain.jl +++ b/src/exchanges/PeriodicDomain.jl @@ -7,7 +7,7 @@ fillval=0.0 ni,nj=Int.(fld.grid.ioSize[:]./fld.grid.fSize[1][:]) s=fld.fSize -FLD=similar(fld;m=fld.meta); +FLD=similar(fld;m=fld.meta) for i=1:ni for j=1:nj @@ -19,24 +19,24 @@ for i=1:ni #step 1 - FLD.f[k]=fill(fillval,s[k].+2N); - @views FLD.f[k][N+1:N+s[k][1],N+1:N+s[k][2]]=fld.f[k]; + FLD.f[k]=fill(fillval,s[k].+2N) + @views FLD.f[k][N+1:N+s[k][1],N+1:N+s[k][2]]=fld.f[k] #step 2 - iW=(s[k][1]-N+1:s[k][1],1:s[k][2]); - iE=(1:N,1:s[k][2]); - jW=(1:N,N+1:N+s[k][2]); - jE=(N+1+s[k][1]:2N+s[k][1],N+1:N+s[k][2]); + iW=(s[k][1]-N+1:s[k][1],1:s[k][2]) + iE=(1:N,1:s[k][2]) + jW=(1:N,N+1:N+s[k][2]) + jE=(N+1+s[k][1]:2N+s[k][1],N+1:N+s[k][2]) FLD.f[k][jW[1],jW[2]]=view(fld.f[kW],iW[1],iW[2]) FLD.f[k][jE[1],jE[2]]=view(fld.f[kE],iE[1],iE[2]) #step 3 - iS=(1:s[k][1],s[k][2]-N+1:s[k][2]); - iN=(1:s[k][1],1:N); - jS=(N+1:N+s[k][1],1:N); - jN=(N+1:N+s[k][1],N+1+s[k][2]:2N+s[k][2]); + iS=(1:s[k][1],s[k][2]-N+1:s[k][2]) + iN=(1:s[k][1],1:N) + jS=(N+1:N+s[k][1],1:N) + jN=(N+1:N+s[k][1],N+1+s[k][2]:2N+s[k][2]) FLD.f[k][jS[1],jS[2]]=view(fld.f[kS],iS[1],iS[2]) FLD.f[k][jN[1],jN[2]]=view(fld.f[kN],iN[1],iN[2]) @@ -53,7 +53,7 @@ function exch_UV_N_PeriodicDomain(fldU,fldV,N) return FLDU,FLDV end -function exch_UV_PeriodicDomain(fldU,fldV); +function exch_UV_PeriodicDomain(fldU,fldV) fillval=0.0 @@ -72,10 +72,10 @@ for i=1:ni #step 1 - FLDU.f[k]=fill(fillval,s[k][1]+1,s[k][2]); - FLDV.f[k]=fill(fillval,s[k][1],s[k][2]+1); - @views FLDU.f[k][1:s[k][1],1:s[k][2]]=fldU.f[k]; - @views FLDV.f[k][1:s[k][1],1:s[k][2]]=fldV.f[k]; + FLDU.f[k]=fill(fillval,s[k][1]+1,s[k][2]) + FLDV.f[k]=fill(fillval,s[k][1],s[k][2]+1) + @views FLDU.f[k][1:s[k][1],1:s[k][2]]=fldU.f[k] + @views FLDV.f[k][1:s[k][1],1:s[k][2]]=fldV.f[k] #step 2 diff --git a/src/grids/NEMO.jl b/src/grids/NEMO.jl index eebd372..69468fe 100644 --- a/src/grids/NEMO.jl +++ b/src/grids/NEMO.jl @@ -100,7 +100,7 @@ variable_in_NEMO(v,vl=variable_list_2d)= #DataFrame version: # filter(p->p.:MITgcm==v,variable_df(variable_list))[1,:NEMO] -function convert_one_grid_variable(grid_data,df_line; +function convert_one_grid_variable(grid_data,df_line; is_3d=false, verbose=false) nam_in=df_line.NEMO nam_out=df_line.MITgcm diff --git a/src/grids/main.jl b/src/grids/main.jl index abd6854..4701e26 100644 --- a/src/grids/main.jl +++ b/src/grids/main.jl @@ -14,7 +14,7 @@ using MeshArrays γ=GridSpec() ``` """ -function GridSpec(category="default", path=tempname(); +function GridSpec(category="default", path=tempname(); np=nothing, ID=:unknown, ioPrec=Float64, verbose=false) if category=="default"&&in(ID,[:unknown, :OISST, :Oscar, :IAP]) verbose ? println("branch 1") : nothing @@ -47,7 +47,7 @@ b = MeshArrays.GridSpec_default(MeshArrays.Grids_simple.xy_OISST(),tile=(1440,10 c = MeshArrays.GridSpec_default(MeshArrays.Grids_simple.xy_OISST(),tile=(720,510)) ``` """ -function GridSpec_default(xy=NamedTuple(); +function GridSpec_default(xy=NamedTuple(); ID=:unknown, ioPrec=Float64, tile=[]) if !isempty(xy) GridSpec_default_xy(xy,ioPrec=ioPrec,tile=tile) @@ -146,7 +146,7 @@ function GridLoad_main(γ=GridSpec(); option=:minimal) if op=="full" list_n=("XC","XG","YC","YG","RAC","RAW","RAS","RAZ","DXC","DXG","DYC","DYG","Depth") if (!isempty(filter(x -> occursin("AngleCS",x), readdir(γ.path)))) - list_n=(list_n...,"AngleCS","AngleSN"); + list_n=(list_n...,"AngleCS","AngleSN") end list_n=(list_n...,"DRC","DRF","RC","RF") list_n=(list_n...,"hFacC","hFacS","hFacW") @@ -196,12 +196,12 @@ true ``` """ function GridLoadVar(nam::String,γ::gcmgrid) - pc=fill(0.5,2); pg=fill(0.0,2); pu=[0.,0.5]; pv=[0.5,0.]; + pc=fill(0.5,2); pg=fill(0.0,2); pu=[0.,0.5]; pv=[0.5,0.] list_n=("XC","XG","YC","YG","RAC","RAW","RAS","RAZ","DXC","DXG","DYC","DYG","Depth","AngleCS","AngleSN") list_u=(u"°",u"°",u"°",u"°",u"m^2",u"m^2",u"m^2",u"m^2",u"m",u"m",u"m",u"m",u"m",1.0,1.0) list_p=(pc,pg,pc,pg,pc,pu,pv,pg,pu,pv,pv,pu,pc,pc,pc) # - list3d_n=("hFacC","hFacS","hFacW"); + list3d_n=("hFacC","hFacS","hFacW") list3d_u=(1.0,1.0,1.0) list3d_p=(fill(0.5,3),[0.,0.5,0.5],[0.5,0.,0.5]) # @@ -226,7 +226,7 @@ function GridLoadVar(nam::String,γ::gcmgrid) n3=Int64(stat(fil).size/reclen) ii=findall(nam.==list3d_n)[1] - m=varmeta(list3d_u[ii],list3d_p[ii],missing,list3d_n[ii],list3d_n[ii]); + m=varmeta(list3d_u[ii],list3d_p[ii],missing,list3d_n[ii],list3d_n[ii]) tmp1=γ.read(joinpath(γ.path,list3d_n[ii]*".data"),MeshArray(γ,γ.ioPrec,n3;meta=m)) else tmp1=missing @@ -269,7 +269,7 @@ function GridAddWS!(Γ::Dict) # tmp1=YC[ff][1:end-2,2:end-1] tmp2=YC[ff][2:end-1,2:end-1] - YW[ff]=(tmp1.+tmp2)./2; + YW[ff]=(tmp1.+tmp2)./2 # tmp1=YC[ff][2:end-1,1:end-2] tmp2=YC[ff][2:end-1,2:end-1] diff --git a/src/grids/simple.jl b/src/grids/simple.jl index 97f1ca2..4a54a74 100644 --- a/src/grids/simple.jl +++ b/src/grids/simple.jl @@ -55,7 +55,7 @@ as specified by `γ` obtained from `GridSpec_ones`. ``` using MeshArrays -γ=MeshArrays.GridSpec_ones("CubeSphere",6,20); +γ=MeshArrays.GridSpec_ones("CubeSphere",6,20) Γ=MeshArrays.GridLoad_ones(γ;option="full") ``` """ @@ -65,25 +65,25 @@ function GridLoad_ones(γ::gcmgrid;option="minimal") Γ=Dict() - pc=fill(0.5,2); pg=fill(0.0,2); pu=[0.,0.5]; pv=[0.5,0.]; + pc=fill(0.5,2); pg=fill(0.0,2); pu=[0.,0.5]; pv=[0.5,0.] if option=="full" - list_n=("XC","XG","YC","YG","RAC","RAW","RAS","RAZ","DXC","DXG","DYC","DYG","Depth","hFacC","hFacS","hFacW"); + list_n=("XC","XG","YC","YG","RAC","RAW","RAS","RAZ","DXC","DXG","DYC","DYG","Depth","hFacC","hFacS","hFacW") list_u=(u"m",u"m",u"m",u"m",u"m^2",u"m^2",u"m^2",u"m^2",u"m",u"m",u"m",u"m",u"m",1.0,1.0,1.0) list_p=(pc,pg,pc,pg,pc,pu,pv,pg,pu,pv,pv,pu,pc,fill(0.5,3),[0.,0.5,0.5],[0.5,0.,0.5]) elseif option=="light" - list_n=("XC","XG","YC","YG","RAC","DXC","DXG","DYC","DYG","Depth"); + list_n=("XC","XG","YC","YG","RAC","DXC","DXG","DYC","DYG","Depth") list_u=(u"m",u"m",u"m",u"m",u"m^2",u"m",u"m",u"m",u"m",u"m") list_p=(pc,pg,pc,pg,pc,pu,pv,pv,pu,pc) else - list_n=("XC","YC"); + list_n=("XC","YC") list_u=(u"°",u"°") list_p=(pc,pc) end - for ii=1:length(list_n); + for ii=1:length(list_n) tmp1=fill(1.,(ioSize[:])) - m=varmeta(list_u[ii],list_p[ii],missing,list_n[ii],list_n[ii]); - tmp1=γ.read(tmp1,MeshArray(γ,Float64;meta=m)); + m=varmeta(list_u[ii],list_p[ii],missing,list_n[ii],list_n[ii]) + tmp1=γ.read(tmp1,MeshArray(γ,Float64;meta=m)) Γ[list_n[ii]]=tmp1 end diff --git a/src/types/gcmarray.jl b/src/types/gcmarray.jl index c0e1fc8..86836cb 100644 --- a/src/types/gcmarray.jl +++ b/src/types/gcmarray.jl @@ -40,9 +40,9 @@ function gcmarray(A::Array{T,N}; ios=[s[1] s[2]] γ=gcmgrid("","PeriodicDomain",1, fs, ios, T, read, write) if N==2 - B=γ.read(A,gcmarray(γ,T;meta=meta)); + B=γ.read(A,gcmarray(γ,T;meta=meta)) else - B=γ.read(A,gcmarray(γ,T,size(A)[3];meta=meta)); + B=γ.read(A,gcmarray(γ,T,size(A)[3];meta=meta)) end B end @@ -56,7 +56,7 @@ function gcmarray(grid::gcmgrid,f::OuterArray{InnerArray{T,N},1}; meta::varmeta=defaultmeta) where {T, N} nFaces=grid.nFaces if N>2 - n3=size(f[1],3); n4=size(f[1],4); + n3=size(f[1],3); n4=size(f[1],4) g=OuterArray{InnerArray{T,2},3}(undef,nFaces,n3,n4) for I in eachindex(view(g,1:nFaces,1:n3,1:n4)) g[I]=view(f[I[1]],:,:,I[2],I[3]) @@ -90,9 +90,9 @@ function gcmarray(grid::gcmgrid,::Type{T}, f=OuterArray{InnerArray{T,2},2}(undef,nFaces,n3) isa(fSize,NTuple) ? fSize=[fSize] : nothing isa(fIndex,Int) ? fIndex=[fIndex] : nothing - for a=1:nFaces; for i3=1:n3; + for a=1:nFaces; for i3=1:n3 f[a,i3]=InnerArray{T,2}(undef,fSize[a]...) - end; end; + end; end gcmarray{T,2,InnerArray{T,2}}(grid,meta,f,fSize,fIndex,thisversion) end @@ -104,9 +104,9 @@ function gcmarray(grid::gcmgrid,::Type{T}, f=OuterArray{InnerArray{T,2},3}(undef,nFaces,n3,n4) isa(fSize,NTuple) ? fSize=[fSize] : nothing isa(fIndex,Int) ? fIndex=[fIndex] : nothing - for a=1:nFaces; for i4=1:n4; for i3=1:n3; + for a=1:nFaces; for i4=1:n4; for i3=1:n3 f[a,i3,i4]=InnerArray{T,2}(undef,fSize[a]...) - end; end; end; + end; end; end gcmarray{T,3,InnerArray{T,2}}(grid,meta,f,fSize,fIndex,thisversion) end diff --git a/src/types/gcmfaces.jl b/src/types/gcmfaces.jl index 21be41b..9bd3eaa 100644 --- a/src/types/gcmfaces.jl +++ b/src/types/gcmfaces.jl @@ -156,11 +156,11 @@ function fijind(A::gcmfaces,ij::Int) tmp11=tmpsize[1]*tmpsize[2] tmp2=tmp1+tmp11 if tmp1Tx,"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)] - plot(LC) - plot(LC[1]) - - x=zeros(γ) - fill!(x,1.0) - y=fill(-1.0,γ) - extrema(y) - @test minimum(y)0) - MeshArrays.nFacesEtc(tmp) - @suppress show(tmp) - - x=tmp[1:10,1,1:2]; y=x[2]; x[3]=1.0 - view(x,1:3,:,1) - MeshArrays.gcmsubset(γ,x.f,x.fSize,x.aSize,x.i,x.iSize) - MeshArrays.fsize(x.f) - MeshArrays.fsize(x.f,2) - size(x) - size(x,3) - - MeshArrays.fijind(tmp,10) - MeshArrays.fsize(tmp) - @test isa(tmp,MeshArrays.gcmfaces) - - tmp=MeshArray(γ) - tmp1=findall(tmp.>0) - tmp[tmp1] - tmp[tmp1].=1.0 - size(tmp1) - tmp1[2] - view(tmp1,:) - @suppress show(tmp1) - similar(tmp1) - #tmp[tmp1] - - @suppress show(tmp) - MeshArrays.getindexetc(tmp,2) - MeshArray(γ,tmp.f,meta=tmp.meta) - MeshArray(γ,meta=tmp.meta) - MeshArray(γ,Float32,3,4) - MeshArray(γ,Float32,tmp.fSize,tmp.fIndex,2,3) - - tmp1=MeshArray(γ,Float32,3) - MeshArray(γ,tmp1.f,meta=tmp1.meta) - MeshArrays.getindexetc(tmp1,2,1) - end -end - -@testset "UnitGrid:" begin - C=MeshArray(randn(20,10)) - D=MeshArray(randn(20,10,3)) - - (Γ,γ)=Grids_simple.UnitGrid( (80,90) , (20,30) ; option="full") - @test isa(γ,gcmgrid) - - γ=Grids_simple.GridSpec_ones("PeriodicDomain",nP=10) - Γ=Grids_simple.GridLoad_ones(γ,option="full") - @test isa(Γ,NamedTuple) - - #various read/write functions - read(write(Γ.XC),γ) - tmp=tempname() - write(tmp,Γ.XC) - MeshArrays.read_tiles(tmp,Γ.XC) - MeshArrays.write_tiles(Γ.XC) - MeshArrays.write_tiles(tmp,Γ.XC) - @test isfile(tmp) - - xy=Grids_simple.xy_OISST() - xy=Grids_simple.xy_Oscar() - - xy=Grids_simple.xy_IAP() - gr=Grids_simple.grid_factors(xy) - dep=[10 100 1000]; msk=ones(gr[:XC].fSize[1]...,3) - gr=Grids_simple.grid_add_z(gr,dep,msk) - - @test haskey(gr,:hFacC) -end - -@testset "NEMO_GRID:" begin - lst=NEMO_GRID.variable_NTA() - nam=NEMO_GRID.variable_in_NEMO(:XC,lst) - @test isa(nam,Symbol) - - grid_data=Dict(:glamt=>zeros(1442,1021)) - XC_a=NEMO_GRID.convert_one_grid_variable(grid_data,lst[14]) - - grid=(XC=XC_a,YC=XC_a,RAC=XC_a,DXG=XC_a,DYG=XC_a) - G=NEMO_GRID.grid_to_MeshArrays(grid) - - XC_e=NEMO_GRID.exchange(G.XC) - @test isa(XC_e,MeshArray_wh) - - G=NEMO_GRID.add_angle_CS_SN(G) - @test haskey(G,:AngleCS) - - z=[1:10]; G=Dict() - grid_data=(gdept_0=z,gdepw_0=z,e3t_0=z,e3w_0=z) - NEMO_GRID.add_one_dim_variables!(G,grid_data) - @test isa(G,Dict) -end - -@testset "GridSpec:" begin - γ = GridSpec(ID=:Oscar) - γ = GridSpec(ID=:IAP) - GridLoad(GridSpec(ID=:OISST)) - GridLoad(GridSpec("ones")) -end - -@testset "Interpolation" begin - λ=interpolation_setup() - @test isa(λ,NamedTuple) - @test all(isfinite.(λ.f)) - # Test LLC90 interpolation coefficients - γ=GridSpec(ID=:LLC90) - λ=interpolation_setup(γ) - @test isa(λ,NamedTuple) - @test all(isfinite.(λ.f)) - # Test LLC270 interpolation coefficients - γ=GridSpec(ID=:LLC270) - λ=interpolation_setup(γ) - @test isa(λ,NamedTuple) - @test all(isfinite.(λ.f)) - # Test outside LLC90/LLC270 - γ=GridSpec(ID=:CS32) - λ=interpolation_setup(γ) - @test isa(λ,NamedTuple) - @test all(isfinite.(λ.f)) - - Γ=GridLoad(γ;option="light") - λ=interpolation_setup(Γ=Γ, - lon=[i for i=-170.:20.0:170., j=-80.:20.0:80.], - lat=[j for i=-170.:20.0:170., j=-80.:20.0:80.]) - @test isa(λ,NamedTuple) - @test all(isfinite.(λ.f)) -end - -@testset "Plotting:" begin - γ=GridSpec(ID=:LLC90) - Γ=GridLoad(γ;option="light") - D=Γ.Depth - λ=interpolation_setup() - - lines(pol_json); lines!(pol_json) - plot(pol_json); plot!(pol_json) - - basins=demo.ocean_basins() - AtlExt=demo.extended_basin(basins,:Atl) - sections,path_sec=demo.ocean_sections(Γ) - my_section=demo.one_section(Γ,[127 127],[-25 -68]) - - fig=MeshArrays.plot_examples(:smoothing_demo,D,D) - (fig1,fig2,fig3)=MeshArrays.plot_examples(:interpolation_demo,Γ) - - fake_ov=40e6*cosd.(360*(1:179)./100)*exp.(-0.1*(-20:29).^2)' - MeshArrays.plot_examples(:meriodional_overturning,Γ,fake_ov) - MeshArrays.plot_examples(:northward_transport,rand(179)) - - MeshArrays.plot_examples(:gradient_EN,λ,D,D) - MeshArrays.plot_examples(:gradient_xy,λ,D,D) - - ## more methods - - scatter(Γ.XC,Γ.YC,color=:black) - heatmap(D,interpolation=λ) - scatter!(current_axis(),Γ.XC,Γ.YC,color=:red) - - heatmap(D) #will display tile by tile - heatmap(D,interpolation=λ,title="ocean depth") #same but w title - - lon0=-160 - proj=Proj.Transformation(MA_preset=2,lon0=lon0) - Dint=reshape(Interpolate(D,λ.f,λ.i,λ.j,λ.w),size(λ.lon)) - - Interpolate(D,λ) - InterpolationFactors(Γ,30.0,30.0) - - ### - - MeshArraysMakieExt = Base.get_extension(MeshArrays, :MeshArraysMakieExt) - pol=MeshArraysMakieExt.pol_to_Makie(pol_shp) - dest="+proj=eqearth +lon_0=$(lon0) +lat_1=0.0 +x_0=0.0 +y_0=0.0 +ellps=GRS80" - MeshArraysMakieExt.split(pol,dest) - MeshArraysMakieExt.split(pol,Observable(dest)) - MeshArraysMakieExt.split(Observable(pol),Observable(dest)) - MeshArraysMakieExt.split(Observable(pol),dest) - - ### - - f = Figure() - ax = f[1, 1] = Axis(f, aspect = DataAspect(), title = "Ocean Depth (m)") - pr_ax=MeshArrays.ProjAxis(ax; proj=proj,lon0=lon0) - for a in [surface! contourf! contour!] - surf = a(pr_ax,λ.lon,λ.lat,0*λ.lat; color=Dint, - colorrange=(0.0,6000.0), colormap=:berlin, shading = NoShading) - end - lines!(pr_ax; polygons=pol_shp,color=:black,linewidth=0.5) - MeshArrays.grid_lines!(pr_ax;color=:lightgreen,linewidth=0.5) - f - - meta=(colorrange=(0.0,6000.0),cmap=:BrBG_10,ttl="Ocean Depth (m)",lon0=lon0) - data=(lon=λ.lon,lat=λ.lat,var=Dint,meta=meta) #,polygons=pol_shp) - plot_examples(:projmap,data,lon0,proj) - plot_examples(:simple_heatmap,data) - - MeshArraysMakieExt.heatmap_globalmap(D) - MeshArraysMakieExt.heatmap_interpolation(D,λ) - MeshArraysMakieExt.heatmap_xy(D,1:10,1:10) - -end - -@testset "nanmath" begin - x=[NaN 1 2] - nansum(x) - nansum(x,1) - nanmax(x,2) - nanmin(x,2) - - nanmean(NaN,1) - nanmean(1,NaN) - nanmean(2,1) - nanmean(NaN,NaN) -end - -@testset "plotting" begin - lon,lat,earth_img=demo.get_basemap() - plot_examples(:basemap,lon,lat,earth_img) - - lon0=-160 - proj=Proj.Transformation(MA_preset=2,lon0=lon0) - plot_examples(:baseproj,proj,lon0,pol=pol_shp) -end - -@testset "polygon operations" begin - fil=MeshArrays.Dataset("countries_shp1",do_read=false) - pol=MeshArrays.read_shp(fil,format=:Shapefile) - name,rule=MeshArrays.within_pol(pol; ID=11) - rule_vec = (x,y) -> rule.(x,y) - - np=10000; lo=-180 .+360*rand(np); la=-90 .+180*rand(np); - np_in=sum(rule_vec(lo,la)) - @test np_in>0 - - path_MITgcm=MITgcm.getdata("mitgcmsmallverif") - path_grid=joinpath(path_MITgcm,"MITgcm","verification","tutorial_held_suarez_cs","input") - pols,pols3D=MeshArrays.Polygons.polygons_demo(path_grid) - Depth=GridLoadVar("Depth",GridSpec(ID=:CS32)) - @test isa(pols[1][1,1],GI.Polygon) - - fig=MeshArrays.plot_examples(:polygons_plot,pols,color=Depth) - MeshArrays.plot_examples(:polygons_plot_dev1,pols,pols3D,sphere_view=true) - MeshArrays.plot_examples(:polygons_plot_dev1,pols,pols3D,sphere_view=false) - @test isa(fig,Makie.FigureAxisPlot) -end - -@testset "doctests" begin - doctest(MeshArrays; manual = false) -end +include("setup.jl") + +# Select testsets via env var, e.g.: JULIA_TESTSETS=transport,gcmfaces julia runtests.jl +# Omit the env var (or set to "all") to run everything. +const _RUN = split(get(ENV, "JULIA_TESTSETS", "all"), ",") +include_maybe(name) = ("all" in _RUN || name in _RUN) && include("testsets/$name.jl") + +include_maybe("mesharray_basic") +include_maybe("vertical_dim") +include_maybe("regional_integration") +include_maybe("transport") +include_maybe("gcmfaces") +include_maybe("unitgrid") +include_maybe("interpolation") +include_maybe("nemo_grid") +include_maybe("gridspec") +include_maybe("plotting_makie") +include_maybe("nanmath") +include_maybe("plotting_basemap") +include_maybe("polygon_ops") +include_maybe("doctests") diff --git a/test/setup.jl b/test/setup.jl new file mode 100644 index 0000000..c9dfe3b --- /dev/null +++ b/test/setup.jl @@ -0,0 +1,15 @@ +using Test, Documenter, Suppressor, MeshArrays, CairoMakie +import DataDeps, JLD2, Shapefile, GeoJSON, Proj, GeometryOps +import MITgcm +import MeshArrays: GI, NEMO_GRID, Grids_simple + +MeshArrays.Dataset("GRID_LL360") +MeshArrays.Dataset("GRID_LLC90") +MeshArrays.Dataset("GRID_LLC270") +MeshArrays.Dataset("GRID_CS32") + +pol_shp=MeshArrays.Dataset("countries_shp1") +pol_json=MeshArrays.Dataset("oceans_geojson1") + +p=dirname(pathof(MeshArrays)) +include(joinpath(p,"../examples/Demos.jl")) diff --git a/test/testsets/doctests.jl b/test/testsets/doctests.jl new file mode 100644 index 0000000..9436ef5 --- /dev/null +++ b/test/testsets/doctests.jl @@ -0,0 +1,3 @@ +@testset "doctests" begin + doctest(MeshArrays; manual = false) +end diff --git a/test/testsets/gcmfaces.jl b/test/testsets/gcmfaces.jl new file mode 100644 index 0000000..2e89430 --- /dev/null +++ b/test/testsets/gcmfaces.jl @@ -0,0 +1,56 @@ +@testset "gcmfaces type:" begin + for ID in (:PeriodicDomain, :CS32, :LLC270) + γ=GridSpec(ID=ID) + MeshArrays.gcmfaces(γ) + MeshArrays.gcmfaces(γ,Float32) + tmp=MeshArrays.gcmfaces(γ,Float32,3) + tmp[3,1,2]=1.0 + view(tmp,:,1,1:2) + MeshArrays.gcmfaces(γ,tmp.f) + MeshArrays.fsize(tmp) + MeshArrays.fsize(tmp,2) + MeshArrays.fsize(tmp.f) + MeshArrays.fsize(tmp.f,2) + size(tmp) + size(tmp,3) + tmp1=similar(tmp) + 2 .*tmp1 + findall(tmp.>0) + MeshArrays.nFacesEtc(tmp) + @suppress show(tmp) + + x=tmp[1:10,1,1:2]; y=x[2]; x[3]=1.0 + view(x,1:3,:,1) + MeshArrays.gcmsubset(γ,x.f,x.fSize,x.aSize,x.i,x.iSize) + MeshArrays.fsize(x.f) + MeshArrays.fsize(x.f,2) + size(x) + size(x,3) + + MeshArrays.fijind(tmp,10) + MeshArrays.fsize(tmp) + @test isa(tmp,MeshArrays.gcmfaces) + + tmp=MeshArray(γ) + tmp1=findall(tmp.>0) + tmp[tmp1] + tmp[tmp1].=1.0 + size(tmp1) + tmp1[2] + view(tmp1,:) + @suppress show(tmp1) + similar(tmp1) + #tmp[tmp1] + + @suppress show(tmp) + MeshArrays.getindexetc(tmp,2) + MeshArray(γ,tmp.f,meta=tmp.meta) + MeshArray(γ,meta=tmp.meta) + MeshArray(γ,Float32,3,4) + MeshArray(γ,Float32,tmp.fSize,tmp.fIndex,2,3) + + tmp1=MeshArray(γ,Float32,3) + MeshArray(γ,tmp1.f,meta=tmp1.meta) + MeshArrays.getindexetc(tmp1,2,1) + end +end diff --git a/test/testsets/gridspec.jl b/test/testsets/gridspec.jl new file mode 100644 index 0000000..0097441 --- /dev/null +++ b/test/testsets/gridspec.jl @@ -0,0 +1,6 @@ +@testset "GridSpec:" begin + γ = GridSpec(ID=:Oscar) + γ = GridSpec(ID=:IAP) + GridLoad(GridSpec(ID=:OISST)) + GridLoad(GridSpec("ones")) +end diff --git a/test/testsets/interpolation.jl b/test/testsets/interpolation.jl new file mode 100644 index 0000000..34b14b2 --- /dev/null +++ b/test/testsets/interpolation.jl @@ -0,0 +1,28 @@ +@testset "Interpolation" begin + λ=interpolation_setup() + @test isa(λ,NamedTuple) + @test all(isfinite.(λ.f)) + # Test LLC90 interpolation coefficients + γ=GridSpec(ID=:LLC90) + λ=interpolation_setup(γ) + @test isa(λ,NamedTuple) + @test all(isfinite.(λ.f)) + # Test LLC270 interpolation coefficients + γ=GridSpec(ID=:LLC270) + λ=interpolation_setup(γ) + @test isa(λ,NamedTuple) + @test all(isfinite.(λ.f)) + # Test outside LLC90/LLC270 + γ=GridSpec(ID=:CS32) + λ=interpolation_setup(γ) + @test isa(λ,NamedTuple) + @test all(isfinite.(λ.f)) + + Γ=GridLoad(γ;option="light") + λ=interpolation_setup(Γ=Γ, + lon=[i for i=-170.:20.0:170., j=-80.:20.0:80.], + lat=[j for i=-170.:20.0:170., j=-80.:20.0:80.]) + @test isa(λ,NamedTuple) + @test all(isfinite.(λ.f)) +end + diff --git a/test/testsets/mesharray_basic.jl b/test/testsets/mesharray_basic.jl new file mode 100644 index 0000000..00a70f0 --- /dev/null +++ b/test/testsets/mesharray_basic.jl @@ -0,0 +1,27 @@ +@testset "MeshArrays tests:" begin + for nTopo=1:4 + if nTopo==1; grTopo="CubeSphere"; nFaces=6; N=200; + elseif nTopo==2; grTopo="LatLonCap"; nFaces=5; N=200; + elseif nTopo==3; grTopo="PeriodicChannel"; nFaces=1; N=400; + elseif nTopo==4; grTopo="PeriodicDomain"; nFaces=1; N=400; + end; + Npt=nFaces*N*N + γ=MeshArrays.GridSpec_ones(grTopo,nP=N) + Γ=MeshArrays.GridLoad_ones(γ;option="full") + @test γ.class == grTopo + Rini= 0.; Rend= 0.; + (Rini,Rend,DXCsm,DYCsm)=demo2(Γ); + @test isa(Rend,MeshArray) + @test sum(isfinite.(Rend)) == Npt + Sini=sqrt(sum(Rini*Rini)/(Npt-1.0)) + Send=sqrt(sum(Rend*Rend)/(Npt-1.0)) + #println([Sini Send]) + @test isapprox(Sini,1.000; atol=1e-2) + @test isapprox(Send,0.093; atol=1e-2) + (dRdx,dRdy)=gradient(Rend,Γ) + MeshArrays.exchange_main(Rend,2) + MeshArrays.exchange_main(dRdx,dRdy,1) + (dRdx_e,dRdy_e)=MeshArrays.exchange_main(dRdx,dRdy,1) + @test isa(dRdx_e,MeshArray_wh) + end +end diff --git a/test/testsets/nanmath.jl b/test/testsets/nanmath.jl new file mode 100644 index 0000000..80da571 --- /dev/null +++ b/test/testsets/nanmath.jl @@ -0,0 +1,12 @@ +@testset "nanmath" begin + x=[NaN 1 2] + nansum(x) + nansum(x,1) + nanmax(x,2) + nanmin(x,2) + + nanmean(NaN,1) + nanmean(1,NaN) + nanmean(2,1) + nanmean(NaN,NaN) +end diff --git a/test/testsets/nemo_grid.jl b/test/testsets/nemo_grid.jl new file mode 100644 index 0000000..59e778e --- /dev/null +++ b/test/testsets/nemo_grid.jl @@ -0,0 +1,22 @@ +@testset "NEMO_GRID:" begin + lst=NEMO_GRID.variable_NTA() + nam=NEMO_GRID.variable_in_NEMO(:XC,lst) + @test isa(nam,Symbol) + + grid_data=Dict(:glamt=>zeros(1442,1021)) + XC_a=NEMO_GRID.convert_one_grid_variable(grid_data,lst[14]) + + grid=(XC=XC_a,YC=XC_a,RAC=XC_a,DXG=XC_a,DYG=XC_a) + G=NEMO_GRID.grid_to_MeshArrays(grid) + + XC_e=NEMO_GRID.exchange(G.XC) + @test isa(XC_e,MeshArray_wh) + + G=NEMO_GRID.add_angle_CS_SN(G) + @test haskey(G,:AngleCS) + + z=[1:10]; G=Dict() + grid_data=(gdept_0=z,gdepw_0=z,e3t_0=z,e3w_0=z) + NEMO_GRID.add_one_dim_variables!(G,grid_data) + @test isa(G,Dict) +end diff --git a/test/testsets/plotting_basemap.jl b/test/testsets/plotting_basemap.jl new file mode 100644 index 0000000..1564422 --- /dev/null +++ b/test/testsets/plotting_basemap.jl @@ -0,0 +1,8 @@ +@testset "plotting" begin + lon,lat,earth_img=demo.get_basemap() + plot_examples(:basemap,lon,lat,earth_img) + + lon0=-160 + proj=Proj.Transformation(MA_preset=2,lon0=lon0) + plot_examples(:baseproj,proj,lon0,pol=pol_shp) +end diff --git a/test/testsets/plotting_makie.jl b/test/testsets/plotting_makie.jl new file mode 100644 index 0000000..d1e71ea --- /dev/null +++ b/test/testsets/plotting_makie.jl @@ -0,0 +1,76 @@ +@testset "Plotting:" begin + γ=GridSpec(ID=:LLC90) + Γ=GridLoad(γ;option="light") + D=Γ.Depth + λ=interpolation_setup(Γ=Γ, + lon=[i for i=-170.:20.0:170., j=-80.:20.0:80.], + lat=[j for i=-170.:20.0:170., j=-80.:20.0:80.]) + λ=interpolation_setup() + + lines(pol_json); lines!(pol_json) + plot(pol_json); plot!(pol_json) + + basins=demo.ocean_basins() + AtlExt=demo.extended_basin(basins,:Atl) + sections,path_sec=demo.ocean_sections(Γ) + my_section=demo.one_section(Γ,[127 127],[-25 -68]) + + fig=MeshArrays.plot_examples(:smoothing_demo,D,D) + (fig1,fig2,fig3)=MeshArrays.plot_examples(:interpolation_demo,Γ) + + fake_ov=40e6*cosd.(360*(1:179)./100)*exp.(-0.1*(-20:29).^2)' + MeshArrays.plot_examples(:meriodional_overturning,Γ,fake_ov) + MeshArrays.plot_examples(:northward_transport,rand(179)) + + MeshArrays.plot_examples(:gradient_EN,λ,D,D) + MeshArrays.plot_examples(:gradient_xy,λ,D,D) + + ## more methods + + scatter(Γ.XC,Γ.YC,color=:black) + heatmap(D,interpolation=λ) + scatter!(current_axis(),Γ.XC,Γ.YC,color=:red) + + heatmap(D) #will display tile by tile + heatmap(D,interpolation=λ,title="ocean depth") #same but w title + + lon0=-160 + proj=Proj.Transformation(MA_preset=2,lon0=lon0) + Dint=reshape(Interpolate(D,λ.f,λ.i,λ.j,λ.w),size(λ.lon)) + + Interpolate(D,λ) + InterpolationFactors(Γ,30.0,30.0) + + ### + + MeshArraysMakieExt = Base.get_extension(MeshArrays, :MeshArraysMakieExt) + pol=MeshArraysMakieExt.pol_to_Makie(pol_shp) + dest="+proj=eqearth +lon_0=$(lon0) +lat_1=0.0 +x_0=0.0 +y_0=0.0 +ellps=GRS80" + MeshArraysMakieExt.split(pol,dest) + MeshArraysMakieExt.split(pol,Observable(dest)) + MeshArraysMakieExt.split(Observable(pol),Observable(dest)) + MeshArraysMakieExt.split(Observable(pol),dest) + + ### + + f = Figure() + ax = f[1, 1] = Axis(f, aspect = DataAspect(), title = "Ocean Depth (m)") + pr_ax=MeshArrays.ProjAxis(ax; proj=proj,lon0=lon0) + for a in [surface! contourf! contour!] + surf = a(pr_ax,λ.lon,λ.lat,0*λ.lat; color=Dint, + colorrange=(0.0,6000.0), colormap=:berlin, shading = NoShading) + end + lines!(pr_ax; polygons=pol_shp,color=:black,linewidth=0.5) + MeshArrays.grid_lines!(pr_ax;color=:lightgreen,linewidth=0.5) + f + + meta=(colorrange=(0.0,6000.0),cmap=:BrBG_10,ttl="Ocean Depth (m)",lon0=lon0) + data=(lon=λ.lon,lat=λ.lat,var=Dint,meta=meta) #,polygons=pol_shp) + plot_examples(:projmap,data,lon0,proj) + plot_examples(:simple_heatmap,data) + + MeshArraysMakieExt.heatmap_globalmap(D) + MeshArraysMakieExt.heatmap_interpolation(D,λ) + MeshArraysMakieExt.heatmap_xy(D,1:10,1:10) + +end diff --git a/test/testsets/polygon_ops.jl b/test/testsets/polygon_ops.jl new file mode 100644 index 0000000..c0b9ee2 --- /dev/null +++ b/test/testsets/polygon_ops.jl @@ -0,0 +1,21 @@ +@testset "polygon operations" begin + fil=MeshArrays.Dataset("countries_shp1",do_read=false) + pol=MeshArrays.read_shp(fil,format=:Shapefile) + name,rule=MeshArrays.within_pol(pol; ID=11) + rule_vec = (x,y) -> rule.(x,y) + + np=10000; lo=-180 .+360*rand(np); la=-90 .+180*rand(np); + np_in=sum(rule_vec(lo,la)) + @test np_in>0 + + path_MITgcm=MITgcm.getdata("mitgcmsmallverif") + path_grid=joinpath(path_MITgcm,"MITgcm","verification","tutorial_held_suarez_cs","input") + pols,pols3D=MeshArrays.Polygons.polygons_demo(path_grid) + Depth=GridLoadVar("Depth",GridSpec(ID=:CS32)) + @test isa(pols[1][1,1],GI.Polygon) + + fig=MeshArrays.plot_examples(:polygons_plot,pols,color=Depth) + MeshArrays.plot_examples(:polygons_plot_dev1,pols,pols3D,sphere_view=true) + MeshArrays.plot_examples(:polygons_plot_dev1,pols,pols3D,sphere_view=false) + @test isa(fig,Makie.FigureAxisPlot) +end diff --git a/test/testsets/regional_integration.jl b/test/testsets/regional_integration.jl new file mode 100644 index 0000000..c039d43 --- /dev/null +++ b/test/testsets/regional_integration.jl @@ -0,0 +1,30 @@ +@testset "Regional Integration:" begin + G,M,files=Integration.example() + @suppress show(M) + + allones=1.0 .+0*G.hFacC + vol0=sum(G.RAC*G.DRF*G.hFacC) + + vol=Integration.volumes(M,G) + @test isapprox(sum(vol),vol0) + + G,M,files=Integration.example(option=:streamlined_loop) + vol=[b(allones) for b in M.h_sum] + @test isapprox(sum(vol),vol0) + + rgns=Integration.define_regions(option=:basins,grid=G) + rgns=Integration.define_regions(option=:dlat_10,grid=G) + rgns=Integration.define_regions(option=(30,10),grid=G) + @test isa(rgns,NamedTuple) + + G,M,files=Integration.example() + files=fill("?",3) + rd0(F,var,tim,tmp)=tim*ones(tmp) + H=Integration.loops(M,files=files,rd=rd0) + @test isa(H,Array) + + G,M,files=Integration.example(option=:streamlined_loop) + files=fill("?",3) + H=Integration.streamlined_loop(M,files=files,rd=rd0) + @test isa(H,Array) +end diff --git a/test/testsets/transport.jl b/test/testsets/transport.jl new file mode 100644 index 0000000..ef8b890 --- /dev/null +++ b/test/testsets/transport.jl @@ -0,0 +1,92 @@ +@testset "Transport computations:" begin + #Load grid and transport / vector field + γ=GridSpec(ID=:LLC90) + @suppress show(γ) + Γ=GridLoad(γ,option=:full) + @suppress show(Γ.XC) + + path=MeshArrays.Dataset("GRID_LLC90") + Tx=γ.read(joinpath(path,"TrspX.bin"),MeshArray(γ,Float32)) + Ty=γ.read(joinpath(path,"TrspY.bin"),MeshArray(γ,Float32)) + plot(Γ.XC) + + hFacC=GridLoadVar("hFacC",γ) + μ=land_mask(hFacC[:,1]) + + lons=[-68 -63]; lats=[-54 -66]; name="Drake Passage" + Trsct=Transect(name,lons,lats,Γ,segment=:long,format=:NamedTuple) + Trsct=Transect(name,lons,lats,Γ) + + mask=demo.extended_basin(demo.ocean_basins(),:Pac) + edge=edge_path("Pacific Ocean Edge",mask,Γ) + + #Various vector operations + hFacW=GridLoadVar("hFacW",γ) + hFacS=GridLoadVar("hFacS",γ) + RAZ=GridLoadVar("RAZ",γ) + + U=0*hFacW; V=0*hFacS; + UVtoTransport(U,V,Γ) + UVtoUEVN(U[:,1],V[:,1],Γ) + curl(U[:,1],V[:,1], merge(Γ,(hFacW=hFacW,hFacS=hFacS,RAZ=RAZ)) ) + dD=zeros(γ) + MeshArrays.UVtoSpeed!(U[:,1],V[:,1],Γ,dD) + + #Meridional transport integral + 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)] + plot(LC) + plot(LC[1]) + + LC_nt=LatitudeCircles(L,Γ,format=:NamedTuple) + @test isa(LC_nt,Array) && isa(LC_nt[1],NamedTuple) + @test LC_nt[1].lat == L[1] + @test isapprox(ThroughFlow(uv,LC_nt[1],Γ), ThroughFlow(uv,LC[1],Γ)) + + x=zeros(γ) + fill!(x,1.0) + y=fill(-1.0,γ) + extrema(y) + @test minimum(y)