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

Merge branch 'master' of https://github.com/simpeg/simpegem

Browse Source
This commit is contained in:
rowanc1
2014-02-18 17:57:04 -08:00
parent bf2021d8e8 915667ad4a
commit 7195adc6c5
3 changed files with 194 additions and 72 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/TDEM/BaseTDEM.py
+2 -7
View File
@@ -135,15 +135,10 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
solveOpts = {'factorize':True,'backend':'scipy'}
def fields(self, m, useThisRhs=None, useThisCalcFields=None):
RHS = useThisRhs or self.getRHS
CalcFields = useThisCalcFields or self.calcFields
def fields(self, m):
self.makeMassMatrices(m)
F = self.getInitialFields()
return self.forward(m, RHS, CalcFields, F=F)
return self.forward(m, self.getRHS, self.calcFields, F=F)
def forward(self, m, RHS, CalcFields, F=None):
simpegEM/TDEM/TDEM_b.py
+42 -6
View File
@@ -64,7 +64,7 @@ class ProblemTDEM_b(ProblemBaseTDEM):
def solveAh(self, m, p):
def AhRHS(tInd, u):
rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p.get_e(tInd) + self.MfMui*p.get_b(tInd)
rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p.get_e(tInd) + p.get_b(tInd)
if tInd == 0:
return rhs
dt = self.getDt(tInd)
@@ -75,8 +75,25 @@ 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}
Y = self.fields(m, useThisRhs=AhRHS, useThisCalcFields=AhCalcFields)
return Y
self.makeMassMatrices(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)
if tInd == self.nTimes-1:
return rhs
dt = self.getDt(tInd+1)
return rhs + 1./dt*self.MfMui*u.get_b(tInd+1)
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}
self.makeMassMatrices(m)
return self.adjoint(m, AhtRHS, AhtCalcFields)
####################################################
# Functions for tests
@@ -87,19 +104,38 @@ class ProblemTDEM_b(ProblemBaseTDEM):
u = self.fields(m)
self.makeMassMatrices(m)
dt = self.getDt(0)
b = 1/dt*u.get_b(0) + self.mesh.edgeCurl*u.get_e(0)
b = 1/dt*self.MfMui*u.get_b(0) + self.MfMui*self.mesh.edgeCurl*u.get_e(0)
e = self.mesh.edgeCurl.T*self.MfMui*u.get_b(0) - self.MeSigma*u.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.times.size):
for i in range(1,self.nTimes):
dt = self.getDt(i)
b = 1/dt*u.get_b(i) + self.mesh.edgeCurl*u.get_e(i) - 1/dt*u.get_b(i-1)
b = 1/dt*self.MfMui*u.get_b(i) + self.MfMui*self.mesh.edgeCurl*u.get_e(i) - 1/dt*self.MfMui*u.get_b(i-1)
e = self.mesh.edgeCurl.T*self.MfMui*u.get_b(i) - self.MeSigma*u.get_e(i)
f.set_b(b, i)
f.set_e(e, i)
return f
def AhtVec(self, m, u=None):
if u is None:
u = self.fields(m)
self.makeMassMatrices(m)
f = FieldsTDEM(self.mesh, 1, self.times.size, 'b')
for i in range(self.nTimes-1):
b = 1/self.getDt(i)*self.MfMui*u.get_b(i) + self.MfMui*self.mesh.edgeCurl*u.get_e(i) - 1/self.getDt(i+1)*self.MfMui*u.get_b(i+1)
e = self.mesh.edgeCurl.T*self.MfMui*u.get_b(i) - self.MeSigma*u.get_e(i)
f.set_b(b, i)
f.set_e(e, i)
N = self.nTimes - 1
b = 1/self.getDt(N)*self.MfMui*u.get_b(N) + self.MfMui*self.mesh.edgeCurl*u.get_e(N)
e = self.mesh.edgeCurl.T*self.MfMui*u.get_b(N) - self.MeSigma*u.get_e(N)
f.set_b(b, N)
f.set_e(e, N)
return f
if __name__ == '__main__':
from SimPEG import *
import simpegEM as EM
simpegEM/Tests/test_forward_EMproblem.py
+150 -59
View File
@@ -10,35 +10,35 @@ class TDEM_bTests(unittest.TestCase):
def setUp(self):
cs = 5.
ncx = 20
ncy = 6
npad = 20
hx = Utils.meshTensors(((0,cs), (ncx,cs), (npad,cs)))
hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs)))
mesh = Mesh.Cyl1DMesh([hx,hy], -hy.sum()/2)
model = Model.Vertical1DModel(mesh)
cs = 5.
ncx = 20
ncy = 6
npad = 20
hx = Utils.meshTensors(((0,cs), (ncx,cs), (npad,cs)))
hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs)))
mesh = Mesh.Cyl1DMesh([hx,hy], -hy.sum()/2)
model = Model.Vertical1DModel(mesh)
opts = {'txLoc':0.,
'txType':'VMD_MVP',
'rxLoc':np.r_[150., 0.],
'rxType':'bz',
'timeCh':np.logspace(-4,-2,20),
}
self.dat = EM.TDEM.DataTDEM1D(**opts)
opts = {'txLoc':0.,
'txType':'VMD_MVP',
'rxLoc':np.r_[150., 0.],
'rxType':'bz',
'timeCh':np.logspace(-4,-2,20),
}
self.dat = EM.TDEM.DataTDEM1D(**opts)
self.prb = EM.TDEM.ProblemTDEM_b(mesh, model)
self.prb.setTimes([1e-5, 5e-5, 2.5e-4], [150, 150, 150])
self.sigma = np.ones(mesh.nCz)*1e-8
self.sigma[mesh.vectorCCz<0] = 0.1
self.prb.pair(self.dat)
self.prb = EM.TDEM.ProblemTDEM_b(mesh, model)
self.prb.setTimes([1e-5, 5e-5, 2.5e-4], [150, 150, 150])
self.sigma = np.ones(mesh.nCz)*1e-8
self.sigma[mesh.vectorCCz<0] = 0.1
self.prb.pair(self.dat)
def test_analitic_b(self):
bz_calc = self.dat.dpred(self.sigma)
bz_ana = mu_0*hzAnalyticDipoleT(self.dat.rxLoc[0], self.prb.times, self.sigma[0])
bz_calc = self.dat.dpred(self.sigma)
bz_ana = mu_0*hzAnalyticDipoleT(self.dat.rxLoc[0], self.prb.times, self.sigma[0])
diff = np.linalg.norm(bz_calc.flatten() - bz_ana.flatten())/np.linalg.norm(bz_ana.flatten())
self.assertTrue(diff<0.05)
diff = np.linalg.norm(bz_calc.flatten() - bz_ana.flatten())/np.linalg.norm(bz_ana.flatten())
self.assertTrue(diff<0.05)
class TDEM_bDerivTests(unittest.TestCase):
@@ -69,62 +69,80 @@ class TDEM_bDerivTests(unittest.TestCase):
self.prb.pair(self.dat)
self.mesh = mesh
def test_AhVec(self):
"""
Test that fields and AhVec produce consistent results
"""
prb = self.prb
sigma = np.ones(self.prb.mesh.nCz)*1e-8
sigma[self.prb.mesh.vectorCCz<0] = 0.1
u = self.prb.fields(sigma)
Ahu = self.prb.AhVec(sigma, u)
self.assertTrue(np.linalg.norm(Ahu.get_b(0)-1/self.prb.getDt(0)*u.get_b(-1))/np.linalg.norm(u.get_b(0)) < 1.e-2)
self.assertTrue(np.linalg.norm(Ahu.get_e(0))/np.linalg.norm(u.get_e(0)) < 1.e-2)
sigma[prb.mesh.vectorCCz<0] = 0.1
u = prb.fields(sigma)
Ahu = prb.AhVec(sigma, u)
V1 = Ahu.get_b(0)
V2 = 1/prb.getDt(0)*prb.MfMui*u.get_b(-1)
self.assertTrue(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)
for i in range(1,u.nTimes):
self.assertTrue(np.linalg.norm(Ahu.get_b(i))/np.linalg.norm(u.get_b(i)) < 1.e-2)
self.assertTrue(np.linalg.norm(Ahu.get_e(i))/np.linalg.norm(u.get_e(i)) < 1.e-2)
dt = prb.getDt(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.get_e(i)
V2 = prb.MeSigma*u.get_e(i)
self.assertTrue(np.linalg.norm(V1)/np.linalg.norm(V2) < 1.e-6)
def test_AhVecVSMat_OneTS(self):
self.prb.setTimes([1e-5], [1])
prb = self.prb
prb.setTimes([1e-5], [1])
sigma = np.ones(self.prb.mesh.nCz)*1e-8
sigma[self.prb.mesh.vectorCCz<0] = 0.1
self.prb.makeMassMatrices(sigma)
sigma = np.ones(prb.mesh.nCz)*1e-8
sigma[prb.mesh.vectorCCz<0] = 0.1
prb.makeMassMatrices(sigma)
dt = self.prb.getDt(0)
a11 = 1/dt*sp.eye(self.prb.mesh.nF)
a12 = self.prb.mesh.edgeCurl
a21 = self.prb.mesh.edgeCurl.T*self.prb.MfMui
a22 = -self.prb.MeSigma
dt = prb.getDt(0)
a11 = 1/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 = self.prb.fields(sigma)
f = prb.fields(sigma)
u1 = A*f.fieldVec()
u2 = self.prb.AhVec(sigma,f).fieldVec()
u2 = prb.AhVec(sigma,f).fieldVec()
self.assertTrue(np.linalg.norm(u1-u2)<1e-12)
self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12)
def test_solveAhVSMat_OneTS(self):
self.prb.setTimes([1e-5], [1])
prb = self.prb
sigma = np.ones(self.prb.mesh.nCz)*1e-8
sigma[self.prb.mesh.vectorCCz<0] = 0.1
self.prb.makeMassMatrices(sigma)
prb.setTimes([1e-5], [1])
dt = self.prb.getDt(0)
a11 = 1/dt*sp.eye(self.prb.mesh.nF)
a12 = self.prb.mesh.edgeCurl
a21 = self.prb.mesh.edgeCurl.T*self.prb.MfMui
a22 = -self.prb.MeSigma
sigma = np.ones(prb.mesh.nCz)*1e-8
sigma[prb.mesh.vectorCCz<0] = 0.1
prb.makeMassMatrices(sigma)
dt = prb.getDt(0)
a11 = 1/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 = self.prb.fields(sigma)
f.set_b(np.zeros((self.prb.mesh.nF,1)),0)
f.set_e(np.random.rand(self.prb.mesh.nE,1),0)
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)
u1 = self.prb.solveAh(sigma,f).fieldVec().flatten()
u1 = prb.solveAh(sigma,f).fieldVec().flatten()
u2 = sp.linalg.spsolve(A.tocsr(),f.fieldVec())
self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
@@ -150,8 +168,6 @@ class TDEM_bDerivTests(unittest.TestCase):
u2 = f_test.fieldVec()
self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
def test_DerivG(self):
"""
Test the derivative of c with respect to sigma
@@ -182,6 +198,7 @@ class TDEM_bDerivTests(unittest.TestCase):
error = np.zeros(num)
order = 0
hv = np.logspace(-1.2,-3, num)
print '\n'
for i, h in enumerate(hv):
f = prb.fields(sigma)
fstep = prb.fields(sigma + h*d_sig)
@@ -211,10 +228,84 @@ class TDEM_bDerivTests(unittest.TestCase):
derChk = lambda m: [prb.data.dpred(m), lambda mx: -prb.Jvec(sigma, mx)]
print '\n'
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=2, eps=1e-20)
self.assertTrue(passed)
def test_projectAdjoint(self):
prb = self.prb
dat = self.dat
mesh = self.mesh
# Generate random fields and data
f = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.times.size, 'b')
for i in range(f.nTimes):
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.nTimes, dat.nTx)
# 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())
self.assertLess((V1-V2)/np.abs(V1), 1e-8)
def test_adjointAhVsAht(self):
prb = self.prb
mesh = self.mesh
sigma = self.sigma
f1 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b')
for i in range(f1.nTimes):
f1.set_b(np.random.rand(mesh.nF, 1), i)
f1.set_e(np.random.rand(mesh.nE, 1), i)
f2 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b')
for i in range(f2.nTimes):
f2.set_b(np.random.rand(mesh.nF, 1), i)
f2.set_e(np.random.rand(mesh.nE, 1), i)
V1 = f2.fieldVec().dot(prb.AhVec(sigma, f1).fieldVec())
V2 = f1.fieldVec().dot(prb.AhtVec(sigma, f2).fieldVec())
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-8)
def test_solveAhtVsAhtVec(self):
prb = self.prb
mesh = self.mesh
sigma = np.random.rand(mesh.nCz)
f1 = EM.TDEM.FieldsTDEM(mesh, 1, prb.nTimes, 'b')
for i in range(f1.nTimes):
f1.set_b(np.random.rand(mesh.nF, 1), i)
f1.set_e(np.random.rand(mesh.nE, 1), i)
f2 = prb.solveAht(sigma, f1)
f3 = prb.AhtVec(sigma, f2)
V1 = np.linalg.norm(f3.fieldVec()-f1.fieldVec())
V2 = np.linalg.norm(f1.fieldVec())
self.assertLess(V1/V2, 1e-8)
def test_adjointsolveAhVssolveAht(self):
prb = self.prb
mesh = self.mesh
sigma = self.sigma
f1 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b')
for i in range(f1.nTimes):
f1.set_b(np.random.rand(mesh.nF, 1), i)
f1.set_e(np.random.rand(mesh.nE, 1), i)
f2 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b')
for i in range(f2.nTimes):
f2.set_b(np.random.rand(mesh.nF, 1), i)
f2.set_e(np.random.rand(mesh.nE, 1), i)
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-8)
if __name__ == '__main__':