Merge pull request #191 from gaelforget/TimeAxis_dev

add TimeAxis (new) to FlowFields, define TimePeriod (mutable)

Author: gaelforget

src/API.jl

@@ -3,6 +3,78 @@
 
 using Dates
 
+
+abstract type AbstractTimeAxis end
+abstract type AbstractTimePeriod end
+
+"""
+    struct TimeAxis{Ty} <: AbstractTimeAxis
+
+Time axis specification for `FlowFields` data structure. This un-mutable specification
+describes the overall range of simulations that can be performed using the `FlowFields`.
+
+- initial : beginning of simulation period
+- final : end of simulation period
+- origin : beginning of the valid time period, when flow fields are available
+- horizon : end of the valid time period, when flow fields are available
+- climatology : false by default ; but can be set to true for time-periodic flow fields.
+
+Rules : 
+
+- In most simple cases `origin,horizon` is just the same as `initial,final`.
+- By default `initial,final` and `origin,horizon` are set to the `FlowField`'s `TimePeriod`.
+- However, any specification where `initial,final` is within `origin,horizon` is correct.
+- This more general approach makes it easy to split simulations into smaller intervals.
+
+```
+using Drifters, Dates
+
+D0=DateTime(2000,1,1)
+D1=DateTime(2000,3,1)
+TA=Drifters.TimeAxis(D0,D1)
+TA=Drifters.TimeAxis(D0,D1,D0,D1,false)
+```
+"""
+struct TimeAxis{Ty} <: AbstractTimeAxis
+    initial::Ty
+    final::Ty
+    origin::Ty
+    horizon::Ty
+    Climatology::Bool
+end
+
+TimeAxis(t0=0.0,t1=1.0) = TimeAxis(t0,t1,t0,t1,false)
+
+"""
+    struct TimePeriod{Ty} <: AbstractTimePeriod
+
+Time period for next simulation of `Individuals` displacement specified 
+as a vector `[initial,final]`. The `value` field is a mutable vector that
+can be changed, e.g. within a time loop when carrying out longer simulations.
+
+Rules : 
+
+- if `final` precedes `initial` then the simulation goes backward in time
+- in most simple cases, `initial,final` are the same in `TimePeriod` and `TimeAxis`.
+- but any `initial,final` within the periods defined in 'TimeAxis` works for `TimePeriod`.
+
+```
+using Drifters, Dates
+
+D0=DateTime(2000,1,1)
+D1=DateTime(2000,3,1)
+TP=Drifters.TimePeriod([D0,D1])
+TA=Drifters.TimeAxis(TP)
+```
+"""
+struct TimePeriod{Ty} <: AbstractTimePeriod
+    value::Array{Ty,1}
+end
+
+TimePeriod(T0,T1) = TimePeriod([T0,T1])
+
+TimeAxis(TP::TimePeriod) = TimeAxis(TP.value...)
+
 """
     abstract type FlowFields
 
@@ -53,17 +125,20 @@ struct uvArrays{Ty} <: FlowFields
     v0::Array{Ty,2}
     v1::Array{Ty,2}
     T::Union{Array{Ty},Array{DateTime}}
+    TimeAxis::AbstractTimeAxis
 end
 
 function FlowFields(u0::Array{Ty,2},u1::Array{Ty,2},
-    v0::Array{Ty,2},v1::Array{Ty,2},T::Union{Array,Tuple}) where Ty
+    v0::Array{Ty,2},v1::Array{Ty,2},T::Union{Array,Tuple}; 
+    time_axis=missing) where Ty
     #ensure T is a vector and enforce type
     T=time_set_type.(collect(T),Ty)
+    TA=(ismissing(time_axis) ? TimeAxis(T...) : time_axis)
     #check array size concistency
     tst=prod([(size(u0)==size(tmp)) for tmp in (u1,v0,v1)])
     !tst ? error("inconsistent array sizes") : nothing
     #call constructor
-    uvArrays(u0,u1,v0,v1,T)
+    uvArrays(u0,u1,v0,v1,T,TA)
 end
 
 struct uvwArrays{Ty} <: FlowFields
@@ -74,6 +149,7 @@ struct uvwArrays{Ty} <: FlowFields
     w0::Array{Ty,3}
     w1::Array{Ty,3}
     T::Union{Array{Ty},Array{DateTime}}
+    TimeAxis::AbstractTimeAxis
 end
 
 """
@@ -90,7 +166,7 @@ F=FlowFields(u=uC,v=vC,period=(0,10.))
 ```
 """
 function FlowFields(; u::Union{Array,Tuple}=[], v::Union{Array,Tuple}=[], w::Union{Array,Tuple}=[], 
-    period::Union{Array,Tuple}=[])
+    period::Union{Array,Tuple}=[], time_axis=missing)
     (isa(u,Tuple)||length(u[:])==2) ? (u0=u[1]; u1=u[2]) : (u0=u; u1=u)
     (isa(v,Tuple)||length(v[:])==2) ? (v0=v[1]; v1=v[2]) : (v0=v; v1=v)
     (isa(w,Tuple)||length(w[:])==2) ? (w0=w[1]; w1=w[2]) : (w0=w; w1=w)
@@ -99,9 +175,9 @@ function FlowFields(; u::Union{Array,Tuple}=[], v::Union{Array,Tuple}=[], w::Uni
     end
     if !isempty(u0) && !isempty(v0)
         if !isempty(w0)
-            FlowFields(u0,u1,v0,v1,w0,w1,period)
+            FlowFields(u0,u1,v0,v1,w0,w1,period,time_axis=time_axis)
         else
-            FlowFields(u0,u1,v0,v1,period)
+            FlowFields(u0,u1,v0,v1,period,time_axis=time_axis)
         end
     else
         []
@@ -139,15 +215,16 @@ to_C_grid!(x;dims=0) = begin
 end
 
 function FlowFields(u0::Array{Ty,3},u1::Array{Ty,3},v0::Array{Ty,3},v1::Array{Ty,3},
-    w0::Array{Ty,3},w1::Array{Ty,3},T::Union{Array,Tuple}) where Ty
+    w0::Array{Ty,3},w1::Array{Ty,3},T::Union{Array,Tuple}; time_axis=missing) where Ty
     #ensure T is a vector and enforce type
     T=time_set_type.(collect(T),Ty)
+    TA=(ismissing(time_axis) ? TimeAxis(T...) : time_axis)
     #check array size concistency
     tst=prod([(size(u0)==size(tmp)) for tmp in (u1,v0,v1)])
     tst=tst*prod([(size(u0)==size(tmp).-(0,0,1)) for tmp in (w0,w1)])
     !tst ? error("inconsistent array sizes") : nothing
     #call constructor
-    uvwArrays(u0,u1,v0,v1,w0,w1,T)
+    uvwArrays(u0,u1,v0,v1,w0,w1,T,TA)
 end
 
 struct uvMeshArrays{Ty} <: FlowFields
@@ -156,26 +233,28 @@ struct uvMeshArrays{Ty} <: FlowFields
     v0::AbstractMeshArray{Ty,1}
     v1::AbstractMeshArray{Ty,1}
     T::Union{Array{Ty},Array{DateTime}}
+    TimeAxis::AbstractTimeAxis
     update_location!::Function
 end
 
 function FlowFields(u0::MeshArrays.MeshArray_wh,u1::MeshArrays.MeshArray_wh,
     v0::MeshArrays.MeshArray_wh,v1::MeshArrays.MeshArray_wh,
-    T::Union{Array,Tuple},update_location!::Function)
+    T::Union{Array,Tuple},update_location!::Function; time_axis=missing)
 
-    FlowFields(u0.MA,u1.MA,v0.MA,v1.MA,T,update_location!)
+    FlowFields(u0.MA,u1.MA,v0.MA,v1.MA,T,update_location!,time_axis=time_axis)
 end
 
 function FlowFields(u0::AbstractMeshArray{Ty,1},u1::AbstractMeshArray{Ty,1},
     v0::AbstractMeshArray{Ty,1},v1::AbstractMeshArray{Ty,1},
-    T::Union{Array,Tuple},update_location!::Function) where Ty
+    T::Union{Array,Tuple},update_location!::Function; time_axis=missing) where Ty
     #ensure T is a vector and enforce type
     T=time_set_type.(collect(T),Ty)
+    TA=(ismissing(time_axis) ? TimeAxis(T...) : time_axis)
     #check array size concistency
     tst=prod([(size(u0)==size(tmp))*(u0.fSize==tmp.fSize) for tmp in (u1,v0,v1)])
     !tst ? error("inconsistent array sizes") : nothing
     #call constructor
-    uvMeshArrays(u0,u1,v0,v1,T,update_location!)
+    uvMeshArrays(u0,u1,v0,v1,T,TA,update_location!)
 end
 
 struct uvwMeshArrays{Ty} <: FlowFields
@@ -186,30 +265,32 @@ struct uvwMeshArrays{Ty} <: FlowFields
     w0::AbstractMeshArray{Ty,2}
     w1::AbstractMeshArray{Ty,2}
     T::Union{Array{Ty},Array{DateTime}}
+    TimeAxis::AbstractTimeAxis
     update_location!::Function
 end
 
 function FlowFields(u0::MeshArrays.MeshArray_wh,u1::MeshArrays.MeshArray_wh,
     v0::MeshArrays.MeshArray_wh,v1::MeshArrays.MeshArray_wh,
     w0::MeshArrays.MeshArray_wh,w1::MeshArrays.MeshArray_wh,
-    T::Union{Array,Tuple},update_location!::Function)
+    T::Union{Array,Tuple},update_location!::Function; time_axis=missing)
 
-    FlowFields(u0.MA,u1.MA,v0.MA,v1.MA,w0.MA,w1.MA,T,update_location!)
+    FlowFields(u0.MA,u1.MA,v0.MA,v1.MA,w0.MA,w1.MA,T,update_location!,time_axis=time_axis)
 end
 
 function FlowFields(u0::AbstractMeshArray{Ty,2},u1::AbstractMeshArray{Ty,2},
     v0::AbstractMeshArray{Ty,2},v1::AbstractMeshArray{Ty,2},
     w0::AbstractMeshArray{Ty,2},w1::AbstractMeshArray{Ty,2},
-    T::Union{Array,Tuple},update_location!::Function) where Ty
+    T::Union{Array,Tuple},update_location!::Function; time_axis=missing) where Ty
     #ensure T is a vector and enforce type
     T=time_set_type.(collect(T),Ty)
+    TA=(ismissing(time_axis) ? TimeAxis(T...) : time_axis)
     #check array size consistency
     tst=prod([(size(u0)==size(tmp))*(u0.fSize==tmp.fSize) for tmp in (u1,v0,v1)])
     tst=tst*prod([(size(u0)==size(tmp).-(0,1))*(u0.fSize==tmp.fSize) for tmp in (w0,w1)])
     !tst ? error("inconsistent array sizes") : nothing
 
     #call constructor
-    uvwMeshArrays(u0,u1,v0,v1,w0,w1,T,update_location!)
+    uvwMeshArrays(u0,u1,v0,v1,w0,w1,T,TA,update_location!)
 end
 
 """
@@ -230,7 +311,7 @@ function ensemble_solver(prob::ODEProblem;solver=Tsit5(),reltol=1e-8,abstol=1e-8
 end
 
 a=fill(0.0,1,1)
-default_flowfields = uvArrays{Float64}(a,a,a,a,[0. 1.])
+default_flowfields = uvArrays{Float64}(a,a,a,a,[0. 1.],TimeAxis(0.,1.))
 default_recorder = DataFrame(ID=Int[], x=Float64[], y=Float64[], t=Float64[])
 default_postproc = (x->x)
 

src/data_wrangling.jl

@@ -30,7 +30,7 @@ function convert_to_FlowFields(U::Array{T,2},V::Array{T,2},t1::T; time_option=:d
         error("unknown time_option")
     end
     )
-    uvMeshArrays{eltype(u.MA)}(u.MA,u.MA,v.MA,v.MA,TT,func)
+    uvMeshArrays{eltype(u.MA)}(u.MA,u.MA,v.MA,v.MA,TT,TimeAxis(TT...),func)
 end
 
 """

src/example_ECCO.jl

@@ -5,6 +5,7 @@ import Dates: DateTime, Year, Month, Day
 
 import Drifters: postprocess_MeshArray, add_lonlat!, OrdinaryDiffEq
 import Drifters: time_in_seconds, time_in_DateTime, monthly_records
+import Drifters: TimeAxis, TimePeriod
 
 import OrdinaryDiffEq: solve, Tsit5, ODEProblem
 import Drifters: update_location!, Individuals, uvMeshArrays, uvwMeshArrays
@@ -67,7 +68,15 @@ function setup_FlowFields(k::Int,Γ::NamedTuple,func::Function,pth::String;
 
     mon=86400.0*365.0/12.0    
     T=(time_unit==:DateTime ? [DateTime(1999,12,15),DateTime(2000,1,15)] : [-mon/2,mon/2])
-    
+    if time_unit==:DateTime
+        D0=DateTime(1992,1,1)
+        D1=DateTime(2011,1,1) #this corresponds to ECCO4 release 1
+        TA=TimeAxis(D0,D1,D0,D1,true) #true corresponds to monthly climatology input
+    else
+        t_final=100*365*86400.0
+        TA=TimeAxis(0.0,t_final,0.0,t_final,true)
+    end
+
     if k==0
         msk=Γ.hFacC
         msk=1.0*(msk .> 0.0)
@@ -76,14 +85,14 @@ function setup_FlowFields(k::Int,Γ::NamedTuple,func::Function,pth::String;
         P=FlowFields(exchange(MeshArray(γ,Float32,nr)),exchange(MeshArray(γ,Float32,nr)),
             exchange(MeshArray(γ,Float32,nr)),exchange(MeshArray(γ,Float32,nr)),
             exchange(MeshArray(γ,Float32,nr+1)),exchange(MeshArray(γ,Float32,nr+1)),
-            T,func)
+            T,func,time_axis=TA)
     else
         msk=Γ.hFacC[:, k]
         msk=1.0*(msk .> 0.0)
         exmsk=exchange(msk).MA
         P=FlowFields(exchange(MeshArray(γ,Float32)),exchange(MeshArray(γ,Float32)),
             exchange(MeshArray(γ,Float32)),exchange(MeshArray(γ,Float32)),
-            T,func)
+            T,func,time_axis=TA)
     end
 
     D = (🔄 = update_FlowFields!, pth=pth, datasets=datasets,
@@ -326,7 +335,7 @@ function init_FlowFields(; k=1, backward_time=false, time_unit=:DateTime, dpth::
 
   #initialize u0,u1 etc arrays
   P,D=setup_FlowFields(k,Γ,func,dpth,
-        backward_time=backward_time, time_unit=time_unit,datasets=datasets)
+        backward_time=backward_time,time_unit=time_unit,datasets=datasets)
 
   #add background map for plotting
   λ=get_interp_coefficients(Γ)

src/example_OCCA.jl

@@ -2,6 +2,7 @@ module OCCA
 
 import Drifters: data_path, uvwMeshArrays, FlowFields
 import Drifters: postprocess_MeshArray, add_lonlat!, interp_to_xy
+import Drifters: TimeAxis, TimePeriod
 
 ##
 
@@ -57,14 +58,12 @@ function setup(;backward_in_time::Bool=false,nmax=Inf)
    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)
    pth=data_path(:OCCA)
 
-   u=s*read(read_data("u",pth,n),MeshArray(γ,Float32,n))
-   v=s*read(read_data("v",pth,n),MeshArray(γ,Float32,n))
+   u=read(read_data("u",pth,n),MeshArray(γ,Float32,n))
+   v=read(read_data("v",pth,n),MeshArray(γ,Float32,n))
 
-   w=s*read_data("w",pth,n)
+   w=read_data("w",pth,n)
    w=-cat(w,zeros(360, 160),dims=3)
    w[:,:,1] .=0.0
    w=read(w,MeshArray(γ,Float32,n+1))
@@ -98,12 +97,16 @@ function setup(;backward_in_time::Bool=false,nmax=Inf)
    Γ.XG[1]=tmpx
    Γ.Depth[1]=circshift(Γ.Depth[1],[-180 0])
 
-   t0=0.0; t1=86400*366*2.0;
-
    (u,v)=exchange(u,v)
    w=exchange(w)
 
-   P=FlowFields(u,u,v,v,w,w,[t0,t1],func)
+   backward_in_time ? s=-1.0 : s=1.0
+   t0=0.0; t1=s*86400*366*2.0;
+
+   t_final=100*365*86400.0
+   TA=TimeAxis(-t_final,t_final,0.0,s*t_final,true)
+
+   P=FlowFields(u,u,v,v,w,w,[t0,t1],func,time_axis=TA)
 
    XC=exchange(Γ.XC)
    YC=exchange(Γ.YC)