forked from simpeg/simpeg
-
Notifications
You must be signed in to change notification settings - Fork 1
/
test_grav_inversion_linear.py
154 lines (117 loc) · 5.06 KB
/
test_grav_inversion_linear.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
from __future__ import print_function
import unittest
import numpy as np
from SimPEG import Mesh
from SimPEG import Utils
from SimPEG import Maps
from SimPEG import Regularization
from SimPEG import DataMisfit
from SimPEG import Optimization
from SimPEG import InvProblem
from SimPEG import Directives
from SimPEG import Inversion
import matplotlib.pyplot as plt
import SimPEG.PF as PF
np.random.seed(43)
class GravInvLinProblemTest(unittest.TestCase):
def setUp(self):
ndv = -100
# Create a self.mesh
dx = 5.
hxind = [(dx, 5, -1.3), (dx, 5), (dx, 5, 1.3)]
hyind = [(dx, 5, -1.3), (dx, 5), (dx, 5, 1.3)]
hzind = [(dx, 5, -1.3), (dx, 6)]
self.mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCC')
# Get index of the center
midx = int(self.mesh.nCx/2)
midy = int(self.mesh.nCy/2)
# Lets create a simple Gaussian topo and set the active cells
[xx, yy] = np.meshgrid(self.mesh.vectorNx, self.mesh.vectorNy)
zz = -np.exp((xx**2 + yy**2) / 75**2) + self.mesh.vectorNz[-1]
# Go from topo to actv cells
topo = np.c_[Utils.mkvc(xx), Utils.mkvc(yy), Utils.mkvc(zz)]
actv = Utils.surface2ind_topo(self.mesh, topo, 'N')
actv = np.asarray([inds for inds, elem in enumerate(actv, 1)
if elem], dtype=int) - 1
# Create active map to go from reduce space to full
self.actvMap = Maps.InjectActiveCells(self.mesh, actv, -100)
nC = len(actv)
# Create and array of observation points
xr = np.linspace(-20., 20., 20)
yr = np.linspace(-20., 20., 20)
X, Y = np.meshgrid(xr, yr)
# Move the observation points 5m above the topo
Z = -np.exp((X**2 + Y**2) / 75**2) + self.mesh.vectorNz[-1] + 5.
# Create a MAGsurvey
locXYZ = np.c_[Utils.mkvc(X.T), Utils.mkvc(Y.T), Utils.mkvc(Z.T)]
rxLoc = PF.BaseGrav.RxObs(locXYZ)
srcField = PF.BaseGrav.SrcField([rxLoc])
survey = PF.BaseGrav.LinearSurvey(srcField)
# We can now create a density model and generate data
# Here a simple block in half-space
model = np.zeros((self.mesh.nCx, self.mesh.nCy, self.mesh.nCz))
model[(midx-2):(midx+2), (midy-2):(midy+2), -6:-2] = 0.5
model = Utils.mkvc(model)
self.model = model[actv]
# Create active map to go from reduce set to full
actvMap = Maps.InjectActiveCells(self.mesh, actv, ndv)
# Create reduced identity map
idenMap = Maps.IdentityMap(nP=nC)
# Create the forward model operator
prob = PF.Gravity.GravityIntegral(
self.mesh,
rhoMap=idenMap,
actInd=actv
)
# Pair the survey and problem
survey.pair(prob)
# Compute linear forward operator and compute some data
d = prob.fields(self.model)
# Add noise and uncertainties (1nT)
data = d + np.random.randn(len(d))*0.001
wd = np.ones(len(data))*.001
survey.dobs = data
survey.std = wd
# PF.Gravity.plot_obs_2D(survey.srcField.rxList[0].locs, d=data)
# Create sensitivity weights from our linear forward operator
wr = PF.Magnetics.get_dist_wgt(self.mesh, locXYZ, actv, 2., 2.)
wr = wr**2.
# Create a regularization
reg = Regularization.Sparse(self.mesh, indActive=actv, mapping=idenMap)
reg.cell_weights = wr
reg.norms = np.c_[0, 0, 0, 0]
reg.gradientType = 'component'
# reg.eps_p, reg.eps_q = 5e-2, 1e-2
# Data misfit function
dmis = DataMisfit.l2_DataMisfit(survey)
dmis.W = 1/wd
# Add directives to the inversion
opt = Optimization.ProjectedGNCG(maxIter=100, lower=-1., upper=1.,
maxIterLS=20, maxIterCG=10,
tolCG=1e-3)
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e+8)
# Here is where the norms are applied
IRLS = Directives.Update_IRLS(f_min_change=1e-4,
minGNiter=1)
update_Jacobi = Directives.UpdatePreconditioner()
self.inv = Inversion.BaseInversion(invProb,
directiveList=[IRLS,
update_Jacobi])
def test_grav_inverse(self):
# Run the inversion
mrec = self.inv.run(self.model)
residual = np.linalg.norm(mrec-self.model) / np.linalg.norm(self.model)
print(residual)
# plt.figure()
# ax = plt.subplot(1, 2, 1)
# midx = int(self.mesh.nCx/2)
# self.mesh.plotSlice(self.actvMap*mrec, ax=ax, normal='Y', ind=midx,
# grid=True, clim=(0, 0.5))
# ax = plt.subplot(1, 2, 2)
# midx = int(self.mesh.nCx/2)
# self.mesh.plotSlice(self.actvMap*self.model, ax=ax, normal='Y', ind=midx,
# grid=True, clim=(0, 0.5))
# plt.show()
self.assertTrue(residual < 0.05)
if __name__ == '__main__':
unittest.main()