Skip to content

Commit 7ea36f7

Browse files
thastthast
committed
init forward
1 parent 35ab478 commit 7ea36f7

31 files changed

Lines changed: 446288 additions & 0 deletions

‎Forward/GRAV_noisydata.obs‎

Lines changed: 963 additions & 0 deletions
Large diffs are not rendered by default.

‎Forward/Gravity_Forming_InverseMesh_matrix.py‎

Lines changed: 175 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,175 @@
1+
"""
2+
3+
"""
4+
import SimPEG.PF as PF
5+
from SimPEG import *
6+
from SimPEG.Utils import io_utils
7+
import matplotlib
8+
import time as tm
9+
import mpl_toolkits.mplot3d as a3
10+
import matplotlib.colors as colors
11+
import matplotlib.pyplot as plt
12+
import scipy as sp
13+
from scipy.interpolate import NearestNDInterpolator
14+
from sklearn.mixture import GaussianMixture
15+
import numpy as np
16+
import copy
17+
from pymatsolver import PardisoSolver
18+
19+
matplotlib.rcParams['font.size'] = 14
20+
21+
# Reproducible Science
22+
np.random.seed(518936)
23+
24+
# We first need to create a susceptibility model.
25+
# Based on a set of parametric surfaces representing TKC,
26+
# we use VTK to discretize the 3-D space.
27+
model_dir = '../Geology_Surfaces/'
28+
29+
# Create our own mesh!
30+
mult_factor = 1.
31+
csx, csy, csz = 10. / mult_factor, 10. / mult_factor, 10. / mult_factor
32+
ncx, ncy, ncz = mult_factor * 60+1, mult_factor * 60+1, mult_factor * 50
33+
npad = 10
34+
hx = [(csx, npad, -1.25), (csx, ncx), (csx, npad, 1.25)]
35+
hy = [(csy, npad, -1.25), (csy, ncy), (csy, npad, 1.25)]
36+
hz = [(csz, npad, -1.25), (csz, ncz)]
37+
mesh = Mesh.TensorMesh([hx, hy, hz], x0="CCN")
38+
xc = 300 + 5.57e5
39+
yc = 600 + 7.133e6
40+
zc = 450.
41+
x0_new = np.r_[mesh.x0[0] + xc, mesh.x0[1] + yc, mesh.x0[2] + zc]
42+
mesh.x0 = x0_new
43+
44+
# Define no-data-value
45+
ndv = -100
46+
47+
# Define survey flight height
48+
Z_bird = csz/2.
49+
print("Bird height: ", Z_bird)
50+
51+
# Read in topo surface
52+
topofile = model_dir + 'TKCtopo.dat'
53+
geosurf = [
54+
[model_dir + 'Till.ts', True, True, 0],
55+
[model_dir + 'PK1.ts', True, True, 2],
56+
[model_dir + 'PK2.ts', True, True, 3],
57+
[model_dir + 'PK3.ts', True, True, 4],
58+
[model_dir + 'HK1.ts', True, True, 5],
59+
[model_dir + 'VK.ts', True, True, 6]
60+
]
61+
62+
print('mesh nC: ', mesh.nC)
63+
# mesh.save('mesh_inverse')
64+
65+
# Create geological model
66+
modelInd = np.ones(mesh.nC) * ndv
67+
for ii in range(len(geosurf)):
68+
tin = tm.time()
69+
print("Computing indices with VTK: " + geosurf[ii][0])
70+
T, S = io_utils.read_GOCAD_ts(geosurf[ii][0])
71+
indx = io_utils.surface2inds(
72+
T, S, mesh,
73+
boundaries=geosurf[ii][1],
74+
internal=geosurf[ii][2]
75+
)
76+
print("VTK operation completed in " + str(tm.time() - tin) + " sec")
77+
modelInd[indx] = geosurf[ii][3]
78+
79+
# Fix artefact for PK1
80+
idx = np.logical_and(modelInd == 2, mesh.gridCC[:, 1] < 350 + 7.133 * 1e6)
81+
modelInd[idx] = ndv*np.ones_like(modelInd[idx])
82+
83+
# Load topography file in UBC format and find the active cells
84+
# Import Topo
85+
topofile = model_dir + 'TKCtopo.dat'
86+
topo = np.genfromtxt(topofile, skip_header=1)
87+
# Find the active cells
88+
actv = Utils.surface2ind_topo(mesh, topo, gridLoc='N')
89+
# Create active map to go from reduce set to full
90+
actvMap = Maps.InjectActiveCells(mesh, actv, ndv)
91+
print("Active cells created from topography!")
92+
93+
94+
def getModel_grav(
95+
Till=0., XVK=0., PK1=-0.8, PK2=0., PK3=0., HK1=-0.2, VK=-0.8, bkgr=0.
96+
):
97+
vals = [Till, XVK, PK1, PK2, PK3, HK1, VK]
98+
model = np.ones(mesh.nC) * bkgr
99+
100+
for ii, den in zip(range(7), vals):
101+
model[modelInd == ii] = den
102+
return model
103+
104+
model = getModel_grav()
105+
model = model[actv]
106+
107+
108+
# Here you can visualize the current model
109+
m_true = actvMap * model
110+
airc = m_true == ndv
111+
m_true[airc] = np.nan
112+
113+
# **Forward system:**
114+
# We create a synthetic survey with observations in cell center.
115+
X, Y = np.meshgrid(
116+
mesh.vectorCCx[npad:-npad:2],
117+
mesh.vectorCCy[npad:-npad:2]
118+
)
119+
120+
# Using our topography, trape the survey and shift it up by the flight height
121+
Ftopo = NearestNDInterpolator(topo[:, :2], topo[:, 2])
122+
Z = Ftopo(Utils.mkvc(X.T), Utils.mkvc(Y.T)) + Z_bird
123+
124+
rxLoc = np.c_[Utils.mkvc(X.T), Utils.mkvc(Y.T), Utils.mkvc(Z.T)]
125+
rxLoc = PF.BaseGrav.RxObs(rxLoc)
126+
print('number of data: ', rxLoc.locs.shape[0])
127+
128+
srcField = PF.BaseGrav.SrcField([rxLoc])
129+
survey = PF.BaseGrav.LinearSurvey(srcField)
130+
131+
# Now that we have a model and a survey we can build the linear system ...
132+
nactv = np.int(np.sum(actv))
133+
# Creat reduced identity map
134+
idenMap = Maps.IdentityMap(nP=nactv)
135+
# Create the forward model operator
136+
prob = PF.Gravity.GravityIntegral(
137+
mesh, rhoMap=idenMap,
138+
actInd=actv,
139+
Solver=PardisoSolver
140+
)
141+
# G_grav = np.load('./G_Grav_Inverse.npy')
142+
# prob._G = G_grav
143+
# Pair the survey and problem
144+
survey.pair(prob)
145+
146+
# Add noise to the data and assign uncertainties
147+
std = 0.0
148+
survey.eps = 0.
149+
# We add some random Gaussian noise
150+
survey.makeSyntheticData(model, std=std)
151+
survey.dobs = survey.dobs + survey.eps * np.random.randn(survey.dobs.shape)
152+
153+
io_utils.writeUBCgravityObservations(
154+
'GRAV_CoarseForward_DoubleContrast_Synthetic_data.obs', survey, survey.dobs
155+
)
156+
np.save('G_Grav_Inverse', prob.G)
157+
158+
Mesh.TensorMesh.writeModelUBC(
159+
mesh,
160+
'model_grav_inverse_mesh.den',
161+
actvMap * model
162+
)
163+
fig, ax = plt.subplots(1, 1, figsize=(6, 6))
164+
d2D = survey.dobs.reshape(X.shape)
165+
dat = ax.contourf(X - xc, Y - yc, d2D, 40)
166+
ax.set_aspect('equal')
167+
ax.plot(X.flatten() - xc, Y.flatten() - yc, 'k.', ms=2)
168+
plt.colorbar(dat)
169+
ax.set_xlabel("Easting (m)")
170+
ax.set_ylabel("Northing (m)")
171+
ax.set_title("Gz (mGal)")
172+
ax.set_xlim(-500, 500)
173+
ax.set_ylim(-500, 500)
174+
fig.savefig('Grav_Data_Inverse_Mesh.png')
175+
plt.show()

‎Forward/MAG_Forming_InverseMesh_matrix.py‎

Lines changed: 180 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,180 @@
1+
"""
2+
3+
"""
4+
5+
import SimPEG.PF as PF
6+
from SimPEG import *
7+
from SimPEG.Utils import io_utils
8+
import matplotlib
9+
import time as tm
10+
import mpl_toolkits.mplot3d as a3
11+
import matplotlib.colors as colors
12+
import matplotlib.pyplot as plt
13+
import scipy as sp
14+
from scipy.interpolate import NearestNDInterpolator
15+
from sklearn.mixture import GaussianMixture
16+
import numpy as np
17+
import copy
18+
from pymatsolver import PardisoSolver
19+
20+
matplotlib.rcParams['font.size'] = 14
21+
22+
# Reproducible Science
23+
np.random.seed(518936)
24+
25+
# We first need to create a susceptibility model.
26+
# Based on a set of parametric surfaces representing TKC,
27+
# we use VTK to discretize the 3-D space.
28+
# model_dir = '../Geological_model/'
29+
30+
# Create our own mesh!
31+
mult_factor = 1.
32+
csx, csy, csz = 10. / mult_factor, 10. / mult_factor, 10. / mult_factor
33+
ncx, ncy, ncz = mult_factor * 60+1, mult_factor * 60+1, mult_factor * 50
34+
npad = 10
35+
hx = [(csx, npad, -1.25), (csx, ncx), (csx, npad, 1.25)]
36+
hy = [(csy, npad, -1.25), (csy, ncy), (csy, npad, 1.25)]
37+
hz = [(csz, npad, -1.25), (csz, ncz)]
38+
mesh = Mesh.TensorMesh([hx, hy, hz], x0="CCN")
39+
xc = 300 + 5.57e5
40+
yc = 600 + 7.133e6
41+
zc = 450.
42+
x0_new = np.r_[mesh.x0[0] + xc, mesh.x0[1] + yc, mesh.x0[2] + zc]
43+
mesh.x0 = x0_new
44+
45+
mesh.writeUBC('mesh_inverse_ubc.msh')
46+
# Define no-data-value
47+
ndv = -100
48+
49+
# Define survey flight height
50+
Z_bird = 20.
51+
52+
# Read in topo surface
53+
model_dir = '../Geology_Surfaces/'
54+
topofile = model_dir + 'TKCtopo.dat'
55+
geosurf = [
56+
[model_dir + 'Till.ts', True, True, 0],
57+
[model_dir + 'PK1.ts', True, True, 2],
58+
[model_dir + 'PK2.ts', True, True, 3],
59+
[model_dir + 'PK3.ts', True, True, 4],
60+
[model_dir + 'HK1.ts', True, True, 5],
61+
[model_dir + 'VK.ts', True, True, 6]
62+
]
63+
64+
print('mesh nC: ', mesh.nC)
65+
66+
# Create geological model
67+
modelInd = np.ones(mesh.nC) * ndv
68+
for ii in range(len(geosurf)):
69+
tin = tm.time()
70+
print("Computing indices with VTK: " + geosurf[ii][0])
71+
T, S = io_utils.read_GOCAD_ts(geosurf[ii][0])
72+
indx = io_utils.surface2inds(
73+
T, S, mesh,
74+
boundaries=geosurf[ii][1],
75+
internal=geosurf[ii][2]
76+
)
77+
print("VTK operation completed in " + str(tm.time() - tin) + " sec")
78+
modelInd[indx] = geosurf[ii][3]
79+
80+
# Fix artefact for PK1
81+
idx = np.logical_and(modelInd == 2, mesh.gridCC[:, 1] < 350 + 7.133 * 1e6)
82+
modelInd[idx] = ndv * np.ones_like(modelInd[idx])
83+
84+
# # Load topography file in UBC format and find the active cells
85+
# Import Topo
86+
topofile = model_dir + 'TKCtopo.dat'
87+
topo = np.genfromtxt(topofile, skip_header=1)
88+
# Find the active cells
89+
actv = Utils.surface2ind_topo(mesh, topo, gridLoc='N')
90+
# Create active map to go from reduce set to full
91+
actvMap = Maps.InjectActiveCells(mesh, actv, ndv)
92+
print("Active cells created from topography!")
93+
94+
modelInd[~actv] = ndv
95+
mesh.writeModelUBC('geomodel_inverse', modelInd)
96+
print('geomodel written')
97+
98+
def getModel_mag(
99+
Till=0.0, XVK=0., PK1=5e-3, PK2=0., PK3=0., HK1=2e-2, VK=5e-3, bkgr=0.
100+
):
101+
vals = [Till, XVK, PK1, PK2, PK3, HK1, VK]
102+
model = np.ones(mesh.nC) * bkgr
103+
104+
for ii, den in zip(range(7), vals):
105+
model[modelInd == ii] = den
106+
return model
107+
108+
model = getModel_mag()
109+
model = model[actv]
110+
111+
# Here you can visualize the current model
112+
m_true = actvMap * model
113+
airc = m_true == ndv
114+
m_true[airc] = np.nan
115+
116+
# **Forward system:**
117+
# We create a synthetic survey with observations in cell center.
118+
X, Y = np.meshgrid(
119+
mesh.vectorCCx[npad:-npad:2],
120+
mesh.vectorCCy[npad:-npad:2]
121+
)
122+
123+
# Using our topography, trape the survey and shift it up by the flight height
124+
Ftopo = NearestNDInterpolator(topo[:, :2], topo[:, 2])
125+
Z = Ftopo(Utils.mkvc(X.T), Utils.mkvc(Y.T)) + Z_bird
126+
127+
rxLoc = np.c_[Utils.mkvc(X.T), Utils.mkvc(Y.T), Utils.mkvc(Z.T)]
128+
rxLoc = PF.BaseGrav.RxObs(rxLoc)
129+
print('number of data: ', rxLoc.locs.shape[0])
130+
131+
# The field parameters at TKC are [H:60,308 nT, I:83.8 d D:25.4 d ]
132+
H0 = (60308., 83.8, 25.4)
133+
srcField = PF.BaseMag.SrcField([rxLoc])
134+
srcField.param = H0
135+
survey = PF.BaseMag.LinearSurvey(srcField)
136+
137+
# Now that we have a model and a survey we can build the linear system ...
138+
nactv = np.int(np.sum(actv))
139+
# Creat reduced identity map
140+
idenMap = Maps.IdentityMap(nP=nactv)
141+
142+
# Create the forward model operator
143+
prob = PF.Magnetics.MagneticIntegral(
144+
mesh, chiMap=idenMap,
145+
actInd=actv,
146+
Solver=PardisoSolver
147+
)
148+
G_grav = np.load('./G_Mag_Inverse.npy')
149+
prob._G = G_grav
150+
# Pair the survey and problem
151+
survey.pair(prob)
152+
153+
# Add noise to the data and assign uncertainties
154+
std = 0.0
155+
survey.eps = 0.
156+
# We add some random Gaussian noise
157+
survey.makeSyntheticData(model, std=std)
158+
survey.dobs = survey.dobs + survey.eps * np.random.randn(survey.dobs.shape)
159+
160+
PF.Magnetics.writeUBCobs(
161+
'MAG_CoarseForward_Synthetic_data.obs', survey, survey.dobs
162+
)
163+
#np.save('G_Mag_Inverse', prob.G)
164+
165+
fig, ax = plt.subplots(1, 1, figsize=(6, 6))
166+
d2D = survey.dobs.reshape(X.shape)
167+
xc = 300 + 5.57e5
168+
yc = 600 + 7.133e6
169+
zc = 450.
170+
dat = ax.contourf(X - xc, Y - yc, d2D, 40)
171+
ax.set_aspect('equal')
172+
ax.plot(X.flatten() - xc, Y.flatten() - yc, 'k.', ms=2)
173+
plt.colorbar(dat)
174+
ax.set_xlabel("Easting (m)")
175+
ax.set_ylabel("Northing (m)")
176+
ax.set_title("B Amplitude (nT)")
177+
ax.set_xlim(-500, 500)
178+
ax.set_ylim(-500, 500)
179+
fig.savefig('Mag_Data_Inverse_Mesh.png')
180+
plt.show()

0 commit comments

Comments
 (0)