Skip to content

Commit

Permalink
Merge pull request #127 from JuliaClimate/less_unicode
Browse files Browse the repository at this point in the history
Less unicode
  • Loading branch information
gaelforget authored Oct 25, 2024
2 parents 7f2ae7b + ecc7a21 commit 8e2e657
Show file tree
Hide file tree
Showing 9 changed files with 314 additions and 316 deletions.
8 changes: 4 additions & 4 deletions docs/src/examples.md
Original file line number Diff line number Diff line change
Expand Up @@ -13,19 +13,19 @@ To rerun an example yourself, the recommended method is to copy the correspondin

Simulate an ensemble of displacements (and trajectories) in a simple 2D configuration.

The [`FlowFields`](@ref) constructor readily wraps a flow field provided in the standard Array format, adds a time range, and returns a `FlowFields` data structure `𝐹`.
The [`FlowFields`](@ref) constructor readily wraps a flow field provided in the standard Array format, adds a time range, and returns a `FlowFields` data structure `F`.

All that is left to do at this stage is to define initial positions for the individuals, put them together with `𝐹` into the [`Individuals`](@ref) data structure `𝐼`, and call `∫!(𝐼)`.
All that is left to do at this stage is to define initial positions for the individuals, put them together with `F` into the [`Individuals`](@ref) data structure `I`, and call `∫!(I)`.

Exercises include the non-periodic domain case, statistics made easy via `DataFrames.jl`, and replacing the flow field with your own.

## Simple Three-Dimensional Flow

[notebook (html)](solid_body_rotation.html)[notebook (code)](https://github.com/JuliaClimate/IndividualDisplacements.jl/blob/master/examples/basics/solid_body_rotation.jl)

Set up a three-dimensional flow field `u,v,w`, initialize a single particle at position `📌`, and wrap everything up within an `Individuals` data structure `𝐼`.
Set up a three-dimensional flow field `u,v,w`, initialize a single particle at position `📌`, and wrap everything up within an `Individuals` data structure `I`.

`𝐼` is displaced by integrating the individual velocity, [moving along through space](https://en.wikipedia.org/wiki/Lagrangian_and_Eulerian_specification_of_the_flow_field), over time `𝑇`. This is the main computation done in this package -- interpolating `u,v,w` to individual positions `𝐼.📌` on the fly, using `𝐼.🚄`, and integrating through time, using `𝐼.∫`.
`I` is displaced by integrating the individual velocity, [moving along through space](https://en.wikipedia.org/wiki/Lagrangian_and_Eulerian_specification_of_the_flow_field), over time `T`. This is the main computation done in this package -- interpolating `u,v,w` to individual positions `I.📌` on the fly, using `I.🚄`, and integrating through time, using `I.∫`.

The flow field consists of [rigid body rotation](https://en.wikipedia.org/wiki/Rigid_body), plus a convergent term, plus a sinking term in the vertical direction. This flow field generates a downward, converging spiral -- a idealized version of a relevant case in the Ocean.

Expand Down
12 changes: 6 additions & 6 deletions docs/src/workflow.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,10 +2,10 @@

As shown in the [Examples](@ref) section, the typical worflow is:

1. set up a `FlowFields` data structure (`𝐹`)
1. set up `Individuals` (`𝐼`) with initial position `📌` and `𝐹`
1. displace `𝐼` by `solve!(𝐼,𝑇)` following `𝐼.𝐹` over `𝑇`
1. post-process by `𝐼.🔧` and record information in `𝐼.🔴`
1. set up a `FlowFields` data structure (`F`)
1. set up `Individuals` (`I`) with initial position `📌` and `F`
1. displace `I` by `solve!(I,T)` following `I.F` over `T`
1. post-process by `I.🔧` and record information in `I.🔴`
1. go back to `step 2` and continue if needed

The data structures for steps `1` and `2` are documented below. Steps `3` and `4` normally take place as part of `solve!` (i.e. `∫!`) which post-processes results, using 🔧, records them in 🔴, and finally updates the positions of individuals in 📌. Since 🔴 is a [DataFrame](https://juliadata.github.io/DataFrames.jl/latest/), it is easily manipulated, plotted, or saved in step `4` or after the fact.
Expand All @@ -28,7 +28,7 @@ For an overview of the examples, please refer to the **example guide**. The rest

## Data Structures

The `Individuals` struct contains a `FlowFields` struct (incl. e.g. arrays), initial positions for the individuals, and the other elements (e.g. functions) involved in `∫!(𝐼,𝑇)` as documented hereafter.
The `Individuals` struct contains a `FlowFields` struct (incl. e.g. arrays), initial positions for the individuals, and the other elements (e.g. functions) involved in `∫!(I,T)` as documented hereafter.

```@autodocs
Modules = [IndividualDisplacements]
Expand All @@ -37,7 +37,7 @@ Order = [:type]

## Main Functions

`∫!(𝐼,𝑇)` displaces individuals 𝐼 continuously over time period 𝑇 according to velocity function 🚄, temporal integration method ∫, and post-processor 🔧 (all embedded within 𝐼).
`∫!(I,T)` displaces individuals I continuously over time period T according to velocity function 🚄, temporal integration method ∫, and post-processor 🔧 (all embedded within I).

```@docs
∫!
Expand Down
102 changes: 51 additions & 51 deletions examples/worldwide/ECCO_FlowFields.jl
Original file line number Diff line number Diff line change
Expand Up @@ -19,34 +19,34 @@ function init_from_file(np ::Int)
end

"""
init_global_randn(np ::Int , 𝐷::NamedTuple)
init_global_randn(np ::Int , D::NamedTuple)
Randomly distribute `np` points over the Earth, within `𝐷.msk`
Randomly distribute `np` points over the Earth, within `D.msk`
region, and return position in grid index space (`i,j,subdomain`).
"""
function init_global_randn(np ::Int , 𝐷::NamedTuple)
function init_global_randn(np ::Int , D::NamedTuple)
(lon, lat) = randn_lonlat(maximum([2*np 10]))
(_,_,_,_,f,x,y)=InterpolationFactors(𝐷.Γ,lon,lat)
(_,_,_,_,f,x,y)=InterpolationFactors(D.Γ,lon,lat)
m=findall( (f.!==0).*((!isnan).(x)) )
n=findall(nearest_to_xy(𝐷.msk,x[m],y[m],f[m]).==1.0)[1:np]
n=findall(nearest_to_xy(D.msk,x[m],y[m],f[m]).==1.0)[1:np]
xyf=permutedims([x[m[n]] y[m[n]] f[m[n]]])
return DataFrame(x=xyf[1,:],y=xyf[2,:],f=xyf[3,:])
end

"""
reset_📌!(𝐼::Individuals,frac::Number,📌::Array)
reset_📌!(I::Individuals,frac::Number,📌::Array)
Randomly select a fraction (frac) of the particles and reset
their positions (𝐼.📌) to a random subset of the specified 📌.
their positions (I.📌) to a random subset of the specified 📌.
"""
function reset_📌!(𝐼::Individuals,frac::Number,📌::Array)
np=length(𝐼.🆔)
function reset_📌!(I::Individuals,frac::Number,📌::Array)
np=length(I.🆔)
n_reset = Int(round(frac*np))
k_reset = rand(1:np, n_reset)
l_reset = rand(1:np, n_reset)
𝐼.📌[k_reset]=deepcopy(📌[l_reset])
isempty(𝐼.🔴.ID) ? m=maximum(𝐼.🆔) : m=max(maximum(𝐼.🔴.ID),maximum(𝐼.🆔))
𝐼.🆔[k_reset]=collect(1:n_reset) .+ m
I.📌[k_reset]=deepcopy(📌[l_reset])
isempty(I.🔴.ID) ? m=maximum(I.🆔) : m=max(maximum(I.🔴.ID),maximum(I.🆔))
I.🆔[k_reset]=collect(1:n_reset) .+ m
end

"""
Expand All @@ -57,7 +57,7 @@ function `func` (e.g., `(u -> MeshArrays.update_location_llc!(u,Γ)))`,
and file location (`pth`).
_Note: the initial implementation approximates month durations to
365 days / 12 months for simplicity and sets 𝑃.𝑇 to [-mon/2,mon/2]_
365 days / 12 months for simplicity and sets 𝑃.T to [-mon/2,mon/2]_
"""
function setup_FlowFields(k::Int::NamedTuple,func::Function,pth::String)
XC=MeshArrays.exchange.XC) #add 1 lon point at each edge
Expand All @@ -80,27 +80,27 @@ function setup_FlowFields(k::Int,Γ::NamedTuple,func::Function,pth::String)
MeshArray(γ,Float32),MeshArray(γ,Float32),[-mon/2,mon/2],func)
end

𝐷 = (🔄 = update_FlowFields!, pth=pth,
D = (🔄 = update_FlowFields!, pth=pth,
XC=XC, YC=YC, iDXC=iDXC, iDYC=iDYC,
k=k, msk=msk, θ0=similar(msk), θ1=similar(msk))

𝐷 = merge(𝐷 , MeshArrays.NeighborTileIndices_cs(Γ))
D = merge(D , MeshArrays.NeighborTileIndices_cs(Γ))

return 𝑃,𝐷
return 𝑃,D
end

"""
update_FlowFields!(𝑃::𝐹_MeshArray2D,𝐷::NamedTuple,t::Float64)
update_FlowFields!(𝑃::F_MeshArray2D,D::NamedTuple,t::Float64)
Update flow field arrays (in 𝑃), 𝑃.𝑇, and ancillary variables (in 𝐷)
Update flow field arrays (in 𝑃), 𝑃.T, and ancillary variables (in D)
according to the chosen time `t` (in `seconds`).
_Note: for now, it is assumed that (1) the time interval `dt` between
consecutive records is diff(𝑃.𝑇), (2) monthly climatologies are used
consecutive records is diff(𝑃.T), (2) monthly climatologies are used
with a periodicity of 12 months, (3) vertical 𝑃.k is selected_
"""
function update_FlowFields!(𝑃::𝐹_MeshArray2D,𝐷::NamedTuple,t::AbstractFloat)
dt=𝑃.𝑇[2]-𝑃.𝑇[1]
function update_FlowFields!(𝑃::F_MeshArray2D,D::NamedTuple,t::AbstractFloat)
dt=𝑃.T[2]-𝑃.T[1]

m0=Int(floor((t+dt/2.0)/dt))
m1=m0+1
Expand All @@ -112,38 +112,38 @@ function update_FlowFields!(𝑃::𝐹_MeshArray2D,𝐷::NamedTuple,t::AbstractF
m1=mod(m1,12)
m1==0 ? m1=12 : nothing

(U,V)=read_velocities(𝑃.u0.grid,m0,𝐷.pth)
u0=U[:,𝐷.k]; v0=V[:,𝐷.k]
(U,V)=read_velocities(𝑃.u0.grid,m0,D.pth)
u0=U[:,D.k]; v0=V[:,D.k]
u0[findall(isnan.(u0))]=0.0; v0[findall(isnan.(v0))]=0.0 #mask with 0s rather than NaNs
u0=u0.*𝐷.iDXC; v0=v0.*𝐷.iDYC; #normalize to grid units
u0=u0.*D.iDXC; v0=v0.*D.iDYC; #normalize to grid units
(u0,v0)=MeshArrays.exchange(u0,v0,1) #add 1 point at each edge for u and v

(U,V)=read_velocities(𝑃.u0.grid,m1,𝐷.pth)
u1=U[:,𝐷.k]; v1=V[:,𝐷.k]
(U,V)=read_velocities(𝑃.u0.grid,m1,D.pth)
u1=U[:,D.k]; v1=V[:,D.k]
u1[findall(isnan.(u1))]=0.0; v1[findall(isnan.(v1))]=0.0 #mask with 0s rather than NaNs
u1=u1.*𝐷.iDXC; v1=v1.*𝐷.iDYC; #normalize to grid units
u1=u1.*D.iDXC; v1=v1.*D.iDYC; #normalize to grid units
(u1,v1)=MeshArrays.exchange(u1,v1,1) #add 1 point at each edge for u and v

𝑃.u0[:]=Float32.(u0[:])
𝑃.u1[:]=Float32.(u1[:])
𝑃.v0[:]=Float32.(v0[:])
𝑃.v1[:]=Float32.(v1[:])
𝑃.𝑇[:]=[t0,t1]
𝑃.T[:]=[t0,t1]

end

"""
update_FlowFields!(𝑃::𝐹_MeshArray3D,𝐷::NamedTuple,t::Float64)
update_FlowFields!(𝑃::F_MeshArray3D,D::NamedTuple,t::Float64)
Update flow field arrays (in 𝑃), 𝑃.𝑇, and ancillary variables (in 𝐷)
Update flow field arrays (in 𝑃), 𝑃.T, and ancillary variables (in D)
according to the chosen time `t` (in `seconds`).
_Note: for now, it is assumed that (1) the time interval `dt` between
consecutive records is diff(𝑃.𝑇), (2) monthly climatologies are used
consecutive records is diff(𝑃.T), (2) monthly climatologies are used
with a periodicity of 12 months, (3) vertical 𝑃.k is selected_
"""
function update_FlowFields!(𝑃::𝐹_MeshArray3D,𝐷::NamedTuple,t::Float64)
dt=𝑃.𝑇[2]-𝑃.𝑇[1]
function update_FlowFields!(𝑃::F_MeshArray3D,D::NamedTuple,t::Float64)
dt=𝑃.T[2]-𝑃.T[1]

m0=Int(floor((t+dt/2.0)/dt))
m1=m0+1
Expand All @@ -155,30 +155,30 @@ function update_FlowFields!(𝑃::𝐹_MeshArray3D,𝐷::NamedTuple,t::Float64)
m1=mod(m1,12)
m1==0 ? m1=12 : nothing

(_,nr)=size(𝐷.Γ.hFacC)
(_,nr)=size(D.Γ.hFacC)

(U,V)=read_velocities(𝑃.u0.grid,m0,𝐷.pth)
(U,V)=read_velocities(𝑃.u0.grid,m0,D.pth)
u0=U; v0=V
u0[findall(isnan.(u0))]=0.0; v0[findall(isnan.(v0))]=0.0 #mask with 0s rather than NaNs
for k=1:nr
u0[:,k]=u0[:,k].*𝐷.iDXC; v0[:,k]=v0[:,k].*𝐷.iDYC; #normalize to grid units
u0[:,k]=u0[:,k].*D.iDXC; v0[:,k]=v0[:,k].*D.iDYC; #normalize to grid units
(tmpu,tmpv)=exchange(u0[:,k],v0[:,k],1) #add 1 point at each edge for u and v
u0[:,k]=tmpu
v0[:,k]=tmpv
end
w0=IndividualDisplacements.read_nctiles(𝐷.pth*"WVELMASS/WVELMASS","WVELMASS",𝑃.u0.grid,I=(:,:,:,m0))
w0=IndividualDisplacements.read_nctiles(D.pth*"WVELMASS/WVELMASS","WVELMASS",𝑃.u0.grid,I=(:,:,:,m0))
w0[findall(isnan.(w0))]=0.0 #mask with 0s rather than NaNs

(U,V)=read_velocities(𝑃.u0.grid,m1,𝐷.pth)
(U,V)=read_velocities(𝑃.u0.grid,m1,D.pth)
u1=U; v1=V
u1[findall(isnan.(u1))]=0.0; v1[findall(isnan.(v1))]=0.0 #mask with 0s rather than NaNs
for k=1:nr
u1[:,k]=u1[:,k].*𝐷.iDXC; v1[:,k]=v1[:,k].*𝐷.iDYC; #normalize to grid units
u1[:,k]=u1[:,k].*D.iDXC; v1[:,k]=v1[:,k].*D.iDYC; #normalize to grid units
(tmpu,tmpv)=exchange(u1[:,k],v1[:,k],1) #add 1 point at each edge for u and v
u1[:,k]=tmpu
v1[:,k]=tmpv
end
w1=IndividualDisplacements.read_nctiles(𝐷.pth*"WVELMASS/WVELMASS","WVELMASS",𝑃.u0.grid,I=(:,:,:,m1))
w1=IndividualDisplacements.read_nctiles(D.pth*"WVELMASS/WVELMASS","WVELMASS",𝑃.u0.grid,I=(:,:,:,m1))
w1[findall(isnan.(w1))]=0.0 #mask with 0s rather than NaNs

𝑃.u0[:,:]=u0[:,:]
Expand All @@ -187,24 +187,24 @@ function update_FlowFields!(𝑃::𝐹_MeshArray3D,𝐷::NamedTuple,t::Float64)
𝑃.v1[:,:]=v1[:,:]
for k=1:nr
tmpw=exchange(-w0[:,k],1)
𝑃.w0[:,k]=tmpw./𝐷.Γ.DRC[k]
𝑃.w0[:,k]=tmpw./D.Γ.DRC[k]
tmpw=exchange(-w1[:,k],1)
𝑃.w1[:,k]=tmpw./𝐷.Γ.DRC[k]
𝑃.w1[:,k]=tmpw./D.Γ.DRC[k]
end
𝑃.w0[:,1]=0*exchange(-w0[:,1],1)
𝑃.w1[:,1]=0*exchange(-w1[:,1],1)
𝑃.w0[:,nr+1]=0*exchange(-w0[:,1],1)
𝑃.w1[:,nr+1]=0*exchange(-w1[:,1],1)

θ0=IndividualDisplacements.read_nctiles(𝐷.pth*"THETA/THETA","THETA",𝑃.u0.grid,I=(:,:,:,m0))
θ0=IndividualDisplacements.read_nctiles(D.pth*"THETA/THETA","THETA",𝑃.u0.grid,I=(:,:,:,m0))
θ0[findall(isnan.(θ0))]=0.0 #mask with 0s rather than NaNs
𝐷.θ0[:,:]=float32.(θ0[:,:])
D.θ0[:,:]=float32.(θ0[:,:])

θ1=IndividualDisplacements.read_nctiles(𝐷.pth*"THETA/THETA","THETA",𝑃.u0.grid,I=(:,:,:,m1))
θ1=IndividualDisplacements.read_nctiles(D.pth*"THETA/THETA","THETA",𝑃.u0.grid,I=(:,:,:,m1))
θ1[findall(isnan.(θ1))]=0.0 #mask with 0s rather than NaNs
𝐷.θ1[:,:]=float32.(θ1[:,:])
D.θ1[:,:]=float32.(θ1[:,:])

𝑃.𝑇[:]=[t0,t1]
𝑃.T[:]=[t0,t1]
end

"""
Expand Down Expand Up @@ -237,8 +237,8 @@ function global_ocean_circulation(;k=1)
func=(u -> MeshArrays.update_location_llc!(u,Γ))

#initialize u0,u1 etc
𝑃,𝐷=setup_FlowFields(k,Γ,func,ScratchSpaces.ECCO)
𝐷.🔄(𝑃,𝐷,0.0)
𝑃,D=setup_FlowFields(k,Γ,func,ScratchSpaces.ECCO)
D.🔄(𝑃,D,0.0)

#add background map for plotting
λ=ECCO_FlowFields.get_interp_coefficients(Γ)
Expand All @@ -249,9 +249,9 @@ function global_ocean_circulation(;k=1)

#add parameters for use in reset!
tmp=(frac=r_reset, Γ=Γ, ODL=ODL)
𝐷=merge(𝐷,tmp)
D=merge(D,tmp)

return 𝑃,𝐷
return 𝑃,D

end

Expand Down
18 changes: 9 additions & 9 deletions examples/worldwide/OCCA_FlowFields.jl
Original file line number Diff line number Diff line change
Expand Up @@ -95,10 +95,10 @@ function setup(;backward_in_time::Bool=false,nmax=Inf)

𝑃=FlowFields(u,u,v,v,w,w,[t0,t1],func)

𝐷 = (θ0=θ, θ1=θ, XC=MeshArrays.exchange.XC), YC=MeshArrays.exchange.YC),
D = (θ0=θ, θ1=θ, XC=MeshArrays.exchange.XC), YC=MeshArrays.exchange.YC),
RF=Γ.RF, RC=Γ.RC,ioSize=(360,160,n), Γ=Γ)

return 𝑃,𝐷
return 𝑃,D

end

Expand All @@ -125,9 +125,9 @@ custom🔴 = DataFrame(ID=Int[], fid=Int[], x=Float64[], y=Float64[],
lon=Float64[], lat=Float64[], dlon=Float64[], dlat=Float64[],
year=Float64[], col=Symbol[])

function custom🔧(sol,𝑃::𝐹_MeshArray3D,𝐷::NamedTuple;id=missing,𝑇=missing)
df=postprocess_MeshArray(sol,𝑃,𝐷,id=id,𝑇=𝑇)
add_lonlat!(df,𝐷.XC,𝐷.YC)
function custom🔧(sol,𝑃::F_MeshArray3D,D::NamedTuple;id=missing,T=missing)
df=postprocess_MeshArray(sol,𝑃,D,id=id,T=T)
add_lonlat!(df,D.XC,D.YC)
df.dlon=0*df.lon
df.dlat=0*df.lat

Expand All @@ -136,14 +136,14 @@ function custom🔧(sol,𝑃::𝐹_MeshArray3D,𝐷::NamedTuple;id=missing,𝑇=

#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(𝐷.RC)
nz=length(D.RC)
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
df.z=D.RF[1 .+ k].*(1 .- w)+D.RF[2 .+ k].*w #depth

#add one isotherm depth
θ=0.5*(𝐷.θ0+𝐷.θ1)
d=MeshArrays.isosurface(θ,15,𝐷)
θ=0.5*(D.θ0+D.θ1)
d=MeshArrays.isosurface(θ,15,D)
d[findall(isnan.(d))].=0.
df.iso=interp_to_xy(df,MeshArrays.exchange(d));

Expand Down
Loading

0 comments on commit 8e2e657

Please sign in to comment.