Hi all,
As a continuation of my efforts to simulate a 3D TDEM survey with a source loop and a grid of receivers, I have come up with the following code, which I’ll post and annotate.
The code runs without errors, but the dB/dt values in the final output are clearly impossible (~e123). Here is what I am doing.
-I have a “filament-like” feature, depicted by a highly conductive dipping tube in an otherwise high resistivity model.
-The source is a ground based (z=0) loop with a radius of ~28 meters.
-Receivers are a uniform grid also on the surface, and record at 31 logspaced points between 1e-7 and 1e-3 (though I have tried various other time ranges).
-I am using a ramp on ramp off style of waveform (I tried others also, always with anomalous amplitudes as a result).
Any help in resolving this would be much appreciated. Here’s the code in blocks:
I am using a few extra packages, namely plotly to visualize the survey.
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.express as px
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
Set up the mesh. I’m using a tensormesh here, with padding. The plotting is commented out. The last line of this section involves the model mapping, which I find somewhat opaque to understand (it seems to have to do with flagging the active cells in a mesh, but in my case ind_active returns all True, since it is not an air-based TDEM study. Is this correct?)
nc = 40 # number of core mesh cells in x, y and z
ncz = 24
dh = 10 # base cell width in x, y and z
dhz = 5
npad = 5 # number of padding cells: these are cells of varying size OUTSIDE the code cell region (usually less dense)
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, 4, -exp), (dhz, ncz)]
mesh = TensorMesh([hx, hy, hz], x0="CCN") # creates the mesh reference
# 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 = 1e-8
ind_active = mesh.gridCC[:, 2] < 0.0
model_map = maps.InjectActiveCells(mesh, ind_active, air_value)
Define the conductivity model. sig_silver is the conductivity of the dipping tube feature.
sig_silver = 10000
sig_background = 2e-3
model = np.log(sig_background) * np.ones(mesh.nC) #background model
#create line in 3D
xo,yo,zo = (0,0,-20) #passes through this point
x1,y1,z1 = (400,400,-250)
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 = 10 #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)
Define the receiver setup. Simple uniform grid of surface receivers of the dB/dt type.
# Observation times for response (time channels)
n_times = 31
time_channels = np.logspace(-7, -3, n_times)
n_tx = 11
# Define receiver locations
xrx, yrx, zrx = np.meshgrid(
np.linspace(-190, 190, n_tx), np.linspace(-190, 190, n_tx), [0]
)
receiver_locations = np.c_[mkvc(xrx), mkvc(yrx), mkvc(zrx)]
receivers_list = []
ntx = np.size(xrx)
# Each unique location defines a new receiver
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, orientation = "z"
)
receivers_list.append(dbzdt_receiver)
Source waveform is a step on, step off waveform, with a stable on time plateau and a final off time of 0s.
waveform_times = np.linspace(-0.002, 0, 21)
waveform = tdem.sources.TrapezoidWaveform(
ramp_on=np.r_[-0.007, -0.006], ramp_off=np.r_[-0.001, 0.0], offTime=0.0
)
waveform_value = [waveform.eval(t) for t in waveform_times]
For those who want visualize the survey, here it the plotly script I’m using.
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.add_scatter3d(
x=receiver_locations[:,0],
y=receiver_locations[:,1],
z=receiver_locations[:,2],
mode='markers',
)
fig.show()
Here is the source layout as a loop with a radius of 28m and a current of 2.5A centered at [0,0,0].
src_ramp_on = tdem.Src.CircularLoop(
receivers_list,
location=np.r_[0.0, 0.0, 0.0],
orientation = "z",
current = 2.5,
radius=28,
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)
Run the simulation.
time_steps = [(1e-5, 20), (1e-5, 10), (1e-4, 30), (1e-3, 5)]
#check total time
dt = meshTensor(time_steps)
tmax = np.hstack([[0], dt]).cumsum()
#Note: waveform on time was defined as -0.007 to 0
simulation = tdem.simulation.Simulation3DMagneticFluxDensity(
mesh, survey=survey, sigmaMap=model_map,verbose=True, solver=Solver, t0=-0.007
)
# Set the time-stepping for the simulation
simulation.time_steps = time_steps
dpred = simulation.dpred(model)
Plot the first and last times of the recordings.
# Data were organized by location, then by time channel
dpred_plotting = np.reshape(dpred, (n_tx ** 2, n_times))
# Plot
fig = plt.figure(figsize=(10, 4))
# dB/dt at early time
v_max = np.max(np.abs(dpred_plotting[:, 0]))
ax11 = fig.add_axes([0.05, 0.05, 0.35, 0.9])
plot2Ddata(
receiver_locations[:, 0:2],
dpred_plotting[:, 0],
ax=ax11,
ncontour=30,
clim=(-v_max, v_max),
contourOpts={"cmap": "bwr"},
)
ax11.set_title("dBz/dt at 0.0000001 s")
ax12 = fig.add_axes([0.42, 0.05, 0.02, 0.9])
norm1 = mpl.colors.Normalize(vmin=-v_max, vmax=v_max)
cbar1 = mpl.colorbar.ColorbarBase(
ax12, norm=norm1, orientation="vertical", cmap=mpl.cm.bwr
)
cbar1.set_label("$T/s$", rotation=270, labelpad=15, size=12)
# dB/dt at later times
v_max = np.max(np.abs(dpred_plotting[:, -1]))
ax21 = fig.add_axes([0.55, 0.05, 0.35, 0.9])
plot2Ddata(
receiver_locations[:, 0:2],
dpred_plotting[:, -1],
ax=ax21,
ncontour=30,
clim=(-v_max, v_max),
contourOpts={"cmap": "bwr"},
)
ax21.set_title("dBz/dt at 0.001 s")
ax22 = fig.add_axes([0.92, 0.05, 0.02, 0.9])
norm2 = mpl.colors.Normalize(vmin=-v_max, vmax=v_max)
cbar2 = mpl.colorbar.ColorbarBase(
ax22, norm=norm2, orientation="vertical", cmap=mpl.cm.bwr
)
cbar2.set_label("$T/s$", rotation=270, labelpad=15, size=12)
plt.show()
The variable dpred_plotting contains an array of the 121x31 of dB/dt. Plotting these results in clearly broken amplitudes and oscillatory behavior, and I’ve tried numerous ways to fix this.
I think I am not setting up the survey correctly (i.e. something to do with putting the source and receivers directly on the boundary). Does anyone have ideas?
Thanks much,
Julien