wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-03 11:45:39 +08:00
Code Issues Packages Projects Releases Wiki Activity

Mesh and model updates.

Browse Source
This commit is contained in:
rowanc1
2014-05-18 16:39:10 -07:00
parent 6e7d88de77
commit 7a15797fbe
8 changed files with 41 additions and 34 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/Base.py
+4 -16
View File
@@ -1,4 +1,4 @@
from SimPEG import Survey, Problem, Utils, np, sp, Solver as SimpegSolver
from SimPEG import Survey, Problem, Utils, Models, np, sp, Solver as SimpegSolver
from scipy.constants import mu_0
class BaseEMProblem(Problem.BaseProblem):
@@ -38,7 +38,7 @@ class BaseEMProblem(Problem.BaseProblem):
def MeSigma(self):
#TODO: hardcoded to sigma as the model
if getattr(self, '_MeSigma', None) is None:
sigma = self.curTModel
sigma = self.curModel.transform
self._MeSigma = self.mesh.getEdgeInnerProduct(sigma)
return self._MeSigma
@@ -46,23 +46,11 @@ class BaseEMProblem(Problem.BaseProblem):
def MeSigmaI(self):
#TODO: hardcoded to sigma as the model
if getattr(self, '_MeSigmaI', None) is None:
sigma = self.curTModel
sigma = self.curModel.transform
self._MeSigmaI = self.mesh.getEdgeInnerProduct(sigma, invMat=True)
return self._MeSigmaI
curModel = Utils.dependentProperty('_curModel', None, ['_MeSigma', '_MeSigmaI', '_curTModel', '_curTModelDeriv'], 'Sets the current model, and removes dependent mass matrices.')
@property
def curTModel(self):
if getattr(self, '_curTModel', None) is None:
self._curTModel = self.mapping.transform(self.curModel)
return self._curTModel
@property
def curTModelDeriv(self):
if getattr(self, '_curTModelDeriv', None) is None:
self._curTModelDeriv = self.mapping.transformDeriv(self.curModel)
return self._curTModelDeriv
deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI']
def fields(self, m):
self.curModel = m
simpegEM/FDEM/FDEM.py
+6 -6
View File
@@ -142,8 +142,8 @@ class ProblemFDEM_e(BaseFDEMProblem):
return C.T*mui*C + 1j*omega(freq)*sig
def getADeriv(self, freq, u, v, adjoint=False):
sig = self.curTModel
dsig_dm = self.curTModelDeriv
sig = self.curModel.transform
dsig_dm = self.curModel.transformDeriv
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=u)
if adjoint:
@@ -222,8 +222,8 @@ class ProblemFDEM_b(BaseFDEMProblem):
mui = self.MfMui
C = self.mesh.edgeCurl
sig = self.curTModel
dsig_dm = self.curTModelDeriv
sig = self.curModel.transform
dsig_dm = self.curModel.transformDeriv
#TODO: This only works if diagonal (no tensors)...
dMeSigmaI_dI = - self.MeSigmaI**2
@@ -275,8 +275,8 @@ class ProblemFDEM_b(BaseFDEMProblem):
def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
b = sol
if fieldType == 'e':
sig = self.curTModel
dsig_dm = self.curTModelDeriv
sig = self.curModel.transform
dsig_dm = self.curModel.transformDeriv
C = self.mesh.edgeCurl
mui = self.MfMui
simpegEM/TDEM/BaseTDEM.py
+14 -4
View File
@@ -18,14 +18,16 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
_FieldsTDEM_pair = FieldsTDEM #: used for the forward calculation only
def fields(self, m):
if self.verbose: print '%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50)
self.curModel = m
# Create a fields storage object
F = self._FieldsTDEM_pair(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)
F = self.forward(m, self.getRHS, self.calcFields, F=F)
if self.verbose: print '%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50)
return F
def forward(self, m, RHS, CalcFields, F=None):
self.curModel = m
@@ -39,11 +41,13 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
if Ainv is not None:
Ainv.clean()
A = self.getA(tInd)
if self.verbose: print 'Factoring... (dt = ' + str(dt) + ')'
if self.verbose: print 'Factoring... (dt = %e)'%dt
Ainv = self.Solver(A, **self.solverOpts)
if self.verbose: print 'Done'
rhs = RHS(tInd, F)
if self.verbose: print ' Solving... (tInd = %d)'%tInd
sol = Ainv * rhs
if self.verbose: print ' Done...'
if sol.ndim == 1:
sol.shape = (sol.size,1)
F[:,:,tInd+1] = CalcFields(sol, tInd)
@@ -62,11 +66,13 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
if Ainv is not None:
Ainv.clean()
A = self.getA(tInd)
if self.verbose: print 'Factoring... (dt = ' + str(dt) + ')'
if self.verbose: print 'Factoring (Adjoint)... (dt = %e)'%dt
Ainv = self.Solver(A, **self.solverOpts)
if self.verbose: print 'Done'
rhs = RHS(tInd, F)
if self.verbose: print ' Solving (Adjoint)... (tInd = %d)'%tInd
sol = Ainv * rhs
if self.verbose: print ' Done...'
if sol.ndim == 1:
sol.shape = (sol.size,1)
F[:,:,tInd+1] = CalcFields(sol, tInd)
@@ -88,12 +94,14 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
* Compute \\\(\\\\vec{w} = -\\\mathbf{Q} \\\\vec{y}\\\)
"""
if self.verbose: print '%s\nCalculating J(v)\n%s'%('*'*50,'*'*50)
self.curModel = m
if u is None:
u = self.fields(m)
p = self.Gvec(m, v, u)
y = self.solveAh(m, p)
Jv = self.survey.projectFieldsDeriv(u, v=y)
if self.verbose: print '%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50)
return - mkvc(Jv)
def Jtvec(self, m, v, u=None):
@@ -111,6 +119,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
* Compute \\\(\\\\vec{w} = -\\\mathbf{G}^\\\\top y\\\)
"""
if self.verbose: print '%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50)
self.curModel = m
if u is None:
u = self.fields(m)
@@ -121,5 +130,6 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
p = self.survey.projectFieldsDeriv(u, v=v, adjoint=True)
y = self.solveAht(m, p)
w = self.Gtvec(m, y, u)
if self.verbose: print '%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50)
return - mkvc(w)
simpegEM/TDEM/TDEM_b.py
+2 -2
View File
@@ -97,7 +97,7 @@ class ProblemTDEM_b(BaseTDEMProblem):
# fake initial 'e' fields
p[:, 'e', 0] = 0.0
c = self.mesh.getEdgeInnerProductDeriv(self.curTModel)*(self.curTModelDeriv*vec)
c = self.mesh.getEdgeInnerProductDeriv(self.curModel.transform)*(self.curModel.transformDeriv*vec)
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)
@@ -127,7 +127,7 @@ class ProblemTDEM_b(BaseTDEMProblem):
if nTx > 1:
vu = vu.sum(axis=1)
tmp += vu
p = -mkvc(self.curTModelDeriv.T*(self.mesh.getEdgeInnerProductDeriv(self.curTModel).T*tmp))
p = -mkvc(self.curModel.transformDeriv.T*(self.mesh.getEdgeInnerProductDeriv(self.curModel.transform).T*tmp))
return p
def solveAh(self, m, p):
simpegEM/Tests/test_FDEM_analytics.py
+6 -1
View File
@@ -28,7 +28,12 @@ class FDEM_analyticTests(unittest.TestCase):
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
prb.pair(survey)
prb.Solver = SolverLU
try:
from pymatsolver import MumpsSolver
prb.Solver = MumpsSolver
except ImportError, e:
prb.Solver = SolverLU
sig = 1e-1
sigma = np.ones(mesh.nC)*sig
simpegEM/Tests/test_TDEM_b_DerivAdjoint.py
+7 -2
View File
@@ -18,8 +18,7 @@ class TDEM_bDerivTests(unittest.TestCase):
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])
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * activeMap
rxOffset = 40.
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
@@ -32,6 +31,12 @@ class TDEM_bDerivTests(unittest.TestCase):
self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
# self.prb.timeSteps = [(1e-05, 100)]
try:
from pymatsolver import MumpsSolver
self.prb.Solver = MumpsSolver
except ImportError, e:
self.prb.Solver = SolverLU
self.sigma = np.ones(mesh.nCz)*1e-8
self.sigma[mesh.vectorCCz<0] = 1e-1
self.sigma = np.log(self.sigma[active])
simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py
+1 -2
View File
@@ -18,8 +18,7 @@ class TDEM_bDerivTests(unittest.TestCase):
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])
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * activeMap
rxOffset = 40.
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
simpegEM/Tests/test_TDEM_combos.py
+1 -1
View File
@@ -20,7 +20,7 @@ def getProb(meshType='CYL',rxTypes='bx,bz',nTx=1):
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])
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * activeMap
rxOffset = 40.