mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-28 18:45:25 +08:00
Merge pull request #6 from simpeg/MultipleTransmitters
Multiple transmitters
This commit is contained in:
Rowan Cockett
@@ -255,25 +255,8 @@ Multiplying **J** onto a vector can be broken into three steps
|
||||
\vec{p}_e^{(n)} = - \diag{\e^{(n)}} \Ace \diag{V} m
|
||||
\end{align}
|
||||
|
||||
First time step
|
||||
|
||||
.. math::
|
||||
|
||||
\begin{align}
|
||||
\frac{1}{\delta t} \MfMui \vec{y}_{b}^{(1)} + \MfMui \dcurl \vec{y}_{e}^{(1)} = \vec{p}_b^{(1)} \\
|
||||
\dcurl^\top \MfMui \vec{y}_b^{(1)} - \MeSig \vec{y}_e^{(1)} = \vec{p}_e^{(1)}
|
||||
\end{align}
|
||||
|
||||
|
||||
.. math::
|
||||
|
||||
\begin{align}
|
||||
\left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(1)} = \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(1)} + \vec{p}_b^{(1)} \\
|
||||
\vec{y}_e^{(1)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(1)} - \MeSig^{-1} \vec{p}_e^{(1)}
|
||||
\end{align}
|
||||
|
||||
|
||||
Remaining time steps:
|
||||
For all time steps:
|
||||
|
||||
.. math::
|
||||
|
||||
@@ -295,6 +278,11 @@ and
|
||||
\vec{y}_e^{(t+1)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t+1)} - \MeSig^{-1} \vec{p}_e^{(t+1)}
|
||||
\end{align}
|
||||
|
||||
.. note::
|
||||
|
||||
For the first time step, \\\(t=0\\\), the term: \\\(\\frac{1}{\\delta t} \\MfMui \\vec{y}_b^{(0)}\\\) is zero.
|
||||
|
||||
|
||||
|
||||
|
||||
Implementing **J** transpose times a vector
|
||||
@@ -319,38 +307,21 @@ Multiplying \\(\\mathbf{J}^\\top\\) onto a vector can be broken into three steps
|
||||
\end{array}
|
||||
\right]
|
||||
|
||||
For the last time-step \\(t=N\\):
|
||||
For the all time-steps (going backwards in time):
|
||||
|
||||
|
||||
.. math::
|
||||
|
||||
\begin{align}
|
||||
\frac{1}{\delta t} \MfMui \vec{y}_{b}^{(N)} + \MfMui \dcurl \vec{y}_{e}^{(N)} = \vec{p}_b^{(N)} \\
|
||||
\dcurl^\top \MfMui \vec{y}_b^{(N)} - \MeSig \vec{y}_e^{(N)} = \vec{p}_e^{(N)}
|
||||
\end{align}
|
||||
A \vec{y}^{(t)} + B \vec{y}^{(t+1)} = \vec{p}^{(t)}
|
||||
|
||||
|
||||
.. math::
|
||||
|
||||
\begin{align}
|
||||
\left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(N)} = \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(N)} + \vec{p}_b^{(N)} \\
|
||||
\vec{y}_e^{(N)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(N)} - \MeSig^{-1} \vec{p}_e^{(N)}
|
||||
\end{align}
|
||||
|
||||
For the rest of the time-steps (going backwards in time)
|
||||
|
||||
|
||||
.. math::
|
||||
|
||||
A \vec{y}^{(t-1)} + B \vec{y}^{(t)} = \vec{p}^{(t-1)}
|
||||
|
||||
|
||||
.. math::
|
||||
|
||||
\begin{align}
|
||||
\frac{1}{\delta t} \MfMui\vec{y}_{b}^{(t-1)} + \MfMui\dcurl \vec{y}_{e}^{(t-1)}
|
||||
- \frac{1}{\delta t} \MfMui \vec{y}_{b}^{(t)}
|
||||
= \vec{p}_b^{(t-1)} \\
|
||||
\dcurl^\top \MfMui \vec{y}_b^{(t-1)} - \MeSig \vec{y}_e^{(t-1)} = \vec{p}_e^{(t-1)}
|
||||
\frac{1}{\delta t} \MfMui\vec{y}_{b}^{(t)} + \MfMui\dcurl \vec{y}_{e}^{(t)}
|
||||
- \frac{1}{\delta t} \MfMui \vec{y}_{b}^{(t+1)}
|
||||
= \vec{p}_b^{(t)} \\
|
||||
\dcurl^\top \MfMui \vec{y}_b^{(t)} - \MeSig \vec{y}_e^{(t)} = \vec{p}_e^{(t)}
|
||||
\end{align}
|
||||
|
||||
and
|
||||
@@ -358,8 +329,13 @@ and
|
||||
.. math::
|
||||
|
||||
\begin{align}
|
||||
\left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(t-1)} =
|
||||
\frac{1}{\delta t} \MfMui \vec{y}_b^{(t)}
|
||||
+ \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(t-1)} + \vec{p}_b^{(t-1)} \\
|
||||
\vec{y}_e^{(t-1)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t-1)} - \MeSig^{-1} \vec{p}_e^{(t-1)}
|
||||
\left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(t)} =
|
||||
\frac{1}{\delta t} \MfMui \vec{y}_b^{(t+1)}
|
||||
+ \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(t)} + \vec{p}_b^{(t)} \\
|
||||
\vec{y}_e^{(t)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t)} - \MeSig^{-1} \vec{p}_e^{(t)}
|
||||
\end{align}
|
||||
|
||||
|
||||
.. note::
|
||||
|
||||
For the last time step, \\\(t=N\\\), the term: \\\(\\frac{1}{\\delta t} \\MfMui \\vec{y}_b^{(N+1)}\\\) is zero.
|
||||
|
||||
@@ -8,7 +8,7 @@ def omega(freq):
|
||||
"""Change frequency to angular frequency, omega"""
|
||||
return 2.*np.pi*freq
|
||||
|
||||
class BaseProblemFDEM(BaseEMProblem):
|
||||
class BaseFDEMProblem(BaseEMProblem):
|
||||
"""
|
||||
We start by looking at Maxwell's equations in the electric field \\(\\vec{E}\\) and the magnetic flux density \\(\\vec{B}\\):
|
||||
|
||||
@@ -106,7 +106,7 @@ class BaseProblemFDEM(BaseEMProblem):
|
||||
return Jtv
|
||||
|
||||
|
||||
class ProblemFDEM_e(BaseProblemFDEM):
|
||||
class ProblemFDEM_e(BaseFDEMProblem):
|
||||
"""
|
||||
By eliminating the magnetic flux density using
|
||||
|
||||
@@ -127,7 +127,7 @@ class ProblemFDEM_e(BaseProblemFDEM):
|
||||
solType = 'e'
|
||||
|
||||
def __init__(self, model, **kwargs):
|
||||
BaseProblemFDEM.__init__(self, model, **kwargs)
|
||||
BaseFDEMProblem.__init__(self, model, **kwargs)
|
||||
|
||||
def getA(self, freq):
|
||||
"""
|
||||
@@ -197,14 +197,14 @@ class ProblemFDEM_e(BaseProblemFDEM):
|
||||
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
|
||||
class ProblemFDEM_b(BaseProblemFDEM):
|
||||
class ProblemFDEM_b(BaseFDEMProblem):
|
||||
"""
|
||||
Solving for b!
|
||||
"""
|
||||
solType = 'b'
|
||||
|
||||
def __init__(self, model, **kwargs):
|
||||
BaseProblemFDEM.__init__(self, model, **kwargs)
|
||||
BaseFDEMProblem.__init__(self, model, **kwargs)
|
||||
|
||||
def getA(self, freq):
|
||||
"""
|
||||
|
||||
+14
-54
@@ -1,7 +1,7 @@
|
||||
from SimPEG import Solver
|
||||
from SimPEG.Problem import BaseTimeProblem
|
||||
from simpegEM.Utils import Sources
|
||||
from SurveyTDEM import FieldsTDEM
|
||||
from SurveyTDEM import FieldsTDEM, SurveyTDEM
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.Utils import sdiag, mkvc
|
||||
from SimPEG import Utils, Mesh
|
||||
@@ -9,49 +9,12 @@ from simpegEM.Base import BaseEMProblem
|
||||
import numpy as np
|
||||
|
||||
|
||||
class MixinInitialFieldCalc(object):
|
||||
"""docstring for MixinInitialFieldCalc"""
|
||||
|
||||
storeTheseFields = 'b'
|
||||
|
||||
def getInitialFields(self):
|
||||
if self.survey.txType == 'VMD_MVP':
|
||||
# Vertical magnetic dipole, magnetic vector potential
|
||||
F = self._getInitialFields_VMD_MVP()
|
||||
else:
|
||||
exStr = 'Invalid txType: ' + str(self.survey.txType)
|
||||
raise Exception(exStr)
|
||||
return F
|
||||
|
||||
def _getInitialFields_VMD_MVP(self):
|
||||
if self.mesh._meshType is 'CYL':
|
||||
if self.mesh.isSymmetric:
|
||||
MVP = Sources.MagneticDipoleVectorPotential(self.survey.txLoc, self.mesh.gridEy, 'y')
|
||||
# MVP = Sources.MagneticDipoleVectorPotential(self.survey.txLoc, np.c_[np.zeros(self.mesh.nN), self.mesh.gridN], 'x')
|
||||
else:
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
elif self.mesh._meshType is 'TENSOR':
|
||||
MVPx = Sources.MagneticDipoleVectorPotential(self.survey.txLoc, self.mesh.gridEx, 'x')
|
||||
MVPy = Sources.MagneticDipoleVectorPotential(self.survey.txLoc, self.mesh.gridEy, 'y')
|
||||
MVPz = Sources.MagneticDipoleVectorPotential(self.survey.txLoc, self.mesh.gridEz, 'z')
|
||||
MVP = np.concatenate((MVPx, MVPy, MVPz))
|
||||
else:
|
||||
raise Exception('Unknown mesh for VMD')
|
||||
|
||||
# Initialize field object
|
||||
F = FieldsTDEM(self.mesh, 1, self.nT, store=self.storeTheseFields)
|
||||
|
||||
# Set initial B
|
||||
F.b0 = self.mesh.edgeCurl*MVP
|
||||
|
||||
return F
|
||||
|
||||
|
||||
class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem, BaseEMProblem):
|
||||
class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
"""docstring for ProblemTDEM1D"""
|
||||
def __init__(self, mesh, mapping=None, **kwargs):
|
||||
BaseTimeProblem.__init__(self, mesh, mapping=mapping, **kwargs)
|
||||
|
||||
surveyPair = SurveyTDEM
|
||||
|
||||
def calcFields(self, sol, solType, tInd):
|
||||
|
||||
@@ -65,18 +28,18 @@ class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem, BaseEMProblem):
|
||||
|
||||
return {'b':b, 'e':e}
|
||||
|
||||
Solver = Solver
|
||||
solveOpts = {}
|
||||
|
||||
def fields(self, m):
|
||||
self.curModel = m
|
||||
F = self.getInitialFields()
|
||||
# Create a fields storage object
|
||||
F = FieldsTDEM(self.mesh, self.survey)
|
||||
for tx in self.survey.txList:
|
||||
# Set the initial conditions
|
||||
F[tx,:,0] = tx.getInitialFields(self.mesh)
|
||||
return self.forward(m, self.getRHS, self.calcFields, F=F)
|
||||
|
||||
|
||||
def forward(self, m, RHS, CalcFields, F=None):
|
||||
if F is None:
|
||||
F = FieldsTDEM(self.mesh, self.survey.nTx, self.nT, store=self.storeTheseFields)
|
||||
F = F or FieldsTDEM(self.mesh, self.survey)
|
||||
|
||||
dtFact = None
|
||||
for tInd, dt in enumerate(self.timeSteps):
|
||||
@@ -84,19 +47,17 @@ class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem, BaseEMProblem):
|
||||
dtFact = dt
|
||||
A = self.getA(tInd)
|
||||
# print 'Factoring... (dt = ' + str(dt) + ')'
|
||||
Asolve = self.Solver(A, **self.solveOpts)
|
||||
Asolve = self.Solver(A, **self.solverOpts)
|
||||
# print 'Done'
|
||||
rhs = RHS(tInd, F)
|
||||
sol = Asolve.solve(rhs)
|
||||
if sol.ndim == 1:
|
||||
sol.shape = (sol.size,1)
|
||||
newFields = CalcFields(sol, self.solType, tInd)
|
||||
F.update(newFields, tInd)
|
||||
F[:,:,tInd+1] = CalcFields(sol, self.solType, tInd)
|
||||
return F
|
||||
|
||||
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 = F or FieldsTDEM(self.mesh, self.survey)
|
||||
|
||||
dtFact = None
|
||||
for tInd, dt in reversed(list(enumerate(self.timeSteps))):
|
||||
@@ -104,13 +65,12 @@ class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem, BaseEMProblem):
|
||||
dtFact = dt
|
||||
A = self.getA(tInd)
|
||||
# print 'Factoring... (dt = ' + str(dt) + ')'
|
||||
Asolve = Solver(A, options=self.solveOpts)
|
||||
Asolve = self.Solver(A, **self.solverOpts)
|
||||
# print 'Done'
|
||||
rhs = RHS(tInd, F)
|
||||
sol = Asolve.solve(rhs)
|
||||
if sol.ndim == 1:
|
||||
sol.shape = (sol.size,1)
|
||||
newFields = CalcFields(sol, self.solType, tInd)
|
||||
F.update(newFields, tInd)
|
||||
F[:,:,tInd+1] = CalcFields(sol, self.solType, tInd)
|
||||
return F
|
||||
|
||||
|
||||
+244
-115
@@ -1,128 +1,257 @@
|
||||
from SimPEG import Utils, np
|
||||
from SimPEG import Utils, Survey, np
|
||||
from SimPEG.Survey import BaseSurvey
|
||||
from simpegEM.Utils import Sources
|
||||
|
||||
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
|
||||
class RxTDEM(Survey.BaseTimeRx):
|
||||
|
||||
knownRxTypes = {
|
||||
'ex':['e', 'Ex'],
|
||||
'ey':['e', 'Ey'],
|
||||
'ez':['e', 'Ez'],
|
||||
|
||||
'bx':['b', 'Fx'],
|
||||
'by':['b', 'Fy'],
|
||||
'bz':['b', 'Fz'],
|
||||
}
|
||||
|
||||
def __init__(self, locs, times, rxType):
|
||||
Survey.BaseTimeRx.__init__(self, locs, times, rxType)
|
||||
|
||||
@property
|
||||
def nTimeCh(self):
|
||||
"""Number of time channels"""
|
||||
return self.timeCh.size
|
||||
def projField(self):
|
||||
"""Field Type projection (e.g. e b ...)"""
|
||||
return self.knownRxTypes[self.rxType][0]
|
||||
|
||||
def __init__(self, **kwargs):
|
||||
BaseSurvey.__init__(self, **kwargs)
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
@property
|
||||
def projGLoc(self):
|
||||
"""Grid Location projection (e.g. Ex Fy ...)"""
|
||||
return self.knownRxTypes[self.rxType][1]
|
||||
|
||||
def projectFields(self, tx, mesh, timeMesh, u):
|
||||
P = self.getP(mesh, timeMesh)
|
||||
u_part = Utils.mkvc(u[tx, self.projField, :])
|
||||
return P*u_part
|
||||
|
||||
def projectFieldsDeriv(self, tx, mesh, timeMesh, u, v, adjoint=False):
|
||||
P = self.getP(mesh, timeMesh)
|
||||
|
||||
if not adjoint:
|
||||
return P * Utils.mkvc(v[tx, self.projField, :])
|
||||
elif adjoint:
|
||||
return P.T * v[tx, self]
|
||||
|
||||
|
||||
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.concatenate((u, b, e))
|
||||
return Utils.mkvc(u)
|
||||
|
||||
class TxTDEM(Survey.BaseTx):
|
||||
rxPair = RxTDEM
|
||||
knownTxTypes = ['VMD_MVP']
|
||||
|
||||
def getInitialFields(self, mesh):
|
||||
F0 = getattr(self, '_getInitialFields_' + self.txType)(mesh)
|
||||
return F0
|
||||
|
||||
def _getInitialFields_VMD_MVP(self, mesh):
|
||||
"""Vertical magnetic dipole, magnetic vector potential"""
|
||||
if mesh._meshType is 'CYL':
|
||||
if mesh.isSymmetric:
|
||||
MVP = Sources.MagneticDipoleVectorPotential(self.loc, mesh.gridEy, 'y')
|
||||
else:
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
elif mesh._meshType is 'TENSOR':
|
||||
MVPx = Sources.MagneticDipoleVectorPotential(self.loc, mesh.gridEx, 'x')
|
||||
MVPy = Sources.MagneticDipoleVectorPotential(self.loc, mesh.gridEy, 'y')
|
||||
MVPz = Sources.MagneticDipoleVectorPotential(self.loc, mesh.gridEz, 'z')
|
||||
MVP = np.concatenate((MVPx, MVPy, MVPz))
|
||||
else:
|
||||
raise Exception('Unknown mesh for VMD')
|
||||
|
||||
return {"b": mesh.edgeCurl*MVP}
|
||||
|
||||
def getJs(self, time):
|
||||
return None
|
||||
|
||||
class SurveyTDEM(Survey.BaseSurvey):
|
||||
"""
|
||||
docstring for SurveyTDEM
|
||||
"""
|
||||
|
||||
txPair = TxTDEM
|
||||
|
||||
def __init__(self, txList, **kwargs):
|
||||
# Sort these by frequency
|
||||
self.txList = txList
|
||||
Survey.BaseSurvey.__init__(self, **kwargs)
|
||||
|
||||
def projectFields(self, u):
|
||||
#TODO: this is hardcoded to 1Tx
|
||||
return self.Qrx.dot(u.b[:,:,0].T).T
|
||||
data = Survey.Data(self)
|
||||
for tx in self.txList:
|
||||
for rx in tx.rxList:
|
||||
data[tx, rx] = rx.projectFields(tx, self.mesh, self.prob.timeMesh, u)
|
||||
return data
|
||||
|
||||
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
|
||||
def projectFieldsDeriv(self, u, v=None, adjoint=False):
|
||||
assert v is not None, 'v to multiply must be provided.'
|
||||
|
||||
####################################################
|
||||
# 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(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
|
||||
|
||||
+102
-59
@@ -1,9 +1,9 @@
|
||||
from BaseTDEM import ProblemBaseTDEM
|
||||
from BaseTDEM import BaseTDEMProblem
|
||||
from SimPEG.Utils import mkvc
|
||||
import numpy as np
|
||||
from SurveyTDEM import SurveyTDEM1D, FieldsTDEM
|
||||
from SurveyTDEM import SurveyTDEM, FieldsTDEM
|
||||
|
||||
class ProblemTDEM_b(ProblemBaseTDEM):
|
||||
class ProblemTDEM_b(BaseTDEMProblem):
|
||||
"""
|
||||
Time-Domain EM problem - B-formulation
|
||||
|
||||
@@ -16,11 +16,11 @@ class ProblemTDEM_b(ProblemBaseTDEM):
|
||||
with \\\(\\b\\\) defined on cell faces and \\\(\e\\\) defined on edges.
|
||||
"""
|
||||
def __init__(self, mesh, mapping=None, **kwargs):
|
||||
ProblemBaseTDEM.__init__(self, mesh, mapping=mapping, **kwargs)
|
||||
BaseTDEMProblem.__init__(self, mesh, mapping=mapping, **kwargs)
|
||||
|
||||
solType = 'b'
|
||||
|
||||
surveyPair = SurveyTDEM1D
|
||||
surveyPair = SurveyTDEM
|
||||
|
||||
####################################################
|
||||
# Internal Methods
|
||||
@@ -32,13 +32,14 @@ class ProblemTDEM_b(ProblemBaseTDEM):
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: A
|
||||
"""
|
||||
|
||||
dt = self.timeSteps[tInd]
|
||||
return self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui + (1.0/dt)*self.MfMui
|
||||
|
||||
def getRHS(self, tInd, F):
|
||||
dt = self.timeSteps[tInd]
|
||||
return (1.0/dt)*self.MfMui*F.get_b(tInd-1)
|
||||
B_n = np.c_[[F[tx,'b',tInd] for tx in self.survey.txList]].T
|
||||
RHS = (1.0/dt)*self.MfMui*B_n
|
||||
return RHS
|
||||
|
||||
|
||||
####################################################
|
||||
@@ -50,15 +51,20 @@ class ProblemTDEM_b(ProblemBaseTDEM):
|
||||
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
|
||||
return - mkvc(w)
|
||||
|
||||
def Gvec(self, m, vec, u=None):
|
||||
"""
|
||||
@@ -68,63 +74,108 @@ class ProblemTDEM_b(ProblemBaseTDEM):
|
||||
:rtype: simpegEM.TDEM.FieldsTDEM
|
||||
:return: f
|
||||
|
||||
Multiply G by a vector where
|
||||
Multiply G by a vector
|
||||
"""
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
p = FieldsTDEM(self.mesh, 1, self.nT, 'b')
|
||||
|
||||
# Note: Fields has shape (nF/E, nTx, nT+1)
|
||||
# However, p will only really fill (:,:,1:nT+1)
|
||||
# meaning the 'initial fields' are zero (:,:,0)
|
||||
p = FieldsTDEM(self.mesh, self.survey)
|
||||
p[:, 'b', :] = 0.0 # b at all times is zero.
|
||||
p[:, 'e', 0] = 0.0 # fake initial fields
|
||||
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 i in range(1,self.nT+1):
|
||||
# TODO: G[1] may be dependent on the model
|
||||
# for a galvanic source (deriv of the dc problem)
|
||||
for tx in self.survey.txList:
|
||||
p[tx, 'e', i] = -u[tx,'e',i]*c # - diag(e) * MsigDeriv * v
|
||||
return p
|
||||
|
||||
def Gtvec(self, m, v, u=None):
|
||||
def Gtvec(self, m, vec, u=None):
|
||||
"""
|
||||
:param numpy.array m: Conductivity model
|
||||
:param numpy.array vec: vector (like a fields)
|
||||
:param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
|
||||
:rtype: np.ndarray (like a model)
|
||||
:return: p
|
||||
|
||||
Multiply G.T by a vector
|
||||
"""
|
||||
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)
|
||||
# Here we can do internal multiplications of Gt*v and then multiply by MsigDeriv.T in one go.
|
||||
for i in range(1,self.nT+1):
|
||||
vu = vec[:,'e',i]*u[:,'e',i]
|
||||
if nTx > 1:
|
||||
vu = vu.sum(axis=1)
|
||||
tmp += vu
|
||||
|
||||
curModel = self.mapping.transform(m)
|
||||
p = -mkvc(self.mapping.transformDeriv(m).T*self.mesh.getEdgeInnerProductDeriv(curModel).T*tmp)
|
||||
return p
|
||||
|
||||
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)
|
||||
|
||||
def AhRHS(tInd, y):
|
||||
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*y[:,'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)
|
||||
return {'b':b, 'e':e}
|
||||
y_b = sol
|
||||
if self.survey.nTx == 1:
|
||||
y_b = mkvc(y_b)
|
||||
y_e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*y_b - self.MeSigmaI*p[:,'e',tInd+1]
|
||||
return {'b':y_b, 'e':y_e}
|
||||
|
||||
self.curModel = m
|
||||
return self.forward(m, AhRHS, AhCalcFields)
|
||||
|
||||
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)
|
||||
# Mini Example:
|
||||
#
|
||||
# nT = 3, len(times) == 4, fields stored in F[:,:,1:4]
|
||||
#
|
||||
# 0 is held for initial conditions (this shifts the storage by +1)
|
||||
# ^
|
||||
# fLoc 0 1 2 3
|
||||
# |-----|-----|-----|
|
||||
# tInd 0 1 2 / /
|
||||
# / __/
|
||||
# 2 (tInd=2 uses fields 3 and would use 4 but it doesn't exist)
|
||||
# / __/
|
||||
# 1 (tInd=1 uses fields 2 and 3)
|
||||
|
||||
def AhtRHS(tInd, y):
|
||||
nTx, nF = self.survey.nTx, self.mesh.nF
|
||||
rhs = np.zeros(nF if nTx == 1 else (nF, nTx))
|
||||
|
||||
if 'e' in p:
|
||||
rhs += self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd+1]
|
||||
if 'b' in p:
|
||||
rhs += p[:,'b',tInd+1]
|
||||
|
||||
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*y[:,'b',tInd+2]
|
||||
|
||||
def AhtCalcFields(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_b = sol
|
||||
if self.survey.nTx == 1:
|
||||
y_b = mkvc(y_b)
|
||||
y_e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*y_b
|
||||
if 'e' in p:
|
||||
y_e += - self.MeSigmaI*p[:,'e',tInd]
|
||||
return {'b':y_b, 'e':y_e}
|
||||
|
||||
self.curModel = m
|
||||
return self.adjoint(m, AhtRHS, AhtCalcFields)
|
||||
@@ -168,18 +219,14 @@ class ProblemTDEM_b(ProblemBaseTDEM):
|
||||
"""
|
||||
|
||||
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):
|
||||
@@ -216,17 +263,13 @@ class ProblemTDEM_b(ProblemBaseTDEM):
|
||||
\\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
|
||||
|
||||
|
||||
|
||||
@@ -1,3 +1,3 @@
|
||||
from BaseTDEM import ProblemBaseTDEM
|
||||
from SurveyTDEM import SurveyTDEM1D, FieldsTDEM
|
||||
from SurveyTDEM import SurveyTDEM, FieldsTDEM, RxTDEM, TxTDEM
|
||||
from BaseTDEM import BaseTDEMProblem
|
||||
from TDEM_b import ProblemTDEM_b
|
||||
|
||||
@@ -21,23 +21,22 @@ 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-5]
|
||||
self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
|
||||
# self.prb.timeSteps = [(1e-05, 100)]
|
||||
|
||||
self.sigma = np.ones(mesh.nCz)*1e-8
|
||||
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):
|
||||
@@ -51,25 +50,26 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
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]
|
||||
|
||||
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):
|
||||
|
||||
@@ -86,8 +86,8 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
A = sp.bmat([[a11,a12],[a21,a22]])
|
||||
|
||||
f = prb.fields(sigma)
|
||||
u1 = A*f.fieldVec()
|
||||
u2 = prb.AhVec(sigma,f).fieldVec()
|
||||
u1 = A*f.tovec()
|
||||
u2 = prb.AhVec(sigma,f).tovec()
|
||||
|
||||
self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12)
|
||||
|
||||
@@ -100,20 +100,23 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
prb.curModel = sigma
|
||||
|
||||
dt = prb.timeSteps[0]
|
||||
a11 = 1/dt*prb.MfMui*sp.eye(prb.mesh.nF)
|
||||
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[:,:,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())
|
||||
|
||||
self.assertLess(np.linalg.norm(u1-u2),1e-8)
|
||||
|
||||
def test_solveAhVsAhVec(self):
|
||||
|
||||
@@ -122,16 +125,16 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
sigma = self.sigma
|
||||
self.prb.curModel = sigma
|
||||
|
||||
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)
|
||||
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)
|
||||
|
||||
Ahf = prb.AhVec(sigma, f)
|
||||
f_test = prb.solveAh(sigma, Ahf)
|
||||
|
||||
u1 = f.fieldVec()
|
||||
u2 = f_test.fieldVec()
|
||||
u1 = f.tovec()
|
||||
u2 = f_test.tovec()
|
||||
self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
|
||||
|
||||
def test_DerivG(self):
|
||||
@@ -146,7 +149,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
dm = 1000*np.random.rand(self.prb.mapping.nP)
|
||||
h = 0.01
|
||||
|
||||
derChk = lambda m: [self.prb.AhVec(m, f).fieldVec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).fieldVec()]
|
||||
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)
|
||||
@@ -161,7 +164,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
dm = 10*np.random.rand(prb.mapping.nP)
|
||||
f = prb.fields(sigma)
|
||||
|
||||
derChk = lambda m: [self.prb.fields(m).fieldVec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).fieldVec()]
|
||||
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)
|
||||
@@ -178,7 +181,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
d_sig = 10*np.random.rand(prb.mapping.nP)
|
||||
|
||||
|
||||
derChk = lambda m: [prb.survey.dpred(m), lambda mx: -prb.Jvec(sigma, mx)]
|
||||
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)
|
||||
@@ -186,19 +189,20 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
def test_projectAdjoint(self):
|
||||
prb = self.prb
|
||||
dat = self.dat
|
||||
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)
|
||||
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)
|
||||
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())
|
||||
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)
|
||||
|
||||
@@ -207,62 +211,63 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
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())
|
||||
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,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 = 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 True:
|
||||
# import matplotlib.pyplot as plt
|
||||
# plt.plot(f3.tovec(),'b')
|
||||
# plt.plot(f1.tovec(),'r')
|
||||
# plt.show()
|
||||
# V1 = np.linalg.norm(f3.tovec()-f1.tovec())
|
||||
# V2 = np.linalg.norm(f1.tovec())
|
||||
# print 'AhtVsAhtVec', V1, V2, f1.tovec()
|
||||
# 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
|
||||
# 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())
|
||||
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
|
||||
# V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec())
|
||||
# V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec())
|
||||
# print V1, V2
|
||||
# self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
|
||||
|
||||
def test_adjointGvecVsGtvec(self):
|
||||
mesh = self.mesh
|
||||
@@ -271,18 +276,18 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
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(1,prb.nT+1):
|
||||
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(1,prb.nT+1):
|
||||
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())
|
||||
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):
|
||||
@@ -291,10 +296,11 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
sigma = self.sigma
|
||||
|
||||
m = np.random.rand(prb.mapping.nP)
|
||||
d = np.random.rand(prb.nT)
|
||||
d = np.random.rand(prb.survey.nD)
|
||||
|
||||
V1 = d.dot(prb.Jvec(sigma, m))
|
||||
V2 = m.dot(prb.Jtvec(sigma, d))
|
||||
print 'AdjointTest', V1, V2
|
||||
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,150 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
import simpegEM as EM
|
||||
|
||||
plotIt = False
|
||||
|
||||
class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
cs = 5.
|
||||
ncx = 20
|
||||
ncy = 6
|
||||
npad = 20
|
||||
hx = [(cs,ncx), (cs,npad,1.3)]
|
||||
hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
|
||||
mesh = Mesh.CylMesh([hx,1,hy], '00C')
|
||||
|
||||
active = mesh.vectorCCz<0.
|
||||
activeMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
|
||||
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])
|
||||
rx2 = EM.TDEM.RxTDEM(np.array([[rxOffset-10, 0., 0.]]), np.logspace(-5,-4, 25), 'bz')
|
||||
tx2 = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx2])
|
||||
|
||||
survey = EM.TDEM.SurveyTDEM([tx,tx2])
|
||||
|
||||
self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
|
||||
# self.prb.timeSteps = [1e-5]
|
||||
self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
|
||||
# self.prb.timeSteps = [(1e-05, 100)]
|
||||
|
||||
self.sigma = np.ones(mesh.nCz)*1e-8
|
||||
self.sigma[mesh.vectorCCz<0] = 1e-1
|
||||
self.sigma = np.log(self.sigma[active])
|
||||
|
||||
self.prb.pair(survey)
|
||||
self.mesh = mesh
|
||||
|
||||
def test_DerivG(self):
|
||||
"""
|
||||
Test the derivative of c with respect to sigma
|
||||
"""
|
||||
|
||||
# Random model and perturbation
|
||||
sigma = np.random.rand(self.prb.mapping.nP)
|
||||
|
||||
f = self.prb.fields(sigma)
|
||||
dm = 1000*np.random.rand(self.prb.mapping.nP)
|
||||
h = 0.01
|
||||
|
||||
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)
|
||||
|
||||
def test_Deriv_dUdM(self):
|
||||
|
||||
prb = self.prb
|
||||
prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
|
||||
mesh = self.mesh
|
||||
sigma = self.sigma
|
||||
|
||||
dm = 10*np.random.rand(prb.mapping.nP)
|
||||
f = prb.fields(sigma)
|
||||
|
||||
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)
|
||||
|
||||
def test_Deriv_J(self):
|
||||
|
||||
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)
|
||||
|
||||
def test_projectAdjoint(self):
|
||||
prb = self.prb
|
||||
survey = prb.survey
|
||||
mesh = self.mesh
|
||||
|
||||
# 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)
|
||||
d = Survey.Data(survey,v=d_vec)
|
||||
|
||||
# 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)
|
||||
|
||||
def test_adjointGvecVsGtvec(self):
|
||||
mesh = self.mesh
|
||||
prb = self.prb
|
||||
|
||||
m = np.random.rand(prb.mapping.nP)
|
||||
sigma = np.random.rand(prb.mapping.nP)
|
||||
|
||||
u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
for i in range(1,prb.nT+1):
|
||||
u[:,'b',i] = np.random.rand(mesh.nF, 2)
|
||||
u[:,'e',i] = np.random.rand(mesh.nE, 2)
|
||||
|
||||
v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
for i in range(1,prb.nT+1):
|
||||
v[:,'b',i] = np.random.rand(mesh.nF, 2)
|
||||
v[:,'e',i] = np.random.rand(mesh.nE, 2)
|
||||
|
||||
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
|
||||
|
||||
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))
|
||||
print 'AdjointTest', V1, V2
|
||||
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -22,22 +22,12 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,
|
||||
actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
|
||||
mapping = Maps.ComboMap(mesh, [Maps.ExpMap, Maps.Vertical1DMap, actMap])
|
||||
|
||||
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-5,-4, 21), 'bz')
|
||||
tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx])
|
||||
|
||||
opts = {'txLoc':np.array([0., 0., 0.]),
|
||||
'txType':'VMD_MVP',
|
||||
'rxLoc':np.array([rxOffset, 0., 0.]),
|
||||
'rxType':'bz',
|
||||
'timeCh':np.logspace(-5,-4, 21),
|
||||
}
|
||||
|
||||
survey = EM.TDEM.SurveyTDEM1D(**opts)
|
||||
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)]
|
||||
|
||||
@@ -46,16 +36,17 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,
|
||||
sigma = np.log(sigma[active])
|
||||
prb.pair(survey)
|
||||
|
||||
bz_ana = mu_0*EM.Utils.Ana.hzAnalyticDipoleT(survey.rxLoc[0]+1e-3, prb.times[1:], sig_half)
|
||||
bz_ana = mu_0*EM.Utils.Ana.hzAnalyticDipoleT(rx.locs[0][0]+1e-3, rx.times, sig_half)
|
||||
|
||||
bz_calc = survey.dpred(sigma)
|
||||
ind = np.logical_and(prb.times[1:] > bounds[0],prb.times[1:] < bounds[1])
|
||||
|
||||
ind = np.logical_and(rx.times > bounds[0],rx.times < bounds[1])
|
||||
log10diff = np.linalg.norm(np.log10(np.abs(bz_calc[ind])) - np.log10(np.abs(bz_ana[ind])))/np.linalg.norm(np.log10(np.abs(bz_ana[ind])))
|
||||
print 'Difference: ', log10diff
|
||||
|
||||
if showIt == True:
|
||||
plt.loglog(prb.times[1:][bz_calc>0], bz_calc[bz_calc>0], 'r', prb.times[1:][bz_calc<0], -bz_calc[bz_calc<0], 'r--')
|
||||
plt.loglog(prb.times[1:], abs(bz_ana), 'b*')
|
||||
plt.loglog(rx.times[bz_calc>0], bz_calc[bz_calc>0], 'r', rx.times[bz_calc<0], -bz_calc[bz_calc<0], 'r--')
|
||||
plt.loglog(rx.times, abs(bz_ana), 'b*')
|
||||
plt.title('sig_half = %e'%sig_half)
|
||||
plt.show()
|
||||
|
||||
|
||||
Reference in New Issue
Block a user