From b62e39a541294ff9864ba5c460d6b47880e4b9d8 Mon Sep 17 00:00:00 2001 From: gaelforget Date: Mon, 10 Feb 2020 13:31:02 -0500 Subject: [PATCH 1/9] in `ex_2_more,jl` explain `ODEProblem(comp_vel..`; use OrdinaryDiffEq --- examples/body_rotation_test.jl | 2 +- examples/ex_2_more.jl | 20 +++++++++++++------- examples/flt_xmpl_single.jl | 2 +- 3 files changed, 15 insertions(+), 9 deletions(-) diff --git a/examples/body_rotation_test.jl b/examples/body_rotation_test.jl index e82d3f3c..ee9a4648 100644 --- a/examples/body_rotation_test.jl +++ b/examples/body_rotation_test.jl @@ -21,7 +21,7 @@ # ## 1. import software and `pkg/flt` trajectories -using IndividualDisplacements, MeshArrays, Plots, DifferentialEquations +using IndividualDisplacements, MeshArrays, OrdinaryDiffEq, Plots # ## 2. Define gridded variables as `MeshArray`s diff --git a/examples/ex_2_more.jl b/examples/ex_2_more.jl index b1260801..f7dc2df4 100644 --- a/examples/ex_2_more.jl +++ b/examples/ex_2_more.jl @@ -18,13 +18,13 @@ # - [x] read gridded output + `float_trajectory*data` => `uvetc` dictionnary. # - [x] recompute `u,v` from gridded output # - [x] compare with `u,v` from `float_traj*data` -# - [x] solve for float trajectory with `DifferentialEquations.jl` +# - [x] solve for trajectory with `OrdinaryDiffEq.jl` (part of `DifferentialEquations.jl`) # # _Notes:_ For documentation see , and # ## 1. import software -using IndividualDisplacements, MeshArrays, DifferentialEquations, Plots +using IndividualDisplacements, MeshArrays, OrdinaryDiffEq, Plots p=dirname(pathof(IndividualDisplacements)) include(joinpath(p,"plot_pyplot.jl")) @@ -32,7 +32,7 @@ include(joinpath(p,"plot_pyplot.jl")) dirIn="flt_example/" prec=Float32 -df=IndividualDisplacements.ReadDisplacements(dirIn,prec) +df=ReadDisplacements(dirIn,prec) #function exported by IndividualDisplacements PyPlot.figure(); PlotBasic(df,300,100000.0) # ## 3. Read gridded variables via `MeshArrays.jl` @@ -120,8 +120,8 @@ colorbar() # ## 6. Recompute displacements from gridded flow fields # + -comp_vel=IndividualDisplacements.VelComp -get_vel=IndividualDisplacements.VelCopy +comp_vel=VelComp #function exported by IndividualDisplacements +get_vel=VelCopy #function exported by IndividualDisplacements uInit=[tmp[1,:lon];tmp[1,:lat]]./uvetc["dx"] nSteps=Int32(tmp[end,:time]/3600)-2 @@ -180,9 +180,15 @@ Plots.plot!(refv) # ## 6. Recompute trajectories from gridded flow fields # -# Solve through time using `DifferentialEquations.jl` +# Solve through time using `OrdinaryDiffEq.jl` with +# +# - `comp_vel` is the function computing `du/dt` +# - `uInit` is the initial condition `u @ tspan[1]` +# - `tspan` is the time interval +# - `uvetc` are parameters for `comp_vel` +# - `Tsit5` is the time-stepping scheme +# - `reltol` and `abstol` are tolerance parameters -using DifferentialEquations tspan = (0.0,nSteps*3600.0) #prob = ODEProblem(get_vel,uInit,tspan,tmp) prob = ODEProblem(comp_vel,uInit,tspan,uvetc) diff --git a/examples/flt_xmpl_single.jl b/examples/flt_xmpl_single.jl index 572cf504..c2709af0 100644 --- a/examples/flt_xmpl_single.jl +++ b/examples/flt_xmpl_single.jl @@ -19,7 +19,7 @@ # ## 1. import software -using IndividualDisplacements, MeshArrays, DifferentialEquations, Plots, Statistics +using IndividualDisplacements, MeshArrays, OrdinaryDiffEq, Plots, Statistics p=dirname(pathof(IndividualDisplacements)); include(joinpath(p,"plot_pyplot.jl")) # ## 2. Read gridded variables as `MeshArray`s From d05040e1972078d1f877a825f5c67dd3fcf45aaa Mon Sep 17 00:00:00 2001 From: gaelforget Date: Mon, 10 Feb 2020 19:36:08 -0500 Subject: [PATCH 2/9] streamline includes --- examples/body_rotation_test.jl | 2 -- examples/ex_2_more.jl | 1 - examples/flt_xmpl_fleet.jl | 4 ++-- examples/flt_xmpl_single.jl | 1 - examples/global_ocean_fleet.jl | 5 ++--- examples/global_ocean_single.jl | 5 ++--- 6 files changed, 6 insertions(+), 12 deletions(-) diff --git a/examples/body_rotation_test.jl b/examples/body_rotation_test.jl index ee9a4648..a7db7585 100644 --- a/examples/body_rotation_test.jl +++ b/examples/body_rotation_test.jl @@ -65,7 +65,6 @@ du=fill(0.0,2); # ## 4. solve through time using DifferentialEquations.jl -using DifferentialEquations tspan = (0.0,nSteps*3600.0) prob = ODEProblem(IndividualDisplacements.VelComp,uInit,tspan,uvetc) sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8) @@ -73,7 +72,6 @@ sol[1:4] # ## 5. visualize trajectory -using Plots Plots.plot(sol[1,:],sol[2,:],linewidth=5,title="Solid body rotation example", xaxis="lon",yaxis="lat",label="Julia Solution") # legend=false #Plots.savefig("SolidBodyRotation.png") diff --git a/examples/ex_2_more.jl b/examples/ex_2_more.jl index f7dc2df4..0d89d123 100644 --- a/examples/ex_2_more.jl +++ b/examples/ex_2_more.jl @@ -210,7 +210,6 @@ for i=1:nSteps-1 end ref=ref./uvetc["dx"] -using Plots Plots.plot(sol[1,:],sol[2,:],linewidth=5,title="Using Recomputed Velocities", xaxis="lon",yaxis="lat",label="Julia Solution") # legend=false Plots.plot!(ref[1,:],ref[2,:],lw=3,ls=:dash,label="MITgcm Solution") diff --git a/examples/flt_xmpl_fleet.jl b/examples/flt_xmpl_fleet.jl index 50ad1c6b..b36139dc 100644 --- a/examples/flt_xmpl_fleet.jl +++ b/examples/flt_xmpl_fleet.jl @@ -19,7 +19,8 @@ # ## 1. import software -using IndividualDisplacements, MeshArrays, DifferentialEquations, Plots, Statistics +using IndividualDisplacements, MeshArrays, OrdinaryDiffEq +using Plots, Statistics, MITgcmTools, DataFrames p=dirname(pathof(IndividualDisplacements)) include(joinpath(p,"plot_pyplot.jl")) @@ -120,7 +121,6 @@ Plots.plot!(tmpv) # ## 5. Solve through time using `DifferentialEquations.jl` -using DifferentialEquations tspan = (0.0,nSteps*3600.0) #prob = ODEProblem(get_vel,uInit,tspan,tmp) prob = ODEProblem(comp_vel,uInit,tspan,uvetc) diff --git a/examples/flt_xmpl_single.jl b/examples/flt_xmpl_single.jl index c2709af0..05a48d27 100644 --- a/examples/flt_xmpl_single.jl +++ b/examples/flt_xmpl_single.jl @@ -119,7 +119,6 @@ Plots.plot!(tmpv) # ## 5. Solve through time using `DifferentialEquations.jl` -using DifferentialEquations tspan = (0.0,nSteps*3600.0) #prob = ODEProblem(get_vel,uInit,tspan,tmp) prob = ODEProblem(comp_vel,uInit,tspan,uvetc) diff --git a/examples/global_ocean_fleet.jl b/examples/global_ocean_fleet.jl index 9cc835a4..82a58cc9 100644 --- a/examples/global_ocean_fleet.jl +++ b/examples/global_ocean_fleet.jl @@ -19,8 +19,8 @@ # ## 1. import software -using IndividualDisplacements, MeshArrays, DifferentialEquations -using Plots, Statistics, Dierckx +using IndividualDisplacements, MeshArrays, OrdinaryDiffEq +using Plots, Statistics, Dierckx, MAT p=dirname(pathof(IndividualDisplacements)) include(joinpath(p,"plot_pyplot.jl")) @@ -42,7 +42,6 @@ GridVariables=Dict("XC" => read(mygrid.path*"XC.latlon.data",MeshArray(mygrid,Fl # Read velocity fields as `MeshArray`s. # + -using MAT file = matopen(mygrid.path*"uv_lonlat.mat") u=read(file, "u") v=read(file, "v") diff --git a/examples/global_ocean_single.jl b/examples/global_ocean_single.jl index 8d6ec3f8..0895795c 100644 --- a/examples/global_ocean_single.jl +++ b/examples/global_ocean_single.jl @@ -19,7 +19,8 @@ # ## 1. import software -using IndividualDisplacements, MeshArrays, DifferentialEquations, Plots, Statistics +using IndividualDisplacements, MeshArrays, OrdinaryDiffEq +using Plots, Statistics, MAT p=dirname(pathof(IndividualDisplacements)) include(joinpath(p,"plot_pyplot.jl")) @@ -41,7 +42,6 @@ GridVariables=Dict("XC" => read(mygrid.path*"XC.latlon.data",MeshArray(mygrid,Fl # Read velocity fields as `MeshArray`s. # + -using MAT file = matopen(mygrid.path*"uv_lonlat.mat") u=read(file, "u") v=read(file, "v") @@ -135,7 +135,6 @@ Plots.plot!(tmpv) # ## 5. Solve through time using `DifferentialEquations.jl` -using DifferentialEquations tspan = (0.0,nSteps*3600.0) #prob = ODEProblem(get_vel,uInit,tspan,tmp) prob = ODEProblem(comp_vel,uInit,tspan,uvetc) From 56a273214bc0ecf4f93ff5af5fe567a23f114dec Mon Sep 17 00:00:00 2001 From: gaelforget Date: Mon, 10 Feb 2020 19:48:49 -0500 Subject: [PATCH 3/9] streamline test --- test/runtests.jl | 36 ++++++++++-------------------------- test/test1_setup.jl | 21 +++++++++++++++++++++ 2 files changed, 31 insertions(+), 26 deletions(-) create mode 100644 test/test1_setup.jl diff --git a/test/runtests.jl b/test/runtests.jl index 1054f04b..9f45939e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,35 +1,19 @@ using Test using IndividualDisplacements, MeshArrays, OrdinaryDiffEq +include("test1_setup.jl") @testset "MeshArrays tests:" begin -mygrid=gcmgrid("flt_example/","ll",1,[(80,42)], [80 42], Float32, read, write) -XC=MeshArray(mygrid,Float32); XC[1]=vec(2500.:5000.:397500.0)*ones(1,42); -XG=MeshArray(mygrid,Float32); XG[1]=vec(0.:5000.:395000.0)*ones(1,42); -YC=MeshArray(mygrid,Float32); YC[1]=ones(80,1)*transpose(vec(2500.:5000.:207500.0)); -YG=MeshArray(mygrid,Float32); YG[1]=ones(80,1)*transpose(vec(0.:5000.:205000.0)); -GridVariables=Dict("XC" => XC,"YC" => YC,"XG" => XG,"YG" => YG,"dx" => 5000.0); + uvetc=test1_setup() -t0=0.0; t1=18001.0*3600.0 -u0=-(YG.-YC[1][40,21])/2000000.; u1=u0 -v0=(XG.-XC[1][40,21])/2000000.; v1=v0 + uInit=[200000.0;0.0]./uvetc["dx"] + nSteps=3000-2 + du=fill(0.0,2); + tspan = (0.0,nSteps*3600.0) + prob = ODEProblem(VelComp,uInit,tspan,uvetc) + sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8) -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); - -uInit=[200000.0;0.0]./uvetc["dx"] -nSteps=3000-2 -du=fill(0.0,2); -tspan = (0.0,nSteps*3600.0) -prob = ODEProblem(IndividualDisplacements.VelComp,uInit,tspan,uvetc) -sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8) - -@test isapprox(sol[1,end],117237.0./uvetc["dx"]; atol=100.) -@test isapprox(sol[2,end],40448.0./uvetc["dx"]; atol=100.) + @test isapprox(sol[1,end],117237.0./uvetc["dx"]; atol=100.) + @test isapprox(sol[2,end],40448.0./uvetc["dx"]; atol=100.) end diff --git a/test/test1_setup.jl b/test/test1_setup.jl new file mode 100644 index 00000000..dc060ce7 --- /dev/null +++ b/test/test1_setup.jl @@ -0,0 +1,21 @@ +function test1_setup() + + mygrid=gcmgrid("flt_example/","ll",1,[(80,42)], [80 42], Float32, read, write) + XC=MeshArray(mygrid,Float32); XC[1]=vec(2500.:5000.:397500.0)*ones(1,42); + XG=MeshArray(mygrid,Float32); XG[1]=vec(0.:5000.:395000.0)*ones(1,42); + YC=MeshArray(mygrid,Float32); YC[1]=ones(80,1)*transpose(vec(2500.:5000.:207500.0)); + YG=MeshArray(mygrid,Float32); YG[1]=ones(80,1)*transpose(vec(0.:5000.:205000.0)); + GridVariables=Dict("XC" => XC,"YC" => YC,"XG" => XG,"YG" => YG,"dx" => 5000.0); + + t0=0.0; t1=18001.0*3600.0 + u0=-(YG.-YC[1][40,21])/2000000.; u1=u0 + v0=(XG.-XC[1][40,21])/2000000.; v1=v0 + + 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) + return merge(uvt,GridVariables) +end From 0257f477ed4addb5ec9c0759834353c948330026 Mon Sep 17 00:00:00 2001 From: gaelforget Date: Mon, 10 Feb 2020 20:13:55 -0500 Subject: [PATCH 4/9] cleanup examples 1 & 2 --- src/IndividualDisplacements.jl | 5 +- src/ReadIndDisp.jl | 93 --------------------------- src/examples.jl | 112 ++++++++++++++++++++++++++++----- 3 files changed, 101 insertions(+), 109 deletions(-) diff --git a/src/IndividualDisplacements.jl b/src/IndividualDisplacements.jl index 6826ed39..0318356e 100644 --- a/src/IndividualDisplacements.jl +++ b/src/IndividualDisplacements.jl @@ -2,15 +2,16 @@ module IndividualDisplacements greet() = print("Get ready for IndividualDisplacements!") +using MeshArrays, OrdinaryDiffEq, StatsBase, DataFrames, Random + include("ReadIndDisp.jl") include("VelocityIndDisp.jl") include("examples.jl") -export VelComp, VelComp!, VelCopy, ReadGriddedFields, ReadDisplacements +export VelComp, VelComp!, VelCopy, ReadDisplacements #include("plot_pyplot.jl") #include("plot_makie.jl") - #export PlotBasic, PlotMapProj, PlotMakie end # module diff --git a/src/ReadIndDisp.jl b/src/ReadIndDisp.jl index 76901adb..ab2c99e5 100644 --- a/src/ReadIndDisp.jl +++ b/src/ReadIndDisp.jl @@ -1,68 +1,4 @@ -using StatsBase -using MeshArrays -using DataFrames - -""" - ReadGriddedFields() - -Read gridded variables from file using MeshArrays and -return result in uvetc Dictionary. -""" -function ReadGriddedFields() - -###### 1) Get gridded variables via MeshArrays.jl - - -mygrid=gcmgrid("flt_example/","ll",1,[(80,42)], [80 42], Float32, read, write) -nr=8 - -## Put grid variables in a dictionary: - -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 a dictionary: - -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)) - -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"] - -## Merge the two dictionaries: - -uvetc=Dict("u0" => u0, "u1" => u1, "v0" => v0, "v1" => v1, "t0" => t0, "t1" => t1) - -uvetc=merge(uvetc,GridVariables) - -## 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) - -end - """ ReadDisplacements(dirIn::String,prec::DataType) @@ -138,32 +74,3 @@ function ReadDisplacements(dirIn::String,prec::DataType) return df end - -""" - myread() - -Read a gridded variable from 2x2 tile files. This is used -in `ReadGriddedFields()` with `flt_example/` -""" -function myread(filRoot::String,x::MeshArray) - prec=eltype(x) - prec==Float64 ? reclen=8 : reclen=4; - - (n1,n2)=Int64.(x.grid.ioSize ./ 2); - fil=filRoot*".001.001.data" - tmp1=stat(fil); - n3=Int64(tmp1.size/n1/n2/reclen); - - v00=x.grid.write(x) - for ii=1:2; for jj=1:2; - fid = open(filRoot*".00$ii.00$jj.data") - fld = Array{prec,1}(undef,(n1*n2*n3)) - read!(fid,fld) - fld = hton.(fld) - - n3>1 ? s=(n1,n2,n3) : s=(n1,n2) - v00[1+(ii-1)*n1:ii*n1,1+(jj-1)*n2:jj*n2,:]=reshape(fld,s) - end; end; - - return x.grid.read(v00,x) -end diff --git a/src/examples.jl b/src/examples.jl index 0937c60f..d7f2cafa 100644 --- a/src/examples.jl +++ b/src/examples.jl @@ -1,36 +1,31 @@ -#using Revise, IndividualDisplacements -#include(joinpath(dirname(pathof(IndividualDisplacements)),"examples.jl")) - -using MeshArrays, OrdinaryDiffEq - """ - ex_1() + example1() Global ocean case -- just reading from file for now. ``` -df=IndividualDisplacements.ex_1() +df=IndividualDisplacements.example1() p=dirname(pathof(IndividualDisplacements)) include(joinpath(p,"plot_pyplot.jl")) PyPlot.figure(); PlotMapProj(df,300); gcf() ``` """ -function ex_1() +function example1() dirIn="run_offflt/" prec=Float32 - df=IndividualDisplacements.ReadDisplacements(dirIn,prec) + df=ReadDisplacements(dirIn,prec) end """ - ex_2() + example2() Reproducing `MITgcm/verification/flt_example/` case. This is based on an extended and modified configuration of the standard MITgcm test case. ``` -(df,ref,sol)=IndividualDisplacements.ex_2(); +(df,ref,sol)=IndividualDisplacements.example2(); p=dirname(pathof(IndividualDisplacements)) include(joinpath(p,"plot_pyplot.jl")) @@ -42,11 +37,11 @@ xaxis="x",yaxis="y",label="Julia Solution") # legend=false pl=Plots.plot!(ref[1,:],ref[2,:],lw=3,ls=:dash,lc=:red,label="MITgcm Solution") ``` """ -function ex_2() +function example2() dirIn="flt_example/" prec=Float32 - df=IndividualDisplacements.ReadDisplacements(dirIn,prec) - uvetc=IndividualDisplacements.ReadGriddedFields() + df=ReadDisplacements(dirIn,prec) + uvetc=IndividualDisplacements.example2_setup() # tmp=df[df.ID .== 200, :] nSteps=Int32(tmp[end,:time]/3600)-2 @@ -73,3 +68,92 @@ function ex_2() # return df,ref,sol end + +""" +example2_setup() + +Read gridded variables from file using MeshArrays and +return result in uvetc Dictionary. +""" +function example2_setup() + + ###### 1) Get gridded variables via MeshArrays.jl + + + mygrid=gcmgrid("flt_example/","ll",1,[(80,42)], [80 42], Float32, read, write) + nr=8 + + ## Put grid variables in a dictionary: + + 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 a dictionary: + + 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)) + + 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"] + + ## Merge the two dictionaries: + + uvetc=Dict("u0" => u0, "u1" => u1, "v0" => v0, "v1" => v1, "t0" => t0, "t1" => t1) + + uvetc=merge(uvetc,GridVariables) + + ## 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) + +end + +""" +myread() + +Read a gridded variable from 2x2 tile files. This is used +in `example2_setup()` with `flt_example/` +""" +function myread(filRoot::String,x::MeshArray) + prec=eltype(x) + prec==Float64 ? reclen=8 : reclen=4; + + (n1,n2)=Int64.(x.grid.ioSize ./ 2); + fil=filRoot*".001.001.data" + tmp1=stat(fil); + n3=Int64(tmp1.size/n1/n2/reclen); + + v00=x.grid.write(x) + for ii=1:2; for jj=1:2; + fid = open(filRoot*".00$ii.00$jj.data") + fld = Array{prec,1}(undef,(n1*n2*n3)) + read!(fid,fld) + fld = hton.(fld) + + n3>1 ? s=(n1,n2,n3) : s=(n1,n2) + v00[1+(ii-1)*n1:ii*n1,1+(jj-1)*n2:jj*n2,:]=reshape(fld,s) + end; end; + + return x.grid.read(v00,x) +end From fb2c251b8be1d3f2a151529ec76631b7d877e404 Mon Sep 17 00:00:00 2001 From: gaelforget Date: Mon, 10 Feb 2020 20:17:49 -0500 Subject: [PATCH 5/9] rename files --- src/IndividualDisplacements.jl | 4 ++-- src/{VelocityIndDisp.jl => compute.jl} | 0 src/{ReadIndDisp.jl => read.jl} | 0 3 files changed, 2 insertions(+), 2 deletions(-) rename src/{VelocityIndDisp.jl => compute.jl} (100%) rename src/{ReadIndDisp.jl => read.jl} (100%) diff --git a/src/IndividualDisplacements.jl b/src/IndividualDisplacements.jl index 0318356e..6a0b5cd0 100644 --- a/src/IndividualDisplacements.jl +++ b/src/IndividualDisplacements.jl @@ -4,8 +4,8 @@ greet() = print("Get ready for IndividualDisplacements!") using MeshArrays, OrdinaryDiffEq, StatsBase, DataFrames, Random -include("ReadIndDisp.jl") -include("VelocityIndDisp.jl") +include("compute.jl") +include("read.jl") include("examples.jl") export VelComp, VelComp!, VelCopy, ReadDisplacements diff --git a/src/VelocityIndDisp.jl b/src/compute.jl similarity index 100% rename from src/VelocityIndDisp.jl rename to src/compute.jl diff --git a/src/ReadIndDisp.jl b/src/read.jl similarity index 100% rename from src/ReadIndDisp.jl rename to src/read.jl From 2c6ee218b88ea84856328618bb4cb75562b6c138 Mon Sep 17 00:00:00 2001 From: gaelforget Date: Mon, 10 Feb 2020 20:28:06 -0500 Subject: [PATCH 6/9] further re-organize functions --- src/compute.jl | 253 ---------------------------------------- src/read.jl | 14 +++ src/update_locations.jl | 234 +++++++++++++++++++++++++++++++++++++ 3 files changed, 248 insertions(+), 253 deletions(-) create mode 100644 src/update_locations.jl diff --git a/src/compute.jl b/src/compute.jl index 1969626c..0327ef18 100644 --- a/src/compute.jl +++ b/src/compute.jl @@ -1,241 +1,4 @@ -using DataFrames - -""" - NeighborTileIndices_dpdo(ni::Int,nj::Int) - -List of W, E, S, N neighbor tile IDs in the case of a doubly -periodic domain with ni x nj tiles. -""" -function NeighborTileIndices_dpdo(ni::Int,nj::Int) - tmp=fill(0,ni*nj,4) - for i=1:ni - for j=1:nj - k=i+ni*(j-1) - kS=j-1; kS==0 ? kS=nj : nothing; kS=i+ni*(kS-1) - kN=j+1; kN==nj+1 ? kN=1 : nothing; kN=i+ni*(kN-1) - kW=i-1; kW==0 ? kW=ni : nothing; kW=kW+ni*(j-1) - kE=i+1; kE==ni+1 ? kE=1 : nothing; kE=kE+ni*(j-1) - tmp[k,1]=kW - tmp[k,2]=kE - tmp[k,3]=kS - tmp[k,4]=kN - end - end - return tmp -end - -""" - UpdateLocation_cs! - -Update location (x,y,fIndex) when out of domain. Note: initially, this -only works for the `dpdo` grid type provided by `MeshArrays.jl`. -""" -function UpdateLocation_cs!(u::Array{Float64,1},grid::Dict) - x,y = u[1:2] - fIndex = Int(u[3]) - nx,ny=grid["XC"].fSize[fIndex] - if x<0||x>nx||y<0||y>ny - j = 0 - x<0 ? j=grid["aW"][fIndex] : nothing - x>nx ? j=grid["aE"][fIndex] : nothing - y<0 ? j=grid["aS"][fIndex] : nothing - y>ny ? j=grid["aN"][fIndex] : nothing - (x,y)=grid["RelocFunctions"][j,fIndex](x,y) - u[1]=x - u[2]=y - u[3]=j - end - # - return u -end - -""" - UpdateLocation_dpdo! - -Update location (x,y,fIndex) when out of domain. Note: initially, this -only works for the `dpdo` grid type provided by `MeshArrays.jl`. -""" -function UpdateLocation_dpdo!(u::Array{Float64,1},grid::gcmgrid) - x,y = u[1:2] - fIndex = Int(u[3]) - # - nx,ny=grid.fSize[fIndex] - ni,nj=Int.(transpose(grid.ioSize)./grid.fSize[1]) - WESN=NeighborTileIndices_dpdo(ni,nj) - # - if x<0 - x=x+nx - u[1]=x - fIndex=WESN[fIndex,1] - u[3]=fIndex - elseif x>=nx - x=x-nx - u[1]=x - fIndex=WESN[fIndex,2] - u[3]=fIndex - end - # - if y<0 - y=y+ny - u[2]=y - fIndex=WESN[fIndex,3] - u[3]=fIndex - elseif y>=ny - y=y-ny - u[2]=y - fIndex=WESN[fIndex,4] - u[3]=fIndex - end - # - return u -end - - -""" - NeighborTileIndices_cs(grid::Dict) - -Derive list of neighboring tile indices for a cs or llc grid + functions that -convert indices from one tile to another. Returns a Dict to merge later. -""" -function NeighborTileIndices_cs(grid::Dict) - s = grid["XC"].fSize - nFaces = length(s) - nFaces == 5 ? s = vcat(s, s[3]) : nothing - aW=Array{Int,1}(undef,nFaces) - aE=similar(aW); aS=similar(aW); aN=similar(aW); - for i = 1:nFaces - (aW[i], aE[i], aS[i], aN[i], _, _, _, _) = MeshArrays.exch_cs_sources(i, s, 1) - end - RelocFunctions=RelocationFunctions_cs(grid["XC"]) - return Dict("aW" => aW, "aE" => aE, "aS" => aS, "aN" => aN, "RelocFunctions" => RelocFunctions) -end - -""" - RelocationFunctions_cs(xmpl) - -Define matrix of functions to convert indices across neighboring tiles -""" -function RelocationFunctions_cs(xmpl::MeshArray) - -# f1 : 0-n,0-n => -n-0,0-n for 1->2, 3->4, 5->6 -# f2 : 0-n,0-n => n-0,-n-0 for 2->4, 4->6, 6->2 -# f3 : 0-n,0-n => 0-n,-n-0 for 2->3, 4->5, 6->1 -# f4 : 0-n,0-n => -n-0,n-0 for 1->3, 3->5, 5->1 -# g1, g2, g3, g4 : the reverse connections - - f1(x, y, nx, ny) = (x .- Float64(nx), y) - f2(x, y, nx, ny) = (Float64(ny) .- y .- 1.0, x .- Float64(nx)) - f3(x, y, nx, ny) = (x, y .- Float64(ny)) - f4(x, y, nx, ny) = (y .- Float64(ny), Float64(nx) .- x .- 1.0) - - g1(x, y, nx, ny) = (x .+ Float64(nx), y) - g2(x, y, nx, ny) = (y .+ Float64(ny), Float64(nx) .- x .- 1.0) - g3(x, y, nx, ny) = (x, y .+ Float64(ny)) - g4(x, y, nx, ny) = (Float64(ny) .- y .- 1.0, x .+ Float64(nx)) - -# - - s = size.(xmpl.f) - nFaces = length(s) - tmp = Array{Function,2}(undef, 6, 6) - -# f1, f2, f3, f4 : always get nx & ny from the source tile - - tmp[2, 1] = (x, y) -> f1(x, y, s[1][1], s[1][2]) - tmp[4, 3] = (x, y) -> f1(x, y, s[3][1], s[3][2]) - tmp[6, 5] = (x, y) -> f1(x, y, s[5][1], s[5][2]) - - tmp[4, 2] = (x, y) -> f2(x, y, s[2][1], s[2][2]) - tmp[6, 4] = (x, y) -> f2(x, y, s[4][1], s[4][2]) - tmp[2, 6] = (x, y) -> f2(x, y, s[6][1], s[6][2]) - - tmp[3, 2] = (x, y) -> f3(x, y, s[2][1], s[2][2]) - tmp[5, 4] = (x, y) -> f3(x, y, s[4][1], s[4][2]) - tmp[1, 6] = (x, y) -> f3(x, y, s[6][1], s[6][2]) - - tmp[3, 1] = (x, y) -> f4(x, y, s[1][1], s[1][2]) - tmp[5, 3] = (x, y) -> f4(x, y, s[3][1], s[3][2]) - tmp[1, 5] = (x, y) -> f4(x, y, s[5][1], s[5][2]) - -# g1, g2, g3, g4 : nx or ny can come from source or target + notice nx/ny flips - - tmp[1, 2] = (x, y) -> g1(x, y, s[1][1], s[2][2]) - tmp[3, 4] = (x, y) -> g1(x, y, s[3][1], s[4][2]) - tmp[5, 6] = (x, y) -> g1(x, y, s[5][1], s[6][2]) - - tmp[2, 4] = (x, y) -> g2(x, y, s[4][1], s[2][1]) - tmp[4, 6] = (x, y) -> g2(x, y, s[6][1], s[4][1]) - tmp[6, 2] = (x, y) -> g2(x, y, s[2][1], s[6][1]) - - tmp[2, 3] = (x, y) -> g3(x, y, s[3][1], s[2][2]) - tmp[4, 5] = (x, y) -> g3(x, y, s[5][1], s[4][2]) - tmp[6, 1] = (x, y) -> g3(x, y, s[1][1], s[6][2]) - - tmp[1, 3] = (x, y) -> g4(x, y, s[1][2], s[3][2]) - tmp[3, 5] = (x, y) -> g4(x, y, s[3][2], s[5][2]) - tmp[5, 1] = (x, y) -> g4(x, y, s[5][2], s[1][2]) - - return tmp - -end - -""" - RelocationFunctions_cs_check(xmpl,RF,trgt) - -Visualize that RelocationFunctions_cs behaves as expected -""" -function RelocationFunctions_cs_check( - xmpl::MeshArray, - RF::Array{Function,2}, - trgt::Int, -) - - s = size.(xmpl.f) - nFaces = length(s) - nFaces == 5 ? s = vcat(s, s[3]) : nothing - - (aW, aE, aS, aN, iW, iE, iS, iN) = MeshArrays.exch_cs_sources(trgt, s, 1) - nx, ny = s[trgt] - p = plot([0.0, nx], [0.0, ny], color = :black) - plot!([0.0, nx], [ny, 0.0], color = :black) - for i = 1:nFaces - (nx, ny) = s[i] - x = [i - 0.5 for i = 1:nx, j = 1:ny] - y = [j - 0.5 for i = 1:nx, j = 1:ny] - c = missing - if aW == i - println("source West=$i") - c = :red - end - if aE == i - println("source East=$i") - c = :orange - end - if aS == i - println("source South=$i") - c = :blue - end - if aN == i - println("source North=$i") - c = :cyan - end - if !ismissing(c) - (x, y) = RF[trgt, i](x, y) - p = scatter!( - x, - y, - color = c, - legend = false, - marker = :rect, - markerstrokewidth = 0.0, - markersize = 1.0, - ) - end - end - return p -end - """ VelComp!(du,u,p::Dict,tim) @@ -344,19 +107,3 @@ function VelComp(du::Array{Float64,2},u::Array{Float64,2},p::Dict,tim) end return du end - -# VelCopy.jl - -""" - VelCopy(du,u,p::DataFrame,t) - -Interpolate velocity from MITgcm float_trajectories output and return -position increment `du`. -""" -function VelCopy(du,u,p::DataFrame,t) - tt=t/3600.0 - tt0=Int32(floor(tt)) - w=tt-tt0 - du[1]=(1.0-w)*p[tt0+1,:uVel]+w*p[tt0+2,:uVel] - du[2]=(1.0-w)*p[tt0+1,:vVel]+w*p[tt0+2,:vVel] -end diff --git a/src/read.jl b/src/read.jl index ab2c99e5..f17d2071 100644 --- a/src/read.jl +++ b/src/read.jl @@ -1,4 +1,18 @@ +""" + VelCopy(du,u,p::DataFrame,t) + +Interpolate velocity from MITgcm float_trajectories output and return +position increment `du`. +""" +function VelCopy(du,u,p::DataFrame,t) + tt=t/3600.0 + tt0=Int32(floor(tt)) + w=tt-tt0 + du[1]=(1.0-w)*p[tt0+1,:uVel]+w*p[tt0+2,:uVel] + du[2]=(1.0-w)*p[tt0+1,:vVel]+w*p[tt0+2,:vVel] +end + """ ReadDisplacements(dirIn::String,prec::DataType) diff --git a/src/update_locations.jl b/src/update_locations.jl new file mode 100644 index 00000000..95a83225 --- /dev/null +++ b/src/update_locations.jl @@ -0,0 +1,234 @@ +""" + NeighborTileIndices_dpdo(ni::Int,nj::Int) + +List of W, E, S, N neighbor tile IDs in the case of a doubly +periodic domain with ni x nj tiles. +""" +function NeighborTileIndices_dpdo(ni::Int,nj::Int) + tmp=fill(0,ni*nj,4) + for i=1:ni + for j=1:nj + k=i+ni*(j-1) + kS=j-1; kS==0 ? kS=nj : nothing; kS=i+ni*(kS-1) + kN=j+1; kN==nj+1 ? kN=1 : nothing; kN=i+ni*(kN-1) + kW=i-1; kW==0 ? kW=ni : nothing; kW=kW+ni*(j-1) + kE=i+1; kE==ni+1 ? kE=1 : nothing; kE=kE+ni*(j-1) + tmp[k,1]=kW + tmp[k,2]=kE + tmp[k,3]=kS + tmp[k,4]=kN + end + end + return tmp +end + +""" + update_location_cs! + +Update location (x,y,fIndex) when out of domain. Note: initially, this +only works for the `dpdo` grid type provided by `MeshArrays.jl`. +""" +function update_location_cs!(u::Array{Float64,1},grid::Dict) + x,y = u[1:2] + fIndex = Int(u[3]) + nx,ny=grid["XC"].fSize[fIndex] + if x<0||x>nx||y<0||y>ny + j = 0 + x<0 ? j=grid["aW"][fIndex] : nothing + x>nx ? j=grid["aE"][fIndex] : nothing + y<0 ? j=grid["aS"][fIndex] : nothing + y>ny ? j=grid["aN"][fIndex] : nothing + (x,y)=grid["RelocFunctions"][j,fIndex](x,y) + u[1]=x + u[2]=y + u[3]=j + end + # + return u +end + +""" + update_location_dpdo! + +Update location (x,y,fIndex) when out of domain. Note: initially, this +only works for the `dpdo` grid type provided by `MeshArrays.jl`. +""" +function update_location_dpdo!(u::Array{Float64,1},grid::gcmgrid) + x,y = u[1:2] + fIndex = Int(u[3]) + # + nx,ny=grid.fSize[fIndex] + ni,nj=Int.(transpose(grid.ioSize)./grid.fSize[1]) + WESN=NeighborTileIndices_dpdo(ni,nj) + # + if x<0 + x=x+nx + u[1]=x + fIndex=WESN[fIndex,1] + u[3]=fIndex + elseif x>=nx + x=x-nx + u[1]=x + fIndex=WESN[fIndex,2] + u[3]=fIndex + end + # + if y<0 + y=y+ny + u[2]=y + fIndex=WESN[fIndex,3] + u[3]=fIndex + elseif y>=ny + y=y-ny + u[2]=y + fIndex=WESN[fIndex,4] + u[3]=fIndex + end + # + return u +end + + +""" + NeighborTileIndices_cs(grid::Dict) + +Derive list of neighboring tile indices for a cs or llc grid + functions that +convert indices from one tile to another. Returns a Dict to merge later. +""" +function NeighborTileIndices_cs(grid::Dict) + s = grid["XC"].fSize + nFaces = length(s) + nFaces == 5 ? s = vcat(s, s[3]) : nothing + aW=Array{Int,1}(undef,nFaces) + aE=similar(aW); aS=similar(aW); aN=similar(aW); + for i = 1:nFaces + (aW[i], aE[i], aS[i], aN[i], _, _, _, _) = MeshArrays.exch_cs_sources(i, s, 1) + end + RelocFunctions=RelocationFunctions_cs(grid["XC"]) + return Dict("aW" => aW, "aE" => aE, "aS" => aS, "aN" => aN, "RelocFunctions" => RelocFunctions) +end + +""" + RelocationFunctions_cs(xmpl) + +Define matrix of functions to convert indices across neighboring tiles +""" +function RelocationFunctions_cs(xmpl::MeshArray) + +# f1 : 0-n,0-n => -n-0,0-n for 1->2, 3->4, 5->6 +# f2 : 0-n,0-n => n-0,-n-0 for 2->4, 4->6, 6->2 +# f3 : 0-n,0-n => 0-n,-n-0 for 2->3, 4->5, 6->1 +# f4 : 0-n,0-n => -n-0,n-0 for 1->3, 3->5, 5->1 +# g1, g2, g3, g4 : the reverse connections + + f1(x, y, nx, ny) = (x .- Float64(nx), y) + f2(x, y, nx, ny) = (Float64(ny) .- y .- 1.0, x .- Float64(nx)) + f3(x, y, nx, ny) = (x, y .- Float64(ny)) + f4(x, y, nx, ny) = (y .- Float64(ny), Float64(nx) .- x .- 1.0) + + g1(x, y, nx, ny) = (x .+ Float64(nx), y) + g2(x, y, nx, ny) = (y .+ Float64(ny), Float64(nx) .- x .- 1.0) + g3(x, y, nx, ny) = (x, y .+ Float64(ny)) + g4(x, y, nx, ny) = (Float64(ny) .- y .- 1.0, x .+ Float64(nx)) + +# + + s = size.(xmpl.f) + nFaces = length(s) + tmp = Array{Function,2}(undef, 6, 6) + +# f1, f2, f3, f4 : always get nx & ny from the source tile + + tmp[2, 1] = (x, y) -> f1(x, y, s[1][1], s[1][2]) + tmp[4, 3] = (x, y) -> f1(x, y, s[3][1], s[3][2]) + tmp[6, 5] = (x, y) -> f1(x, y, s[5][1], s[5][2]) + + tmp[4, 2] = (x, y) -> f2(x, y, s[2][1], s[2][2]) + tmp[6, 4] = (x, y) -> f2(x, y, s[4][1], s[4][2]) + tmp[2, 6] = (x, y) -> f2(x, y, s[6][1], s[6][2]) + + tmp[3, 2] = (x, y) -> f3(x, y, s[2][1], s[2][2]) + tmp[5, 4] = (x, y) -> f3(x, y, s[4][1], s[4][2]) + tmp[1, 6] = (x, y) -> f3(x, y, s[6][1], s[6][2]) + + tmp[3, 1] = (x, y) -> f4(x, y, s[1][1], s[1][2]) + tmp[5, 3] = (x, y) -> f4(x, y, s[3][1], s[3][2]) + tmp[1, 5] = (x, y) -> f4(x, y, s[5][1], s[5][2]) + +# g1, g2, g3, g4 : nx or ny can come from source or target + notice nx/ny flips + + tmp[1, 2] = (x, y) -> g1(x, y, s[1][1], s[2][2]) + tmp[3, 4] = (x, y) -> g1(x, y, s[3][1], s[4][2]) + tmp[5, 6] = (x, y) -> g1(x, y, s[5][1], s[6][2]) + + tmp[2, 4] = (x, y) -> g2(x, y, s[4][1], s[2][1]) + tmp[4, 6] = (x, y) -> g2(x, y, s[6][1], s[4][1]) + tmp[6, 2] = (x, y) -> g2(x, y, s[2][1], s[6][1]) + + tmp[2, 3] = (x, y) -> g3(x, y, s[3][1], s[2][2]) + tmp[4, 5] = (x, y) -> g3(x, y, s[5][1], s[4][2]) + tmp[6, 1] = (x, y) -> g3(x, y, s[1][1], s[6][2]) + + tmp[1, 3] = (x, y) -> g4(x, y, s[1][2], s[3][2]) + tmp[3, 5] = (x, y) -> g4(x, y, s[3][2], s[5][2]) + tmp[5, 1] = (x, y) -> g4(x, y, s[5][2], s[1][2]) + + return tmp + +end + +""" + RelocationFunctions_cs_check(xmpl,RF,trgt) + +Visualize that RelocationFunctions_cs behaves as expected +""" +function RelocationFunctions_cs_check( + xmpl::MeshArray, + RF::Array{Function,2}, + trgt::Int, +) + + s = size.(xmpl.f) + nFaces = length(s) + nFaces == 5 ? s = vcat(s, s[3]) : nothing + + (aW, aE, aS, aN, iW, iE, iS, iN) = MeshArrays.exch_cs_sources(trgt, s, 1) + nx, ny = s[trgt] + p = plot([0.0, nx], [0.0, ny], color = :black) + plot!([0.0, nx], [ny, 0.0], color = :black) + for i = 1:nFaces + (nx, ny) = s[i] + x = [i - 0.5 for i = 1:nx, j = 1:ny] + y = [j - 0.5 for i = 1:nx, j = 1:ny] + c = missing + if aW == i + println("source West=$i") + c = :red + end + if aE == i + println("source East=$i") + c = :orange + end + if aS == i + println("source South=$i") + c = :blue + end + if aN == i + println("source North=$i") + c = :cyan + end + if !ismissing(c) + (x, y) = RF[trgt, i](x, y) + p = scatter!( + x, + y, + color = c, + legend = false, + marker = :rect, + markerstrokewidth = 0.0, + markersize = 1.0, + ) + end + end + return p +end From f4f96c50fee1be4f449ab1b4577b1c68615e3d14 Mon Sep 17 00:00:00 2001 From: gaelforget Date: Mon, 10 Feb 2020 21:16:30 -0500 Subject: [PATCH 7/9] improve documentation --- README.md | 17 ++++------------- docs/src/index.md | 8 +++++++- 2 files changed, 11 insertions(+), 14 deletions(-) diff --git a/README.md b/README.md index 69754fdc..1374a7c0 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ [![](https://img.shields.io/badge/docs-dev-blue.svg)](https://JuliaClimate.github.io/IndividualDisplacements.jl/dev) [![DOI](https://zenodo.org/badge/208676176.svg)](https://zenodo.org/badge/latestdoi/208676176) -**IndividualDisplacements.jl** computes point displacements over a gridded domain. It is geared towards the analysis of Climate, Ocean, etc models (`Arakawa C-grids` are natively supported) and material transports within the Earth System (e.g. simulation of plastics or planktons over the Global Ocean; dusts or chemicals in the Atmosphere). +**IndividualDisplacements.jl** computes point displacements over a gridded domain. It is geared towards the analysis of Climate, Ocean, etc models (`Arakawa C-grids` are natively supported) and the simulation of material transports within the Earth System (e.g. plastics or planktons in the Ocean; dusts or chemicals in the Atmosphere). Inter-operability with popular climate model grids via [MeshArrays.jl](https://github.com/JuliaClimate/MeshArrays.jl) is an important aspect. The package can read and write individual displacement collection files, including those generated by the [MIT general circulation model](https://mitgcm.readthedocs.io/en/latest/?badge=latest). `IndividualDisplacements`'s initial test suite is based on global ocean model simulations called [ECCO (v4r2)](https://eccov4.readthedocs.io/en/latest/) and [CBIOMES (alpha)](https://cbiomes.readthedocs.io/en/latest/) (see [Forget et al. 2015](https://doi.org/10.5194/gmd-8-3071-2015)). @@ -17,19 +17,15 @@ Pkg.add("IndividualDisplacements") Pkg.test("IndividualDisplacements") ``` -![alt-text-1](examples/ex_2_mitgcm.png "From MITgcm") | ![alt-text-2](examples/ex_1_mitgcm.png "From MITgcm") -:------------------------------:|:---------------------------------: -![alt-text-3](examples/PeriodicDomainRandomFlow.png "Computed in Julia") | ![alt-text-4](examples/LatLonCap300mDepth.png "Computed in Julia") + ### Example -This example reproduces an individual trajectory originally computed by the `MITgcm` fortran implementation: +This example reproduces an individual trajectory computed by the `MITgcm` (in fortran): ``` #import software using IndividualDisplacements -p=dirname(pathof(IndividualDisplacements)); -cd(p*"/../examples") #download data if !isdir("flt_example") @@ -37,12 +33,7 @@ if !isdir("flt_example") end #run example -(df,pl)=IndividualDisplacements.ex_2() - -#display result -display(pl) +(df,ref,sol)=IndividualDisplacements.example2() ``` ![alt-text-5](examples/ex_2.png "Computed in Julia") - -![alt-text-6](examples/SolidBodyRotation.png "Unit test") diff --git a/docs/src/index.md b/docs/src/index.md index 99d9bb13..c2db9429 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,6 +1,12 @@ # IndividualDisplacements.jl -**IndividualDisplacements.jl** computes elementary point displacements over a gridded model domain. It can also read / write them from / to file. A typical application is the simulation and analysis of materials drifting or moving over the Global Ocean (e.g. plastics or planktons) or Atmosphere (e.g. dust or chemicals). Inter-operability with popular climate model grids and [MeshArrays.jl](https://github.com/JuliaClimate/MeshArrays.jl) is an important prospect. `IndividualDisplacements.jl ` was initially designed in relation to [MITgcm](https://mitgcm.readthedocs.io/en/latest/?badge=latest), [ECCOv4](https://eccov4.readthedocs.io/en/latest/) ([Forget et al. 2015](https://doi.org/10.5194/gmd-8-3071-2015)), and [CBIOMES](https://cbiomes.readthedocs.io/en/latest/) model simulations. +**IndividualDisplacements.jl** computes elementary point displacements over a gridded Earth domain (e.g. a climate model C-grid). A typical application is the simulation and analysis of materials moving through atmospheric flows (e.g. dust or chemicals) or oceanic flows (e.g. plastics or planktons). + +Inter-operability with popular climate model grids and their representation in [MeshArrays.jl](https://github.com/JuliaClimate/MeshArrays.jl) is a central element. The package can also read and plot trajectory simulation output from e.g. the [MITgcm](https://mitgcm.readthedocs.io/en/latest/?badge=latest). It was originally developed using [ECCOv4](https://eccov4.readthedocs.io/en/latest/) and [CBIOMES](https://cbiomes.readthedocs.io/en/latest/) ocean model simulations ([Forget et al. 2015](https://doi.org/10.5194/gmd-8-3071-2015)). + +The `VelComp!` and `VelComp` functions compute the velocity of tracked points. `tests/runtests.jl` uses solid body rotation as a benchmark (see below). + + ## API Guide From 8221ac3ca37190ec1e591003566277e9d7fe81eb Mon Sep 17 00:00:00 2001 From: gaelforget Date: Mon, 10 Feb 2020 21:17:02 -0500 Subject: [PATCH 8/9] bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 43ef9086..0188a980 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "IndividualDisplacements" uuid = "b92f0c32-5b7e-11e9-1d7b-238b2da8b0e6" authors = ["gaelforget "] -version = "0.1.5" +version = "0.1.6" [deps] DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" From e944e08ac660ee01fb5c88e7f96fe4b95cc55d88 Mon Sep 17 00:00:00 2001 From: gaelforget Date: Mon, 10 Feb 2020 21:29:09 -0500 Subject: [PATCH 9/9] rm un-necessary comments in code example --- README.md | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/README.md b/README.md index 1374a7c0..65fa1f5c 100644 --- a/README.md +++ b/README.md @@ -21,18 +21,15 @@ Pkg.test("IndividualDisplacements") ### Example -This example reproduces an individual trajectory computed by the `MITgcm` (in fortran): +This example reproduces an individual trajectory computed by `MITgcm` but using `Julia`: ``` -#import software using IndividualDisplacements -#download data if !isdir("flt_example") run(`git clone https://github.com/gaelforget/flt_example`) end -#run example (df,ref,sol)=IndividualDisplacements.example2() ```