SPDE with variable anisotropy¶
In this tutorial, we show how the use of SPDE for Varying Anisotropy in the Simulation process
In [1]:
import gstlearn as gl
import gstlearn.plot as gp
import gstlearn.document as gdoc
import numpy as np
import matplotlib.pyplot as plt
import numpy as np
import random
import math
gdoc.setNoScroll()
Defining some global parameters
In [2]:
#Extension of the simulation grid
nxmax = 500
nymax = 200
#Well Definition
nwell = 6
vlag = 3
#Anisotropy ratio
ratio=1.5
range=150
#Some seeds
seed1 = 34556643
seed2 = 244212
seednc = 432432
seedw = 2432145
# Color Scale
zlim = [-1.6, 2.5]
nostatKeys = ["A","R2"]
# Graphic scale
gp.setDefault(dims=[10,7])
Internal function
In [3]:
def make_well(db,res,nwell=nwell,vlag=vlag,seed=seedw):
nxmax = res.getNX(0)
indexes = gl.VectorHelper.sampleRanks(ntotal=nxmax, number=nwell, seed=seed, optSort=1)
x1 = np.ones(0)
x2 = np.ones(0)
for i in indexes:
bot = math.ceil(db.getValue("W1",i))
top = math.floor(db.getValue("W2",i))
temp = np.arange(bot, top, step=vlag)
x1 = np.concatenate((x1, i * np.ones(len(temp), dtype=int) + 0.2))
x2 = np.concatenate((x2, temp))
db_sample = gl.Db.createFromSamples(len(x1),gl.ELoadBy.COLUMN,np.concatenate((x1,x2)))
db_sample.setName("New.1","x1")
db_sample.setName("New.2","x2")
db_sample.setLocator("x1",gl.ELoc.X,0)
db_sample.setLocator("x2",gl.ELoc.X,1)
err = gl.migrate(res, db_sample, "*simu*")
return db_sample
Simulating the layer boundaries
In [4]:
db = gl.DbGrid.create(nx=nxmax)
model = gl.Model.createFromParam(gl.ECov.GAUSSIAN,range=200,space=gl.SpaceRN(1))
err = gl.simtub(None,dbout=db,model=model,nbtuba=1000,seed=seed1,namconv=gl.NamingConvention("W1"))
err = gl.simtub(None,dbout=db,model=model,nbtuba=1000,seed=seed2,namconv=gl.NamingConvention("W2"))
db["W1"]=db["W1"]-min(db["W1"])
db["W2"]=db["W2"]-min(db["W2"])
db["W2"]=db["W1"]+db["W2"]+1
db["W1"]=nymax*db["W1"]/max(db["W2"])
db["W2"]=nymax*db["W2"]/max(db["W2"])
db.display()
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 1 Number of Columns = 4 Total number of samples = 500 Grid characteristics: --------------------- Origin : 0.000 Mesh : 1.000 Number : 500 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = x1 - Locator = x1 Column = 2 - Name = W1 - Locator = NA Column = 3 - Name = W2 - Locator = z1
Plotting the limits of the layer
In [5]:
ax = gp.grid1D(db,name="W1", color="blue")
ax = gp.grid1D(db,name="W2", color="green")
ax.decoration(title="Layer limits")
Creation of the varying anisotropy ("directed" by the two layers)
In [6]:
model = gl.Model.createFromParam(gl.ECov.BESSEL_K,range=range,param=1,space=gl.SpaceRN(2))
dbgrid = gl.DbGrid.create([nxmax,nymax])
ind = (dbgrid["x1"]).reshape(1,-1)[0].astype(int)
dbgrid["sel"] = (dbgrid["x2"] > db[ind,"W1"]) & (dbgrid["x2"] < db[ind,"W2"])
anglesi = np.arctan(db["W1"][1:]-db["W1"][:-1])/np.pi*180
angless = np.arctan(db["W2"][1:]-db["W2"][:-1])/np.pi*180
anglesi = np.insert(anglesi, 0, anglesi[0])
angless = np.insert(angless, 0, angless[0])
aniso = (dbgrid["x2"]-db[ind,"W1"]) / (db[ind,"W2"]-db[ind,"W1"])
aniso = anglesi[ind] + aniso * (angless[ind]-anglesi[ind])
ratio = ratio*(db[ind,"W2"]-db[ind,"W1"])/max(db["W2"]-db["W1"])
dbgrid.addColumns(aniso,"aniso")
dbgrid.addColumns(ratio,"ratio")
dbgrid.setLocator("sel",gl.ELoc.SEL)
dbgrid.setLocators(["aniso","ratio"],gl.ELoc.NOSTAT)
dbgrid.display()
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 6 Total number of samples = 100000 Number of active samples = 43700 Grid characteristics: --------------------- Origin : 0.000 0.000 Mesh : 1.000 1.000 Number : 500 200 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = x1 - Locator = x1 Column = 2 - Name = x2 - Locator = x2 Column = 3 - Name = sel - Locator = sel Column = 4 - Name = aniso - Locator = nostat1 Column = 5 - Name = ratio - Locator = nostat2
In [7]:
ax = dbgrid.plot("aniso")
ax.decoration(title="Anisotropy Angle")
In [8]:
ax = dbgrid.plot("ratio")
ax.decoration(title="Anisotropy Ratio")
Display the anisotropy maps (on a coarser grid)
Creating the Meshing
In [9]:
mesh = gl.MeshETurbo.createFromGrid(dbgrid, mode=0)
Assigning non-stationarity to the Model
In [10]:
dbcoarse = dbgrid.coarsify([20,20])
nostat = gl.NoStatArray(nostatKeys, dbcoarse)
err = model.addNoStat(nostat)
model.display()
Model characteristics ===================== Space dimension = 2 Number of variable(s) = 1 Number of basic structure(s) = 1 Number of drift function(s) = 0 Number of drift equation(s) = 0 Covariance Part --------------- K-Bessel (Third Parameter = 1) - Sill = 1.000 - Range = 150.000 - Theo. Range = 43.301 Total Sill = 1.000 Non-Stationary Parameters ------------------------- Angle : GRF=1 Str=1 V#1=1 Range : GRF=1 Str=1 V#1=2 Known Mean(s) 0.000
In [11]:
ax = gp.modelOnGrid(model=model, db=dbcoarse, scale=50, flagOrtho=False)
ax = gp.grid1D(db,name="W1", color="blue")
ax = gp.grid1D(db,name="W2", color="green")
ax.decoration(title="Layer limits")
In [12]:
nostat = gl.NoStatArray(nostatKeys, dbgrid)
err = model.addNoStat(nostat)
Non conditional simulation¶
In [13]:
spdeRes = gl.SPDE(model=model,domain=dbgrid,calcul=gl.ESPDECalcMode.SIMUNONCOND, mesh=mesh)
Process the simulation on the output grid and visualize the result (the resulting variable is stored in rank 'iuid' and name 'spde.simu'.
In [14]:
iuid = spdeRes.compute(dbgrid, namconv=gl.NamingConvention("spde.simu",False))
ax = dbgrid.plot("spde.simu")
ax.decoration(title="Non conditional Simulation")
Kriging¶
Creating a set of fictitious wells (extracted from the non-conditional simulation)
In [15]:
db_sample = make_well(db,dbgrid,nwell=nwell,vlag=vlag,seed=seedw)
In [16]:
ax = gp.grid1D(db,name="W1", color="blue")
ax = gp.grid1D(db,name="W2", color="green")
ax = gp.point(db_sample, nameCoorY="x2")
ax.decoration(title="Layer limits")
In [17]:
cgParam = gl.CGParam(5000)
spdeParam = gl.SPDEParam()
spdeParam.setCGparams(cgParam)
spdeRes = gl.SPDE(model=model,domain=dbgrid,data=db_sample,calcul=gl.ESPDECalcMode.KRIGING,
mesh=mesh,useCholesky=-1,params=spdeParam)
In [18]:
iuid = spdeRes.compute(dbgrid, namconv=gl.NamingConvention("spde",False))
dbgrid.display()
ax = dbgrid.plot("spde.estim")
ax.decoration(title="Kriging")
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 8 Total number of samples = 100000 Number of active samples = 43700 Grid characteristics: --------------------- Origin : 0.000 0.000 Mesh : 1.000 1.000 Number : 500 200 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = x1 - Locator = x1 Column = 2 - Name = x2 - Locator = x2 Column = 3 - Name = sel - Locator = sel Column = 4 - Name = aniso - Locator = nostat1 Column = 5 - Name = ratio - Locator = nostat2 Column = 6 - Name = spde.simu - Locator = NA Column = 7 - Name = spde.estim - Locator = z1
Conditional Simulation¶
In [19]:
spdeRes = gl.SPDE(model=model,domain=dbgrid,data=db_sample,calcul=gl.ESPDECalcMode.SIMUCOND,mesh=mesh)
In [20]:
iuid = spdeRes.compute(dbgrid, namconv=gl.NamingConvention("spde.condsimu",False))
dbgrid.display()
ax = dbgrid.plot("spde.condsimu")
ax.decoration(title="Conditional Simulation")
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 9 Total number of samples = 100000 Number of active samples = 43700 Grid characteristics: --------------------- Origin : 0.000 0.000 Mesh : 1.000 1.000 Number : 500 200 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = x1 - Locator = x1 Column = 2 - Name = x2 - Locator = x2 Column = 3 - Name = sel - Locator = sel Column = 4 - Name = aniso - Locator = nostat1 Column = 5 - Name = ratio - Locator = nostat2 Column = 6 - Name = spde.simu - Locator = NA Column = 7 - Name = spde.estim - Locator = NA Column = 8 - Name = spde.condsimu - Locator = z1
