gher-uliege / DIVAnd.jl

DIVAnd performs an n-dimensional variational analysis of arbitrarily located observations
GNU General Public License v2.0
70 stars 11 forks source link

[DOC] documentation of `diffusion` function #120

Closed ctroupin closed 1 year ago

ctroupin commented 1 year ago

Request from users: could we describe the use of the function diffusion (which they use to create a field of variable correlation length)

https://github.com/gher-uliege/DIVAnd.jl/blob/master/src/varanalysis.jl#L23 AeEBl131PEQFVI06 in particular the parameters:

Alexander-Barth commented 1 year ago

Is this clear with this information: https://github.com/gher-uliege/DIVAnd.jl/blob/9cc4f2dcd1da3092602edbb4509d1e457c613921/src/varanalysis.jl#L22

jmbeckers commented 1 year ago

I suppose when no boundary condition is given, zero gradient is applied naturally?

qcazenave commented 1 year ago

Thanks for the added comments ! I wonder though how come it works with pmn defined in degrees-1 and slen in meters as it seems to be the case in the code "combined_correlation_len.jl" ?

ctroupin commented 1 year ago

https://github.com/gher-uliege/EMODnet-Chemistry/blob/master/julia/european-climatologies/combined_correlation_len.jl#L67

slen = (100e3,100e3) made me think it's expressed in meters but that the invert of pmn to that depends.

qcazenave commented 1 year ago

If pmn is defined with: mask,pmn = DIVAnd.domain(bathname,bathisglobal,lonr,latr) and lonr and latr are given in degrees, doesn't that mean that pmn is in degrees-1 ? And yet, after that, we have slen = (800e3,800e3) lenf = DIVAnd.diffusion(mask,pmn,slen,len) the result is good but I wonder whether it is by chance or there is something I do not understand ?

jmbeckers commented 1 year ago

I think indeed that the call is not appropriate, MWE:


using DIVAnd
using PyPlot

testsize=100

xi,yi= ndgrid(range(9.5,stop=30.9,length=testsize),range(53.,stop=60.,length=testsize));

# reference field

# all points are valid points
mask = trues(size(xi));

# this problem has a simple cartesian metric
# pm is the inverse of the resolution along the 1st dimension
# pn is the inverse of the resolution along the 2nd dimension

pm = ones(size(xi)) / (xi[2,1]-xi[1,1]);
pn = ones(size(xi)) / (yi[1,2]-yi[1,1]);

len=zeros(Float64,size(xi))

len[50,50]=1.0

@show size(xi)

slen = (800e3,800e3) 

 # decrease to avoid too long calculations ...
slen=(100.,100.)
lenf = DIVAnd.diffusion(mask,(pm,pn),slen,len)

pcolor(xi,yi,len),colorbar()
figure()
pcolor(xi,yi,lenf),colorbar()

extrema(lenf)

slen=(1.,1.)
lenf = DIVAnd.diffusion(mask,(pm,pn),slen,len)

figure()
pcolor(xi,yi,lenf),colorbar()

extrema(lenf)

It will just smooth out everything until a constant value is reached (I wonder how long the execution actually takes ?)

Now why does the final result look fine ? It think it is because

https://github.com/gher-uliege/EMODnet-Chemistry/blob/master/julia/european-climatologies/combined_correlation_len.jl#L89

just decreased the value to a fifth of the constant value when approaching a coast ?

Alexander-Barth commented 1 year ago

If pmn is defined with: mask,pmn = DIVAnd.domain(bathname,bathisglobal,lonr,latr) and lonr and latr are given in degrees, doesn't that mean that pmn is in degrees-1 ?

Actually, DIVAnd.domain will always return pm and pn in m⁻¹. You can verify it by looking at the order of magnitude (typicall 1e-5 m⁻¹ for a resolution of 100 km).

This has been added to the documentation to avoid a misunderstand in future: https://github.com/gher-uliege/DIVAnd.jl/commit/232c039e5eb8241529df5f5fe86096509d956328

jmbeckers commented 1 year ago

That was unexpected for me, as almost all functions in DIVAnd are unaware of units and dimensions so I assumed it was just doing what I wrote by hand in the MWE). Maybe we should list the functions that are specifically related to analyses on Earth and/or check all function descriptions for units ?

qcazenave commented 1 year ago

Ah sorry, I was so sure I had already checked the order of magnitude of pmn that I didn't look at it before sending my message. Thanks for the documentation.

Alexander-Barth commented 1 year ago

@jmbeckers I think the natural place to document this would be in the function doc string. Also the main function diva3d makes this assumption (as correlation length are given in meters). These are the function who implicetly depend on deg2m (convertion of degree to meters) that I have seen:

DIVAnd_metric
domain
diva3d
jmbeckers commented 1 year ago

Yes, documenting them clearly should be enough.