Merge pull request #59 from gaelforget/v0p2p20

V0p2p20

Project.toml

@@ -1,7 +1,7 @@
 name = "MeshArrays"
 uuid = "cb8c808f-1acf-59a3-9d2b-6e38d009f683"
 authors = ["gaelforget <[email protected]>"]
-version = "0.2.20"
+version = "0.2.21"
 
 [deps]
 CatViews = "81a5f4ea-a946-549a-aa7e-2a7f63a27d31"

README.md

@@ -16,25 +16,30 @@
 ```
 using Pkg
 Pkg.add("MeshArrays")
-Pkg.test("MeshArrays")
 ```
 
 ### Workflow Example
 
-The diffusive smoother presented here as an example uses `MeshArrays.jl` to compute partial derivatives over a global domain / grid, which involves transfering data between neighboring subdomain arrays. In this workflow example, we 
+In this example, we use a _doubly periodic_ domain with _16 subdomains_ of _20x20 grid points_ each. The chosen workflow applies a diffusion-based smoother to a noise field. Thus it demonstrate the use of `MeshArrays.jl` to compute partial derivatives and budgets over a global domain, which notably requires transfering data between neighboring subdomain arrays. 
+
+The workflow example proceeds as follows:
 
 1. generate a global grid decomposition
 2. seed random noise across global domain
 3. smooth out noise by applying diffusion globally
-4. plots the `outer` array of subdomain / `inner` arrays
+4. plots the `outer` array of subdomains / `inner` arrays
 
 ```
-using MeshArrays; p=dirname(pathof(MeshArrays))
-γ,Γ=GridOfOnes("PeriodicDomain",16,20)
+using MeshArrays
+
+(nP,nQ,nF)=(20,20,16)
+γ=gcmgrid("","PeriodicDomain",nF,fill((nP,nQ),nF), 
+			[nP nQ*nF],Float32, read, write)
+Γ=UnitGrid(γ) #step 1
 
+p=dirname(pathof(MeshArrays))
 include(joinpath(p,"../examples/Demos.jl"))
-(xi,xo,_,_)=demo2(Γ);
-show(xo)
+(xi,xo,_,_)=demo2(Γ) #steps 2 & 3
 
 using Plots
 include(joinpath(p,"../examples/Plots.jl"))
@@ -45,7 +50,7 @@ heatmap(xo,clims=(-0.25,0.25),colorbar=false,tickfont = (4, :black))
 
 ### Global Grids
 
-In the previous example we used a basic _doubly periodic_  domain with _16 subdomains_ of _40x40 grid points_ each. However, `MeshArrays.jl` also readily supports more elaborate global grid configurations, such as the ones shown below, which are commonly used in modeling climate.
+`MeshArrays.jl` also readily supports more elaborate global grid configurations, such as the ones shown below, which are commonly used in modeling climate.
 
 <img src="docs/images/sphere_all.png" width="40%">
 

docs/src/index.md

@@ -1,9 +1,9 @@
 # MeshArrays.jl
 
-`MeshArrays.jl` defines an array type that can contain / organize / distribute collections of inter-connected arrays as generally done in climate models (see [Earth Grids](@ref)). `MeshArrays.jl`'s data structures can be used to simulate and analyze key variables of the climate system such as global [transports](https://doi.org/10.1038/s41561-019-0333-7) and particle [displacements](https://doi.org/10.21105/joss.02813).
+`MeshArrays.jl` defines an array type that can contain / organize / distribute collections of inter-connected arrays as generally done in climate models (see [Earth Grids](@ref)). `MeshArrays.jl`'s data structures can be used to simulate and analyze key variables of the climate system such as [particles](https://doi.org/10.21105/joss.02813) and [transports](https://doi.org/10.1038/s41561-019-0333-7).
 
 !!! note
-    [Global Ocean Notebooks](https://github.com/JuliaClimate/GlobalOceanNotebooks.git), [IndividualDisplacements.jl](https://juliaclimate.github.io/IndividualDisplacements.jl/dev/), and [MITgcmUtils.jl]() illustrate various use cases for example.
+    [These Notebooks](https://github.com/JuliaClimate/GlobalOceanNotebooks.git), [IndividualDisplacements.jl](https://juliaclimate.github.io/IndividualDisplacements.jl/dev/), and [MITgcmUtils.jl]() illustrate various use cases for example.
 
 
 ## Contents

docs/src/main.md

@@ -2,11 +2,11 @@
 
 ## Summary
 
-The `MeshArray` type is a sub-type of `AbstractArray` with an `outer array` where each element is itself a 2D `inner array`. This setup potentially allows different choices for the outer and inner arrays – for example `DistributedArrays` and `AxisArrays`, respectively, could be an option. `MeshArrays.jl` thus provides a simple but general solution to analyze or e.g. simulate climate system variables. 
+The `MeshArray` type is a sub-type of `AbstractArray` with an `outer array` where each element is itself a 2D `inner array`. By default, outer and inner arrays are of all of the standard `Array` type. However, this setup potentially allows different choices for the outer and inner arrays – for example `DistributedArrays` and `AxisArrays`, respectively, could be an option. `MeshArrays.jl` thus provides a simple but general solution to analyze or e.g. simulate climate system variables. 
 
-The internals of a `MeshArray` are regulated by its `gcmgrid` -- a struct containing just a few index ranges, array size specifications, and inter-connection rules. A second  lightweight struct, `varmeta`, contains the `MeshArray` variable name, its unit, and position in grid space. A general approach like this is useful because climate models often involve advanced domain decompositions (see *Earth Model Grids*), and many variables, which puts a burden on users. 
+The internals of a `MeshArray` are regulated by its `gcmgrid` -- a struct containing just a few index ranges, array size specifications, and inter-connection rules. A second  lightweight struct, `varmeta`, contains the `MeshArray` variable name, its unit, time, and position in grid space. A general approach like this is useful because climate models often involve advanced domain decompositions (see, e.g., [Earth Grids](@ref)), and many variables, which can put a burden on users. 
 
-Encoding the grid specification inside the `MeshArray` data type allows user to manipulate `MeshArray`s just like they would manipulate `Array`s without having to keep track of model grid details. In addition, the provided `exchange` methods readily transfer data between connected subdomains to extend them at the sides. This makes it easy to compute e.g. partial derivatives and related operators like gradients, curl, or divergences over subdomain edges as often needed for precise computation of transports, budgets, etc using climate model output.
+Encoding the grid specification inside the `MeshArray` data type allows user to manipulate `MeshArray`s just like they would manipulate `Array`s without having to invoke model grid details explicitely. In addition, the provided `exchange` methods readily transfer data between connected subdomains to extend them at the sides. This makes it easy to compute e.g. partial derivatives and related operators like gradients, curl, or divergences over subdomain edges as often needed for precise computation of transports, budgets, etc using climate model output (see, e.g., [Tutorial](@ref)).
 
 ## Data Structures
 
@@ -83,27 +83,37 @@ julia> show(exchange(D))
 A `MeshArray` includes a `gcmgrid` specification which can be constructed as outlined below.
 
 ```
-gcmgrid(path::String, class::String, 
-        nFaces::Int, fSize::Array{NTuple{2, Int},1}, 
-        ioSize::Array{Int64,2}, ioPrec::Type, 
-        read::Function, write::Function)
+struct gcmgrid
+  path::String
+  class::String
+  nFaces::Int
+  fSize::Array{NTuple{2, Int},1}
+  ioSize::Union{NTuple{2, Int},Array{Int64,2}}
+  ioPrec::Type
+  read::Function
+  write::Function
+end
 ```
 
 Importantly, a `gcmgrid` does **not** contain any actual grid data -- hence its memory footprint is minimal. Grid variables are instead read to memory only when needed e.g. as shown below. To make this easy, each `gcmgrid` includes a pair of `read` / `write` methods to allow for basic `I/O` at any time. These methods are typically specified by the user although defaults are provided. 
 
 ```
 using MeshArrays, Unitful
 γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
-m=MeshArrays.varmeta(u"m",fill(0.5,2),"Depth","Depth")
+m=MeshArrays.varmeta(u"m",fill(0.5,2),missing,"Depth","Depth")
 D=γ.read(γ.path*"Depth.data",MeshArray(γ,Float64;meta=m))
 ```
 
 The above commands define a `MeshArray` called `D` which is the one displayed at the top of this section. A definition of the `varmeta` structure is reported below. The `position` of a `D` point within its grid cell is given as `x ∈ [0. 1.]` in each direction.
 
 ```
-varmeta(unit::Union{Unitful.AbstractQuantity,Number},
-        position::Array{Float64,1},
-        name::String,long_name::String)
+struct varmeta
+  unit::Union{Unitful.Units,Number,Missing}
+  position::Array{Float64,1}
+  time::Union{DateTime,Missing,Array{DateTime,1}}
+  name::String
+  long_name::String
+end
 ```
 
 ## Interpolation, Plots

docs/src/start.md

@@ -2,41 +2,76 @@
 
 ## Install
 
+To install `MeshArrays.jl` and verify it works as expected, open the `Julia` REPL and type:
+
 ```
 using Pkg
 Pkg.add("MeshArrays")
 Pkg.test("MeshArrays")
 ```
 
-!!! note
-    The `Julia` package manager, `Pkg.jl`, is documented [here](https://docs.julialang.org/en/v1/) and [further here](https://julialang.github.io/Pkg.jl/v1/).
+## Use
+
+To create your first `MeshArray`, open the `Julia` REPL and type:
+
+```
+using MeshArrays
+tmp=MeshArray(randn(20,10))
+show(tmp)
+```
 
 ## Tutorial
 
-The sequence of examples in this tutorial is as follows :
+In this tutorial, the same example is repeated three times -- for three different grid configurations commonly used in numerical models (see [Earth Grids](@ref)). The example proceeds as follows in each grid case:
 
-1. generate a grid configuration
-1. create a map of random noise
-1. apply diffusion (a transport process)
-1. plot results for each subdomain array
+1. generate grid configuration
+1. initialize map with random noise
+1. apply smoothing across whole domain
+1. plot results for the subdomain arrays
 
-These rely on grid configurations commonly used in global models (see [Earth Grids](@ref)). In step 3, the smoothing is based on integrating a lateral diffusion equation through time over the global domain. This illustrates how `MeshArrays.jl` computes partial derivatives between neighboring subdomains. 
+Step 3 illustrates how `MeshArrays.jl` computes partial derivatives between neighboring subdomains. The `MeshArrays.smooth()` method is indeed based on lateral diffusion, a transport process, which is integrated through time over the global domain. Other processes like e.g. advection by ocean currents would work similarly.
 
 ### 1. Doubly Periodic Domain
 
-Let's setup a doubly periodic domain with 16 subdomains. Each one contains an array of 20 by 20 grid points.
+Let's start with a doubly periodic domain split into `nF=16` subdomains. Each subdomain corresponds to an array of `nP=20` by `nQ=20` grid points. 
+
+```
+(nP,nQ,nF)=(20,20,16)
+facesSize=fill((nP,nQ),nF)
+ioSize=[nP nQ*nF]
+```
+
+The `UnitGrid()` function is used to generate such a grid configuration.
 
 ```
-using MeshArrays; p=dirname(pathof(MeshArrays))
-γ,Γ=GridOfOnes("PeriodicDomain",16,20)
+using MeshArrays
+γ=gcmgrid("","PeriodicDomain",
+			nF,facesSize, ioSize,
+			Float32, read, write)
+Γ=UnitGrid(γ)
+```
 
-include(joinpath(p,"../examples/Demos.jl"))
-(xi,xo,_,_)=demo2(Γ);
-show(xo)
+Then we do steps 2 (`zin`) and 3 (`zout`) as follows.
+
+```
+#initialize 2D field of random numbers
+tmp1=randn(Float32,Tuple(γ.ioSize))
+zin =γ.read(tmp1,MeshArray(γ,Float32))
+
+#smoothing length scales in x, y directions
+Lx=3*Γ.DXC; Ly=3*Γ.DYC
+
+#apply smoother
+zout=smooth(zin,Lx,Ly,Γ)
+```
+
+The `heatmap` function allows you to visualize that `zout` is indeed smoother (and therefore muted) than is `zin`.
 
+```
 using Plots
+p=dirname(pathof(MeshArrays))
 include(joinpath(p,"../examples/Plots.jl"))
-heatmap(xo,clims=(-0.25,0.25),colorbar=false,tickfont = (4, :black))
+heatmap(zout,clims=(-0.25,0.25),tickfont = (4, :black))
 ```
 
 Grid scale noise           |  Smoothed noise
@@ -45,21 +80,38 @@ Grid scale noise           |  Smoothed noise
 
 ### 2. Cube Sphere
 
-This grid has 6 subdomains, 100x100 points each, covering the six faces of a cube.
+Now we instead use a grid that has 6 subdomains, 100x100 points each, covering the six faces of a cube. We note that this _cube sphere_ topology involves connections between subdomain that are slightly more complicated than in the first example.
 
 ```
-γ,Γ=GridOfOnes("CubeSphere",6,100)
+(nP,nQ,nF)=(32,32,6)
+facesSize=fill((nP,nQ),nF)
+ioSize=[nP nQ*nF]
+
+using MeshArrays
+γ=gcmgrid("","CubeSphere",
+			nF,facesSize, ioSize,
+			Float32, read, write)
+Γ=UnitGrid(γ)
+```
+
+Next we call `demo2()` which combines steps 2 and 3. 
+
+```
+p=dirname(pathof(MeshArrays))
+include(joinpath(p,"../examples/Demos.jl"))
 D=demo2(Γ)
 ```
 
 ### 3. LLC90 grid
 
 The [Lat-Lon-Cap grid](http://www.geosci-model-dev.net/8/3071/2015/) (or LLC) is a global ocean model grid which is widely used in the [MITgcm user community](https://mitgcm.readthedocs.io/en/latest/). It has 5 uneven subdomains, variable grid spacing, and continents [(Forget et al 2015)](http://www.geosci-model-dev.net/8/3071/2015/). LLC90's resolution is one degree albeit with modications in the Arctic and along the Equator.
 
+In this case, the grid variables are read from files found in the `MeshArrays.GRID_LLC90` folder by the `GridLoad` function. The `GridSpec` function provides the corresponding domain sizes for this commonly used, `LLC90`, grid. 
+
 ```
-Γ=GridLoad(GridSpec("LatLonCap",MeshArrays.GRID_LLC90))
+γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
+Γ=GridLoad(γ)
 D=demo2(Γ)
-heatmap(D[2],clims=(-0.25,0.25))
 ```
 
 ## Earth Grids

src/Grids.jl

@@ -1,6 +1,40 @@
 
 Dict_to_NamedTuple(tmp::Dict) = (; zip(Symbol.(keys(tmp)), values(tmp))...)
 
+"""
+    UnitGrid(γ::gcmgrid)
+
+Generate a unit grid, where every grid spacing and area is 1, according to γ. 
+"""
+function UnitGrid(γ::gcmgrid)
+    nFaces=γ.nFaces
+    ioSize=(γ.ioSize[1],γ.ioSize[2])
+
+    Γ=Dict()
+    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)
+    pc=fill(0.5,2); pg=fill(0.0,2); pu=[0.,0.5]; pv=[0.5,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])
+
+    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));
+        Γ[list_n[ii]]=tmp1
+    end
+
+    for i in 1:nFaces
+        (np,nq)=γ.fSize[i]
+        Γ["XC"][i]=vec(0.5:1.0:np-0.5)*ones(1,nq)
+        Γ["XG"][i]=vec(0.0:1.0:np-1.0)*ones(1,nq)
+        Γ["YC"][i]=ones(np,1)*transpose(vec(0.5:1.0:nq-0.5))
+        Γ["YG"][i]=ones(np,1)*transpose(vec(0.0:1.0:nq-1.0))
+    end
+
+    return Dict_to_NamedTuple(Γ)
+
+end
+
 """
     simple_periodic_domain(np::Integer,nq=missing)
 
@@ -28,26 +62,7 @@ function simple_periodic_domain(np::Integer,nq=missing)
     ioPrec=Float32
     γ=gcmgrid("","PeriodicDomain",1,facesSize, ioSize, ioPrec, read, write)
 
-    Γ=Dict()
-    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)
-    pc=fill(0.5,2); pg=fill(0.0,2); pu=[0.,0.5]; pv=[0.5,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])
-    for ii=1:length(list_n);
-        tmp1=fill(1.,np,nq*nFaces);
-        m=varmeta(list_u[ii],list_p[ii],missing,list_n[ii],list_n[ii]);
-        tmp1=γ.read(tmp1,MeshArray(γ,Float64;meta=m));
-        tmp2=Symbol(list_n[ii]);
-        @eval (($tmp2) = ($tmp1))
-        Γ[list_n[ii]]=tmp1
-    end
-
-    Γ["XC"][1]=vec(0.5:1.0:np-0.5)*ones(1,nq)
-    Γ["XG"][1]=vec(0.0:1.0:np-1.0)*ones(1,nq)
-    Γ["YC"][1]=ones(np,1)*transpose(vec(0.5:1.0:nq-0.5))
-    Γ["YG"][1]=ones(np,1)*transpose(vec(0.0:1.0:nq-1.0))
-
-    return Dict_to_NamedTuple(Γ)
+    return UnitGrid(γ)
 end
 
 ## GridSpec function with default GridName argument:
@@ -240,22 +255,7 @@ function GridOfOnes(grTp,nF,nP)
 
     γ=gcmgrid(grDir,grTopo,nFaces,facesSize, ioSize, ioPrec, read, write)
 
-    Γ=Dict()
-    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)
-    pc=fill(0.5,2); pg=fill(0.0,2); pu=[0.,0.5]; pv=[0.5,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])
-
-    for ii=1:length(list_n);
-        tmp1=fill(1.,nP,nP*nF); m=varmeta(list_u[ii],list_p[ii],missing,list_n[ii],list_n[ii]);
-        tmp1=γ.read(tmp1,MeshArray(γ,Float64;meta=m));
-        tmp2=Symbol(list_n[ii]);
-        @eval (($tmp2) = ($tmp1))
-        Γ[list_n[ii]]=tmp1
-    end
-
-    return γ, Dict_to_NamedTuple(Γ)
-
+    return γ, UnitGrid(γ)
 end
 
 """

src/MeshArrays.jl

@@ -27,7 +27,7 @@ include("Interpolation.jl")
 
 export AbstractMeshArray, MeshArray, InnerArray, OuterArray, varmeta
 export gcmgrid, exchange, gradient, convergence, smooth, mask
-export simple_periodic_domain, GridSpec, GridLoad, GridOfOnes, GridAddWS!
+export UnitGrid, simple_periodic_domain, GridSpec, GridLoad, GridOfOnes, GridAddWS!
 export Tiles, Interpolate, InterpolationFactors, knn
 export ScalarPotential, VectorPotential, ThroughFlow
 export StereographicProjection, LatitudeCircles