Merge pull request #103 from gaelforget/v0p3p11c

streamline examples and configurations

examples/jupyter/flow_fields.jl

@@ -1,140 +1,4 @@
-using MeshArrays, OceanStateEstimation, MITgcmTools
-import IndividualDisplacements.NetCDF as NetCDF
-
-"""
-    solid_body_rotation(np,nz)
-
-Set up an idealized flow field which consists of 
-[rigid body rotation](https://en.wikipedia.org/wiki/Rigid_body), 
-plus a convergent term, plus a sinking term.
-
-```
-u,v,w=solid_body_rotation(12,4)
-```
-"""
-function solid_body_rotation(np,nz)
-    Γ=simple_periodic_domain(np)
-    Γ = UnitGrid(Γ.XC.grid;option="full")
-    γ=Γ.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])
-
-    #... 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])
-
-    #Replicate u,v in vertical dimension
-    uu=MeshArray(γ,γ.ioPrec,nz)
-    [uu[k]=u[1] for k=1:nz]
-    vv=MeshArray(γ,γ.ioPrec,nz)
-    [vv[k]=v[1] for k=1:nz]
-
-    #Vertical velocity component w    
-    w=fill(-0.01,MeshArray(γ,γ.ioPrec,nz+1))
-
-    return write(uu),write(vv),write(w)
-end
-
-"""
-    OCCA_FlowFields(;backward_in_time::Bool=false,nmax=Inf)
-
-Define gridded variables and return result as NamedTuple
-"""
-function OCCA_FlowFields(;backward_in_time::Bool=false,nmax=Inf)
-
-   γ=GridSpec("PeriodicChannel",MeshArrays.GRID_LL360)
-   Γ=GridLoad(γ;option="full")
-   n=length(Γ.RC)
-   isfinite(nmax) ? n=min(n,Int(nmax)) : nothing
-
-   g=Γ.XC.grid
-   func=(u -> MeshArrays.update_location_dpdo!(u,g))
-
-   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)
-
-   function rd(filename, varname,n)
-   fil = NetCDF.open(filename, varname)
-   siz = size(fil)
-   tmp = zeros(siz[1:2]...,n)
-   [tmp .+= fil[:,:,1:n,t] for t=1:12]
-   tmp ./= 12.0
-   tmp[findall(tmp.<-1e22)] .= 0.0
-   return tmp
-   end
-
-   fileIn=OCCAclim_path*"DDuvel.0406clim.nc"
-   u=s*read(rd(fileIn,"u",n),MeshArray(γ,Float32,n))
-
-   fileIn=OCCAclim_path*"DDvvel.0406clim.nc"
-   v=s*read(rd(fileIn,"v",n),MeshArray(γ,Float32,n))
-
-   fileIn=OCCAclim_path*"DDwvel.0406clim.nc"
-   w=s*rd(fileIn,"w",n)
-   w=-cat(w,zeros(360, 160),dims=3)
-   w[:,:,1] .=0.0
-   w=read(w,MeshArray(γ,Float32,n+1))
-
-   fileIn=OCCAclim_path*"DDtheta.0406clim.nc"
-   θ=read(rd(fileIn,"theta",n),MeshArray(γ,Float32,n))
-
-#   fileIn=OCCAclim_path*"DDsalt.0406clim.nc"
-#   𝑆=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]
-   end
-
-   for i in eachindex(u)
-      u[i]=circshift(u[i],[-180 0])
-      v[i]=circshift(v[i],[-180 0])
-      θ[i]=circshift(θ[i],[-180 0])
-#      𝑆[i]=circshift(𝑆[i],[-180 0])
-   end
-
-   for i in eachindex(w)
-      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[1:180,:]=tmpx[1:180,:] .- 360.0
-   Γ.XC[1]=tmpx
-
-   tmpx=circshift(Γ.XG[1],[-180 0])
-   tmpx[1:180,:]=tmpx[1:180,:] .- 360.0
-   Γ.XG[1]=tmpx
-   Γ.Depth[1]=circshift(Γ.Depth[1],[-180 0])
-
-   t0=0.0; t1=86400*366*2.0;
-
-   for k=1:n
-    (tmpu,tmpv)=exchange(u[:,k],v[:,k],1)
-    u[:,k]=tmpu
-    v[:,k]=tmpv
-   end
-   for k=1:n+1
-    tmpw=exchange(w[:,k],1)
-    w[:,k]=tmpw
-   end
-
-   𝑃=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))
-
-   return 𝑃,𝐷,Γ
-
-end
+using MeshArrays, MITgcmTools
 
 """
     test1_setup()

examples/jupyter/solid_body_rotation.jl

@@ -26,8 +26,6 @@
 
 using IndividualDisplacements
 import IndividualDisplacements.DataFrames: DataFrame
-p=dirname(pathof(IndividualDisplacements))
-include(joinpath(p,"../examples/jupyter/flow_fields.jl"))
 
 #nb # %% {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
 # ### 1.2  Flow Fields
@@ -38,7 +36,7 @@ include(joinpath(p,"../examples/jupyter/flow_fields.jl"))
 
 np,nz=16,4 #gridded domain size (horizontal and vertical)
 
-u,v,w=solid_body_rotation(np,nz) #staggered velocity arrays
+u,v,w=IndividualDisplacements.solid_body_rotation(np,nz) #staggered velocity arrays
 
 𝐹=FlowFields(u,u,v,v,0*w,1*w,[0,19.95*2*pi]); #FlowFields data structure
 

examples/jupyter/three_dimensional_ocean.jl

@@ -17,90 +17,32 @@
 # ## 1. Load Software
 #
 
-using IndividualDisplacements, OceanStateEstimation
+using IndividualDisplacements
 import IndividualDisplacements.DataFrames: DataFrame
 
 p=dirname(pathof(IndividualDisplacements))
-include(joinpath(p,"../examples/jupyter/helper_functions.jl"))
-
-OceanStateEstimation.get_occa_velocity_if_needed();
+include(joinpath(p,"../examples/worldwide/OCCA_FlowFields.jl"))
 
 #nb # %% {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
 # ## 2.1 Ocean Circulation Setup
 #
 
-𝑃,𝐷,Γ=OCCA_FlowFields(nmax=5);
-
-#nb # %% {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
-# ## 2.3 Initialize Individual Positions
-#
-
-"""
-    initial_positions(Γ; nf=10000, lon_rng=(-160.0,-159.0), lat_rng=(30.0,31.0))
-
-Randomly assign initial positions in longitude,latitude ranges. Positions are 
-expressed in, normalized, grid point units (x,y in the 0,nx and 0,ny range). 
-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(Γ::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])
-   return x,y
-end
+𝑃,𝐷=OCCA_FlowFields.setup(nmax=5);
 
 #nb # %% {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
-# ## 2.3 Diagnostic / Post-Processing Setup
+# ## 2.2 Initialize Positions
 #
 
-custom🔴 = DataFrame(ID=Int[], fid=Int[], x=Float64[], y=Float64[], 
-   k=Float64[], z=Float64[], iso=Float64[], t=Float64[], 
-   lon=Float64[], lat=Float64[], year=Float64[], col=Symbol[])
-
-function custom🔧(sol,𝑃::𝐹_MeshArray3D,𝐷::NamedTuple;id=missing,𝑇=missing)
-   df=postprocess_MeshArray(sol,𝑃,𝐷,id=id,𝑇=𝑇)
-   add_lonlat!(df,𝐷.XC,𝐷.YC)
-
-   #add year (convenience time axis for plotting)
-   df.year=df.t ./86400/365
-
-   #add depth (i.e. the 3rd, vertical, coordinate)
-   k=[[sol[i][3,1] for i in 1:size(sol,3)];[sol[i][3,end] for i in 1:size(sol,3)]]
-
-   nz=length(𝐼.𝑃.u1)
-   df.k=min.(max.(k[:],Ref(0.0)),Ref(nz)) #level
-   k=Int.(floor.(df.k)); w=(df.k-k); 
-   df.z=𝐷.RF[1 .+ k].*(1 .- w)+𝐷.RF[2 .+ k].*w #depth
-
-   #add one isotherm depth
-   θ=0.5*(𝐷.θ0+𝐷.θ1)
-   d=isosurface(θ,15,𝐷)
-   d[findall(isnan.(d))].=0.
-   df.iso=interp_to_xy(df,exchange(d));
-
-   #add color = f(iso-z)
-   c=fill(:gold,length(df.iso))
-   c[findall(df.iso.<df.z)].=:violet
-   df.col=c
-
-   #to plot e.g. Pacific Ocean transports, shift longitude convention?
-   df.lon[findall(df.lon .< 0.0 )] = df.lon[findall(df.lon .< 0.0 )] .+360.0
-   return df
-end
+nf=100; lo=(-160.0,-150.0); la=(30.0,40.0); level=2.5; 
+df=OCCA_FlowFields.initial_positions(𝐷.Γ, nf, lo, la, level)
 
 # ## 2.4 Individuals Data Structure
 #
-# Set up `Individuals` data structure with `nf` particles moving in 3D 
+# Set up `Individuals` data structure with `nf` iindividuals moving in 3D 
 # on a regular 1 degree resolution grid covering most of the Globe.
 
-nf=100; lo=(-160.0,-150.0); la=(30.0,40.0); kk=2.5; 
-df=DataFrame(:z => fill(kk,nf),:f => fill(1,nf))
-(df.x,df.y)=initial_positions(Γ, nf, lo, la)
-
-𝐼=Individuals(𝑃,df.x,df.y,df.z,df.f,(🔴=custom🔴,🔧=custom🔧, 𝐷=𝐷))
+𝐼=Individuals(𝑃,df.x,df.y,df.z,df.f,
+   (🔴=OCCA_FlowFields.custom🔴,🔧=OCCA_FlowFields.custom🔧, 𝐷=𝐷))
 
 #nb # %% {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
 # ## 3.1 Compute Displacements
@@ -110,30 +52,7 @@ df=DataFrame(:z => fill(kk,nf),:f => fill(1,nf))
 
 ∫!(𝐼,𝑇)
 
-#nb # %% {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
-# ## 3.2 Analyze Results
-#
-# The recorded simulation output, 🔴, is a in the [DataFrames](https://juliadata.github.io/DataFrames.jl/latest/) tabular format, which is easily manipulated or plotted.
-
-# - either `Plots.jl`:
-#
-# ```
-# include(joinpath(p,"../examples/recipes_plots.jl"))
-# p=plot(𝐼)
-# #p=map(𝐼,OceanDepthLog(Γ))
-# display(p)
-# ```
-
-# - or `Makie.jl`:
-#
-# ```
-# include(joinpath(p,"../examples/recipes_Makie.jl"))
-# #p=plot(𝐼)
-# p=plot_paths(𝐼.🔴,100,180.);
-# display(p)
-# ```
-
-# ## 3.3 Alternatives (optional / unit testing)
+# ## 3.2 Alternatives (optional / unit testing)
 
 (x,y,z,f)=𝐼.📌[1]
 𝐽=Individuals(𝐼.𝑃,x,y,z,f)