# Crosswell EM inversion in 2d cylindrical coordinate

Hi All,

Using SimPEG, my colleagues and I perform crosswell EM modeling in 2d cylindrical coordinate system. We were able to generate crosswell EM forward modeling data and confirmed good agreement between SimPEG and reference solutions.

Unfortunately, we were not successful in inverting the same synthetic data. We assume that we made a mistake in our script and are looking for an EM inversion example in 2d cylindrical coordinate system. If anyone can share an example with us, please let us know. Thanks!

Best,
Evan

1 Like

Hi Evan,

Glad to hear that you are using SimPEG!

If you can share the script that you are using, I could take a look at that as well.

All the best,

Seogi.

Hi Seogi,

Thanks for your kind help. Yes, we use SimPEG in our various projects. Thanks to all SimPEG developers!

After I ran a crosswell EM forward model in the 2d cylindrical coordinate system, I tried to invert it but got stuck to an error. Below is my script. Could you take a look? I assume that I made a mistake but I couldn’t make it run.

Best regards,
Evan

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from scipy.constants import mu_0, inch, foot
import os
import ipywidgets
import time
from string import ascii_lowercase
from matplotlib import rcParams
from matplotlib import gridspec

import discretize
from discretize import utils
from SimPEG.electromagnetics import frequency_domain as fdem
from SimPEG import maps, utils
from pymatsolver import Pardiso

from SimPEG import (
optimization,
data_misfit,
regularization,
inverse_problem,
inversion,
directives,
Report,
)

import casingSimulations as casing_sim

get_ipython().run_line_magic(‘matplotlib’, ‘inline’)

from matplotlib import rcParams
rcParams[“font.size”] = 14

freq=500.0
sigma_bgd1=1.0/3.0
sigma_bgd2=1.0/4.0
sigma_bgd3=1.0/6.0
sigma_reservoir=1.0/40.0
mur_slug=np.r_[1.25]
sigma_slug=1.0/7.6
slug_h=20
slug_w=np.r_[30.0, 60.0, 85.0]

skin_depth=503*np.sqrt(1/sigma_reservoir/freq)
print (skin_depth)

ncr = 100
ncz = 100
dr = 10
dz = 10
exp_r = 1.15
exp_z = 1.15

hr = [(dr, ncr), (dr, npad_r, exp_r)]

mesh = discretize.CylMesh([hr, 1, hz], x0=“00C”)

x0 = mesh.x0

nC = mesh.nC

cc = mesh.gridCC

fig, ax = plt.subplots(1,1)
mesh.plotGrid(ax=ax)

print(
f"The mesh has {mesh.nC} cells \n"
f" * x-extent: {mesh.vectorNx[-1]:1.1e} m\n"
f" * z-extent: [{mesh.vectorNz[0]:1.1e} m, {mesh.vectorNz[-1]:1.1e} m]"
)

fig, ax = plt.subplots(1,1)
mesh.plotGrid(ax=ax)
ax.set_xlim([0, 500])
ax.set_ylim([-300, 300])

print(“the mesh’s lower left corner {}”.format(mesh.x0))
print(“the mesh has {} cells”.format(mesh.nC))

def change_local_property (mesh, array_model, local_value, z1, z2, r1, r2):
inds = (
(mesh.gridCC[:, 2] >= z1) & (mesh.gridCC[:, 2] <= z2) & (mesh.gridCC[:, 0] >= r1) & (mesh.gridCC[:, 0] <= r2)
)
array_model[inds] = local_value
return array_model

m_sigma_bgd = sigma_bgd1 * np.ones(mesh.nC) # background model
m_sigma_x1 = sigma_bgd1 * np.ones(mesh.nC) # x=30
m_sigma_x2 = sigma_bgd1 * np.ones(mesh.nC) # x=60
m_sigma_x3 = sigma_bgd1 * np.ones(mesh.nC) # x=85

m_sigma_bgd = change_local_property (mesh, m_sigma_bgd, 1/13/3, -10, 10, 0, 1000)
m_sigma_x1 = change_local_property (mesh, m_sigma_x1, sigma_reservoir, -10, 10, 0, 1000)
m_sigma_x2 = change_local_property (mesh, m_sigma_x2, sigma_reservoir, -10, 10, 0, 1000)
m_sigma_x3 = change_local_property (mesh, m_sigma_x3, sigma_reservoir, -10, 10, 0, 1000)

m_sigma_bgd = change_local_property (mesh, m_sigma_bgd, sigma_bgd2, -210, -10, 0, 1000)
m_sigma_x1 = change_local_property (mesh, m_sigma_x1, sigma_bgd2, -210, -10, 0, 1000)
m_sigma_x2 = change_local_property (mesh, m_sigma_x2, sigma_bgd2, -210, -10, 0, 1000)
m_sigma_x3 = change_local_property (mesh, m_sigma_x3, sigma_bgd2, -210, -10, 0, 1000)

m_sigma_bgd = change_local_property (mesh, m_sigma_bgd, sigma_bgd3, -10000, -210, 0, 1000)
m_sigma_x1 = change_local_property (mesh, m_sigma_x1, sigma_bgd3, -10000, -210, 0, 1000)
m_sigma_x2 = change_local_property (mesh, m_sigma_x2, sigma_bgd3, -10000, -210, 0, 1000)
m_sigma_x3 = change_local_property (mesh, m_sigma_x3, sigma_bgd3, -10000, -210, 0, 1000)

m_sigma_x1 = change_local_property (mesh, m_sigma_x1, sigma_slug, -slug_h/2, slug_h/2, 0, slug_w[0])
m_sigma_x2 = change_local_property (mesh, m_sigma_x2, sigma_slug, -slug_h/2, slug_h/2, 0, slug_w[1])
m_sigma_x3 = change_local_property (mesh, m_sigma_x3, sigma_slug, -slug_h/2, slug_h/2, 0, slug_w[2])

fig, ax = plt.subplots(1, 3, figsize=(15,5))
cb = plt.colorbar(mesh.plotImage(np.log10(m_sigma_x1), ax=ax[0], grid=False)[0], ax=ax[0])
ax[0].set_title(“x=30m”)
ax[0].set_xlim([0, 250])
ax[0].set_ylim([-300, 300])

cb = plt.colorbar(mesh.plotImage(np.log10(m_sigma_x2), ax=ax[1], grid=False)[0], ax=ax[1])
ax[1].set_title(“x=60m”)
ax[1].set_xlim([0, 250])
ax[1].set_ylim([-300, 300])

cb = plt.colorbar(mesh.plotImage(np.log10(m_sigma_x3), ax=ax[2], grid=False)[0], ax=ax[2])
cb.set_label(“Log10(sigma)”)
ax[2].set_title(“x=85m”)
ax[2].set_xlim([0, 250])
ax[2].set_ylim([-300, 300])

fig, ax = plt.subplots(1, 1, figsize=(5,5))
cb = plt.colorbar(mesh.plotImage(m_sigma_bgd, ax=ax, grid=False)[0], ax=ax)
cb.set_label(“Sigma”)
ax.set_title(“BGD”)
ax.set_xlim([0, 300])
ax.set_ylim([-200, 200])

m_mur_one = 1.0 * np.ones(mesh.nC)

m_mur1p25_x1 = 1.0 * np.ones(mesh.nC)
m_mur1p25_x2 = 1.0 * np.ones(mesh.nC)
m_mur1p25_x3 = 1.0 * np.ones(mesh.nC)

m_mur1p25_x1 = change_local_property (mesh, m_mur1p25_x1, 1, -slug_h/2, slug_h/2, 30, 1000000)
m_mur1p25_x2 = change_local_property (mesh, m_mur1p25_x2, 1, -slug_h/2, slug_h/2, 60, 1000000)
m_mur1p25_x3 = change_local_property (mesh, m_mur1p25_x3, 1, -slug_h/2, slug_h/2, 85, 1000000)

fig, ax = plt.subplots(1, 3, figsize=(15,5))
cb = plt.colorbar(mesh.plotImage(m_mur1p25_x1, ax=ax[0], grid=False)[0], ax=ax[0])
ax[0].set_title(“x=30m”)
ax[0].set_xlim([0, 250])
ax[0].set_ylim([-200, 200])

cb = plt.colorbar(mesh.plotImage(m_mur1p25_x2, ax=ax[1], grid=False)[0], ax=ax[1])
ax[1].set_title(“x=60m”)
ax[1].set_xlim([0, 250])
ax[1].set_ylim([-200, 200])

cb = plt.colorbar(mesh.plotImage(m_mur1p25_x3, ax=ax[2], grid=False)[0], ax=ax[2])
cb.set_label(“Mur”)
ax[2].set_title(“x=85m”)
ax[2].set_xlim([0, 250])
ax[2].set_ylim([-200, 200])

orientation = “z”

src_z=np.r_[-150:160:10]
nsrc=len(src_z)

rx_z =np.r_[-150:160:10]
rx_loc=np.zeros((len(rx_z),3))
rx_loc[:,0]=250
rx_loc[:,2]=rx_z

srcList = []
for isrc_z in src_z:
src_loc=np.r_[0.0, 0.0, isrc_z]
rx_real = fdem.Rx.PointMagneticFluxDensity(locations=rx_loc, orientation=orientation, component=“real”)
rx_imag = fdem.Rx.PointMagneticFluxDensity(locations=rx_loc, orientation=orientation, component=“imag”)

``````src = fdem.Src.MagDipole(
loc=src_loc,
orientation=orientation,
freq=freq,
)

srcList.append(src)
``````

survey = fdem.Survey(srcList)
t = time.time()

prob=fdem.Simulation3DMagneticFluxDensity(mesh,survey=survey,sigmaMap=maps.IdentityMap(mesh),Solver=Pardiso) #, mu=mu_0*m_mur1p25_x3)
data0 = prob.make_synthetic_data(m_sigma_x3, relative_error=0.05, noise_floor=1e-15, add_noise=False)
sol3=data0.dclean

print(“Done forward simulation. Elapsed time = {:1.2f} s”.format(time.time() - t))

nsrc=len(src_z)
nrx=len(rx_z)
sol3_amp=np.zeros(nrx*nsrc)

for isrc in range(nsrc):
i_sol3_real=sol3[(2*isrc)nrx:(2isrc+1)nrx]
i_sol3_imag=sol3[(2
isrc+1)nrx:(2isrc+2)nrx]
i_sol3_amp=np.sqrt(i_sol3_real2+i_sol3_imag2)
sol3_amp[isrc
nrx:(isrc+1)*nrx]=i_sol3_amp

img3=np.reshape(sol3_amp, (nsrc, nrx))
plt.imshow(np.log10(img3))
plt.clim(-17.0, -14.7)
plt.colorbar()
plt.title(‘x=85m, mur=1.25’)
plt.xlabel(‘Rx location(m)’)
plt.ylabel(‘Tx location(m)’)
ax0 = plt.gca()
ax0.set_xticks(np.arange(0, 31, 10))
ax0.set_yticks(np.arange(0,31,10))
ax0.set_xticklabels(np.arange(-150, 160, 100))
ax0.set_yticklabels(np.arange(-150, 160, 100))
plt.tight_layout()

dmisfit = data_misfit.L2DataMisfit(simulation=prob, data=data0)
reg = regularization.Simple(mesh)
opt = optimization.InexactGaussNewton(maxIterCG=10, remember=“xc”)
invProb = inverse_problem.BaseInvProblem(dmisfit, reg, opt)

betaest = directives.BetaEstimate_ByEig(beta0_ratio=0.25)
target = directives.TargetMisfit()

directiveList = [betaest, target]
inv = inversion.BaseInversion(invProb, directiveList=directiveList)

print(“The target misfit is {:1.2f}”.format(target.target))

m0 = m_sigma_bgd
t = time.time()
mrec = inv.run(m0)
print("\n Inversion Complete. Elapsed Time = {:1.2f} s".format(time.time() - t))

Hi Evan,

May be you were not able to run inversion I guess.
Had to fix couple of things. Have a look what I have modified in the below curvenote:

https://curvenote.com/@sgkang09/researches/!C6tVMWngm9vU5xfPcQcC.3

I compared signals with and without the slug, and it seems they are almost coincident, so you reach the target misfit with a single iteration… anyway inversion machinery works…

Cheers,

Seogi.

1 Like

Hi Seogi,

I will compare mine with your modified version and learn. Once again, thanks for your kind help.
Yes, this configuration does not have enough sensitivity to the target. If it converges after one iteration, it is what we should expect.

Best regards,
Evan

Hi Seogi and SimPEG community people,

Thank Seogi. I was able to run the example. I am wondering if I could ask a few more questions. In advance, thanks for your kind help.

Best regards,
Evan

1. How can I enforce the lower and upper limit of conductivity values in the inversion? In some examples, I see a script like
opt = optimization.ProjectedGNCG(
maxIter=15, lower=0.0, upper=0.25, maxIterLS=5, maxIterCG=10, tolCG=1e-4
)

In contrast, optimization.InexactGaussNewton cannot take such an parameter.
opt = optimization.InexactGaussNewton(maxIterCG=10, remember=“xc”)

1. Active cells in the inversion
I modified the code to constrain the volume of the inversion domain as shown below but ended up with dimension mismatch error again: The problem would come from # of entire mesh cells != # of cells in the constrained volume.

ind_active1 = (mesh_reg.gridCC[:, 1] <= 10.1)
ind_active2 = (mesh_reg.gridCC[:, 1] >= -10.1)

ind_active3= ind_active1 & ind_active2

nC = int(ind_active3.sum())
model_map = maps.IdentityMap(nP=nC) # model consists of a value for each cell

starting_model = sigma_reservoir * np.ones(nC)

dmisfit = data_misfit.L2DataMisfit(simulation=prob, data=data1)

reg = regularization.Sparse(
mesh_reg,
indActive=ind_active3,
mapping=model_map,
mref=np.log(starting_model),
alpha_s=1,
alpha_x=1,
alpha_y=1,
alpha_z=1,
)

opt = optimization.InexactGaussNewton(maxIterCG=10, remember=“xc”)

invProb = inverse_problem.BaseInvProblem(dmisfit, reg, opt)

betaest = directives.BetaEstimate_ByEig(beta0_ratio=0.25)
target = directives.TargetMisfit()

directiveList = [betaest, target]
inv = inversion.BaseInversion(invProb, directiveList=directiveList)

print(“The target misfit is {:1.2f}”.format(target.target))

m0 = np.log(starting_model)

t = time.time()
mrec = inv.run(m0)
print("\n Inversion Complete. Elapsed Time = {:1.2f} s".format(time.time() - t))