mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-10 00:18:24 +08:00
Work on the time domain inverse problem.
This commit is contained in:
Dave Marchant
+74
-15
@@ -2,6 +2,8 @@ from SimPEG import Utils, Solver
|
||||
from SimPEG.Data import BaseData
|
||||
from SimPEG.Problem import BaseProblem
|
||||
from simpegEM.Utils import Sources
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.Utils import sdiag, mkvc
|
||||
import numpy as np
|
||||
|
||||
class DataTDEM1D(BaseData):
|
||||
@@ -19,9 +21,8 @@ class DataTDEM1D(BaseData):
|
||||
BaseData.__init__(self, **kwargs)
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
|
||||
def dpred(self, sigma, F=None):
|
||||
if F is None: F = self.prob.field(sigma)
|
||||
return self.Qrx.dot(F.b[:,:,0].T)
|
||||
def projectField(self, u):
|
||||
return self.Qrx.dot(u.b[:,:,0].T)
|
||||
|
||||
####################################################
|
||||
# Interpolation Matrices
|
||||
@@ -110,10 +111,56 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
|
||||
"""docstring for ProblemTDEM1D"""
|
||||
def __init__(self, mesh, model, **kwargs):
|
||||
BaseProblem.__init__(self, mesh, model, **kwargs)
|
||||
|
||||
|
||||
####################################################
|
||||
# Physical Properties
|
||||
####################################################
|
||||
|
||||
@property
|
||||
def sigma(self):
|
||||
return self._sigma
|
||||
@sigma.setter
|
||||
def sigma(self, value):
|
||||
self._sigma = value
|
||||
_sigma = None
|
||||
|
||||
####################################################
|
||||
# Mass Matrices
|
||||
####################################################
|
||||
|
||||
@property
|
||||
def MfMui(self): return self._MfMui
|
||||
|
||||
@property
|
||||
def MeSigmaI(self): return self._MeSigmaI
|
||||
|
||||
def makeMassMatrices(self, m):
|
||||
MeSigma = self.mesh.getMass(m, loc='e')
|
||||
self._MeSigmaI = sdiag(1/MeSigma.diagonal())
|
||||
self._MfMui = self.mesh.getMass(1/mu_0, loc='f')
|
||||
|
||||
|
||||
def calcFields(self, sol, solType, tInd):
|
||||
|
||||
if solType == 'b':
|
||||
b = sol
|
||||
e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b
|
||||
# Todo: implement non-zero js
|
||||
else:
|
||||
errStr = 'solType: ' + solType
|
||||
raise NotImplementedError(errStr)
|
||||
|
||||
return {'b':b, 'e':e}
|
||||
|
||||
solveOpts = {'factorize':True,'backend':'scipy'}
|
||||
|
||||
def field(self, m):
|
||||
def field(self, m, useThisRhs=None, useThisCalcFields=None):
|
||||
RHS = useThisRhs or self.getRHS
|
||||
CalcFields = useThisCalcFields or self.calcFields
|
||||
|
||||
self.makeMassMatrices(m)
|
||||
|
||||
F = self.getInitialFields()
|
||||
dtFact = None
|
||||
for tInd, t in enumerate(self.times):
|
||||
@@ -121,15 +168,18 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
|
||||
if dt!=dtFact:
|
||||
dtFact = dt
|
||||
A = self.getA(tInd)
|
||||
print 'Factoring... (dt = ' + str(dt) + ')'
|
||||
# print 'Factoring... (dt = ' + str(dt) + ')'
|
||||
Asolve = Solver(A,options=self.solveOpts)
|
||||
print 'Done'
|
||||
rhs = self.getRHS(tInd, F)
|
||||
# print 'Done'
|
||||
rhs = RHS(tInd, F)
|
||||
sol = Asolve.solve(rhs)
|
||||
if sol.ndim == 1:
|
||||
sol.shape = (sol.size,1)
|
||||
F.update(sol, tInd, self.solType)
|
||||
newFields = CalcFields(sol, self.solType, tInd)
|
||||
F.update(newFields, tInd)
|
||||
return F
|
||||
|
||||
|
||||
|
||||
class FieldsTDEM(object):
|
||||
"""docstring for FieldsTDEM"""
|
||||
@@ -154,12 +204,9 @@ class FieldsTDEM(object):
|
||||
self.nTx = nTx #: Number of transmitters
|
||||
self.mesh = mesh
|
||||
|
||||
def update(self, sol, tInd, solType):
|
||||
if solType == 'b':
|
||||
self.set_b(sol, tInd)
|
||||
else:
|
||||
errStr = 'solType: ' + solType
|
||||
raise NotImplementedError(errStr)
|
||||
def update(self, newFields, tInd):
|
||||
self.set_b(newFields['b'], tInd)
|
||||
self.set_e(newFields['e'], tInd)
|
||||
|
||||
####################################################
|
||||
# Get Methods
|
||||
@@ -171,6 +218,12 @@ class FieldsTDEM(object):
|
||||
else:
|
||||
return self.b[ind,:,:]
|
||||
|
||||
def get_e(self, ind):
|
||||
if ind == -1:
|
||||
return self.e0
|
||||
else:
|
||||
return self.e[ind,:,:]
|
||||
|
||||
####################################################
|
||||
# Set Methods
|
||||
####################################################
|
||||
@@ -179,4 +232,10 @@ class FieldsTDEM(object):
|
||||
if self.b is None:
|
||||
self.b = np.zeros((self.nTimes, np.sum(self.mesh.nF), self.nTx))
|
||||
self.b[:] = np.nan
|
||||
self.b[ind, :] = b
|
||||
self.b[ind,:,:] = b
|
||||
|
||||
def set_e(self, e, ind):
|
||||
if self.e is None:
|
||||
self.e = np.zeros((self.nTimes, np.sum(self.mesh.nE), self.nTx))
|
||||
self.e[:] = np.nan
|
||||
self.e[ind,:,:] = e
|
||||
|
||||
+71
-54
@@ -1,6 +1,6 @@
|
||||
from BaseTDEM import ProblemBaseTDEM
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.Utils import sdiag
|
||||
from BaseTDEM import FieldsTDEM
|
||||
import numpy as np
|
||||
|
||||
class ProblemTDEM_b(ProblemBaseTDEM):
|
||||
"""
|
||||
@@ -10,43 +10,6 @@ class ProblemTDEM_b(ProblemBaseTDEM):
|
||||
ProblemBaseTDEM.__init__(self, mesh, model, **kwargs)
|
||||
|
||||
solType = 'b'
|
||||
|
||||
####################################################
|
||||
# Physical Properties
|
||||
####################################################
|
||||
|
||||
@property
|
||||
def sigma(self):
|
||||
return self._sigma
|
||||
@sigma.setter
|
||||
def sigma(self, value):
|
||||
self._sigma = value
|
||||
_sigma = None
|
||||
|
||||
####################################################
|
||||
# Mass Matrices
|
||||
####################################################
|
||||
|
||||
@property
|
||||
def MfMui(self):
|
||||
if self._MfMui is None:
|
||||
self._MfMui = self.mesh.getMass(1/mu_0, loc='f')
|
||||
return self._MfMui
|
||||
@MfMui.setter
|
||||
def MfMui(self, value):
|
||||
self._MfMui = value
|
||||
_MfMui = None
|
||||
|
||||
@property
|
||||
def MeSigmaI(self):
|
||||
if self._MeSigmaI is None:
|
||||
MeSigma = self.mesh.getMass(self.sigma, loc='e')
|
||||
self._MeSigmaI = sdiag(1/MeSigma.diagonal())
|
||||
return self._MeSigmaI
|
||||
@MeSigmaI.setter
|
||||
def MeSigmaI(self, value):
|
||||
self._MeSigmaI = value
|
||||
_MeSigmaI = None
|
||||
|
||||
####################################################
|
||||
# Internal Methods
|
||||
@@ -58,7 +21,48 @@ class ProblemTDEM_b(ProblemBaseTDEM):
|
||||
|
||||
def getRHS(self, tInd, F):
|
||||
dt = self.getDt(tInd)
|
||||
return (1/dt)*self.MfMui*F.get_b(tInd-1)
|
||||
return (1/dt)*self.MfMui*F.get_b(tInd-1)
|
||||
|
||||
|
||||
####################################################
|
||||
# Derivatives
|
||||
####################################################
|
||||
|
||||
def J(self, m, v, u=None):
|
||||
if u is None:
|
||||
u = self.field(m)
|
||||
p = self.G(m, v, u)
|
||||
y = self.solveAh(m, p)
|
||||
return self.data.projectField(y)
|
||||
|
||||
def G(self, m, v, u=None):
|
||||
if u is None:
|
||||
u = self.field(m)
|
||||
p = FieldsTDEM(self.mesh, 1, self.times.size, 'b')
|
||||
c = self.mesh.getEdgeMassDeriv()*self.model.transformDeriv(m)*v
|
||||
for i in range(self.times.size):
|
||||
ei = u.get_e(i)
|
||||
pVal = np.empty_like(ei)
|
||||
for j in range(ei.shape[1]):
|
||||
pVal[:,j] = -ei[:,j]*c
|
||||
p.set_e(pVal,i)
|
||||
return p
|
||||
|
||||
def solveAh(self, m, p):
|
||||
def AhRHS(tInd, u):
|
||||
if tInd == 0:
|
||||
return self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p.get_e(tInd)
|
||||
else:
|
||||
dt = self.getDt(tInd)
|
||||
return self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p.get_e(tInd) + 1./dt*self.MfMui*u.get_b(tInd-1)
|
||||
|
||||
def AhCalcFields(sol, solType, tInd):
|
||||
b = sol
|
||||
e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b - self.MeSigmaI*p.get_e(tInd)
|
||||
return {'b':b, 'e':e}
|
||||
|
||||
Y = self.field(m, useThisRhs=AhRHS, useThisCalcFields=AhCalcFields)
|
||||
return Y
|
||||
|
||||
if __name__ == '__main__':
|
||||
from SimPEG import *
|
||||
@@ -76,24 +80,37 @@ if __name__ == '__main__':
|
||||
mesh = Mesh.Cyl1DMesh([hx,hy], -hy.sum()/2)
|
||||
model = Model.Vertical1DModel(mesh)
|
||||
|
||||
txLoc = 0.
|
||||
txType = 'VMD_MVP'
|
||||
rxLoc = np.r_[150., 0.]
|
||||
rxType = 'bz'
|
||||
timeCh = np.logspace(-4,-2,20)
|
||||
dat = EM.TDEM.DataTDEM1D(txLoc=txLoc, txType=txType, rxLoc=rxLoc, rxType=rxType, timeCh=timeCh)
|
||||
opts = {'txLoc':0.,
|
||||
'txType':'VMD_MVP',
|
||||
'rxLoc':np.r_[150., 0.],
|
||||
'rxType':'bz',
|
||||
'timeCh':np.logspace(-4,-2,20),
|
||||
}
|
||||
dat = EM.TDEM.DataTDEM1D(**opts)
|
||||
|
||||
prb = EM.TDEM.ProblemTDEM_b(mesh, model)
|
||||
prb.setTimes([1e-5, 5e-5, 2.5e-4], [150, 150, 150])
|
||||
prb.sigma = np.ones(mesh.nCz)*1e-8
|
||||
prb.sigma[mesh.vectorCCz<0] = 0.1
|
||||
sigma = np.ones(mesh.nCz)*1e-8
|
||||
sigma[mesh.vectorCCz<0] = 0.1
|
||||
prb.pair(dat)
|
||||
f = prb.field(sigma)
|
||||
# prb.G(prb.sigma, prb.sigma)
|
||||
# prb.solveAh(prb.sigma, f)
|
||||
# prb.J(prb.sigma, prb.sigma, f)
|
||||
|
||||
bz_calc = dat.dpred(-999999)
|
||||
bz_ana = mu_0*hzAnalyticDipoleT(dat.rxLoc[0], prb.times, prb.sigma[0])
|
||||
from SimPEG.Tests import checkDerivative
|
||||
m0 = sigma
|
||||
dx = np.zeros_like(sigma)
|
||||
dx[prb.mesh.vectorCCz<0] = 1e-4
|
||||
derChk = lambda m: [dat.dpred(m), lambda mx: prb.J(m0, mx, u=f)]
|
||||
passed = checkDerivative(derChk, m0, dx=dx, plotIt=False)
|
||||
|
||||
# bz_calc = dat.dpred(sigma)
|
||||
# bz_ana = mu_0*hzAnalyticDipoleT(dat.rxLoc[0], prb.times, sigma[0])
|
||||
|
||||
# plt.loglog(prb.times, np.abs(bz_calc.flatten()), label='TDEM_b')
|
||||
# plt.loglog(prb.times, np.abs(bz_ana), 'r', label='Analytic')
|
||||
# plt.legend()
|
||||
# plt.show()
|
||||
|
||||
|
||||
plt.loglog(prb.times, np.abs(bz_calc.flatten()), label='TDEM_b')
|
||||
plt.loglog(prb.times, np.abs(bz_ana), 'r', label='Analytic')
|
||||
plt.legend()
|
||||
plt.show()
|
||||
|
||||
@@ -1,15 +1,52 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
import simpegEM as EM
|
||||
from scipy.constants import mu_0
|
||||
from simpegEM.Utils.Ana import hzAnalyticDipoleT
|
||||
|
||||
|
||||
class EMProblemTests(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
pass
|
||||
|
||||
cs = 5.
|
||||
ncx = 20
|
||||
ncy = 6
|
||||
npad = 20
|
||||
hx = Utils.meshTensors(((0,cs), (ncx,cs), (npad,cs)))
|
||||
hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs)))
|
||||
mesh = Mesh.Cyl1DMesh([hx,hy], -hy.sum()/2)
|
||||
model = Model.Vertical1DModel(mesh)
|
||||
|
||||
opts = {'txLoc':0.,
|
||||
'txType':'VMD_MVP',
|
||||
'rxLoc':np.r_[150., 0.],
|
||||
'rxType':'bz',
|
||||
'timeCh':np.logspace(-4,-2,20),
|
||||
}
|
||||
self.dat = EM.TDEM.DataTDEM1D(**opts)
|
||||
|
||||
self.prb = EM.TDEM.ProblemTDEM_b(mesh, model)
|
||||
self.prb.setTimes([1e-5, 5e-5, 2.5e-4], [150, 150, 150])
|
||||
self.sigma = np.ones(mesh.nCz)*1e-8
|
||||
self.sigma[mesh.vectorCCz<0] = 0.1
|
||||
self.prb.pair(self.dat)
|
||||
|
||||
|
||||
def test_analitic_b(self):
|
||||
bz_calc = self.dat.dpred(self.sigma)
|
||||
bz_ana = mu_0*hzAnalyticDipoleT(self.dat.rxLoc[0], self.prb.times, self.sigma[0])
|
||||
|
||||
diff = np.linalg.norm(bz_calc.flatten() - bz_ana.flatten())/np.linalg.norm(bz_ana.flatten())
|
||||
print diff
|
||||
self.assertTrue(diff<0.05)
|
||||
|
||||
def test_awesome(self):
|
||||
self.assertTrue(True)
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
Reference in New Issue
Block a user