Merge pull request JuliaClimate#11 from JuliaClimate/docsEtc-gf02

Docs etc gf02

Project.toml

@@ -1,7 +1,7 @@
 name = "IndividualDisplacements"
 uuid = "b92f0c32-5b7e-11e9-1d7b-238b2da8b0e6"
 authors = ["gaelforget <gforget@mit.edu>"]
-version = "0.1.6"
+version = "0.1.7"
 
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"

docs/src/index.md

@@ -6,7 +6,35 @@ Inter-operability with popular climate model grids and their representation in [
 
 The `VelComp!` and `VelComp` functions compute the velocity of tracked points. `tests/runtests.jl` uses solid body rotation as a benchmark (see below).
 
-<img src="https://github.com/JuliaClimate/IndividualDisplacements.jl/blob/master/examples/SolidBodyRotation.png" width="50%">
+![SolidBodyRotation](https://github.com/JuliaClimate/IndividualDisplacements.jl/raw/master/examples/SolidBodyRotation.png)
+
+## List Of Examples
+
+Solid body rotation is used for unit testing:
+
+```
+test/runtests.jl
+examples/SolidBodyRotation.jl
+```
+
+The two examples documented under `API` + extenstions to `example2`:
+
+```
+src/examples.jl
+examples/example2fleet.jl
+examples/example2more.jl
+```
+
+An intermediate example: `examples/PeriodicChannel_fleet.jl`
+
+Two examples use `VelComp!` and `update_locations.jl`:
+
+```
+examples/PeriodicDomainRandom_fleet.jl
+examples/GlobalDomain_fleet.jl
+```
+
+![SolidBodyRotation](https://github.com/JuliaClimate/IndividualDisplacements.jl/raw/master/examples/LatLonCap300mDepth.png)
 
 ## API Guide
 

examples/llc90_global_fleet.jl → examples/GlobalDomain_fleet.jl

@@ -27,7 +27,8 @@ include(joinpath(dirname(pathof(MeshArrays)),"../examples/Plots.jl"))
 
 # Put grid variables in a dictionary.
 
-mygrid=GridSpec("LatLonCap","GRID_LLC90/"); GridVariables=GridLoad(mygrid);
+mygrid=GridSpec("LatLonCap","GRID_LLC90/"); 
+GridVariables=GridLoad(mygrid);
 GridVariables=merge(GridVariables,
     IndividualDisplacements.NeighborTileIndices_cs(GridVariables));
 
@@ -117,7 +118,7 @@ sol_one = solve(prob,Tsit5(),reltol=1e-4,abstol=1e-4)
 sol_two = solve(prob,Euler(),dt=1e6)
 size(sol_one)
 
-# Define initial condition array.
+# Define initial condition array (**retired code**).
 
 if false
         fIndex = 1
@@ -139,6 +140,8 @@ if false
         du = fill(0.0, size(uInitS))
 end
 
+# Define initial condition array
+
 # +
 uInitS = Array{Float64,2}(undef, 3, prod(XC.grid.ioSize))
 kk = 0

examples/global_ocean_fleet.jl → examples/PeriodicChannel_fleet.jl

@@ -29,7 +29,7 @@ include(joinpath(p,"plot_pyplot.jl"))
 # Put grid variables in a dictionary.
 
 # +
-mygrid=gcmgrid("llc90_latlon/","ll",1,[(360,178)], [360 178], Float32, read, write)
+mygrid=gcmgrid("llc90_latlon/","Periodic",1,[(360,178)], [360 178], Float32, read, write)
 
 GridVariables=Dict("XC" => read(mygrid.path*"XC.latlon.data",MeshArray(mygrid,Float32)),
 "YC" => read(mygrid.path*"YC.latlon.data",MeshArray(mygrid,Float32)),

examples/body_rotation_test.jl → examples/SolidBodyRotation.jl

@@ -28,7 +28,7 @@ using IndividualDisplacements, MeshArrays, OrdinaryDiffEq, Plots
 # Put grid variables in a dictionary:
 
 # +
-mygrid=gcmgrid("flt_example/","ll",1,[(80,42)], [80 42], Float32, read, write)
+mygrid=gcmgrid("flt_example/","PeriodicChanel",1,[(80,42)], [80 42], Float32, read, write)
 nr=8
 
 XC=MeshArray(mygrid,Float32); XC[1]=vec(2500.:5000.:397500.0)*ones(1,42);

examples/example2fleet.jl

@@ -0,0 +1,83 @@
+# ---
+# jupyter:
+#   jupytext:
+#     formats: ipynb,jl:light
+#     text_representation:
+#       extension: .jl
+#       format_name: light
+#       format_version: '1.4'
+#       jupytext_version: 1.2.4
+#   kernelspec:
+#     display_name: Julia 1.3.0-rc4
+#     language: julia
+#     name: julia-1.3
+# ---
+
+# # This notebook
+#
+# _Notes:_ For documentation see <https://gaelforget.github.io/MeshArrays.jl/stable/>, <https://docs.juliadiffeq.org/latest/solvers/ode_solve.html> and <https://en.wikipedia.org/wiki/Displacement_(vector)>
+
+# ## 1. import software
+
+using IndividualDisplacements, MeshArrays, OrdinaryDiffEq
+using Plots, Statistics, MITgcmTools, DataFrames
+p=dirname(pathof(IndividualDisplacements))
+include(joinpath(p,"plot_pyplot.jl"))
+
+# ## 2. Read gridded variables as `MeshArray`s
+
+dirIn="flt_example/"
+prec=Float32
+df=ReadDisplacements(dirIn,prec)
+uvetc=IndividualDisplacements.example2_setup()
+
+# ## 2. Recompute displacements from gridded flow fields
+
+# Define the ODE problem
+
+# +
+nSteps=2998
+tspan = (0.0,nSteps*3600.0)
+
+comp_vel=IndividualDisplacements.VelComp
+get_vel=IndividualDisplacements.VelCopy
+# -
+
+# Set up initial conditions
+#
+# _Note how the size of sol (ie nb of steps) depends on initial location:_
+
+# +
+#ii1=1:10:80; ii2=1:10:42; #->sol is (2, 40, 40065)
+#ii1=30:37; ii2=16:20; #->sol is (2, 40, 9674)
+#ii1=10:17; ii2=16:20; #->sol is (2, 40, 51709)
+ii1=5:5:40; ii2=5:5:25; #->sol is (2, 40, 51709)
+
+n1=length(ii1); n2=length(ii2);
+uInitS=Array{Float64,2}(undef,(2,n1*n2))
+for i1 in eachindex(ii1); for i2 in eachindex(ii2);
+        i=i1+(i2-1)*n1
+        uInitS[1,i]=ii1[i1]-0.5
+        uInitS[2,i]=ii2[i2]-0.5       
+end; end;
+
+prob = ODEProblem(comp_vel,uInitS,tspan,uvetc)
+# -
+
+# Compute solution
+
+sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
+size(sol)
+
+# Display results
+
+# +
+ID=collect(1:size(sol,2))*ones(1,size(sol,3))
+lon=5000* mod.(sol[1,:,:],80)
+lat=5000* mod.(sol[2,:,:],42)
+df = DataFrame(ID=Int.(ID[:]), lon=lon[:], lat=lat[:])
+
+PyPlot.figure(); PlotBasic(df,size(sol,2),100000.0)
+# -
+
+

examples/ex_2_more.jl → examples/example2more.jl

@@ -37,66 +37,15 @@ PyPlot.figure(); PlotBasic(df,300,100000.0)
 
 # ## 3. Read gridded variables via `MeshArrays.jl`
 #
-# Put grid variables in a dictionary.
-#
-# _Note:_ `myread` function deals with tiled files from `flt_example/`.
-
-# +
-import IndividualDisplacements: myread
-
-mygrid=gcmgrid("flt_example/","ll",1,[(80,42)], [80 42], Float32, read, write)
-nr=8
-
-GridVariables=Dict("XC" => myread(mygrid.path*"XC",MeshArray(mygrid,Float32)),
-"YC" => myread(mygrid.path*"YC",MeshArray(mygrid,Float32)),
-"XG" => myread(mygrid.path*"XG",MeshArray(mygrid,Float32)),
-"YG" => myread(mygrid.path*"YG",MeshArray(mygrid,Float32)),
-"dx" => 5000.0);
-# -
-
-# Put velocity fields in in another dictionary and merge the two dictionaries.
-
-# +
-t0=0.0 #approximation / simplification
-t1=18001.0*3600.0
-
-u0=myread(mygrid.path*"U.0000000001",MeshArray(mygrid,Float32,nr))
-u1=myread(mygrid.path*"U.0000018001",MeshArray(mygrid,Float32,nr))
-v0=myread(mygrid.path*"V.0000000001",MeshArray(mygrid,Float32,nr))
-v1=myread(mygrid.path*"V.0000018001",MeshArray(mygrid,Float32,nr))
+# Put gridded variables in a dictionary.
 
-kk=3 #3 to match -1406.25 in pkg/flt output
-u0=u0[:,kk]; u1=u1[:,kk];
-v0=v0[:,kk]; v1=v1[:,kk];
-
-u0=u0./GridVariables["dx"]
-u1=u1./GridVariables["dx"]
-v0=v0./GridVariables["dx"]
-v1=v1./GridVariables["dx"]
-
-uvt = Dict("u0" => u0, "u1" => u1, "v0" => v0, "v1" => v1, "t0" => t0, "t1" => t1)
-
-uvetc=merge(uvt,GridVariables);
-# -
+uvetc=IndividualDisplacements.example2_setup()
 
 # ## 4. Visualize velocity fields
 
-# +
-mskW=myread(mygrid.path*"hFacW",MeshArray(mygrid,Float32,nr))
-mskW=1.0 .+ 0.0 * mask(mskW[:,kk],NaN,0.0)
-mskS=myread(mygrid.path*"hFacS",MeshArray(mygrid,Float32,nr))
-mskS=1.0 .+ 0.0 * mask(mskS[:,kk],NaN,0.0)
-
-msk=Dict("mskW" => mskW, "mskS" => mskS)
-
-uvetc=merge(uvetc,msk);
-# -
-
-heatmap(mskW[1,1].*u0[1,1],title="U at the start")
-
-heatmap(mskS[1,1].*v0[1,1],title="V at the start")
+heatmap(uvetc["mskW"][1,1].*uvetc["u0"][1,1],title="U at the start")
 
-heatmap(mskW[1,1].*u1[1,1]-u0[1,1],title="U end - U start")
+heatmap(uvetc["mskW"][1,1].*uvetc["u1"][1,1]-uvetc["u0"][1,1],title="U end - U start")
 
 # ## 5. Visualize trajectories from `MITgcm/pkg/flt`
 #
@@ -107,13 +56,15 @@ tmp[1:4,:]
 
 # Super-impose trajectory over velocity field (first for u ...)
 
-PyPlot.contourf(GridVariables["XG"].f[1], GridVariables["YC"].f[1], mskW.f[1].*u0.f[1])
+PyPlot.contourf(uvetc["XG"].f[1], uvetc["YC"].f[1], 
+    uvetc["mskW"].f[1].*uvetc["u0"].f[1])
 PyPlot.plot(tmp[:,:lon],tmp[:,:lat],color="r")
 colorbar()
 
 # Super-impose trajectory over velocity field (... then for v)
 
-PyPlot.contourf(GridVariables["XG"].f[1], GridVariables["YC"].f[1], mskS.f[1].*v0.f[1])
+PyPlot.contourf(uvetc["XG"].f[1], uvetc["YC"].f[1], 
+    uvetc["mskS"].f[1].*uvetc["v0"].f[1])
 PyPlot.plot(tmp[:,:lon],tmp[:,:lat],color="r")
 colorbar()