Diffusion Advection¶
In [1]:
import gstlearn as gl
import numpy as np
import matplotlib.pyplot as plt
In [2]:
cov = gl.CovDiffusionAdvection()
vel = [0.,0]
scaleTime = 0.2
sigma2 = 3
p1 = 1
p2 = 1
scale = 2
ctxt = gl.CovContext()
cova1 = gl.CovAniso.createIsotropic(type=gl.ECov.MATERN,
range = scale,
param =p1, flagRange= False,
ctxt=ctxt)
cova2 = gl.CovAniso.createIsotropic(type=gl.ECov.MATERN,
range = scale,
param =p2,
flagRange= False,ctxt=ctxt)
cova1 = None
cov = gl.CovDiffusionAdvection.create(cova1,cova2,scaleTime,vel,sigma2)
covprod = cov.getSpatialTrace()
covprod.display()
Markov (Third Parameter = 1) - Sill = 3.000 - Range = 12.000 - Theo. Range = 2.000
In [3]:
cov.getSpatialTrace().getMarkovCoeffs()
Out[3]:
array([1., 2., 1.])
In [4]:
cov.getGlobalCorrec()
Out[4]:
0.24989756768264187
Spatial Trace¶
In [5]:
mH1 = 50
mH2 = 50
N=2**8
result = np.array(cov.evalCovFFT([mH1,mH2],0,N).getValues()).reshape(N,N)
trace= np.array(covprod.evalCovFFT([mH1,mH2],N).getValues()).reshape(N,N)
In [6]:
plt.imshow(result,extent = [-mH1,mH1,-mH2,mH2])
plt.colorbar()
plt.show()
plt.imshow(trace,extent = [-mH1,mH1,-mH2,mH2])
c = plt.colorbar()
In [7]:
X1 = np.linspace(-mH1,mH1,N)
X2 = np.linspace(-mH2,mH2,N)
plt.plot(X1,trace[int(N/2),:],label = "Trace")
plt.plot(X1,result[int(N/2),:],"--",label="FFT")
l = plt.legend()
In [8]:
N=2**7
ntimes = 40
restime = np.zeros(shape=ntimes)
result = np.zeros(shape = (N*N,ntimes))
times = np.arange(ntimes) * 0.5
for it,time in enumerate(times):
result[:,it] = np.array(cov.evalCovFFT([mH1,mH2],time,N).getValues())
restime[it] = result[:,it].reshape(N,N)[int(N/2),int(N/2)]
In [9]:
fig, axs = plt.subplots(nrows=4,ncols=5, figsize=(10,5))
k = 0
for i in range(4) :
for j in range(5):
pt = axs[i,j].imshow(result[:,k].reshape(N,N),origin = "lower",extent = [-mH1,mH1,-mH2,mH2])
pt.set_clim([result.min(),result.max()])
k+=1
Temporal covariance¶
In [10]:
plt.plot(times,restime,label = "FFT")
if cova1 is None and np.sum(vel) == 0:
plt.plot(times,sigma2 * np.exp(-times*scaleTime),"--",label = "Theoretical")
#plt.ylim((0,restime.max()))
ax = plt.legend()
