wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-06 05:16:51 +08:00
Code Issues Packages Projects Releases Wiki Activity

Forward problem working. Not yet tested with multi Tx.

Browse Source
This commit is contained in:
rowanc1
2014-04-26 23:15:31 -07:00
parent 78b3a281f4
commit 9ebbe7613d
6 changed files with 134 additions and 78 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/FDEM/FDEM.py
+5 -5
View File
@@ -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):
"""
simpegEM/TDEM/BaseTDEM.py
+12 -50
View File
@@ -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,14 +47,13 @@ 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):
@@ -104,7 +66,7 @@ class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem, BaseEMProblem):
dtFact = dt
A = self.getA(tInd)
# print 'Factoring... (dt = ' + str(dt) + ')'
Asolve = Solver(A, options=self.solveOpts)
Asolve = Solver(A, options=self.solverOpts)
# print 'Done'
rhs = RHS(tInd, F)
sol = Asolve.solve(rhs)
simpegEM/TDEM/SurveyTDEM.py
+99 -2
View File
@@ -1,5 +1,102 @@
from SimPEG import Utils, np
from SimPEG import Utils, Survey, np
from SimPEG.Survey import BaseSurvey
from simpegEM.Utils import Sources
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 projField(self):
"""Field Type projection (e.g. e b ...)"""
return self.knownRxTypes[self.rxType][0]
@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 * v
elif adjoint:
return P.T * v
class FieldsTDEM(Survey.TimeFields):
"""Fancy Field Storage for a TDEM survey."""
knownFields = {'b': 'F', 'e': 'E'}
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):
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 projectFieldsDeriv(self, u):
raise Exception('Use Transmitters to project fields deriv.')
class SurveyTDEM1D(BaseSurvey):
"""
@@ -52,7 +149,7 @@ class SurveyTDEM1D(BaseSurvey):
_Qrx = None
class FieldsTDEM(object):
class FieldsTDEM_OLD(object):
"""docstring for FieldsTDEM"""
phi0 = None #: Initial electric potential
simpegEM/TDEM/TDEM_b.py
+8 -7
View File
@@ -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.concatenate([F[tx,'b',tInd] for tx in self.survey.txList], axis=1)
RHS = (1.0/dt)*self.MfMui*B_n
return RHS
####################################################
simpegEM/TDEM/__init__.py
+2 -2
View File
@@ -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
simpegEM/Tests/test_TDEM_forward_Analytic.py
+8 -12
View File
@@ -22,15 +22,10 @@ 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:
@@ -46,16 +41,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()