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

Switch to BaseTimeProblem

Browse Source
This commit is contained in:
rowanc1
2014-04-24 16:50:58 -07:00
parent d9e91a5a1d
commit 200a1680ca
6 changed files with 71 additions and 115 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/TDEM/BaseTDEM.py
+11 -55
View File
@@ -1,9 +1,10 @@
from SimPEG import Solver
from SimPEG.Problem import BaseProblem
from SimPEG.Problem import BaseTimeProblem
from simpegEM.Utils import Sources
from FieldsTDEM import FieldsTDEM
from scipy.constants import mu_0
from SimPEG.Utils import sdiag, mkvc
from SimPEG import Utils, Mesh
import numpy as np
@@ -37,63 +38,18 @@ class MixinInitialFieldCalc(object):
raise Exception('Unknown mesh for VMD')
# Initialize field object
F = FieldsTDEM(self.mesh, 1, self.times.size, store=self.storeTheseFields)
F = FieldsTDEM(self.mesh, 1, self.times[1:].size, store=self.storeTheseFields)
# Set initial B
F.b0 = self.mesh.edgeCurl*MVP
return F
class MixinTimeStuff(object):
"""docstring for MixinTimeStuff"""
def dt():
doc = "Size of time steps"
def fget(self):
return self._dt
def fdel(self):
del self._dt
return locals()
dt = property(**dt())
def nsteps():
doc = "Number of steps to take"
def fget(self):
return self._nsteps
def fdel(self):
del self._nsteps
return locals()
nsteps = property(**nsteps())
def times():
doc = "Modeling times"
def fget(self):
t = np.r_[1:self.nsteps[0]+1]*self.dt[0]
for i in range(1,self.dt.size):
t = np.r_[t, np.r_[1:self.nsteps[i]+1]*self.dt[i]+t[-1]]
return t
return locals()
times = property(**times())
def getDt(self, tInd):
return np.concatenate([self.dt[i].repeat(self.nsteps[i]) for i in range(self.dt.size)])[tInd]
def setTimes(self, dt, nsteps):
dt = np.array(dt)
nsteps = np.array(nsteps)
assert dt.size==nsteps.size, "dt, nsteps must be same length"
self._dt = dt
self._nsteps = nsteps
@property
def nTimes(self):
return self.times.size
class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem):
"""docstring for ProblemTDEM1D"""
def __init__(self, mesh, mapping=None, **kwargs):
BaseProblem.__init__(self, mesh, mapping=mapping, **kwargs)
BaseTimeProblem.__init__(self, mesh, mapping=mapping, **kwargs)
####################################################
@@ -151,11 +107,11 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
def forward(self, m, RHS, CalcFields, F=None):
if F is None:
F = FieldsTDEM(self.mesh, self.survey.nTx, self.nTimes, store=self.storeTheseFields)
F = FieldsTDEM(self.mesh, self.survey.nTx, self.nT, store=self.storeTheseFields)
dtFact = None
for tInd, t in enumerate(self.times):
dt = self.getDt(tInd)
for tInd, t in enumerate(self.times[1:]):
dt = self.timeSteps[tInd]
if dt!=dtFact:
dtFact = dt
A = self.getA(tInd)
@@ -172,11 +128,11 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
def adjoint(self, m, RHS, CalcFields, F=None):
if F is None:
F = FieldsTDEM(self.mesh, self.survey.nTx, self.nTimes, store=self.storeTheseFields)
F = FieldsTDEM(self.mesh, self.survey.nTx, self.nT, store=self.storeTheseFields)
dtFact = None
for tInd, t in reversed(list(enumerate(self.times))):
dt = self.getDt(tInd)
for tInd, t in reversed(list(enumerate(self.times[1:]))):
dt = self.timeSteps[tInd]
if dt!=dtFact:
dtFact = dt
A = self.getA(tInd)
simpegEM/TDEM/FieldsTDEM.py
+5 -5
View File
@@ -18,9 +18,9 @@ class FieldsTDEM(object):
j = None #: Current density
h = None #: Magnetic field
def __init__(self, mesh, nTx, nTimes, store='b'):
def __init__(self, mesh, nTx, nT, store='b'):
self.nTimes = nTimes #: Number of times
self.nT = nT #: Number of times
self.nTx = nTx #: Number of transmitters
self.mesh = mesh
@@ -30,7 +30,7 @@ class FieldsTDEM(object):
def fieldVec(self):
u = np.ndarray((0, self.nTx))
for i in range(self.nTimes):
for i in range(self.nT):
u = np.r_[u, self.get_b(i), self.get_e(i)]
if self.nTx == 1:
u = u.flatten()
@@ -58,7 +58,7 @@ class FieldsTDEM(object):
def set_b(self, b, ind):
if self.b is None:
self.b = np.zeros((self.nTimes, np.sum(self.mesh.nF), self.nTx))
self.b = np.zeros((self.nT, np.sum(self.mesh.nF), self.nTx))
self.b[:] = np.nan
if len(b.shape) == 1:
b = b[:, np.newaxis]
@@ -66,7 +66,7 @@ class FieldsTDEM(object):
def set_e(self, e, ind):
if self.e is None:
self.e = np.zeros((self.nTimes, np.sum(self.mesh.nE), self.nTx))
self.e = np.zeros((self.nT, np.sum(self.mesh.nE), self.nTx))
self.e[:] = np.nan
if len(e.shape) == 1:
e = e[:, np.newaxis]
simpegEM/TDEM/SurveyTDEM.py
+3 -3
View File
@@ -29,11 +29,11 @@ class SurveyTDEM1D(BaseSurvey):
def projectFieldsAdjoint(self, d):
# TODO: make the following self.nTimeCh
d = d.reshape((self.prob.nTimes, self.nTx), order='F')
d = d.reshape((self.prob.nT, self.nTx), order='F')
#TODO: *Qtime.T need to multiply by a time projection. (outside for loop??)
ii = 0
F = FieldsTDEM(self.prob.mesh, self.nTx, self.prob.nTimes, 'b')
for ii in range(self.prob.nTimes):
F = FieldsTDEM(self.prob.mesh, self.nTx, self.prob.nT, 'b')
for ii in range(self.prob.nT):
b = self.Qrx.T*d[ii,:]
F.set_b(b, ii)
F.set_e(np.zeros((self.prob.mesh.nE,self.nTx)), ii)
simpegEM/TDEM/TDEM_b.py
+17 -17
View File
@@ -34,11 +34,11 @@ class ProblemTDEM_b(ProblemBaseTDEM):
:return: A
"""
dt = self.getDt(tInd)
dt = self.timeSteps[tInd]
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)
dt = self.timeSteps[tInd]
return (1.0/dt)*self.MfMui*F.get_b(tInd-1)
@@ -73,10 +73,10 @@ class ProblemTDEM_b(ProblemBaseTDEM):
"""
if u is None:
u = self.fields(m)
p = FieldsTDEM(self.mesh, 1, self.times.size, 'b')
p = FieldsTDEM(self.mesh, 1, self.nT, 'b')
curModel = self.mapping.transform(m)
c = self.mesh.getEdgeInnerProductDeriv(curModel)*self.mapping.transformDeriv(m)*vec
for i in range(self.times.size):
for i in range(self.nT):
ei = u.get_e(i)
pVal = np.empty_like(ei)
for j in range(ei.shape[1]):
@@ -90,7 +90,7 @@ class ProblemTDEM_b(ProblemBaseTDEM):
if u is None:
u = self.fields(m)
tmp = np.zeros((self.mesh.nE,self.survey.nTx))
for i in range(self.nTimes):
for i in range(self.nT):
tmp += v.get_e(i)*u.get_e(i)
curModel = self.mapping.transform(m)
@@ -102,7 +102,7 @@ class ProblemTDEM_b(ProblemBaseTDEM):
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)
dt = self.timeSteps[tInd]
return rhs + 1.0/dt*self.MfMui*u.get_b(tInd-1)
def AhCalcFields(sol, solType, tInd):
@@ -117,9 +117,9 @@ class ProblemTDEM_b(ProblemBaseTDEM):
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:
if tInd == self.nT-1:
return rhs
dt = self.getDt(tInd+1)
dt = self.timeSteps[tInd+1]
return rhs + 1.0/dt*self.MfMui*u.get_b(tInd+1)
def AhtCalcFields(sol, solType, tInd):
@@ -169,14 +169,14 @@ class ProblemTDEM_b(ProblemBaseTDEM):
"""
self.makeMassMatrices(sigma)
dt = self.getDt(0)
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)
f = FieldsTDEM(self.mesh, 1, self.times.size, 'b')
f = FieldsTDEM(self.mesh, 1, self.nT, 'b')
f.set_b(b, 0)
f.set_e(e, 0)
for i in range(1,self.nTimes):
dt = self.getDt(i)
for i in range(1,self.nT):
dt = self.timeSteps[i]
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)
@@ -217,14 +217,14 @@ class ProblemTDEM_b(ProblemBaseTDEM):
\\right] \\\\
"""
self.makeMassMatrices(sigma)
f = FieldsTDEM(self.mesh, 1, self.times.size, 'b')
for i in range(self.nTimes-1):
b = 1/self.getDt(i)*self.MfMui*vec.get_b(i) + self.MfMui*self.mesh.edgeCurl*vec.get_e(i) - 1/self.getDt(i+1)*self.MfMui*vec.get_b(i+1)
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)
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)
N = self.nTimes - 1
b = 1/self.getDt(N)*self.MfMui*vec.get_b(N) + self.MfMui*self.mesh.edgeCurl*vec.get_e(N)
N = self.nT - 1
b = 1.0/self.timeSteps[N]*self.MfMui*vec.get_b(N) + self.MfMui*self.mesh.edgeCurl*vec.get_e(N)
e = self.mesh.edgeCurl.T*self.MfMui*vec.get_b(N) - self.MeSigma*vec.get_e(N)
f.set_b(b, N)
f.set_e(e, N)
simpegEM/Tests/test_TDEM_b_DerivAdjoint.py
+30 -30
View File
@@ -29,7 +29,7 @@ class TDEM_bDerivTests(unittest.TestCase):
self.dat = EM.TDEM.SurveyTDEM1D(**opts)
self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
self.prb.setTimes([1e-5, 5e-5, 2.5e-4], [10, 10, 10])
self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (0.00025, 10)]
self.sigma = np.ones(mesh.nCz)*1e-8
self.sigma[mesh.vectorCCz<0] = 1e-1
@@ -50,16 +50,16 @@ 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.timeSteps[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)
self.assertTrue(np.linalg.norm(V1) < 1.e-6)
for i in range(1,u.nTimes):
for i in range(1,u.nT):
dt = prb.getDt(i)
dt = prb.timeSteps[i]
V1 = Ahu.get_b(i)
V2 = 1/dt*prb.MfMui*u.get_b(i-1)
@@ -72,11 +72,11 @@ class TDEM_bDerivTests(unittest.TestCase):
def test_AhVecVSMat_OneTS(self):
prb = self.prb
prb.setTimes([1e-5], [1])
prb.timeSteps = [(1e-05, 1)]
sigma = self.sigma
prb.makeMassMatrices(sigma)
dt = prb.getDt(0)
dt = prb.timeSteps[0]
a11 = 1/dt*prb.MfMui*sp.eye(prb.mesh.nF)
a12 = prb.MfMui*prb.mesh.edgeCurl
a21 = prb.mesh.edgeCurl.T*prb.MfMui
@@ -92,12 +92,12 @@ class TDEM_bDerivTests(unittest.TestCase):
def test_solveAhVSMat_OneTS(self):
prb = self.prb
prb.setTimes([1e-5], [1])
prb.timeSteps = [(1e-05, 1)]
sigma = self.sigma
prb.makeMassMatrices(sigma)
dt = prb.getDt(0)
dt = prb.timeSteps[0]
a11 = 1/dt*prb.MfMui*sp.eye(prb.mesh.nF)
a12 = prb.MfMui*prb.mesh.edgeCurl
a21 = prb.mesh.edgeCurl.T*prb.MfMui
@@ -120,8 +120,8 @@ class TDEM_bDerivTests(unittest.TestCase):
sigma = self.sigma
self.prb.makeMassMatrices(sigma)
f = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.times.size, 'b')
for i in range(f.nTimes):
f = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b')
for i in range(f.nT):
f.set_b(np.zeros((mesh.nF, 1)), i)
f.set_e(np.random.rand(mesh.nE, 1), i)
@@ -152,7 +152,7 @@ class TDEM_bDerivTests(unittest.TestCase):
def test_Deriv_dUdM(self):
prb = self.prb
prb.setTimes([1e-5, 1e-4, 1e-3], [10, 10, 10])
prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
mesh = self.mesh
sigma = self.sigma
@@ -168,7 +168,7 @@ class TDEM_bDerivTests(unittest.TestCase):
def test_Deriv_J(self):
prb = self.prb
prb.setTimes([1e-5, 1e-4, 1e-3], [10, 10, 10])
prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
mesh = self.mesh
sigma = self.sigma
@@ -188,11 +188,11 @@ class TDEM_bDerivTests(unittest.TestCase):
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 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b')
for i in range(f.nT):
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)
d = np.random.rand(dat.prob.nT, dat.nTx)
# Check that d.T*Q*f = f.T*Q.T*d
V1 = d.T.dot(dat.projectFields(f))
@@ -205,13 +205,13 @@ class TDEM_bDerivTests(unittest.TestCase):
mesh = self.mesh
sigma = self.sigma
f1 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b')
for i in range(f1.nTimes):
f1 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b')
for i in range(f1.nT):
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 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b')
for i in range(f2.nT):
f2.set_b(np.random.rand(mesh.nF, 1), i)
f2.set_e(np.random.rand(mesh.nE, 1), i)
@@ -224,8 +224,8 @@ class TDEM_bDerivTests(unittest.TestCase):
mesh = self.mesh
sigma = np.random.rand(prb.mapping.nP)
f1 = EM.TDEM.FieldsTDEM(mesh, 1, prb.nTimes, 'b')
for i in range(f1.nTimes):
f1 = EM.TDEM.FieldsTDEM(mesh, 1, prb.nT, 'b')
for i in range(f1.nT):
f1.set_b(np.random.rand(mesh.nF, 1), i)
f1.set_e(np.random.rand(mesh.nE, 1), i)
@@ -241,13 +241,13 @@ class TDEM_bDerivTests(unittest.TestCase):
mesh = self.mesh
sigma = self.sigma
f1 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b')
for i in range(f1.nTimes):
f1 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b')
for i in range(f1.nT):
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 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b')
for i in range(f2.nT):
f2.set_b(np.random.rand(mesh.nF, 1), i)
f2.set_e(np.random.rand(mesh.nE, 1), i)
@@ -262,13 +262,13 @@ class TDEM_bDerivTests(unittest.TestCase):
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):
u = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b')
for i in range(u.nT):
u.set_b(np.random.rand(mesh.nF, 1), i)
u.set_e(np.random.rand(mesh.nE, 1), i)
v = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b')
for i in range(v.nTimes):
v = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b')
for i in range(v.nT):
v.set_b(np.random.rand(mesh.nF, 1), i)
v.set_e(np.random.rand(mesh.nE, 1), i)
@@ -282,7 +282,7 @@ class TDEM_bDerivTests(unittest.TestCase):
sigma = self.sigma
m = np.random.rand(prb.mapping.nP)
d = np.random.rand(prb.nTimes)
d = np.random.rand(prb.nT)
V1 = d.dot(prb.Jvec(sigma, m))
V2 = m.dot(prb.Jtvec(sigma, d))
simpegEM/Tests/test_TDEM_forward_Analytic.py
+5 -5
View File
@@ -39,23 +39,23 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,
# except ImportError, e:
# pass
prb.setTimes([1e-6, 5e-6, 1e-5, 5e-5, 1e-4, 5e-4], [40, 40, 40, 40, 40, 40])
prb.timeSteps = [(1e-06, 40), (5e-06, 40), (1e-05, 40), (5e-05, 40), (0.0001, 40), (0.0005, 40)]
sigma = np.ones(mesh.nCz)*1e-8
sigma[active] = sig_half
sigma = np.log(sigma[active])
prb.pair(survey)
bz_ana = mu_0*EM.Utils.Ana.hzAnalyticDipoleT(survey.rxLoc[0], prb.times, sig_half)
bz_ana = mu_0*EM.Utils.Ana.hzAnalyticDipoleT(survey.rxLoc[0]+1e-3, prb.times[1:], sig_half)
bz_calc = survey.dpred(sigma)
ind = np.logical_and(prb.times > bounds[0],prb.times < bounds[1])
ind = np.logical_and(prb.times[1:] > bounds[0],prb.times[1:] < bounds[1])
log10diff = np.linalg.norm(np.log10(np.abs(bz_calc[ind])) - np.log10(np.abs(bz_ana[ind])))/np.linalg.norm(np.log10(np.abs(bz_ana[ind])))
print 'Difference: ', log10diff
if showIt == True:
plt.loglog(prb.times[bz_calc>0], bz_calc[bz_calc>0], 'r', prb.times[bz_calc<0], -bz_calc[bz_calc<0], 'r--')
plt.loglog(prb.times, abs(bz_ana), 'b*')
plt.loglog(prb.times[1:][bz_calc>0], bz_calc[bz_calc>0], 'r', prb.times[1:][bz_calc<0], -bz_calc[bz_calc<0], 'r--')
plt.loglog(prb.times[1:], abs(bz_ana), 'b*')
plt.title('sig_half = %e'%sig_half)
plt.show()