SPDE¶
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 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]
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[i, "W1"][0])
top = math.floor(db[i, "W2"][0])
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.create(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]:
gp.grid1D(db, name="W1", color="blue")
gp.grid1D(db, name="W2", color="green")
gp.decoration(title="Layer limits")
Creation of the varying anisotropy ("directed" by the two layers)
In [6]:
model = gl.Model.createFromParam(
gl.ECov.MATERN, range=range, param=1, space=gl.SpaceRN.create(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.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 = NA Column = 5 - Name = ratio - Locator = NA
In [7]:
gp.plot(dbgrid, "aniso")
gp.decoration(title="Anisotropy Angle")
In [8]:
gp.plot(dbgrid, "ratio")
gp.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([10, 10])
err = model.getCovAniso(0).makeAngleNoStatDb("aniso", 0, dbgrid)
err = model.getCovAniso(0).makeRangeNoStatDb("ratio", 1)
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 --------------- Matern (Third Parameter = 1) - Sill = 1.000 - Range = 150.000 - Theo. Range = 43.301 List of Non-Stationary Parameters 1 - Angle : IdAngle=1 2 - Range : IDir=2 Total Sill = 1.000 Known Mean(s) 0.000
In [11]:
gp.covaOnGrid(model.getCovAniso(0), db=dbcoarse, scale=50, flagOrtho=False)
gp.grid1D(db, name="W1", color="blue")
gp.grid1D(db, name="W2", color="green")
gp.decoration(title="Layer limits")
Non conditional simulation¶
Process the simulation on the output grid and visualize the result (the resulting variable is stored in rank 'iuid' and name 'spde.simu'.
In [12]:
meshes = gl.VectorMeshes([mesh])
In [13]:
err = gl.simulateSPDE(
None, dbgrid, model, meshesS=meshes, namconv=gl.NamingConvention("spde.simu", False)
)
In [14]:
gp.plot(dbgrid, "spde.simu", flagLegend=True)
gp.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]:
gp.grid1D(db, name="W1", color="blue")
gp.grid1D(db, name="W2", color="green")
gp.symbol(db_sample, nameCoorY="x2")
gp.decoration(title="Layer limits")
In [17]:
err = gl.krigingSPDE(
db_sample,
dbgrid,
model,
meshesK=meshes,
namconv=gl.NamingConvention("spde.kriging", False, False),
)
In [18]:
gp.grid1D(db, name="W1", color="blue")
gp.grid1D(db, name="W2", color="green")
gp.plot(dbgrid, "spde.kriging")
gp.symbol(db_sample, nameCoorY="x2")
gp.decoration(title="Kriging")
Conditional Simulation¶
In [19]:
err = gl.simulateSPDE(
db_sample,
dbgrid,
model,
meshesK=meshes,
namconv=gl.NamingConvention("spde.condSimu", False, False),
)
In [20]:
gp.grid1D(db, name="W1", color="blue")
gp.grid1D(db, name="W2", color="green")
gp.plot(dbgrid, "spde.condSimu")
gp.symbol(db_sample, nameCoorY="x2")
gp.decoration(title="Conditional Simulation")
