Hi all,
I’m attempting to figure out if simpeg is the way to go to simulate ground based TDEM (so dB/dt) for situations where I’ve run a large source loop and have receivers logging dB/dt measurements on points along lines that are both inside and outside the loop. I’m quite sure 3D effects are significant given that I’m looking for highly conductive linear bodies at various depths and dips. Happy to take the dive but would be grateful for any leads.
Cheers
Julien
I have since “taken the dive” and am modifying several examples to simulate a 3D TDEM system based on a single transmitting loop and multiple surface receivers.
I am currently running to the following error:
ValueError: Points outside of mesh
A similar issue was raised in this thread:
I have checked to make sure that my receivers are not attempting to record samples outside the time stepping range (they do not). Using verbose mode, the simulation does indeed make it through all the time steps, and crashes afterward. Here is the code I used (I’m inputing a highly conductive dipping “filament” type body into a background high resistivity body). Any ideas?
from discretize import TensorMesh
from discretize import TreeMesh
from discretize.utils import mkvc, refine_tree_xyz
import matplotlib.pyplot as plt
from SimPEG.utils import plot2Ddata, surface2ind_topo
from SimPEG import maps
import SimPEG.electromagnetics.time_domain as tdem
from SimPEG.utils import meshTensor
import numpy as np
import numpy.matlib as matlib
import matplotlib as mpl
import matplotlib.pyplot as plt
import os
import plotly.offline as pyo
import plotly.graph_objs as go
import plotly.io as pio
pio.renderers.default='browser'
%matplotlib inline
try:
from pymatsolver import Pardiso as Solver
except ImportError:
from SimPEG import SolverLU as Solver
save_file = False
'''
Create a mesh
'''
nc = 20 # number of core mesh cells in x, y and z
ncz = 30
dh = 100 # base cell width in x, y and z
dhz = 5
npad = 5 # number of padding cells
exp = 1.25 # expansion rate of padding cells
hx = [(dh, nc)] #reads like this: dilation left of zero, core cells, dilation positive from zero
hy = [(dh, nc)]
hz = [(dhz, 7, -exp), (dhz, ncz)]
mesh = TensorMesh([hx, hy, hz], x0="CCN") # creates the mesh
# The bottom southwest corner
x0 = mesh.x0 #indicates the southwest deepest cell
# The total number of cells
nC = mesh.nC
# An (nC, 3) array containing the cell-center locations
cc = mesh.gridCC
# A boolean array specifying which cells lie on the boundary
bInd = mesh.cellBoundaryInd
# The cell volumes
v = mesh.vol
# Plot all cells volumes or plot cell volumes for a particular horizontal slice
# fig = plt.figure(figsize=(9, 4))
# ax1 = fig.add_subplot(121)
# ax2 = fig.add_subplot(122)
# mesh.plotImage(np.log10(v), grid=True, ax=ax1)
# ax1.set_title("All Cell Log-Volumes")
# cplot = mesh.plotSlice(np.log10(v), grid=True, ax=ax2, normal="Z", ind=2)
# cplot[0].set_clim(np.min(np.log10(v)), np.max(np.log10(v)))
# ax2.set_title("Cell Log-Volumes #2")
air_value = 0.0
ind_active = mesh.gridCC[:, 2] < 0.0
model_map = maps.InjectActiveCells(mesh, ind_active, air_value)
'''
Define the conductivity model
'''
sig_silver = 1/1e-4
sig_background = 1/2000
model = np.log(sig_background) * np.ones(mesh.nC) #background model
#create line in 3D
xo,yo,zo = (0,0,-50) #passes through this point
x1,y1,z1 = (400,400,-180)
Lx,Ly,Lz = [np.linspace(xo,x1,100),np.linspace(yo,y1,100),np.linspace(zo,z1,100)]
#function to calculate the distance of every cell to a point
def dist(x,msh):
dis = np.zeros(len(msh))
p1 = np.array([mesh.gridCC[:,0], mesh.gridCC[:,1], mesh.gridCC[:,2]])
p2 = np.transpose(matlib.repmat(np.transpose(x),len(p1[0]),1))
dist = (((p2[0,:]-p1[0,:])**2)+((p2[1,:]-p1[1,:])**2)+((p2[2,:]-p1[2,:])**2))**(1/2)
return dist
#Note: mesh.gridCC[x1,xyz] returns one element (x) and its 3d x,y,z position in the volume
thk = 20 #vein thickness
for i in range(len(Lx)):
ds = dist(np.array([Lx[i],Ly[i],Lz[i]]),mesh)
model[ds < thk] = np.log(sig_silver)
#display model as an isosurface plot
fig= go.Figure(data=go.Isosurface(
x= mesh.gridCC[:,0],
y= mesh.gridCC[:,1],
z= mesh.gridCC[:,2],
value=model,
isomin=-6,
isomax=10,
))
fig.show()
'''
Define the receivers
'''
# Observation times for response (time channels)
n_times = 31
time_channels = np.logspace(-4, -3, n_times)
n_tx = 11
# Define receiver locations
xrx, yrx, zrx = np.meshgrid(
np.linspace(-500, 500, n_tx), np.linspace(-190, 190, n_tx), [30]
)
receiver_locations = np.c_[mkvc(xrx), mkvc(yrx), mkvc(zrx)]
receivers_list = []
ntx = np.size(xrx)
# Each unique location defines a new transmitter
for ii in range(ntx):
# Here we define receivers that measure the h-field in A/m
dbzdt_receiver = tdem.receivers.PointMagneticFluxTimeDerivative(
receiver_locations[ii, :], time_channels, "z"
)
receivers_list = [
dbzdt_receiver
] # Make a list containing all receivers even if just one
'''
Define the source waveforms
'''
waveform_times = np.linspace(-7.81e-3, 0, 51)
# Set this up to emulate nanotem
waveform = tdem.sources.TrapezoidWaveform(
ramp_on=np.r_[-7.81e-3, -7.8096e-3], ramp_off=np.r_[-0.0004e-3, 0.0], off_time=0.0
)
waveform_value = [waveform.eval(t) for t in waveform_times]
'''
Define the source
'''
src_ramp_on = tdem.Src.CircularLoop(
receivers_list,
location=np.r_[0.0, 0.0, 0.0],
orientation="z",
radius=100,
waveform=waveform,
)
src_list_ramp_on = [src_ramp_on]
survey_ramp_on = tdem.Survey(src_list_ramp_on)
survey = tdem.Survey(src_list_ramp_on)
'''
Simulation
'''
# dt and how many dt will result (so here we have 40 time steps in total)
time_steps = [(1e-5, 20), (1e-5, 10), (1e-4, 30), (1e-3, 20)]
#check total time
dt = meshTensor(time_steps)
tmax = np.hstack([[0], dt]).cumsum()
#Note: waveform on time was defined as -0.002 to 0
simulation = tdem.simulation.Simulation3DMagneticFluxDensity(
mesh, survey=survey, sigmaMap=model_map,verbose=True, solver=Solver, t0=-7.81e-3
)
# Set the time-stepping for the simulation
simulation.time_steps = time_steps
dpred = simulation.dpred(model)Solved; stupid error (receivers were above the surface of the mesh, which was not defined).