forked from simpeg/simpeg
-
Notifications
You must be signed in to change notification settings - Fork 1
/
plot_inductive_src_permeable_target.py
269 lines (207 loc) · 8.03 KB
/
plot_inductive_src_permeable_target.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
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
"""
EM: TDEM: Permeable Target, Inductive Source
============================================
In this example, we demonstrate 2 approaches for simulating TDEM data when
a permeable target is present in the simulation domain. In the first, we
use a step-on waveform (QuarterSineRampOnWaveform) and look at the magnetic
flux at a late on-time. In the second, we solve the magnetostatic problem
to compute the initial magnetic flux so that a step-off waveform may be used.
A cylindrically symmetric mesh is employed and a circular loop source is used
"""
import discretize
from discretize import utils
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from scipy.constants import mu_0
from pymatsolver import Pardiso
import time
from SimPEG.EM import TDEM
from SimPEG import Utils, Maps, versions
###############################################################################
# Model Parameters
# ----------------
#
# Here, we define our simulation parameters. The target has a relative
# permeability of 100 :math:`\mu_0`
target_mur = 100 # permeability of the target
target_l = 500 # length of target
target_r = 50 # radius of the target
sigma_back = 1e-5 # conductivity of the background
radius_loop = 100 # radius of the transmitter loop
###############################################################################
# Mesh
# ----
#
# Next, we create a cylindrically symmteric tensor mesh
csx = 5. # core cell size in the x-direction
csz = 5. # core cell size in the z-direction
domainx = 100 # use a uniform cell size out to a radius of 100m
# padding parameters
npadx, npadz = 15, 15 # number of padding cells
pfx = 1.4 # expansion factor for the padding to infinity in the x-direction
pfz = 1.4 # expansion factor for the padding to infinity in the z-direction
ncz = int(target_l/csz) # number of z cells in the core region
# create the cyl mesh
mesh = discretize.CylMesh([
[(csx, int(domainx/csx)), (csx, npadx, pfx)],
1,
[(csz, npadz, -pfz), (csz, ncz), (csz, npadz, pfz)]
])
# put the origin at the top of the target
mesh.x0 = [0, 0, -mesh.hz[:npadz + ncz].sum()]
# plot the mesh
mesh.plotGrid()
###############################################################################
# Assign physical properties on the mesh
mur_model = np.ones(mesh.nC)
# find the indices of the target
x_inds = mesh.gridCC[:, 0] < target_r
z_inds = (mesh.gridCC[:, 2] <= 0) & (mesh.gridCC[:, 2] >= -target_l)
mur_model[x_inds & z_inds] = target_mur
mu_model = mu_0 * mur_model
sigma = np.ones(mesh.nC) * sigma_back
###############################################################################
# Plot the models
xlim = np.r_[-200, 200] # x-limits in meters
zlim = np.r_[-1.5*target_l, 10.] # z-limits in meters. (z-positive up)
fig, ax = plt.subplots(1, 1, figsize=(6, 5))
# plot the permeability
plt.colorbar(mesh.plotImage(
mur_model, ax=ax,
pcolorOpts={'norm': LogNorm()}, # plot on a log-scale
mirror=True
)[0], ax=ax)
ax.plot(np.r_[radius_loop], np.r_[0.], 'wo', markersize=8)
ax.plot(np.r_[-radius_loop], np.r_[0.], 'wx', markersize=8)
ax.set_title("Relative permeability", fontsize=13)
ax.set_xlim(xlim)
ax.set_ylim(zlim)
###############################################################################
# Waveform for the Long On-Time Simulation
# ----------------------------------------
#
# Here, we define our time-steps for the simulation where we will use a
# waveform with a long on-time to reach a steady-state magnetic field and
# define a quarter-sine ramp-on waveform as our transmitter waveform
ramp = [
(1e-5, 20), (1e-4, 20), (3e-4, 20), (1e-3, 20), (3e-3, 20), (1e-2, 20),
(3e-2, 20), (1e-1, 20), (3e-1, 20), (1, 50)
]
time_mesh = discretize.TensorMesh([ramp])
# define an off time past when we will simulate to keep the transmitter on
offTime = 100
quarter_sine = TDEM.Src.QuarterSineRampOnWaveform(
ramp_on=np.r_[0., 3], ramp_off= offTime - np.r_[1., 0]
)
# evaluate the waveform at each time in the simulation
quarter_sine_plt = [quarter_sine.eval(t) for t in time_mesh.gridN]
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.plot(time_mesh.gridN, quarter_sine_plt)
ax.plot(time_mesh.gridN, np.zeros(time_mesh.nN), 'k|', markersize=2)
ax.set_title('quarter sine waveform')
###############################################################################
# Sources for the 2 simulations
# -----------------------------
#
# We use two sources, one for the magnetostatic simulation and one for the
# ramp on simulation.
# For the magnetostatic simulation. The default waveform is a step-off
src_magnetostatic = TDEM.Src.CircularLoop(
[], loc=np.r_[0., 0., 0.], orientation="z", radius=100,
)
# For the long on-time simulation. We use the ramp-on waveform
src_ramp_on = TDEM.Src.CircularLoop(
[], loc=np.r_[0., 0., 0.], orientation="z", radius=100,
waveform=quarter_sine
)
src_list_magnetostatic = [src_magnetostatic]
src_list_ramp_on = [src_ramp_on]
###############################################################################
# Create the simulations
# ----------------------
#
# To simulate magnetic flux data, we use the b-formulation of Maxwell's
# equations
prob_magnetostatic = TDEM.Problem3D_b(
mesh=mesh, sigmaMap=Maps.IdentityMap(mesh), timeSteps=ramp,
Solver=Pardiso
)
prob_ramp_on = TDEM.Problem3D_b(
mesh=mesh, sigmaMap=Maps.IdentityMap(mesh), timeSteps=ramp,
Solver=Pardiso
)
survey_magnetostatic = TDEM.Survey(srcList=src_list_magnetostatic)
survey_ramp_on = TDEM.Survey(src_list_ramp_on)
prob_magnetostatic.pair(survey_magnetostatic)
prob_ramp_on.pair(survey_ramp_on)
###############################################################################
# Run the long on-time simulation
# -------------------------------
t = time.time()
print('--- Running Long On-Time Simulation ---')
prob_ramp_on.mu = mu_model
fields = prob_ramp_on.fields(sigma)
print(" ... done. Elapsed time {}".format(time.time() - t))
print('\n')
# grab the last time-step in the simulation
b_ramp_on = utils.mkvc(fields[:, 'b', -1])
###############################################################################
# Compute Magnetostatic Fields from the step-off source
# -----------------------------------------------------
prob_magnetostatic.mu = mu_model
prob_magnetostatic.model = sigma
b_magnetostatic = src_magnetostatic.bInitial(prob_magnetostatic)
###############################################################################
# Plot the results
# -----------------------------------------------------
def plotBFieldResults(
ax=None, clim_min=None, clim_max=None,
max_depth=1.5*target_l, max_r=100,
top=10., view="magnetostatic"
):
if ax is None:
plt.subplots(1, 1, figsize=(6, 7))
assert view.lower() in ["magnetostatic", "late_ontime", 'diff']
xlim = max_r*np.r_[-1, 1] # x-limits in meters
zlim = np.r_[-max_depth, top] # z-limits in meters. (z-positive up)
clim = None
if clim_max is not None and clim_max != 0.:
clim = clim_max * np.r_[-1, 1]
if clim_min is not None and clim_min != 0.:
clim[0] = clim_min
if view == "magnetostatic":
plotme = b_magnetostatic
elif view == "late_ontime":
plotme = b_ramp_on
elif view == "diff":
plotme = b_magnetostatic-b_ramp_on
cb = plt.colorbar(mesh.plotImage(
plotme,
view='vec', vType='F',
ax=ax, range_x=xlim, range_y=zlim,
sample_grid=np.r_[np.diff(xlim)/100., np.diff(zlim)/100.],
mirror=True,
pcolorOpts={'norm': LogNorm()}
)[0], ax=ax)
cb.set_clim(clim)
ax.set_title('{}'.format(view), fontsize=13)
ax.set_xlim(xlim)
ax.set_ylim(zlim)
cb.update_ticks()
return ax
fig, ax = plt.subplots(1, 3, figsize=(12, 5))
for a, v in zip(ax, ["magnetostatic", "late_ontime", "diff"]):
a = plotBFieldResults(
ax=a,
clim_min=1e-15,
clim_max=1e-7,
view=v,
max_r=200
)
plt.tight_layout()
###############################################################################
# Print the version of SimPEG and dependencies
# --------------------------------------------
#
versions()