# Mise a la masse experiment¶

Here, we deomstrate setting up a 3D mise-a-la-masse experiment. A point electrode is connected to a conductive target at depth and potential measurements are taken on the surface

In [1]:
import numpy as np
from SimPEG import Mesh, Maps, Utils
from SimPEG.EM.Static import DC

try:
from pymatsolver import Pardiso as Solver
except ImportError:
from SimPEG import SolverLU as Solver
print("Using LU solver. Will Be slow")

import matplotlib.pyplot as plt
plt.rcParams['image.cmap'] = 'viridis'

%matplotlib inline


## Mesh¶

In [2]:
cs = 10.
core = 600.
ncore = np.ceil(core/cs)

hx = Utils.meshTensor([(cs, npad, -1.3), (cs, ncore), (cs, npad, 1.3)])
hz = Utils.meshTensor([(cs, npad, -1.3), (cs, ncore)])

mesh = Mesh.TensorMesh([hx, hx, hz], x0=-np.r_[hx.sum()/2, hx.sum()/2, hz.sum()])

In [3]:
mesh.plotGrid()
print(mesh.nC)

643264


## Conductivity Model¶

In [4]:
sigma_background = 1e-5
sigma_target = 1e-1

sigma = np.ones(mesh.nC)*sigma_background

In [5]:
target_thickness = 30.
target_width = 100.
target_length = 200.
target_center = -450.
target_dip = -80.

In [6]:
# find the coordinates of the top and bottom of the place
dz = 0.5 * target_length * np.sin(np.deg2rad(target_dip))
dx = 0.5 * target_length * np.cos(np.deg2rad(target_dip))

In [7]:
target_line = lambda x: np.tan(np.deg2rad(target_dip)) * x + target_center

In [8]:
target_right = target_line(mesh.gridCC[:,0] - target_thickness/2.)
target_left = target_line(mesh.gridCC[:,0] + target_thickness/2.)
target_top = target_center - dz
target_bottom = target_center + dz

target_inds = (
(mesh.gridCC[:, 2] <= target_right) & (mesh.gridCC[:, 2] >= target_left) &
(mesh.gridCC[:, 1] >= -target_width/2.) & (mesh.gridCC[:, 1] <= target_width/2.) &
(mesh.gridCC[:, 2] <= target_top) & (mesh.gridCC[:, 2] >= target_bottom)
)

In [9]:
sigma[target_inds] = sigma_target

In [10]:
import ipywidgets

In [11]:
def plot_model(i=mesh.vnC[1]/2):
i = int(i)
fig, ax = plt.subplots(1, 1)
plt.colorbar(mesh.plotSlice(np.log10(sigma), normal='Y', ax=ax, ind=i)[0], ax=ax)
ax.set_xlim([-250, 250])
ax.set_ylim([-800, 0])

In [ ]:
ipywidgets.interact(
plot_model,
i = ipywidgets.IntSlider(min=0, max=mesh.vnC[1], value= int(mesh.vnC[1]/2.))
)

Out[ ]:
<function __main__.plot_model(i=46.0)>
In [ ]:
src_a = np.r_[0., 0., target_center]
src_b = np.r_[0., 1500., -cs/2. ]

In [ ]:
fig, ax = plt.subplots(1,1)
plt.colorbar(mesh.plotSlice(np.log10(sigma), normal='Y', ax=ax)[0], ax=ax)
ax.plot(src_a[0], src_a[2], 'ro')
ax.set_xlim([-250, 250])
ax.set_ylim([-800, 0])

Out[ ]:
(-800, 0)

## Survey parameters¶

Specify the receiver locations

In [ ]:
rx_x = mesh.vectorCCx[(mesh.vectorCCx >= -300) & (mesh.vectorCCx <= 300)][::2]
rx_y = mesh.vectorCCy[(mesh.vectorCCy >= -300) & (mesh.vectorCCy <= 300)][::2]
rx_z = np.r_[mesh.vectorCCz.max()]

rx_locs = Utils.ndgrid([rx_x, rx_y, rx_z])
plt.plot(rx_locs[:,0], rx_locs[:,1], '.')

Out[ ]:
[<matplotlib.lines.Line2D at 0x124d7cd30>]
In [ ]:
rx = DC.Rx.Pole(rx_locs)
src = DC.Src.Dipole([rx], src_a, src_b)
survey = DC.Survey([src])

In [ ]:
prob = DC.Problem3D_CC(mesh, sigmaMap=Maps.IdentityMap(mesh), Solver=Solver)
prob.pair(survey)


## Compute Data¶

In [ ]:
%%time
dpred = survey.dpred(sigma)

In [ ]:
plt.colorbar(
plt.contourf(
rx_x, rx_y,
dpred.reshape(len(rx_x), len(rx_y), order='F'), 30,
),
)

In [ ]: