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

Make a base EM problem that is shared between TD and FD

Browse Source
This commit is contained in:
rowanc1
2014-04-26 18:55:02 -07:00
parent 9a801a87d5
commit 05bb79f710
7 changed files with 94 additions and 116 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/Base.py
+68
View File
@@ -0,0 +1,68 @@
from SimPEG import Survey, Problem, Utils, np, sp, Solver as SimpegSolver
from scipy.constants import mu_0
class BaseEMProblem(Problem.BaseProblem):
def __init__(self, mesh, **kwargs):
Problem.BaseProblem.__init__(self, mesh, **kwargs)
solType = None
storeTheseFields = ['e', 'b']
surveyPair = Survey.BaseSurvey
dataPair = Survey.Data
Solver = SimpegSolver
solverOpts = {}
####################################################
# Mass Matrices
####################################################
@property
def MfMui(self):
#TODO: assuming constant mu
if getattr(self, '_MfMui', None) is None:
self._MfMui = self.mesh.getFaceInnerProduct(1/mu_0)
return self._MfMui
@property
def Me(self):
if getattr(self, '_Me', None) is None:
self._Me = self.mesh.getEdgeInnerProduct()
return self._Me
@property
def MeSigma(self):
#TODO: hardcoded to sigma as the model
if getattr(self, '_MeSigma', None) is None:
sigma = self.curTModel
self._MeSigma = self.mesh.getEdgeInnerProduct(sigma)
return self._MeSigma
@property
def MeSigmaI(self):
#TODO: hardcoded to sigma as the model
if getattr(self, '_MeSigmaI', None) is None:
sigma = self.curTModel
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
def fields(self, m):
self.curModel = m
F = self.forward(m, self.getRHS, self.calcFields)
return F
simpegEM/FDEM/FDEM.py
+2 -62
View File
@@ -2,12 +2,13 @@ from SimPEG import Survey, Problem, Utils, np, sp, Solver as SimpegSolver
from scipy.constants import mu_0
from SurveyFDEM import SurveyFDEM, FieldsFDEM
from simpegEM.Utils import Sources
from simpegEM.Base import BaseEMProblem
def omega(freq):
"""Change frequency to angular frequency, omega"""
return 2.*np.pi*freq
class BaseProblemFDEM(Problem.BaseProblem):
class BaseProblemFDEM(BaseEMProblem):
"""
We start by looking at Maxwell's equations in the electric field \\(\\vec{E}\\) and the magnetic flux density \\(\\vec{B}\\):
@@ -17,69 +18,8 @@ class BaseProblemFDEM(Problem.BaseProblem):
\\nabla \\times \\mu^{-1} \\vec{B} - \\sigma \\vec{E} = \\vec{J_s}
"""
def __init__(self, model, **kwargs):
Problem.BaseProblem.__init__(self, model, **kwargs)
solType = None
storeTheseFields = ['e', 'b']
surveyPair = SurveyFDEM
dataPair = Survey.Data
Solver = SimpegSolver
solverOpts = {}
####################################################
# Mass Matrices
####################################################
@property
def MfMui(self):
#TODO: assuming constant mu
if getattr(self, '_MfMui', None) is None:
self._MfMui = self.mesh.getFaceInnerProduct(1/mu_0)
return self._MfMui
@property
def Me(self):
if getattr(self, '_Me', None) is None:
self._Me = self.mesh.getEdgeInnerProduct()
return self._Me
@property
def MeSigma(self):
#TODO: hardcoded to sigma as the model
if getattr(self, '_MeSigma', None) is None:
sigma = self.curTModel
self._MeSigma = self.mesh.getEdgeInnerProduct(sigma)
return self._MeSigma
@property
def MeSigmaI(self):
#TODO: hardcoded to sigma as the model
if getattr(self, '_MeSigmaI', None) is None:
sigma = self.curTModel
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
def fields(self, m):
self.curModel = m
F = self.forward(m, self.getRHS, self.calcFields)
return F
def forward(self, m, RHS, CalcFields):
simpegEM/FDEM/RHSem.py
+3 -4
View File
@@ -50,10 +50,6 @@ def path2edgeModel(mesh, pts):
edgeModel = np.r_[edm_x, edm_y, edm_z]
return edgeModel
def rho(x1, y1, x, y):
r = np.sqrt((x-x1)**2+(y-y1)**2)
return r
def MMRhalf(loc1, loc2, x, y):
""" Anaytic function for MMR response (B^{1D})
- loc1=(x1,y1): x, y location for (+) charge
@@ -67,6 +63,9 @@ def MMRhalf(loc1, loc2, x, y):
y2=loc2[1]
mu0 = 4*np.pi*1e-7
I = 1
rho = lambda x1, y1, x, y: np.sqrt((x-x1)**2+(y-y1)**2)
By =mu0*I/(4*np.pi)*np.array((x-x1)/rho(x1,y1,x,y)**2-(x-x2)/rho(x2,y2,x,y)**2)
Bx =mu0*I/(4*np.pi)*np.array(-(y-y1)/rho(x1,y1,x,y)**2+(y-y2)/rho(x2,y2,x,y)**2)
simpegEM/TDEM/BaseTDEM.py
+5 -36
View File
@@ -5,6 +5,7 @@ from SurveyTDEM import FieldsTDEM
from scipy.constants import mu_0
from SimPEG.Utils import sdiag, mkvc
from SimPEG import Utils, Mesh
from simpegEM.Base import BaseEMProblem
import numpy as np
@@ -46,44 +47,12 @@ class MixinInitialFieldCalc(object):
return F
class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem):
class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem, BaseEMProblem):
"""docstring for ProblemTDEM1D"""
def __init__(self, mesh, mapping=None, **kwargs):
BaseTimeProblem.__init__(self, mesh, mapping=mapping, **kwargs)
####################################################
# Physical Properties
####################################################
@property
def sigma(self):
return self._sigma
@sigma.setter
def sigma(self, value):
self._sigma = value
_sigma = None
####################################################
# Mass Matrices
####################################################
@property
def MfMui(self): return self._MfMui
@property
def MeSigma(self): return self._MeSigma
@property
def MeSigmaI(self): return self._MeSigmaI
def makeMassMatrices(self, m):
sig = self.mapping.transform(m)
self._MeSigma = self.mesh.getEdgeInnerProduct(sig)
self._MeSigmaI = Utils.sdInv(self.MeSigma)
self._MfMui = self.mesh.getFaceInnerProduct(1.0/mu_0)
def calcFields(self, sol, solType, tInd):
if solType == 'b':
@@ -100,7 +69,7 @@ class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem):
solveOpts = {}
def fields(self, m):
self.makeMassMatrices(m)
self.curModel = m
F = self.getInitialFields()
return self.forward(m, self.getRHS, self.calcFields, F=F)
@@ -111,7 +80,7 @@ class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem):
dtFact = None
for tInd, dt in enumerate(self.timeSteps):
if dt!=dtFact:
if dt != dtFact:
dtFact = dt
A = self.getA(tInd)
# print 'Factoring... (dt = ' + str(dt) + ')'
@@ -131,7 +100,7 @@ class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem):
dtFact = None
for tInd, dt in reversed(list(enumerate(self.timeSteps))):
if dt!=dtFact:
if dt != dtFact:
dtFact = dt
A = self.getA(tInd)
# print 'Factoring... (dt = ' + str(dt) + ')'
simpegEM/TDEM/TDEM_b.py
+10 -10
View File
@@ -109,7 +109,7 @@ class ProblemTDEM_b(ProblemBaseTDEM):
e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b - self.MeSigmaI*p.get_e(tInd)
return {'b':b, 'e':e}
self.makeMassMatrices(m)
self.curModel = m
return self.forward(m, AhRHS, AhCalcFields)
def solveAht(self, m, p):
@@ -126,16 +126,16 @@ class ProblemTDEM_b(ProblemBaseTDEM):
e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b - self.MeSigmaI*p.get_e(tInd)
return {'b':b, 'e':e}
self.makeMassMatrices(m)
self.curModel = m
return self.adjoint(m, AhtRHS, AhtCalcFields)
####################################################
# Functions for tests
####################################################
def AhVec(self, sigma, vec):
def AhVec(self, m, vec):
"""
:param numpy.array sigma: Conductivity model
:param numpy.array m: Conductivity model
:param simpegEM.TDEM.FieldsTDEM vec: Fields object
:rtype: simpegEM.TDEM.FieldsTDEM
:return: f
@@ -167,7 +167,7 @@ class ProblemTDEM_b(ProblemBaseTDEM):
\\right] \\\\
"""
self.makeMassMatrices(sigma)
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)
@@ -182,9 +182,9 @@ class ProblemTDEM_b(ProblemBaseTDEM):
f.set_e(e, i)
return f
def AhtVec(self, sigma, vec):
def AhtVec(self, m, vec):
"""
:param numpy.array sigma: Conductivity model
:param numpy.array m: Conductivity model
:param simpegEM.TDEM.FieldsTDEM vec: Fields object
:rtype: simpegEM.TDEM.FieldsTDEM
:return: f
@@ -215,7 +215,7 @@ class ProblemTDEM_b(ProblemBaseTDEM):
\end{array}
\\right] \\\\
"""
self.makeMassMatrices(sigma)
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)
@@ -257,9 +257,9 @@ if __name__ == '__main__':
# prb.setTimes([1e-5, 5e-5, 2.5e-4], [10, 10, 10])
prb.timeSteps = [(1e-5, 10)]
prb.pair(survey)
sigma = np.random.rand(mesh.nCz)
m = np.random.rand(mesh.nCz)
print survey.dpred(sigma)
print survey.dpred(m)
simpegEM/Tests/test_TDEM_b_DerivAdjoint.py
+5 -4
View File
@@ -17,7 +17,7 @@ class TDEM_bDerivTests(unittest.TestCase):
mesh = Mesh.CylMesh([hx,1,hy], '00C')
active = mesh.vectorCCz<0.
activeMap = Maps.ActiveCells(mesh, active, -8, nC=mesh.nCz)
activeMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ComboMap(mesh,
[Maps.ExpMap, Maps.Vertical1DMap, activeMap])
@@ -76,7 +76,7 @@ class TDEM_bDerivTests(unittest.TestCase):
prb = self.prb
prb.timeSteps = [1e-05]
sigma = self.sigma
prb.makeMassMatrices(sigma)
prb.curModel = sigma
dt = prb.timeSteps[0]
a11 = 1/dt*prb.MfMui*sp.eye(prb.mesh.nF)
@@ -97,7 +97,7 @@ class TDEM_bDerivTests(unittest.TestCase):
prb.timeSteps = [1e-05]
sigma = self.sigma
prb.makeMassMatrices(sigma)
prb.curModel = sigma
dt = prb.timeSteps[0]
a11 = 1/dt*prb.MfMui*sp.eye(prb.mesh.nF)
@@ -120,7 +120,7 @@ class TDEM_bDerivTests(unittest.TestCase):
prb = self.prb
mesh = self.prb.mesh
sigma = self.sigma
self.prb.makeMassMatrices(sigma)
self.prb.curModel = sigma
f = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b')
for i in range(f.nT):
@@ -242,6 +242,7 @@ class TDEM_bDerivTests(unittest.TestCase):
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)
def test_adjointsolveAhVssolveAht(self):
simpegEM/__init__.py
+1
View File
@@ -2,3 +2,4 @@
import Utils
import TDEM
import FDEM
import Base