Crosswell EM inversion in 2d cylindrical coordinate

evan · March 15, 2021, 3:20pm

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

sgkang09 · March 17, 2021, 6:40pm

Hi Evan,

Glad to hear that you are using SimPEG!

How about this example?

github.com

simpeg/simpeg/blob/master/examples/06-tdem/plot_inv_tdem_1D.py

"""
EM: TDEM: 1D: Inversion
=======================

Here we will create and run a TDEM 1D inversion.
"""
import numpy as np
from SimPEG.electromagnetics import time_domain
from SimPEG import (
    optimization,
    discretize,
    maps,
    data_misfit,
    regularization,
    inverse_problem,
    inversion,
    directives,
    utils,
)
import matplotlib.pyplot as plt

This file has been truncated. show original

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

All the best,

Seogi.

evan · March 17, 2021, 7:47pm

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
npad_r = 40
npad_z = 40
exp_r = 1.15
exp_z = 1.15

hr = [(dr, ncr), (dr, npad_r, exp_r)]
hz = [(dz, npad_z, -exp_z), (dz, ncz), (dz, npad_z, exp_z)]

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(
    receiver_list=[rx_real, rx_imag],
    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[(2isrc+1)nrx:(2isrc+2)nrx]
i_sol3_amp=np.sqrt(i_sol3_real2+i_sol3_imag2)
sol3_amp[isrcnrx:(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))

sgkang09 · March 19, 2021, 4:50am

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.

evan · March 19, 2021, 5:18am

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

evan · March 25, 2021, 2:08am

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

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”)

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))

Topic		Replies	Views
Trying to model and invert the data set Electromagnetics	4	1272	February 12, 2020
FDEM HLEM coplanar loop inversion Electromagnetics	0	119	December 5, 2023
About 2d tdem forward and inversion Electromagnetics	0	278	September 22, 2023
3d effects in TDEM Electromagnetics	4	1001	July 20, 2020
Ground based loop TDEM amplitude anomalies Electromagnetics	2	606	June 10, 2022

Crosswell EM inversion in 2d cylindrical coordinate

Related topics