-
Notifications
You must be signed in to change notification settings - Fork 0
/
grnded3raw.py
95 lines (86 loc) · 3.24 KB
/
grnded3raw.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
from SimPEG import *
from SimPEG import EM
from scipy.constants import mu_0
import numpy as np
import scipy.sparse as sp
from simpegEMIP.TDEM import geteref, Problem3DIP_Linear, SurveyLinear
from simpegEMIP.TDEM import Survey, Rx
from simpegEMIP.TDEM import Problem3DEM_e, Problem3D_e
import matplotlib.pyplot as plt
from pymatsolver import PardisoSolver
eta, tau, c = 0.1, 0.01, 0.5
cs, ncx, ncz, npad = 10., 25, 20, 18
hx = [(cs,ncx), (cs,npad,1.3)]
hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
mesh = Mesh.CylMesh([hx,1,hz], '00C')
sigmaInf = np.ones(mesh.nC) * 0.001
actinds = mesh.gridCC[:,2]<0.
layerind = (np.logical_and(mesh.gridCC[:,2]<0, mesh.gridCC[:,2]>-50.)) & (mesh.gridCC[:,0]<100.)
sigmaInf[~actinds] = 1e-8
sigmaInf[layerind] = 0.1
eta = np.zeros(mesh.nC)
eta[layerind] = 0.5
tau = np.ones(mesh.nC) * 1.
c = np.ones(mesh.nC) * 0.5
def get_em_data(sigma, eta=None, tau=None, c=None, data_type='em'):
rxloc = np.array([[0., 0., 30.]])
srcloc = np.array([[0., 0., 30.]])
dt = 1.47e-3
tpeak = 2.73e-3
t0 = tpeak + dt
rx_vtem = Rx.Point_dbdt(rxloc, np.logspace(np.log10(2e-5), np.log10(0.009), 51)+t0, orientation='z')
src_vtem = EM.TDEM.Src.CircularLoop([rx_vtem], waveform=EM.TDEM.Src.VTEMWaveform(offTime=t0, peakTime=tpeak, a=3.), loc=srcloc)
survey = Survey([src_vtem])
if data_type == 'em':
prb = Problem3DEM_e(mesh, sigma=sigma)
elif data_type == 'emip':
prb = Problem3D_e(mesh, sigmaInf=sigma, eta=eta, tau=tau, c=c)
prb.timeSteps = [(tpeak/10, 10), ((t0-tpeak)/10, 10), (1e-06, 5), (2.5e-06, 5), (5e-06, 5), (1e-05, 10), (2e-05, 10), (4e-05, 10), (8e-05, 10), (1.6e-04, 10), (3.2e-04, 20)]
prb.Solver = PardisoSolver
prb.pair(survey)
e = prb.fields(sigmaInf)
data = survey.dpred(sigmaInf, f=e)
# waveform
cur = []
for t in prb.times:
cur.append(src_vtem.waveform.eval(t))
cur = np.hstack(cur)
return e, data, cur
def get_ip_data(sigma, eref, eta, tau, c):
rxloc = np.array([[0., 0., 30.]])
srcloc = np.array([[0., 0., 30.]])
rx_ip = Rx.Point_dbdt(rxloc, np.logspace(np.log10(2e-5), np.log10(0.009), 51), 'z')
src_ip = EM.TDEM.Src.CircularLoop([rx_ip], loc=srcloc)
dt = 1.47e-3
tpeak = 2.73e-3
t0 = tpeak + dt
survey_ip = SurveyLinear([src_ip])
t1, t2, t3 = dt, t0-0.001365, t0
prb_ip = Problem3DIP_Linear(
mesh,
sigmaInf=sigmaInf,
eta=eta,
tau=tau,
c=c,
actinds = actinds,
tlags = [0., t1, t2, t3]
)
prb_ip.Solver = PardisoSolver
prb_ip.pair(survey_ip)
prb_ip.set_eref(eref)
ip_approx = survey_ip.dpred([])
return ip_approx
mesh.plotImage(np.log10(sigmaInf))
plt.xlim(0,200)
plt.ylim(-200,200)
e_emip, data_emip, cur = get_em_data(sigmaInf, eta, tau, c, data_type='emip')
e_em, data_em, cur = get_em_data(sigmaInf)
eref = geteref(e_em[:,0,:], mesh, option=None, tInd=20)
ip = get_ip_data(sigmaInf, eref, eta, tau, c)
data = ip + data_em
times = np.logspace(np.log10(2e-5), np.log10(0.009), 51)
plt.loglog(times, -data, 'k')
plt.loglog(times, data, '--k')
plt.loglog(times, -data_emip, 'kx')
plt.loglog(times, data_emip, 'kx')
plt.loglog(times, -data_em, 'b')