mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-24 13:20:06 +08:00
Updates to TDEM (Jtvec still not working.)
This commit is contained in:
rowanc1
@@ -58,7 +58,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
|
||||
def adjoint(self, m, RHS, CalcFields, F=None):
|
||||
if F is None:
|
||||
F = FieldsTDEM(self.mesh, self.survey.nTx, self.nT, store=self.storeTheseFields)
|
||||
F = FieldsTDEM(self.mesh, self.survey)
|
||||
|
||||
dtFact = None
|
||||
for tInd, dt in reversed(list(enumerate(self.timeSteps))):
|
||||
@@ -73,6 +73,6 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
if sol.ndim == 1:
|
||||
sol.shape = (sol.size,1)
|
||||
newFields = CalcFields(sol, self.solType, tInd)
|
||||
F.update(newFields, tInd)
|
||||
F[:,:,tInd] = newFields
|
||||
return F
|
||||
|
||||
|
||||
+160
-127
@@ -37,15 +37,32 @@ class RxTDEM(Survey.BaseTimeRx):
|
||||
P = self.getP(mesh, timeMesh)
|
||||
|
||||
if not adjoint:
|
||||
return P * v
|
||||
return P * Utils.mkvc(v[tx, self.projField, :])
|
||||
elif adjoint:
|
||||
return P.T * v
|
||||
Ptv = P.T * v[tx, self]
|
||||
return Ptv
|
||||
|
||||
|
||||
class FieldsTDEM(Survey.TimeFields):
|
||||
"""Fancy Field Storage for a TDEM survey."""
|
||||
knownFields = {'b': 'F', 'e': 'E'}
|
||||
|
||||
def tovec(self):
|
||||
nTx, nF, nE = self.survey.nTx, self.mesh.nF, self.mesh.nE
|
||||
u = np.empty(0 if nTx == 1 else (0, nTx))
|
||||
|
||||
for i in range(self.survey.prob.nT):
|
||||
if 'b' in self:
|
||||
b = self[:,'b',i+1]
|
||||
else:
|
||||
b = np.zeros(nF if nTx == 1 else (nF, nTx))
|
||||
|
||||
if 'e' in self:
|
||||
e = self[:,'e',i+1]
|
||||
else:
|
||||
e = np.zeros(nE if nTx == 1 else (nE, nTx))
|
||||
u = np.r_[u, b, e]
|
||||
return u
|
||||
|
||||
class TxTDEM(Survey.BaseTx):
|
||||
rxPair = RxTDEM
|
||||
@@ -94,132 +111,148 @@ class SurveyTDEM(Survey.BaseSurvey):
|
||||
data[tx, rx] = rx.projectFields(tx, self.mesh, self.prob.timeMesh, u)
|
||||
return data
|
||||
|
||||
def projectFieldsDeriv(self, u):
|
||||
raise Exception('Use Transmitters to project fields deriv.')
|
||||
def projectFieldsDeriv(self, u, v=None, adjoint=False):
|
||||
assert v is not None, 'v to multiply must be provided.'
|
||||
|
||||
|
||||
class SurveyTDEM1D(BaseSurvey):
|
||||
"""
|
||||
docstring for SurveyTDEM1D
|
||||
"""
|
||||
|
||||
txLoc = None #: txLoc
|
||||
txType = None #: txType
|
||||
rxLoc = None #: rxLoc
|
||||
rxType = None #: rxType
|
||||
timeCh = None #: timeCh
|
||||
nTx = 1 #: Number of transmitters
|
||||
|
||||
@property
|
||||
def nTimeCh(self):
|
||||
"""Number of time channels"""
|
||||
return self.timeCh.size
|
||||
|
||||
def __init__(self, **kwargs):
|
||||
BaseSurvey.__init__(self, **kwargs)
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
|
||||
def projectFields(self, u):
|
||||
#TODO: this is hardcoded to 1Tx
|
||||
return self.Qrx.dot(u.b[:,:,0].T).T
|
||||
|
||||
def projectFieldsAdjoint(self, d):
|
||||
# TODO: make the following self.nTimeCh
|
||||
d = d.reshape((self.prob.nT, self.nTx), order='F')
|
||||
#TODO: *Qtime.T need to multiply by a time projection. (outside for loop??)
|
||||
ii = 0
|
||||
F = FieldsTDEM(self.prob.mesh, self.nTx, self.prob.nT, 'b')
|
||||
for ii in range(self.prob.nT):
|
||||
b = self.Qrx.T*d[ii,:]
|
||||
F.set_b(b, ii)
|
||||
F.set_e(np.zeros((self.prob.mesh.nE,self.nTx)), ii)
|
||||
return F
|
||||
|
||||
####################################################
|
||||
# Interpolation Matrices
|
||||
####################################################
|
||||
|
||||
@property
|
||||
def Qrx(self):
|
||||
if self._Qrx is None:
|
||||
if self.rxType == 'bz':
|
||||
locType = 'Fz'
|
||||
self._Qrx = self.prob.mesh.getInterpolationMat(self.rxLoc, locType=locType)
|
||||
return self._Qrx
|
||||
_Qrx = None
|
||||
|
||||
|
||||
class FieldsTDEM_OLD(object):
|
||||
"""docstring for FieldsTDEM"""
|
||||
|
||||
phi0 = None #: Initial electric potential
|
||||
A0 = None #: Initial magnetic vector potential
|
||||
e0 = None #: Initial electric field
|
||||
b0 = None #: Initial magnetic flux density
|
||||
j0 = None #: Initial current density
|
||||
h0 = None #: Initial magnetic field
|
||||
|
||||
phi = None #: Electric potential
|
||||
A = None #: Magnetic vector potential
|
||||
e = None #: Electric field
|
||||
b = None #: Magnetic flux density
|
||||
j = None #: Current density
|
||||
h = None #: Magnetic field
|
||||
|
||||
def __init__(self, mesh, nTx, nT, store='b'):
|
||||
|
||||
self.nT = nT #: Number of times
|
||||
self.nTx = nTx #: Number of transmitters
|
||||
self.mesh = mesh
|
||||
|
||||
def update(self, newFields, tInd):
|
||||
self.set_b(newFields['b'], tInd)
|
||||
self.set_e(newFields['e'], tInd)
|
||||
|
||||
def fieldVec(self):
|
||||
u = np.ndarray((0, self.nTx))
|
||||
for i in range(self.nT):
|
||||
u = np.r_[u, self.get_b(i), self.get_e(i)]
|
||||
if self.nTx == 1:
|
||||
u = u.flatten()
|
||||
return u
|
||||
|
||||
####################################################
|
||||
# Get Methods
|
||||
####################################################
|
||||
|
||||
def get_b(self, ind):
|
||||
if ind == -1:
|
||||
return self.b0
|
||||
if not adjoint:
|
||||
data = Survey.Data(self)
|
||||
for tx in self.txList:
|
||||
for rx in tx.rxList:
|
||||
data[tx, rx] = rx.projectFieldsDeriv(tx, self.mesh, self.prob.timeMesh, u, v)
|
||||
return data
|
||||
else:
|
||||
return self.b[ind,:,:]
|
||||
|
||||
def get_e(self, ind):
|
||||
if ind == -1:
|
||||
return self.e0
|
||||
else:
|
||||
return self.e[ind,:,:]
|
||||
|
||||
####################################################
|
||||
# Set Methods
|
||||
####################################################
|
||||
|
||||
def set_b(self, b, ind):
|
||||
if self.b is None:
|
||||
self.b = np.zeros((self.nT, np.sum(self.mesh.nF), self.nTx))
|
||||
self.b[:] = np.nan
|
||||
if len(b.shape) == 1:
|
||||
b = b[:, np.newaxis]
|
||||
self.b[ind,:,:] = b
|
||||
|
||||
def set_e(self, e, ind):
|
||||
if self.e is None:
|
||||
self.e = np.zeros((self.nT, np.sum(self.mesh.nE), self.nTx))
|
||||
self.e[:] = np.nan
|
||||
if len(e.shape) == 1:
|
||||
e = e[:, np.newaxis]
|
||||
self.e[ind,:,:] = e
|
||||
f = FieldsTDEM(self.mesh, self)
|
||||
for tx in self.txList:
|
||||
for rx in tx.rxList:
|
||||
Ptv = rx.projectFieldsDeriv(tx, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
|
||||
Ptv = Ptv.reshape((-1, 1, self.prob.timeMesh.nN), order='F')
|
||||
f[tx, rx.projField, :] = Ptv
|
||||
return f
|
||||
|
||||
|
||||
def __contains__(self, key):
|
||||
return key in self.children
|
||||
|
||||
# class SurveyTDEM1D(BaseSurvey):
|
||||
# """
|
||||
# docstring for SurveyTDEM1D
|
||||
# """
|
||||
|
||||
# txLoc = None #: txLoc
|
||||
# txType = None #: txType
|
||||
# rxLoc = None #: rxLoc
|
||||
# rxType = None #: rxType
|
||||
# timeCh = None #: timeCh
|
||||
# nTx = 1 #: Number of transmitters
|
||||
|
||||
# @property
|
||||
# def nTimeCh(self):
|
||||
# """Number of time channels"""
|
||||
# return self.timeCh.size
|
||||
|
||||
# def __init__(self, **kwargs):
|
||||
# BaseSurvey.__init__(self, **kwargs)
|
||||
# Utils.setKwargs(self, **kwargs)
|
||||
|
||||
# def projectFields(self, u):
|
||||
# #TODO: this is hardcoded to 1Tx
|
||||
# return self.Qrx.dot(u.b[:,:,0].T).T
|
||||
|
||||
# def projectFieldsAdjoint(self, d):
|
||||
# # TODO: make the following self.nTimeCh
|
||||
# d = d.reshape((self.prob.nT, self.nTx), order='F')
|
||||
# #TODO: *Qtime.T need to multiply by a time projection. (outside for loop??)
|
||||
# ii = 0
|
||||
# F = FieldsTDEM(self.prob.mesh, self.nTx, self.prob.nT, 'b')
|
||||
# for ii in range(self.prob.nT):
|
||||
# b = self.Qrx.T*d[ii,:]
|
||||
# F.set_b(b, ii)
|
||||
# F.set_e(np.zeros((self.prob.mesh.nE,self.nTx)), ii)
|
||||
# return F
|
||||
|
||||
# ####################################################
|
||||
# # Interpolation Matrices
|
||||
# ####################################################
|
||||
|
||||
# @property
|
||||
# def Qrx(self):
|
||||
# if self._Qrx is None:
|
||||
# if self.rxType == 'bz':
|
||||
# locType = 'Fz'
|
||||
# self._Qrx = self.prob.mesh.getInterpolationMat(self.rxLoc, locType=locType)
|
||||
# return self._Qrx
|
||||
# _Qrx = None
|
||||
|
||||
|
||||
# class FieldsTDEM_OLD(object):
|
||||
# """docstring for FieldsTDEM"""
|
||||
|
||||
# phi0 = None #: Initial electric potential
|
||||
# A0 = None #: Initial magnetic vector potential
|
||||
# e0 = None #: Initial electric field
|
||||
# b0 = None #: Initial magnetic flux density
|
||||
# j0 = None #: Initial current density
|
||||
# h0 = None #: Initial magnetic field
|
||||
|
||||
# phi = None #: Electric potential
|
||||
# A = None #: Magnetic vector potential
|
||||
# e = None #: Electric field
|
||||
# b = None #: Magnetic flux density
|
||||
# j = None #: Current density
|
||||
# h = None #: Magnetic field
|
||||
|
||||
# def __init__(self, mesh, nTx, nT, store='b'):
|
||||
|
||||
# self.nT = nT #: Number of times
|
||||
# self.nTx = nTx #: Number of transmitters
|
||||
# self.mesh = mesh
|
||||
|
||||
# def update(self, newFields, tInd):
|
||||
# self.set_b(newFields['b'], tInd)
|
||||
# self.set_e(newFields['e'], tInd)
|
||||
|
||||
# def fieldVec(self):
|
||||
# u = np.ndarray((0, self.nTx))
|
||||
# for i in range(self.nT):
|
||||
# u = np.r_[u, self.get_b(i), self.get_e(i)]
|
||||
# if self.nTx == 1:
|
||||
# u = u.flatten()
|
||||
# return u
|
||||
|
||||
# ####################################################
|
||||
# # Get Methods
|
||||
# ####################################################
|
||||
|
||||
# def get_b(self, ind):
|
||||
# if ind == -1:
|
||||
# return self.b0
|
||||
# else:
|
||||
# return self.b[ind,:,:]
|
||||
|
||||
# def get_e(self, ind):
|
||||
# if ind == -1:
|
||||
# return self.e0
|
||||
# else:
|
||||
# return self.e[ind,:,:]
|
||||
|
||||
# ####################################################
|
||||
# # Set Methods
|
||||
# ####################################################
|
||||
|
||||
# def set_b(self, b, ind):
|
||||
# if self.b is None:
|
||||
# self.b = np.zeros((self.nT, np.sum(self.mesh.nF), self.nTx))
|
||||
# self.b[:] = np.nan
|
||||
# if len(b.shape) == 1:
|
||||
# b = b[:, np.newaxis]
|
||||
# self.b[ind,:,:] = b
|
||||
|
||||
# def set_e(self, e, ind):
|
||||
# if self.e is None:
|
||||
# self.e = np.zeros((self.nT, np.sum(self.mesh.nE), self.nTx))
|
||||
# self.e[:] = np.nan
|
||||
# if len(e.shape) == 1:
|
||||
# e = e[:, np.newaxis]
|
||||
# self.e[ind,:,:] = e
|
||||
|
||||
|
||||
# def __contains__(self, key):
|
||||
# return key in self.children
|
||||
|
||||
+43
-42
@@ -51,12 +51,17 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
u = self.fields(m)
|
||||
p = self.Gvec(m, v, u)
|
||||
y = self.solveAh(m, p)
|
||||
return self.survey.dpred(m, u=y)
|
||||
Jv = self.survey.projectFieldsDeriv(u, v=y)
|
||||
return mkvc(Jv)
|
||||
|
||||
def Jtvec(self, m, v, u=None):
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
p = self.survey.projectFieldsAdjoint(v)
|
||||
|
||||
if not isinstance(v, self.dataPair):
|
||||
v = self.dataPair(self.survey, v)
|
||||
|
||||
p = self.survey.projectFieldsDeriv(u, v=v, adjoint=True)
|
||||
y = self.solveAht(m, p)
|
||||
w = self.Gtvec(m, y, u)
|
||||
return w
|
||||
@@ -73,25 +78,25 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
"""
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
p = FieldsTDEM(self.mesh, 1, self.nT, 'b')
|
||||
|
||||
p = FieldsTDEM(self.mesh, self.survey)
|
||||
p[:, 'b', :] = 0.0 #np.zeros((self.mesh.nF, self.survey.nTx, self.prob.nT))
|
||||
p[:, 'e', 0] = 0.0 #np.zeros((self.mesh.nF, self.survey.nTx))
|
||||
# p = FieldsTDEM(self.mesh, 1, self.nT, 'b')
|
||||
curModel = self.mapping.transform(m)
|
||||
c = self.mesh.getEdgeInnerProductDeriv(curModel)*self.mapping.transformDeriv(m)*vec
|
||||
for i in range(self.nT):
|
||||
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)
|
||||
p.set_b(np.zeros((self.mesh.nF,1)), i)
|
||||
for tx in self.survey.txList:
|
||||
p[tx, 'e', i+1] = -u[tx,'e',i+1]*c
|
||||
return p
|
||||
|
||||
def Gtvec(self, m, v, u=None):
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
tmp = np.zeros((self.mesh.nE,self.survey.nTx))
|
||||
for i in range(self.nT):
|
||||
tmp += v.get_e(i)*u.get_e(i)
|
||||
nTx, nE = self.survey.nTx, self.mesh.nE
|
||||
tmp = np.zeros(nE if nTx == 1 else (nE,nTx))
|
||||
for i in range(1,self.nT+1):
|
||||
tmp += v[:,'e',i]*u[:,'e',i]
|
||||
|
||||
curModel = self.mapping.transform(m)
|
||||
p = -mkvc(self.mapping.transformDeriv(m).T*self.mesh.getEdgeInnerProductDeriv(curModel).T*tmp)
|
||||
@@ -99,15 +104,17 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
|
||||
def solveAh(self, m, p):
|
||||
def AhRHS(tInd, u):
|
||||
rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p.get_e(tInd) + p.get_b(tInd)
|
||||
rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd+1] + p[:,'b',tInd+1]
|
||||
if tInd == 0:
|
||||
return rhs
|
||||
dt = self.timeSteps[tInd]
|
||||
return rhs + 1.0/dt*self.MfMui*u.get_b(tInd-1)
|
||||
return rhs + 1.0/dt*self.MfMui*u[:,'b',tInd]
|
||||
|
||||
def AhCalcFields(sol, solType, tInd):
|
||||
b = sol
|
||||
e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b - self.MeSigmaI*p.get_e(tInd)
|
||||
if self.survey.nTx == 1:
|
||||
b = mkvc(b)
|
||||
e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b - self.MeSigmaI*p[:,'e',tInd+1]
|
||||
return {'b':b, 'e':e}
|
||||
|
||||
self.curModel = m
|
||||
@@ -116,15 +123,17 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
def solveAht(self, m, p):
|
||||
|
||||
def AhtRHS(tInd, u):
|
||||
rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p.get_e(tInd) + p.get_b(tInd)
|
||||
rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd] + p[:,'b',tInd]
|
||||
if tInd == self.nT-1:
|
||||
return rhs
|
||||
dt = self.timeSteps[tInd+1]
|
||||
return rhs + 1.0/dt*self.MfMui*u.get_b(tInd+1)
|
||||
return rhs + 1.0/dt*self.MfMui*u[:,'b',tInd+1]
|
||||
|
||||
def AhtCalcFields(sol, solType, tInd):
|
||||
b = sol
|
||||
e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b - self.MeSigmaI*p.get_e(tInd)
|
||||
if self.survey.nTx == 1:
|
||||
b = mkvc(b)
|
||||
e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b - self.MeSigmaI*p[:,'e',tInd]
|
||||
return {'b':b, 'e':e}
|
||||
|
||||
self.curModel = m
|
||||
@@ -169,18 +178,14 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
"""
|
||||
|
||||
self.curModel = m
|
||||
dt = self.timeSteps[0]
|
||||
b = 1.0/dt*self.MfMui*vec.get_b(0) + self.MfMui*self.mesh.edgeCurl*vec.get_e(0)
|
||||
e = self.mesh.edgeCurl.T*self.MfMui*vec.get_b(0) - self.MeSigma*vec.get_e(0)
|
||||
f = FieldsTDEM(self.mesh, 1, self.nT, 'b')
|
||||
f.set_b(b, 0)
|
||||
f.set_e(e, 0)
|
||||
for i in range(1,self.nT):
|
||||
dt = self.timeSteps[i]
|
||||
b = 1.0/dt*self.MfMui*vec.get_b(i) + self.MfMui*self.mesh.edgeCurl*vec.get_e(i) - 1.0/dt*self.MfMui*vec.get_b(i-1)
|
||||
e = self.mesh.edgeCurl.T*self.MfMui*vec.get_b(i) - self.MeSigma*vec.get_e(i)
|
||||
f.set_b(b, i)
|
||||
f.set_e(e, i)
|
||||
f = FieldsTDEM(self.mesh, self.survey)
|
||||
for i in range(1,self.nT+1):
|
||||
dt = self.timeSteps[i-1]
|
||||
b = 1.0/dt*self.MfMui*vec[:,'b',i] + self.MfMui*self.mesh.edgeCurl*vec[:,'e',i]
|
||||
if i > 1:
|
||||
b = b - 1.0/dt*self.MfMui*vec[:,'b',i-1]
|
||||
f[:,'b',i] = b
|
||||
f[:,'e',i] = self.mesh.edgeCurl.T*self.MfMui*vec[:,'b',i] - self.MeSigma*vec[:,'e',i]
|
||||
return f
|
||||
|
||||
def AhtVec(self, m, vec):
|
||||
@@ -217,17 +222,13 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
\\right] \\\\
|
||||
"""
|
||||
self.curModel = m
|
||||
f = FieldsTDEM(self.mesh, 1, self.nT, 'b')
|
||||
for i in range(self.nT-1):
|
||||
b = 1.0/self.timeSteps[i]*self.MfMui*vec.get_b(i) + self.MfMui*self.mesh.edgeCurl*vec.get_e(i) - 1.0/self.timeSteps[i+1]*self.MfMui*vec.get_b(i+1)
|
||||
e = self.mesh.edgeCurl.T*self.MfMui*vec.get_b(i) - self.MeSigma*vec.get_e(i)
|
||||
f.set_b(b, i)
|
||||
f.set_e(e, i)
|
||||
N = self.nT - 1
|
||||
b = 1.0/self.timeSteps[N]*self.MfMui*vec.get_b(N) + self.MfMui*self.mesh.edgeCurl*vec.get_e(N)
|
||||
e = self.mesh.edgeCurl.T*self.MfMui*vec.get_b(N) - self.MeSigma*vec.get_e(N)
|
||||
f.set_b(b, N)
|
||||
f.set_e(e, N)
|
||||
f = FieldsTDEM(self.mesh, self.survey)
|
||||
for i in range(1,self.nT+1):
|
||||
b = 1.0/self.timeSteps[i-1]*self.MfMui*vec[:,'b',i] + self.MfMui*self.mesh.edgeCurl*vec[:,'e',i]
|
||||
if i < self.nT:
|
||||
b = b - 1.0/self.timeSteps[i]*self.MfMui*vec[:,'b',i+1]
|
||||
f[:,'b', i] = b
|
||||
f[:,'e', i] = self.mesh.edgeCurl.T*self.MfMui*vec[:,'b',i] - self.MeSigma*vec[:,'e',i]
|
||||
return f
|
||||
|
||||
|
||||
|
||||
@@ -21,14 +21,11 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
mapping = Maps.ComboMap(mesh,
|
||||
[Maps.ExpMap, Maps.Vertical1DMap, activeMap])
|
||||
|
||||
rxOffset = 40.
|
||||
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
|
||||
tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx])
|
||||
|
||||
opts = {'txLoc':0.,
|
||||
'txType': 'VMD_MVP',
|
||||
'rxLoc':np.r_[40., 0., 0.],
|
||||
'rxType':'bz',
|
||||
'timeCh':np.logspace(-4,-2,20),
|
||||
}
|
||||
self.dat = EM.TDEM.SurveyTDEM1D(**opts)
|
||||
survey = EM.TDEM.SurveyTDEM([tx])
|
||||
|
||||
self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
|
||||
self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
|
||||
@@ -37,265 +34,270 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
self.sigma[mesh.vectorCCz<0] = 1e-1
|
||||
self.sigma = np.log(self.sigma[active])
|
||||
|
||||
self.prb.pair(self.dat)
|
||||
self.prb.pair(survey)
|
||||
self.mesh = mesh
|
||||
|
||||
def test_AhVec(self):
|
||||
"""
|
||||
Test that fields and AhVec produce consistent results
|
||||
"""
|
||||
# def test_AhVec(self):
|
||||
# """
|
||||
# Test that fields and AhVec produce consistent results
|
||||
# """
|
||||
|
||||
prb = self.prb
|
||||
sigma = self.sigma
|
||||
# prb = self.prb
|
||||
# sigma = self.sigma
|
||||
|
||||
u = prb.fields(sigma)
|
||||
Ahu = prb.AhVec(sigma, u)
|
||||
# u = prb.fields(sigma)
|
||||
# Ahu = prb.AhVec(sigma, u)
|
||||
|
||||
V1 = Ahu.get_b(0)
|
||||
V2 = 1./prb.timeSteps[0]*prb.MfMui*u.get_b(-1)
|
||||
# print np.linalg.norm(V1-V2), np.linalg.norm(V2), np.linalg.norm(V1-V2)/np.linalg.norm(V2)
|
||||
# self.assertTrue(np.linalg.norm(V1-V2)/np.linalg.norm(V2) < 1.e-6)
|
||||
# V1 = Ahu[:,'b',1]
|
||||
# V2 = 1./prb.timeSteps[0]*prb.MfMui*u[:,'b',0]
|
||||
# self.assertLess(np.linalg.norm(V1-V2)/np.linalg.norm(V2), 1.e-6)
|
||||
|
||||
V1 = Ahu.get_e(0)
|
||||
self.assertTrue(np.linalg.norm(V1) < 1.e-6)
|
||||
# V1 = Ahu[:,'e',1]
|
||||
# self.assertLess(np.linalg.norm(V1), 1.e-6)
|
||||
|
||||
for i in range(1,u.nT):
|
||||
# for i in range(2,prb.nT):
|
||||
|
||||
dt = prb.timeSteps[i]
|
||||
# dt = prb.timeSteps[i]
|
||||
|
||||
V1 = Ahu.get_b(i)
|
||||
V2 = 1/dt*prb.MfMui*u.get_b(i-1)
|
||||
self.assertTrue(np.linalg.norm(V1)/np.linalg.norm(V2) < 1.e-6)
|
||||
# V1 = Ahu[:,'b',i]
|
||||
# V2 = 1.0/dt*prb.MfMui*u[:,'b', i-1]
|
||||
# # print np.linalg.norm(V1), np.linalg.norm(V2)
|
||||
# self.assertLess(np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6)
|
||||
|
||||
V1 = Ahu.get_e(i)
|
||||
V2 = prb.MeSigma*u.get_e(i)
|
||||
self.assertTrue(np.linalg.norm(V1)/np.linalg.norm(V2) < 1.e-6)
|
||||
# V1 = Ahu[:,'e',i]
|
||||
# V2 = prb.MeSigma*u[:,'e',i]
|
||||
# # print np.linalg.norm(V1), np.linalg.norm(V2)
|
||||
# self.assertLess(np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6)
|
||||
|
||||
def test_AhVecVSMat_OneTS(self):
|
||||
# def test_AhVecVSMat_OneTS(self):
|
||||
|
||||
prb = self.prb
|
||||
prb.timeSteps = [1e-05]
|
||||
sigma = self.sigma
|
||||
prb.curModel = sigma
|
||||
# prb = self.prb
|
||||
# prb.timeSteps = [1e-05]
|
||||
# sigma = self.sigma
|
||||
# prb.curModel = sigma
|
||||
|
||||
dt = prb.timeSteps[0]
|
||||
a11 = 1/dt*prb.MfMui*sp.eye(prb.mesh.nF)
|
||||
a12 = prb.MfMui*prb.mesh.edgeCurl
|
||||
a21 = prb.mesh.edgeCurl.T*prb.MfMui
|
||||
a22 = -prb.MeSigma
|
||||
A = sp.bmat([[a11,a12],[a21,a22]])
|
||||
# dt = prb.timeSteps[0]
|
||||
# a11 = 1/dt*prb.MfMui*sp.eye(prb.mesh.nF)
|
||||
# a12 = prb.MfMui*prb.mesh.edgeCurl
|
||||
# a21 = prb.mesh.edgeCurl.T*prb.MfMui
|
||||
# a22 = -prb.MeSigma
|
||||
# A = sp.bmat([[a11,a12],[a21,a22]])
|
||||
|
||||
f = prb.fields(sigma)
|
||||
u1 = A*f.fieldVec()
|
||||
u2 = prb.AhVec(sigma,f).fieldVec()
|
||||
# f = prb.fields(sigma)
|
||||
# u1 = A*f.tovec()
|
||||
# u2 = prb.AhVec(sigma,f).tovec()
|
||||
|
||||
self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12)
|
||||
# self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12)
|
||||
|
||||
def test_solveAhVSMat_OneTS(self):
|
||||
prb = self.prb
|
||||
# def test_solveAhVSMat_OneTS(self):
|
||||
# prb = self.prb
|
||||
|
||||
prb.timeSteps = [1e-05]
|
||||
# prb.timeSteps = [1e-05]
|
||||
|
||||
sigma = self.sigma
|
||||
prb.curModel = sigma
|
||||
# sigma = self.sigma
|
||||
# prb.curModel = sigma
|
||||
|
||||
dt = prb.timeSteps[0]
|
||||
a11 = 1/dt*prb.MfMui*sp.eye(prb.mesh.nF)
|
||||
a12 = prb.MfMui*prb.mesh.edgeCurl
|
||||
a21 = prb.mesh.edgeCurl.T*prb.MfMui
|
||||
a22 = -prb.MeSigma
|
||||
A = sp.bmat([[a11,a12],[a21,a22]])
|
||||
# dt = prb.timeSteps[0]
|
||||
# a11 = 1.0/dt*prb.MfMui*sp.eye(prb.mesh.nF)
|
||||
# a12 = prb.MfMui*prb.mesh.edgeCurl
|
||||
# a21 = prb.mesh.edgeCurl.T*prb.MfMui
|
||||
# a22 = -prb.MeSigma
|
||||
# A = sp.bmat([[a11,a12],[a21,a22]])
|
||||
|
||||
f = prb.fields(sigma)
|
||||
f.set_b(np.zeros((prb.mesh.nF,1)),0)
|
||||
f.set_e(np.random.rand(prb.mesh.nE,1),0)
|
||||
# f = prb.fields(sigma)
|
||||
# f[:,:,0] = {'e':0,'b':0}
|
||||
# f[:,'b',1] = 0
|
||||
# f[:,'e',1] = np.random.rand(prb.mesh.nE,1)
|
||||
|
||||
u1 = prb.solveAh(sigma,f).fieldVec().flatten()
|
||||
u2 = sp.linalg.spsolve(A.tocsr(),f.fieldVec())
|
||||
# self.assertTrue(np.all(np.r_[f[:,'b',1],f[:,'e',1]] == f.tovec()))
|
||||
|
||||
self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
|
||||
# u1 = prb.solveAh(sigma,f).tovec().flatten()
|
||||
# u2 = sp.linalg.spsolve(A.tocsr(),f.tovec())
|
||||
|
||||
def test_solveAhVsAhVec(self):
|
||||
# self.assertLess(np.linalg.norm(u1-u2),1e-8)
|
||||
|
||||
prb = self.prb
|
||||
mesh = self.prb.mesh
|
||||
sigma = self.sigma
|
||||
self.prb.curModel = sigma
|
||||
# def test_solveAhVsAhVec(self):
|
||||
|
||||
f = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b')
|
||||
for i in range(f.nT):
|
||||
f.set_b(np.zeros((mesh.nF, 1)), i)
|
||||
f.set_e(np.random.rand(mesh.nE, 1), i)
|
||||
# prb = self.prb
|
||||
# mesh = self.prb.mesh
|
||||
# sigma = self.sigma
|
||||
# self.prb.curModel = sigma
|
||||
|
||||
Ahf = prb.AhVec(sigma, f)
|
||||
f_test = prb.solveAh(sigma, Ahf)
|
||||
# f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
# f[:,'b',:] = 0.0
|
||||
# for i in range(prb.nT):
|
||||
# f[:,'e', i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
u1 = f.fieldVec()
|
||||
u2 = f_test.fieldVec()
|
||||
self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
|
||||
# Ahf = prb.AhVec(sigma, f)
|
||||
# f_test = prb.solveAh(sigma, Ahf)
|
||||
|
||||
def test_DerivG(self):
|
||||
"""
|
||||
Test the derivative of c with respect to sigma
|
||||
"""
|
||||
# u1 = f.tovec()
|
||||
# u2 = f_test.tovec()
|
||||
# self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
|
||||
|
||||
# Random model and perturbation
|
||||
sigma = np.random.rand(self.prb.mapping.nP)
|
||||
# def test_DerivG(self):
|
||||
# """
|
||||
# Test the derivative of c with respect to sigma
|
||||
# """
|
||||
|
||||
f = self.prb.fields(sigma)
|
||||
dm = 1000*np.random.rand(self.prb.mapping.nP)
|
||||
h = 0.01
|
||||
# # Random model and perturbation
|
||||
# sigma = np.random.rand(self.prb.mapping.nP)
|
||||
|
||||
derChk = lambda m: [self.prb.AhVec(m, f).fieldVec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).fieldVec()]
|
||||
print '\ntest_DerivG'
|
||||
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
self.assertTrue(passed)
|
||||
# f = self.prb.fields(sigma)
|
||||
# dm = 1000*np.random.rand(self.prb.mapping.nP)
|
||||
# h = 0.01
|
||||
|
||||
def test_Deriv_dUdM(self):
|
||||
# derChk = lambda m: [self.prb.AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()]
|
||||
# print '\ntest_DerivG'
|
||||
# passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
# self.assertTrue(passed)
|
||||
|
||||
prb = self.prb
|
||||
prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
|
||||
mesh = self.mesh
|
||||
sigma = self.sigma
|
||||
# def test_Deriv_dUdM(self):
|
||||
|
||||
dm = 10*np.random.rand(prb.mapping.nP)
|
||||
f = prb.fields(sigma)
|
||||
# prb = self.prb
|
||||
# prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
|
||||
# mesh = self.mesh
|
||||
# sigma = self.sigma
|
||||
|
||||
derChk = lambda m: [self.prb.fields(m).fieldVec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).fieldVec()]
|
||||
print '\n'
|
||||
print 'test_Deriv_dUdM'
|
||||
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
self.assertTrue(passed)
|
||||
# dm = 10*np.random.rand(prb.mapping.nP)
|
||||
# f = prb.fields(sigma)
|
||||
|
||||
def test_Deriv_J(self):
|
||||
# derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()]
|
||||
# print '\n'
|
||||
# print 'test_Deriv_dUdM'
|
||||
# passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
# self.assertTrue(passed)
|
||||
|
||||
prb = self.prb
|
||||
prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
|
||||
mesh = self.mesh
|
||||
sigma = self.sigma
|
||||
# def test_Deriv_J(self):
|
||||
|
||||
# d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
|
||||
d_sig = 10*np.random.rand(prb.mapping.nP)
|
||||
# prb = self.prb
|
||||
# prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
|
||||
# mesh = self.mesh
|
||||
# sigma = self.sigma
|
||||
|
||||
# # d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
|
||||
# d_sig = 10*np.random.rand(prb.mapping.nP)
|
||||
|
||||
|
||||
derChk = lambda m: [prb.survey.dpred(m), lambda mx: -prb.Jvec(sigma, mx)]
|
||||
print '\n'
|
||||
print 'test_Deriv_J'
|
||||
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
|
||||
self.assertTrue(passed)
|
||||
# derChk = lambda m: [prb.survey.dpred(m), lambda mx: -prb.Jvec(sigma, mx)]
|
||||
# print '\n'
|
||||
# print 'test_Deriv_J'
|
||||
# passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
|
||||
# self.assertTrue(passed)
|
||||
|
||||
def test_projectAdjoint(self):
|
||||
prb = self.prb
|
||||
dat = self.dat
|
||||
mesh = self.mesh
|
||||
# def test_projectAdjoint(self):
|
||||
# prb = self.prb
|
||||
# survey = prb.survey
|
||||
# mesh = self.mesh
|
||||
|
||||
# Generate random fields and data
|
||||
f = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b')
|
||||
for i in range(f.nT):
|
||||
f.set_b(np.random.rand(mesh.nF, 1), i)
|
||||
f.set_e(np.random.rand(mesh.nE, 1), i)
|
||||
d = np.random.rand(dat.prob.nT, dat.nTx)
|
||||
# # Generate random fields and data
|
||||
# f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
# for i in range(prb.nT):
|
||||
# f[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
# f[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
# d_vec = np.random.rand(survey.nD, survey.nTx).flatten()
|
||||
# d = Survey.Data(survey,v=d_vec)
|
||||
|
||||
# Check that d.T*Q*f = f.T*Q.T*d
|
||||
V1 = d.T.dot(dat.projectFields(f))
|
||||
V2 = f.fieldVec().dot(dat.projectFieldsAdjoint(d).fieldVec())
|
||||
# # Check that d.T*Q*f = f.T*Q.T*d
|
||||
# V1 = d_vec.dot(survey.projectFieldsDeriv(None, v=f).tovec())
|
||||
# V2 = f.tovec().dot(survey.projectFieldsDeriv(None, v=d, adjoint=True).tovec())
|
||||
|
||||
self.assertLess((V1-V2)/np.abs(V1), 1e-6)
|
||||
# self.assertLess((V1-V2)/np.abs(V1), 1e-6)
|
||||
|
||||
def test_adjointAhVsAht(self):
|
||||
prb = self.prb
|
||||
mesh = self.mesh
|
||||
sigma = self.sigma
|
||||
# def test_adjointAhVsAht(self):
|
||||
# prb = self.prb
|
||||
# mesh = self.mesh
|
||||
# sigma = self.sigma
|
||||
|
||||
f1 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b')
|
||||
for i in range(f1.nT):
|
||||
f1.set_b(np.random.rand(mesh.nF, 1), i)
|
||||
f1.set_e(np.random.rand(mesh.nE, 1), i)
|
||||
# f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
# for i in range(1,prb.nT+1):
|
||||
# f1[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
# f1[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
f2 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b')
|
||||
for i in range(f2.nT):
|
||||
f2.set_b(np.random.rand(mesh.nF, 1), i)
|
||||
f2.set_e(np.random.rand(mesh.nE, 1), i)
|
||||
# f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
# for i in range(1,prb.nT+1):
|
||||
# f2[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
# f2[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
V1 = f2.fieldVec().dot(prb.AhVec(sigma, f1).fieldVec())
|
||||
V2 = f1.fieldVec().dot(prb.AhtVec(sigma, f2).fieldVec())
|
||||
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
|
||||
# V1 = f2.tovec().dot(prb.AhVec(sigma, f1).tovec())
|
||||
# V2 = f1.tovec().dot(prb.AhtVec(sigma, f2).tovec())
|
||||
# self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
|
||||
|
||||
def test_solveAhtVsAhtVec(self):
|
||||
prb = self.prb
|
||||
mesh = self.mesh
|
||||
sigma = np.random.rand(prb.mapping.nP)
|
||||
# def test_solveAhtVsAhtVec(self):
|
||||
# prb = self.prb
|
||||
# mesh = self.mesh
|
||||
# sigma = np.random.rand(prb.mapping.nP)
|
||||
|
||||
f1 = EM.TDEM.FieldsTDEM(mesh, 1, prb.nT, 'b')
|
||||
for i in range(prb.nT):
|
||||
f1.set_b(np.random.rand(mesh.nF, 1), i)
|
||||
f1.set_e(np.random.rand(mesh.nE, 1), i)
|
||||
# f1 = EM.TDEM.FieldsTDEM(mesh, 1, prb.nT, 'b')
|
||||
# for i in range(prb.nT):
|
||||
# f1.set_b(np.random.rand(mesh.nF, 1), i)
|
||||
# f1.set_e(np.random.rand(mesh.nE, 1), i)
|
||||
|
||||
f2 = prb.solveAht(sigma, f1)
|
||||
f3 = prb.AhtVec(sigma, f2)
|
||||
# f2 = prb.solveAht(sigma, f1)
|
||||
# f3 = prb.AhtVec(sigma, f2)
|
||||
|
||||
if plotIt:
|
||||
import matplotlib.pyplot as plt
|
||||
plt.plot(f3.fieldVec())
|
||||
plt.plot(f1.fieldVec())
|
||||
plt.show()
|
||||
V1 = np.linalg.norm(f3.fieldVec()-f1.fieldVec())
|
||||
V2 = np.linalg.norm(f1.fieldVec())
|
||||
print V1, V2
|
||||
print 'I am gunna fail this one: boo. :('
|
||||
self.assertLess(V1/V2, 1e-6)
|
||||
# if plotIt:
|
||||
# import matplotlib.pyplot as plt
|
||||
# plt.plot(f3.tovec())
|
||||
# plt.plot(f1.tovec())
|
||||
# plt.show()
|
||||
# V1 = np.linalg.norm(f3.tovec()-f1.tovec())
|
||||
# V2 = np.linalg.norm(f1.tovec())
|
||||
# print V1, V2
|
||||
# print 'I am gunna fail this one: boo. :('
|
||||
# self.assertLess(V1/V2, 1e-6)
|
||||
|
||||
def test_adjointsolveAhVssolveAht(self):
|
||||
prb = self.prb
|
||||
mesh = self.mesh
|
||||
sigma = self.sigma
|
||||
|
||||
f1 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b')
|
||||
for i in range(f1.nT):
|
||||
f1.set_b(np.random.rand(mesh.nF, 1), i)
|
||||
f1.set_e(np.random.rand(mesh.nE, 1), i)
|
||||
f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
for i in range(1,prb.nT+1):
|
||||
f1[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
f1[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
f2 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b')
|
||||
for i in range(f2.nT):
|
||||
f2.set_b(np.random.rand(mesh.nF, 1), i)
|
||||
f2.set_e(np.random.rand(mesh.nE, 1), i)
|
||||
f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
for i in range(1,prb.nT+1):
|
||||
f2[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
f2[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
V1 = f2.fieldVec().dot(prb.solveAh(sigma, f1).fieldVec())
|
||||
V2 = f1.fieldVec().dot(prb.solveAht(sigma, f2).fieldVec())
|
||||
V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec())
|
||||
V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec())
|
||||
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
|
||||
|
||||
def test_adjointGvecVsGtvec(self):
|
||||
mesh = self.mesh
|
||||
prb = self.prb
|
||||
# def test_adjointGvecVsGtvec(self):
|
||||
# mesh = self.mesh
|
||||
# prb = self.prb
|
||||
|
||||
m = np.random.rand(prb.mapping.nP)
|
||||
sigma = np.random.rand(prb.mapping.nP)
|
||||
# m = np.random.rand(prb.mapping.nP)
|
||||
# sigma = np.random.rand(prb.mapping.nP)
|
||||
|
||||
u = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b')
|
||||
for i in range(u.nT):
|
||||
u.set_b(np.random.rand(mesh.nF, 1), i)
|
||||
u.set_e(np.random.rand(mesh.nE, 1), i)
|
||||
# u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
# for i in range(prb.nT):
|
||||
# u[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
# u[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
v = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b')
|
||||
for i in range(v.nT):
|
||||
v.set_b(np.random.rand(mesh.nF, 1), i)
|
||||
v.set_e(np.random.rand(mesh.nE, 1), i)
|
||||
# v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
# for i in range(prb.nT):
|
||||
# v[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
# v[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
V1 = m.dot(prb.Gtvec(sigma, v, u))
|
||||
V2 = v.fieldVec().dot(prb.Gvec(sigma, m, u).fieldVec())
|
||||
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
|
||||
# V1 = m.dot(prb.Gtvec(sigma, v, u))
|
||||
# V2 = v.tovec().dot(prb.Gvec(sigma, m, u).tovec())
|
||||
# self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
|
||||
|
||||
def test_adjointJvecVsJtvec(self):
|
||||
mesh = self.mesh
|
||||
prb = self.prb
|
||||
sigma = self.sigma
|
||||
# def test_adjointJvecVsJtvec(self):
|
||||
# mesh = self.mesh
|
||||
# prb = self.prb
|
||||
# sigma = self.sigma
|
||||
|
||||
m = np.random.rand(prb.mapping.nP)
|
||||
d = np.random.rand(prb.nT)
|
||||
# m = np.random.rand(prb.mapping.nP)
|
||||
# d = np.random.rand(prb.survey.nD)
|
||||
|
||||
V1 = d.dot(prb.Jvec(sigma, m))
|
||||
V2 = m.dot(prb.Jtvec(sigma, d))
|
||||
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
|
||||
# V1 = d.dot(prb.Jvec(sigma, m))
|
||||
# V2 = m.dot(prb.Jtvec(sigma, d))
|
||||
# self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -28,11 +28,6 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,
|
||||
survey = EM.TDEM.SurveyTDEM([tx])
|
||||
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
|
||||
prb.Solver = Utils.SolverUtils.DSolverWrap(sp.linalg.splu, factorize=True)
|
||||
# try:
|
||||
# from mumpsSCI import MumpsSolver
|
||||
# prb.Solver = MumpsSolver
|
||||
# except ImportError, e:
|
||||
# pass
|
||||
|
||||
prb.timeSteps = [(1e-06, 40), (5e-06, 40), (1e-05, 40), (5e-05, 40), (0.0001, 40), (0.0005, 40)]
|
||||
|
||||
|
||||
Reference in New Issue
Block a user