wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-08-02 13:00:17 +08:00
Code Issues Packages Projects Releases Wiki Activity

Updates to for mapping changes in SimPEG

Browse Source
This commit is contained in:
rowanc1
2014-04-15 12:07:58 -07:00
parent 113c97ca87
commit eab87c892b
7 changed files with 61 additions and 47 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/FDEM/FDEM.py
+3 -3
View File
@@ -66,13 +66,13 @@ class BaseProblemFDEM(Problem.BaseProblem):
@property
def curTModel(self):
if getattr(self, '_curTModel', None) is None:
self._curTModel = self.model.transform(self.curModel)
self._curTModel = self.mapping.transform(self.curModel)
return self._curTModel
@property
def curTModelDeriv(self):
if getattr(self, '_curTModelDeriv', None) is None:
self._curTModelDeriv = self.model.transformDeriv(self.curModel)
self._curTModelDeriv = self.mapping.transformDeriv(self.curModel)
return self._curTModelDeriv
def fields(self, m):
@@ -132,7 +132,7 @@ class BaseProblemFDEM(Problem.BaseProblem):
if not isinstance(v, self.dataPair):
v = self.dataPair(self.survey, v)
Jtv = np.zeros(self.model.nP)
Jtv = np.zeros(self.mapping.nP)
for freq in self.survey.freqs:
AT = self.getA(freq).T
simpegEM/TDEM/BaseTDEM.py
+3 -3
View File
@@ -92,8 +92,8 @@ class MixinTimeStuff(object):
class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
"""docstring for ProblemTDEM1D"""
def __init__(self, model, **kwargs):
BaseProblem.__init__(self, model, **kwargs)
def __init__(self, mesh, mapping=None, **kwargs):
BaseProblem.__init__(self, mesh, mapping=mapping, **kwargs)
####################################################
@@ -122,7 +122,7 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
def MeSigmaI(self): return self._MeSigmaI
def makeMassMatrices(self, m):
sig = self.model.transform(m)
sig = self.mapping.transform(m)
self._MeSigma = self.mesh.getEdgeInnerProduct(sig)
self._MeSigmaI = sdiag(1/self.MeSigma.diagonal())
self._MfMui = self.mesh.getFaceInnerProduct(1/mu_0)
simpegEM/TDEM/TDEM_b.py
+13 -10
View File
@@ -16,8 +16,8 @@ class ProblemTDEM_b(ProblemBaseTDEM):
with \\\(\\b\\\) defined on cell faces and \\\(\e\\\) defined on edges.
"""
def __init__(self, model, **kwargs):
ProblemBaseTDEM.__init__(self, model, **kwargs)
def __init__(self, mesh, mapping=None, **kwargs):
ProblemBaseTDEM.__init__(self, mesh, mapping=mapping, **kwargs)
solType = 'b'
@@ -35,11 +35,11 @@ class ProblemTDEM_b(ProblemBaseTDEM):
"""
dt = self.getDt(tInd)
return self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui + (1/dt)*self.MfMui
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.getDt(tInd)
return (1/dt)*self.MfMui*F.get_b(tInd-1)
return (1.0/dt)*self.MfMui*F.get_b(tInd-1)
####################################################
@@ -74,7 +74,8 @@ class ProblemTDEM_b(ProblemBaseTDEM):
if u is None:
u = self.fields(m)
p = FieldsTDEM(self.mesh, 1, self.times.size, 'b')
c = self.mesh.getEdgeMassDeriv()*self.model.transformDeriv(m)*vec
curModel = self.mapping.transform(m)
c = self.mesh.getEdgeInnerProductDeriv(curModel)*self.mapping.transformDeriv(m)*vec
for i in range(self.times.size):
ei = u.get_e(i)
pVal = np.empty_like(ei)
@@ -91,7 +92,9 @@ class ProblemTDEM_b(ProblemBaseTDEM):
tmp = np.zeros((self.mesh.nE,self.survey.nTx))
for i in range(self.nTimes):
tmp += v.get_e(i)*u.get_e(i)
p = -mkvc(self.model.transformDeriv(m).T*self.mesh.getEdgeMassDeriv().T*tmp)
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):
@@ -100,7 +103,7 @@ class ProblemTDEM_b(ProblemBaseTDEM):
if tInd == 0:
return rhs
dt = self.getDt(tInd)
return rhs + 1./dt*self.MfMui*u.get_b(tInd-1)
return rhs + 1.0/dt*self.MfMui*u.get_b(tInd-1)
def AhCalcFields(sol, solType, tInd):
b = sol
@@ -117,7 +120,7 @@ class ProblemTDEM_b(ProblemBaseTDEM):
if tInd == self.nTimes-1:
return rhs
dt = self.getDt(tInd+1)
return rhs + 1./dt*self.MfMui*u.get_b(tInd+1)
return rhs + 1.0/dt*self.MfMui*u.get_b(tInd+1)
def AhtCalcFields(sol, solType, tInd):
b = sol
@@ -167,14 +170,14 @@ class ProblemTDEM_b(ProblemBaseTDEM):
self.makeMassMatrices(sigma)
dt = self.getDt(0)
b = 1/dt*self.MfMui*vec.get_b(0) + self.MfMui*self.mesh.edgeCurl*vec.get_e(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.times.size, 'b')
f.set_b(b, 0)
f.set_e(e, 0)
for i in range(1,self.nTimes):
dt = self.getDt(i)
b = 1/dt*self.MfMui*vec.get_b(i) + self.MfMui*self.mesh.edgeCurl*vec.get_e(i) - 1/dt*self.MfMui*vec.get_b(i-1)
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)
simpegEM/Tests/test_FDEM.py
+10 -4
View File
@@ -15,7 +15,7 @@ def getProblem(fdemType, comp):
hz = Utils.meshTensors(((npad,cs), (ncz,cs), (npad,cs)))
mesh = Mesh.TensorMesh([hx,hy,hz],[-hx.sum()/2., -hy.sum()/2., -hz.sum()/2.])
model = Model.LogModel(mesh)
mapping = Maps.ExpMap(mesh)
x = np.linspace(-30,30,6)
XYZ = Utils.ndgrid(x,x,np.r_[0])
@@ -25,13 +25,19 @@ def getProblem(fdemType, comp):
survey = EM.FDEM.SurveyFDEM([Tx0])
if fdemType == 'e':
prb = EM.FDEM.ProblemFDEM_e(model)
prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping)
elif fdemType == 'b':
prb = EM.FDEM.ProblemFDEM_b(model)
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
else:
raise NotImplementedError()
prb.pair(survey)
try:
from mumpsSCI import MumpsSolver
prb.Solver = MumpsSolver
except ImportError, e:
pass
return prb
def adjointTest(fdemType, comp):
@@ -42,7 +48,7 @@ def adjointTest(fdemType, comp):
survey = prb.survey
v = np.random.rand(survey.nD)
w = np.random.rand(prb.model.nP)
w = np.random.rand(prb.mapping.nP)
u = prb.fields(m)
vJw = v.dot(prb.Jvec(m, w, u=u))
simpegEM/Tests/test_FDEM_analytics.py
+2 -2
View File
@@ -17,7 +17,7 @@ class FDEM_analyticTests(unittest.TestCase):
hz = Utils.meshTensors(((npad,cs), (ncz,cs), (npad,cs)))
mesh = Mesh.TensorMesh([hx,hy,hz], x0=[-hx.sum()/2.,-hy.sum()/2.,-hz.sum()/2.,])
model = Model.LogModel(mesh)
mapping = Maps.ExpMap(mesh)
x = np.linspace(-10,10,5)
XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
@@ -26,7 +26,7 @@ class FDEM_analyticTests(unittest.TestCase):
survey = EM.FDEM.SurveyFDEM([Tx0])
prb = EM.FDEM.ProblemFDEM_b(model)
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
prb.pair(survey)
prb.Solver = Utils.SolverUtils.DSolverWrap(sp.linalg.splu, checkAccuracy=False)
simpegEM/Tests/test_TDEM_b_DerivAdjoint.py
+16 -15
View File
@@ -12,23 +12,23 @@ class TDEM_bDerivTests(unittest.TestCase):
npad = 20
hx = Utils.meshTensors(((0,cs), (ncx,cs), (npad,cs)))
hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs)))
mesh = Mesh.CylMesh([hx,hy], -hy.sum()/2)
mesh = Mesh.CylMesh([hx,1,hy], [0,0,-hy.sum()/2])
active = mesh.vectorCCz<0.
model = Model.ActiveModel(mesh, active, -8, nC=mesh.nCz)
model = Model.ComboModel(mesh,
[Model.LogModel, Model.Vertical1DModel, model])
activeMap = Maps.ActiveCells(mesh, active, -8, nC=mesh.nCz)
mapping = Maps.ComboMap(mesh,
[Maps.ExpMap, Maps.Vertical1DMap, activeMap])
opts = {'txLoc':0.,
'txType': 'VMD_MVP',
'rxLoc':np.r_[150., 0.],
'rxLoc':np.r_[40., 0., 0.],
'rxType':'bz',
'timeCh':np.logspace(-4,-2,20),
}
self.dat = EM.TDEM.SurveyTDEM1D(**opts)
self.prb = EM.TDEM.ProblemTDEM_b(model)
self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
self.prb.setTimes([1e-5, 5e-5, 2.5e-4], [10, 10, 10])
self.sigma = np.ones(mesh.nCz)*1e-8
@@ -50,7 +50,8 @@ class TDEM_bDerivTests(unittest.TestCase):
Ahu = prb.AhVec(sigma, u)
V1 = Ahu.get_b(0)
V2 = 1/prb.getDt(0)*prb.MfMui*u.get_b(-1)
V2 = 1./prb.getDt(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.get_e(0)
@@ -137,10 +138,10 @@ class TDEM_bDerivTests(unittest.TestCase):
"""
# Random model and perturbation
sigma = np.random.rand(self.prb.model.nP)
sigma = np.random.rand(self.prb.mapping.nP)
f = self.prb.fields(sigma)
dm = 1000*np.random.rand(self.prb.model.nP)
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()]
@@ -155,7 +156,7 @@ class TDEM_bDerivTests(unittest.TestCase):
mesh = self.mesh
sigma = self.sigma
dm = 10*np.random.rand(prb.model.nP)
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()]
@@ -172,7 +173,7 @@ class TDEM_bDerivTests(unittest.TestCase):
sigma = self.sigma
# d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
d_sig = 10*np.random.rand(prb.model.nP)
d_sig = 10*np.random.rand(prb.mapping.nP)
derChk = lambda m: [prb.survey.dpred(m), lambda mx: -prb.Jvec(sigma, mx)]
@@ -221,7 +222,7 @@ class TDEM_bDerivTests(unittest.TestCase):
def test_solveAhtVsAhtVec(self):
prb = self.prb
mesh = self.mesh
sigma = np.random.rand(prb.model.nP)
sigma = np.random.rand(prb.mapping.nP)
f1 = EM.TDEM.FieldsTDEM(mesh, 1, prb.nTimes, 'b')
for i in range(f1.nTimes):
@@ -258,8 +259,8 @@ class TDEM_bDerivTests(unittest.TestCase):
mesh = self.mesh
prb = self.prb
m = np.random.rand(prb.model.nP)
sigma = np.random.rand(prb.model.nP)
m = np.random.rand(prb.mapping.nP)
sigma = np.random.rand(prb.mapping.nP)
u = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b')
for i in range(u.nTimes):
@@ -280,7 +281,7 @@ class TDEM_bDerivTests(unittest.TestCase):
prb = self.prb
sigma = self.sigma
m = np.random.rand(prb.model.nP)
m = np.random.rand(prb.mapping.nP)
d = np.random.rand(prb.nTimes)
V1 = d.dot(prb.Jvec(sigma, m))
simpegEM/Tests/test_forward_EMproblem.py
+14 -10
View File
@@ -26,9 +26,10 @@ class TDEM_bTests(unittest.TestCase):
mesh = Mesh.TensorMesh([hx,hy,hz], [-hx.sum()/2.,-hy.sum()/2.,-hz.sum()/2.])
active = mesh.vectorCCz<0.
model = Model.ActiveModel(mesh, active, -8, nC=mesh.nCz)
model = Model.ComboModel(mesh,
[Model.LogModel, Model.Vertical1DModel, model])
actMap = Maps.ActiveCells(mesh, active, -8, nC=mesh.nCz)
mapping = Maps.ComboMap(mesh,
[Maps.ExpMap, Maps.Vertical1DMap, actMap])
opts = {'txLoc':np.array([0., 0., 0.]),
'txType':'VMD_MVP',
@@ -40,11 +41,14 @@ class TDEM_bTests(unittest.TestCase):
sig_half = 1e-3
self.sig_half = sig_half
self.dat = EM.TDEM.SurveyTDEM1D(**opts)
self.prb = EM.TDEM.ProblemTDEM_b(model)
from mumpsSCI import MumpsSolver
self.prb.Solver = MumpsSolver
self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
try:
from mumpsSCI import MumpsSolver
self.prb.Solver = MumpsSolver
except ImportError:
pass
self.prb.setTimes([1e-6], [100])
# self.prb.setTimes([1e-6, 5e-6, 1e-5, 5e-5, 1e-4, 5e-4], [40, 40, 40, 40, 40, 40])
self.sigma = np.ones(mesh.nCz)*1e-8
self.sigma[active] = sig_half
@@ -78,15 +82,15 @@ class TDEM_bTests(unittest.TestCase):
chi_half = 0.
Logmodel = EM1D.BaseEM1D.BaseEM1DModel(mesh1D)
modelReal = Model.ComboModel(mesh1D, [Logmodel])
expmap = EM1D.BaseEM1D.BaseEM1DMap(mesh1D)
mappingReal = Maps.ComboMap(mesh1D, [expmap])
m_1D = np.log(np.ones(nlay)*sig_half)
TDsurvey.rxType = 'Bz'
WT0, WT1, YBASE = EM1D.DigFilter.LoadWeights()
options = {'WT0': WT0, 'WT1': WT1, 'YBASE': YBASE}
prob = EM1D.EM1D.EM1D(modelReal, **options)
prob = EM1D.EM1D.EM1D(mesh1D, mapping=mappingReal, **options)
prob.pair(TDsurvey)
prob.chi = np.zeros(TDsurvey.nlay)