From d27e7d3abab38c690859627ff8a5a55733a07e72 Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Sat, 15 Feb 2014 14:16:13 -0800 Subject: [PATCH 1/4] Simplified fields method. Changed solveAh, AhVec & tests to be consistent with notes with regards to MfMui --- simpegEM/TDEM/BaseTDEM.py | 9 +- simpegEM/TDEM/TDEM_b.py | 9 +- simpegEM/Tests/test_forward_EMproblem.py | 135 +++++++++++++---------- 3 files changed, 84 insertions(+), 69 deletions(-) diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py index d0e7f992..75241c37 100644 --- a/simpegEM/TDEM/BaseTDEM.py +++ b/simpegEM/TDEM/BaseTDEM.py @@ -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): diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index a72dda53..7f00485b 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -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,7 +75,8 @@ 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) + self.makeMassMatrices(m) + Y = self.forward(m, AhRHS, AhCalcFields) return Y #################################################### @@ -87,14 +88,14 @@ 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): 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) diff --git a/simpegEM/Tests/test_forward_EMproblem.py b/simpegEM/Tests/test_forward_EMproblem.py index 1334a632..5ec65518 100644 --- a/simpegEM/Tests/test_forward_EMproblem.py +++ b/simpegEM/Tests/test_forward_EMproblem.py @@ -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): @@ -75,56 +75,75 @@ class TDEM_bDerivTests(unittest.TestCase): 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 +169,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 +199,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,6 +229,7 @@ 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) From 1128fbe39b111c202307dbb76129a132e27f6c17 Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Sat, 15 Feb 2014 14:16:48 -0800 Subject: [PATCH 2/4] Added adjoint test for projection. --- simpegEM/Tests/test_forward_EMproblem.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/simpegEM/Tests/test_forward_EMproblem.py b/simpegEM/Tests/test_forward_EMproblem.py index 5ec65518..1f5750d8 100644 --- a/simpegEM/Tests/test_forward_EMproblem.py +++ b/simpegEM/Tests/test_forward_EMproblem.py @@ -234,6 +234,23 @@ class TDEM_bDerivTests(unittest.TestCase): 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.assertTrue((V1-V2)/V1<1e-12) if __name__ == '__main__': From 7bb70103aa5857785740c49f418f391fbccdf15c Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Tue, 18 Feb 2014 15:49:54 -0800 Subject: [PATCH 3/4] Added AhtVec method & test. --- simpegEM/TDEM/TDEM_b.py | 21 ++++++++++++++++++++- simpegEM/Tests/test_forward_EMproblem.py | 21 +++++++++++++++++++-- 2 files changed, 39 insertions(+), 3 deletions(-) diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 7f00485b..2eada637 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -93,7 +93,7 @@ class ProblemTDEM_b(ProblemBaseTDEM): 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*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) @@ -101,6 +101,25 @@ class ProblemTDEM_b(ProblemBaseTDEM): 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 diff --git a/simpegEM/Tests/test_forward_EMproblem.py b/simpegEM/Tests/test_forward_EMproblem.py index 1f5750d8..42ad57dd 100644 --- a/simpegEM/Tests/test_forward_EMproblem.py +++ b/simpegEM/Tests/test_forward_EMproblem.py @@ -69,7 +69,6 @@ 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 @@ -233,7 +232,6 @@ class TDEM_bDerivTests(unittest.TestCase): 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 @@ -252,6 +250,25 @@ class TDEM_bDerivTests(unittest.TestCase): self.assertTrue((V1-V2)/V1<1e-12) + 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)/V1, 1e-12) + if __name__ == '__main__': unittest.main() From 915667ad4ab450fb91ceea77864535bbc863a2be Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Tue, 18 Feb 2014 16:27:09 -0800 Subject: [PATCH 4/4] solvent method and tests. --- simpegEM/TDEM/TDEM_b.py | 20 +++++++++-- simpegEM/Tests/test_forward_EMproblem.py | 42 ++++++++++++++++++++++-- 2 files changed, 58 insertions(+), 4 deletions(-) diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 2eada637..35aaf5c9 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -76,8 +76,24 @@ class ProblemTDEM_b(ProblemBaseTDEM): return {'b':b, 'e':e} self.makeMassMatrices(m) - Y = self.forward(m, AhRHS, AhCalcFields) - return Y + 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 diff --git a/simpegEM/Tests/test_forward_EMproblem.py b/simpegEM/Tests/test_forward_EMproblem.py index 42ad57dd..95dcb72f 100644 --- a/simpegEM/Tests/test_forward_EMproblem.py +++ b/simpegEM/Tests/test_forward_EMproblem.py @@ -248,7 +248,7 @@ class TDEM_bDerivTests(unittest.TestCase): V1 = d.T.dot(dat.projectFields(f)) V2 = f.fieldVec().dot(dat.projectFieldsAdjoint(d).fieldVec()) - self.assertTrue((V1-V2)/V1<1e-12) + self.assertLess((V1-V2)/np.abs(V1), 1e-8) def test_adjointAhVsAht(self): prb = self.prb @@ -267,7 +267,45 @@ class TDEM_bDerivTests(unittest.TestCase): V1 = f2.fieldVec().dot(prb.AhVec(sigma, f1).fieldVec()) V2 = f1.fieldVec().dot(prb.AhtVec(sigma, f2).fieldVec()) - self.assertLess(np.abs(V1-V2)/V1, 1e-12) + 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__':