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

Minor Memory improvements and speedups.

Browse Source 
Speedups are from brackets mat-vec prods vs mat-mat
Memory is from not setting all 'b' fields to zero in Gvec
This commit is contained in:
rowanc1
2014-05-15 08:52:30 -07:00
parent 06b34d8f7b
commit d58cb21123
5 changed files with 41 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/TDEM/BaseTDEM.py
+4 -2
View File
@@ -88,7 +88,8 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
"""
self.curModel = m
u = u or self.fields(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)
@@ -110,7 +111,8 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
"""
self.curModel = m
u = u or self.fields(m)
if u is None:
u = self.fields(m)
if not isinstance(v, self.dataPair):
v = self.dataPair(self.survey, v)
simpegEM/TDEM/TDEM_b.py
+24 -15
View File
@@ -45,7 +45,7 @@ class ProblemTDEM_b(BaseTDEMProblem):
if self.solType == 'b':
b = sol
e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b
e = self.MeSigmaI*(self.mesh.edgeCurl.T*(self.MfMui*b))
# Todo: implement non-zero js
else:
raise NotImplementedError('solType "%s" is not implemented in CalcFields.' % self.solType)
@@ -67,21 +67,27 @@ class ProblemTDEM_b(BaseTDEMProblem):
Multiply G by a vector
"""
u = u or self.fields(m)
if u is None:
u = self.fields(m)
self.curModel = m
# Note: Fields has shape (nF/E, nTx, nT+1)
# However, p will only really fill (:,:,1:nT+1)
# meaning the 'initial fields' are zero (:,:,0)
p = FieldsTDEM(self.mesh, self.survey)
p[:, 'b', :] = 0.0 # b at all times is zero.
p[:, 'e', 0] = 0.0 # fake initial fields
c = self.mesh.getEdgeInnerProductDeriv(self.curTModel)*self.curTModelDeriv*vec
# 'b' at all times is zero.
# However, to save memory we will **not** do:
#
# p[:, 'b', :] = 0.0
# fake initial 'e' fields
p[:, 'e', 0] = 0.0
c = self.mesh.getEdgeInnerProductDeriv(self.curTModel)*(self.curTModelDeriv*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)
for tx in self.survey.txList:
p[tx, 'e', i] = -u[tx,'e',i]*c # - diag(e) * MsigDeriv * v
p[tx, 'e', i] = -u[tx,'e',i]*c # i.e.: - diag(e) * MsigDeriv * v
return p
def Gtvec(self, m, vec, u=None):
@@ -94,7 +100,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
Multiply G.T by a vector
"""
u = u or self.fields(m)
if u is None:
u = self.fields(m)
self.curModel = m
nTx, nE = self.survey.nTx, self.mesh.nE
@@ -143,7 +150,9 @@ class ProblemTDEM_b(BaseTDEMProblem):
"""
def AhRHS(tInd, y):
rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd+1] + p[:,'b',tInd+1]
rhs = self.MfMui*(self.mesh.edgeCurl*(self.MeSigmaI*p[:,'e',tInd+1]))
if 'b' in p:
rhs = rhs + p[:,'b',tInd+1]
if tInd == 0:
return rhs
dt = self.timeSteps[tInd]
@@ -153,7 +162,7 @@ class ProblemTDEM_b(BaseTDEMProblem):
y_b = sol
if self.survey.nTx == 1:
y_b = mkvc(y_b)
y_e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*y_b - self.MeSigmaI*p[:,'e',tInd+1]
y_e = self.MeSigmaI*(self.mesh.edgeCurl.T*(self.MfMui*y_b)) - self.MeSigmaI*p[:,'e',tInd+1]
return {'b':y_b, 'e':y_e}
return self.forward(m, AhRHS, AhCalcFields)
@@ -211,7 +220,7 @@ class ProblemTDEM_b(BaseTDEMProblem):
rhs = np.zeros(nF if nTx == 1 else (nF, nTx))
if 'e' in p:
rhs += self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd+1]
rhs += self.MfMui*(self.mesh.edgeCurl*(self.MeSigmaI*p[:,'e',tInd+1]))
if 'b' in p:
rhs += p[:,'b',tInd+1]
@@ -224,7 +233,7 @@ class ProblemTDEM_b(BaseTDEMProblem):
y_b = sol
if self.survey.nTx == 1:
y_b = mkvc(y_b)
y_e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*y_b
y_e = self.MeSigmaI*(self.mesh.edgeCurl.T*(self.MfMui*y_b))
if 'e' in p:
y_e += - self.MeSigmaI*p[:,'e',tInd+1]
return {'b':y_b, 'e':y_e}
@@ -273,11 +282,11 @@ class ProblemTDEM_b(BaseTDEMProblem):
f = FieldsTDEM(self.mesh, self.survey)
for i in range(1,self.nT+1):
dt = self.timeSteps[i-1]
b = 1.0/dt*self.MfMui*vec[:,'b',i] + self.MfMui*self.mesh.edgeCurl*vec[:,'e',i]
b = 1.0/dt*self.MfMui*vec[:,'b',i] + self.MfMui*(self.mesh.edgeCurl*vec[:,'e',i])
if i > 1:
b = b - 1.0/dt*self.MfMui*vec[:,'b',i-1]
f[:,'b',i] = b
f[:,'e',i] = self.mesh.edgeCurl.T*self.MfMui*vec[:,'b',i] - self.MeSigma*vec[:,'e',i]
f[:,'e',i] = self.mesh.edgeCurl.T*(self.MfMui*vec[:,'b',i]) - self.MeSigma*vec[:,'e',i]
return f
def _AhtVec(self, m, vec):
@@ -316,9 +325,9 @@ class ProblemTDEM_b(BaseTDEMProblem):
self.curModel = m
f = FieldsTDEM(self.mesh, self.survey)
for i in range(self.nT):
b = 1.0/self.timeSteps[i]*self.MfMui*vec[:,'b',i+1] + self.MfMui*self.mesh.edgeCurl*vec[:,'e',i+1]
b = 1.0/self.timeSteps[i]*self.MfMui*vec[:,'b',i+1] + self.MfMui*(self.mesh.edgeCurl*vec[:,'e',i+1])
if i < self.nT-1:
b = b - 1.0/self.timeSteps[i+1]*self.MfMui*vec[:,'b',i+2]
f[:,'b', i+1] = b
f[:,'e', i+1] = self.mesh.edgeCurl.T*self.MfMui*vec[:,'b',i+1] - self.MeSigma*vec[:,'e',i+1]
f[:,'e', i+1] = self.mesh.edgeCurl.T*(self.MfMui*vec[:,'b',i+1]) - self.MeSigma*vec[:,'e',i+1]
return f
simpegEM/Tests/test_FDEM.py
+1 -1
View File
@@ -33,7 +33,7 @@ def getProblem(fdemType, comp):
prb.pair(survey)
try:
from mumpsSCI import MumpsSolver
from pymatsolver import MumpsSolver
prb.Solver = MumpsSolver
except ImportError, e:
pass
simpegEM/Tests/test_TDEM_combos.py
+6
View File
@@ -2,6 +2,11 @@ import unittest
from SimPEG import *
import simpegEM as EM
try:
from pymatsolver import MumpsSolver
except ImportError, e:
MumpsSolver = Utils.SolverUtils.DSolverWrap(sp.linalg.splu, factorize=True)
plotIt = False
def getProb(meshType='CYL',rxTypes='bx,bz',nTx=1):
@@ -30,6 +35,7 @@ def getProb(meshType='CYL',rxTypes='bx,bz',nTx=1):
# prb.timeSteps = [1e-5]
prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
# prb.timeSteps = [(1e-05, 100)]
prb.Solver = MumpsSolver
sigma = np.ones(mesh.nCz)*1e-8
sigma[mesh.vectorCCz<0] = 1e-1
simpegEM/Tests/test_TDEM_forward_Analytic.py
+6 -1
View File
@@ -4,6 +4,11 @@ import simpegEM as EM
from scipy.constants import mu_0
import matplotlib.pyplot as plt
try:
from pymatsolver import MumpsSolver
except ImportError, e:
MumpsSolver = Utils.SolverUtils.DSolverWrap(sp.linalg.splu, factorize=True)
def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,1e-3], showIt=False):
if meshType == 'CYL':
@@ -27,7 +32,7 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,
survey = EM.TDEM.SurveyTDEM([tx])
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
prb.Solver = Utils.SolverUtils.DSolverWrap(sp.linalg.splu, factorize=True)
prb.Solver = MumpsSolver
prb.timeSteps = [(1e-06, 40), (5e-06, 40), (1e-05, 40), (5e-05, 40), (0.0001, 40), (0.0005, 40)]