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

Merge branch 'master' into em/primsec

Browse Source 
# Conflicts:
#	SimPEG/EM/FDEM/SrcFDEM.py
This commit is contained in:
Lindsey Heagy
2016-02-10 16:15:01 -08:00
parent 436bc7af95 8ba9815137
commit 944fb1dfee
105 changed files with 4446 additions and 2910 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/DataMisfit.py
+11 -17
View File
@@ -59,20 +59,6 @@ class BaseDataMisfit(object):
"""
raise NotImplementedError('This method should be overwritten.')
# TODO: implement target misfit as a property, or possibly as an inversion directive.
# def target(self, forward):
# """target(forward)
# Target for data misfit. By default this is the number of data,
# which satisfies the Discrepancy Principle.
# :rtype: float
# :return: data misfit target
# """
# prob, survey = self.splitForward(forward)
# return survey.nD
class l2_DataMisfit(BaseDataMisfit):
@@ -103,10 +89,18 @@ class l2_DataMisfit(BaseDataMisfit):
"""
if getattr(self, '_Wd', None) is None:
print 'SimPEG.l2_DataMisfit is creating default weightings for Wd.'
survey = self.survey
eps = np.linalg.norm(Utils.mkvc(survey.dobs),2)*1e-5
self._Wd = Utils.sdiag(1/(abs(survey.dobs)*survey.std+eps))
if getattr(survey,'std', None) is None:
print 'SimPEG.DataMisfit.l2_DataMisfit assigning default std of 5%'
survey.std = 0.05
if getattr(survey, 'eps', None) is None:
print 'SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||'
survey.eps = np.linalg.norm(Utils.mkvc(survey.dobs),2)*1e-5
self._Wd = Utils.sdiag(1/(abs(survey.dobs)*survey.std+survey.eps))
return self._Wd
@Wd.setter
SimPEG/Directives.py
+30
View File
@@ -206,6 +206,36 @@ class SaveOutputEveryIteration(_SaveEveryIteration):
f.write(' %3d %1.4e %1.4e %1.4e %1.4e\n'%(self.opt.iter, self.invProb.beta, self.invProb.phi_d, self.invProb.phi_m, self.opt.f))
f.close()
class SaveOutputDictEveryIteration(_SaveEveryIteration):
"""SaveOutputDictEveryIteration"""
def initialize(self):
print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName
def endIter(self):
# Save the data.
ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
phi_ms = 0.5*ms.dot(ms)
if self.reg.smoothModel == True:
mref = self.reg.mref
else:
mref = 0
mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
phi_mx = 0.5 * mx.dot(mx)
if self.prob.mesh.dim==2:
my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
phi_my = 0.5 * my.dot(my)
else:
phi_my = 'NaN'
if self.prob.mesh.dim==3:
mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) )
phi_mz = 0.5 * mz.dot(mz)
else:
phi_mz = 'NaN'
# Save the file as a npz
np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
SimPEG/EM/Analytics/FDEM.py
+2 -2
View File
@@ -2,7 +2,6 @@ from __future__ import division
import numpy as np
from scipy.constants import mu_0, pi
from scipy.special import erf
import matplotlib.pyplot as plt
from SimPEG import Utils
@@ -59,8 +58,9 @@ def MagneticDipoleWholeSpace(XYZ, srcLoc, sig, f, moment=1., orientation='X', mu
from SimPEG import EM
import matplotlib.pyplot as plt
from scipy.constants import mu_0
freqs = np.logspace(-2,5,100)
Bx, By, Bz = EM.Analytics.FDEM.AnalyticMagDipoleWholeSpace([0,100,0], [0,0,0], 1e-2, freqs, m=1, orientation='Z')
Bx, By, Bz = EM.Analytics.FDEM.MagneticDipoleWholeSpace([0,100,0], [0,0,0], 1e-2, freqs, moment=1, orientation='Z')
plt.loglog(freqs, np.abs(Bz.real)/mu_0, 'b')
plt.loglog(freqs, np.abs(Bz.imag)/mu_0, 'r')
plt.legend(('real','imag'))
SimPEG/EM/Examples/__init__.py
-1
View File
@@ -1 +0,0 @@
import CylInversion
SimPEG/EM/FDEM/FDEM.py
+248 -113
View File
@@ -15,18 +15,20 @@ class BaseFDEMProblem(BaseEMProblem):
.. math ::
\mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\
{\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{M^e} \mathbf{s_e}}
{\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}}
if using the E-B formulation (:code:`Problem_e`
or :code:`Problem_b`) or the magnetic field
or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity.
If we write Maxwell's equations in terms of
\\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\)
.. math ::
\mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{M^e} \mathbf{s_m} \\\\
\mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{s_m} \\\\
\mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`).
if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
"""
@@ -36,7 +38,11 @@ class BaseFDEMProblem(BaseEMProblem):
def fields(self, m=None):
"""
Solve the forward problem for the fields.
Solve the forward problem for the fields.
:param numpy.array m: inversion model (nP,)
:rtype numpy.array:
:return F: forward solution
"""
self.curModel = m
@@ -50,16 +56,22 @@ class BaseFDEMProblem(BaseEMProblem):
Srcs = self.survey.getSrcByFreq(freq)
ftype = self._fieldType + 'Solution'
F[Srcs, ftype] = sol
Ainv.clean()
return F
def Jvec(self, m, v, f=None):
def Jvec(self, m, v, u=None):
"""
Sensitivity times a vector
Sensitivity times a vector.
:param numpy.array m: inversion model (nP,)
:param numpy.array v: vector which we take sensitivity product with (nP,)
:param SimPEG.EM.FDEM.Fields u: fields object
:rtype numpy.array:
:return: Jv (ndata,)
"""
if f is None:
f = self.fields(m)
if u is None:
u = self.fields(m)
self.curModel = m
@@ -71,33 +83,41 @@ class BaseFDEMProblem(BaseEMProblem):
for src in self.survey.getSrcByFreq(freq):
ftype = self._fieldType + 'Solution'
u_src = f[src, ftype]
u_src = u[src, ftype]
dA_dm = self.getADeriv_m(freq, u_src, v)
dRHS_dm = self.getRHSDeriv_m(freq, src, v)
du_dm = Ainv * ( - dA_dm + dRHS_dm )
for rx in src.rxList:
df_duFun = getattr(f, '_%sDeriv_u'%rx.projField, None)
df_duFun = getattr(u, '_%sDeriv_u'%rx.projField, None)
df_dudu_dm = df_duFun(src, du_dm, adjoint=False)
df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None)
df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None)
df_dm = df_dmFun(src, v, adjoint=False)
Df_Dm = np.array(df_dudu_dm + df_dm,dtype=complex)
P = lambda v: rx.projectFieldsDeriv(src, self.mesh, f, v) # wrt u, also have wrt m
P = lambda v: rx.projectFieldsDeriv(src, self.mesh, u, v) # wrt u, also have wrt m
Jv[src, rx] = P(Df_Dm)
Ainv.clean()
return Utils.mkvc(Jv)
def Jtvec(self, m, v, f=None):
def Jtvec(self, m, v, u=None):
"""
Sensitivity transpose times a vector
Sensitivity transpose times a vector
:param numpy.array m: inversion model (nP,)
:param numpy.array v: vector which we take adjoint product with (nP,)
:param SimPEG.EM.FDEM.Fields u: fields object
:rtype numpy.array:
:return: Jv (ndata,)
"""
if f is None:
f = self.fields(m)
if u is None:
u = self.fields(m)
self.curModel = m
@@ -113,12 +133,12 @@ class BaseFDEMProblem(BaseEMProblem):
for src in self.survey.getSrcByFreq(freq):
ftype = self._fieldType + 'Solution'
u_src = f[src, ftype]
u_src = u[src, ftype]
for rx in src.rxList:
PTv = rx.projectFieldsDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
PTv = rx.projectFieldsDeriv(src, self.mesh, u, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
df_duTFun = getattr(f, '_%sDeriv_u'%rx.projField, None)
df_duTFun = getattr(u, '_%sDeriv_u'%rx.projField, None)
df_duT = df_duTFun(src, PTv, adjoint=True)
ATinvdf_duT = ATinv * df_duT
@@ -127,11 +147,12 @@ class BaseFDEMProblem(BaseEMProblem):
dRHS_dmT = self.getRHSDeriv_m(freq,src, ATinvdf_duT, adjoint=True)
du_dmT = -dA_dmT + dRHS_dmT
df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None)
df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None)
dfT_dm = df_dmFun(src, PTv, adjoint=True)
du_dmT += dfT_dm
# TODO: this should be taken care of by the reciever
real_or_imag = rx.projComp
if real_or_imag is 'real':
Jtv += np.array(du_dmT,dtype=complex).real
@@ -139,16 +160,18 @@ class BaseFDEMProblem(BaseEMProblem):
Jtv += - np.array(du_dmT,dtype=complex).real
else:
raise Exception('Must be real or imag')
ATinv.clean()
return Jtv
return Utils.mkvc(Jtv)
def getSourceTerm(self, freq):
"""
Evaluates the sources for a given frequency and puts them in matrix form
Evaluates the sources for a given frequency and puts them in matrix form
:param float freq: Frequency
:rtype: numpy.ndarray (nE or nF, nSrc)
:return: S_m, S_e
:param float freq: Frequency
:rtype: (numpy.ndarray, numpy.ndarray)
:return: S_m, S_e (nE or nF, nSrc)
"""
Srcs = self.survey.getSrcByFreq(freq)
if self._eqLocs is 'FE':
@@ -172,20 +195,22 @@ class BaseFDEMProblem(BaseEMProblem):
class Problem_e(BaseFDEMProblem):
"""
By eliminating the magnetic flux density using
.. math ::
\mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} + \mathbf{s_m}\\right)
we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
By eliminating the magnetic flux density using
.. math ::
\\left(\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{C}+ i \omega \mathbf{M^e_{\sigma}} \\right)\mathbf{e} = \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M^e}\mathbf{s_e}
\mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} + \mathbf{s_m}\\right)
which we solve for \\\(\\\mathbf{e}\\\).
we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
.. math ::
\\left(\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{C}+ i \omega \mathbf{M^e_{\sigma}} \\right)\mathbf{e} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M^e}\mathbf{s_e}
which we solve for :math:`\mathbf{e}`.
:param SimPEG.Mesh mesh: mesh
"""
_fieldType = 'e'
@@ -197,13 +222,16 @@ class Problem_e(BaseFDEMProblem):
def getA(self, freq):
"""
.. math ::
\mathbf{A} = \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{C} + i \omega \mathbf{M^e_{\sigma}}
System matrix
.. math ::
\mathbf{A} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{C} + i \omega \mathbf{M^e_{\sigma}}
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
:return: A
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
:return: A
"""
MfMui = self.MfMui
MeSigma = self.MeSigma
C = self.mesh.edgeCurl
@@ -212,6 +240,20 @@ class Problem_e(BaseFDEMProblem):
def getADeriv_m(self, freq, u, v, adjoint=False):
"""
Product of the derivative of our system matrix with respect to the model and a vector
.. math ::
\\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = i \omega \\frac{d \mathbf{M^e_{\sigma}}\mathbf{v} }{d\mathbf{m}}
:param float freq: frequency
:param numpy.ndarray u: solution vector (nE,)
:param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
"""
dsig_dm = self.curModel.sigmaDeriv
dMe_dsig = self.MeSigmaDeriv(u)
@@ -222,26 +264,37 @@ class Problem_e(BaseFDEMProblem):
def getRHS(self, freq):
"""
.. math ::
\mathbf{RHS} = \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M_e}\mathbf{s_e}
Right hand side for the system
:param float freq: Frequency
:rtype: numpy.ndarray (nE, nSrc)
:return: RHS
.. math ::
\mathbf{RHS} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M_e}\mathbf{s_e}
:param float freq: Frequency
:rtype: numpy.ndarray
:return: RHS (nE, nSrc)
"""
S_m, S_e = self.getSourceTerm(freq)
C = self.mesh.edgeCurl
MfMui = self.MfMui
RHS = C.T * (MfMui * S_m) -1j * omega(freq) * S_e
return RHS
return C.T * (MfMui * S_m) -1j * omega(freq) * S_e
def getRHSDeriv_m(self, freq, src, v, adjoint=False):
"""
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.Src src: FDEM source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of rhs deriv with a vector
"""
C = self.mesh.edgeCurl
MfMui = self.MfMui
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
if adjoint:
dRHS = MfMui * (C * v)
@@ -252,20 +305,22 @@ class Problem_e(BaseFDEMProblem):
class Problem_b(BaseFDEMProblem):
"""
We eliminate \\\(\\\mathbf{e}\\\) using
We eliminate :math:`\mathbf{e}` using
.. math ::
.. math ::
\mathbf{e} = \mathbf{M^e_{\sigma}}^{-1} \\left(\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{s_e}\\right)
\mathbf{e} = \mathbf{M^e_{\sigma}}^{-1} \\left(\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{s_e}\\right)
and solve for \\\(\\\mathbf{b}\\\) using:
and solve for :math:`\mathbf{b}` using:
.. math ::
.. math ::
\\left(\mathbf{C} \mathbf{M^e_{\sigma}}^{-1} \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} + i \omega \\right)\mathbf{b} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{M^e}\mathbf{s_e}
\\left(\mathbf{C} \mathbf{M^e_{\sigma}}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} + i \omega \\right)\mathbf{b} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{M^e}\mathbf{s_e}
.. note ::
The inverse problem will not work with full anisotropy
.. note ::
The inverse problem will not work with full anisotropy
:param SimPEG.Mesh mesh: mesh
"""
_fieldType = 'b'
@@ -277,12 +332,14 @@ class Problem_b(BaseFDEMProblem):
def getA(self, freq):
"""
.. math ::
\mathbf{A} = \mathbf{C} \mathbf{M^e_{\sigma}}^{-1} \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} + i \omega
System matrix
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
:return: A
.. math ::
\mathbf{A} = \mathbf{C} \mathbf{M^e_{\sigma}}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} + i \omega
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
:return: A
"""
MfMui = self.MfMui
@@ -298,6 +355,20 @@ class Problem_b(BaseFDEMProblem):
def getADeriv_m(self, freq, u, v, adjoint=False):
"""
Product of the derivative of our system matrix with respect to the model and a vector
.. math ::
\\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = \mathbf{C} \\frac{\mathbf{M^e_{\sigma}} \mathbf{v}}{d\mathbf{m}}
:param float freq: frequency
:param numpy.ndarray u: solution vector (nF,)
:param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
"""
MfMui = self.MfMui
C = self.mesh.edgeCurl
MeSigmaIDeriv = self.MeSigmaIDeriv
@@ -317,12 +388,14 @@ class Problem_b(BaseFDEMProblem):
def getRHS(self, freq):
"""
.. math ::
\mathbf{RHS} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{s_e}
Right hand side for the system
:param float freq: Frequency
:rtype: numpy.ndarray (nE, nSrc)
:return: RHS
.. math ::
\mathbf{RHS} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{s_e}
:param float freq: Frequency
:rtype: numpy.ndarray
:return: RHS (nE, nSrc)
"""
S_m, S_e = self.getSourceTerm(freq)
@@ -338,6 +411,17 @@ class Problem_b(BaseFDEMProblem):
return RHS
def getRHSDeriv_m(self, freq, src, v, adjoint=False):
"""
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.Src src: FDEM source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of rhs deriv with a vector
"""
C = self.mesh.edgeCurl
S_m, S_e = src.eval(self)
MfMui = self.MfMui
@@ -346,7 +430,7 @@ class Problem_b(BaseFDEMProblem):
v = self.MfMui * v
MeSigmaIDeriv = self.MeSigmaIDeriv(S_e)
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
if not adjoint:
RHSderiv = C * (MeSigmaIDeriv * v)
@@ -369,21 +453,22 @@ class Problem_b(BaseFDEMProblem):
class Problem_j(BaseFDEMProblem):
"""
We eliminate \\\(\\\mathbf{h}\\\) using
We eliminate \\\(\\\mathbf{h}\\\) using
.. math ::
.. math ::
\mathbf{h} = \\frac{1}{i \omega} \mathbf{M_{\mu}^e}^{-1} \\left(-\mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{j} + \mathbf{M^e} \mathbf{s_m} \\right)
\mathbf{h} = \\frac{1}{i \omega} \mathbf{M_{\mu}^e}^{-1} \\left(-\mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + \mathbf{M^e} \mathbf{s_m} \\right)
and solve for \\\(\\\mathbf{j}\\\) using
and solve for \\\(\\\mathbf{j}\\\) using
.. math ::
.. math ::
\\left(\mathbf{C} \mathbf{M_{\mu}^e}^{-1} \mathbf{C}^T \mathbf{M_{\\rho}^f} + i \omega\\right)\mathbf{j} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1} \mathbf{M^e} \mathbf{s_m} -i\omega\mathbf{s_e}
\\left(\mathbf{C} \mathbf{M_{\mu}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} + i \omega\\right)\mathbf{j} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1} \mathbf{M^e} \mathbf{s_m} -i\omega\mathbf{s_e}
.. note::
This implementation does not yet work with full anisotropy!!
.. note::
This implementation does not yet work with full anisotropy!!
:param SimPEG.Mesh mesh: mesh
"""
_fieldType = 'j'
@@ -395,12 +480,14 @@ class Problem_j(BaseFDEMProblem):
def getA(self, freq):
"""
.. math ::
\\mathbf{A} = \\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C}^T \\mathbf{M^f_{\\sigma^{-1}}} + i\\omega
System matrix
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
:return: A
.. math ::
\\mathbf{A} = \\mathbf{C} \\mathbf{M^e_{\\mu^{-1}}} \\mathbf{C}^{\\top} \\mathbf{M^f_{\\sigma^{-1}}} + i\\omega
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
:return: A
"""
MeMuI = self.MeMuI
@@ -417,12 +504,20 @@ class Problem_j(BaseFDEMProblem):
def getADeriv_m(self, freq, u, v, adjoint=False):
"""
In this case, we assume that electrical conductivity, \\\(\\\sigma\\\) is the physical property of interest (i.e. \\\(\\\sigma\\\) = model.transform). Then we want
Product of the derivative of our system matrix with respect to the model and a vector
.. math ::
In this case, we assume that electrical conductivity, :math:`\sigma` is the physical property of interest (i.e. :math:`\sigma` = model.transform). Then we want
\\frac{\mathbf{A(\sigma)} \mathbf{v}}{d \\mathbf{m}} &= \\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C^T} \\frac{d \\mathbf{M^f_{\\sigma^{-1}}}}{d \\mathbf{m}}
&= \\mathbf{C} \\mathbf{M^e_{mu}^{-1}} \\mathbf{C^T} \\frac{d \\mathbf{M^f_{\\sigma^{-1}}}}{d \\mathbf{\\sigma^{-1}}} \\frac{d \\mathbf{\\sigma^{-1}}}{d \\mathbf{\\sigma}} \\frac{d \\mathbf{\\sigma}}{d \\mathbf{m}}
.. math ::
\\frac{\mathbf{A(\sigma)} \mathbf{v}}{d \mathbf{m}} = \mathbf{C} \mathbf{M^e_{mu^{-1}}} \mathbf{C^{\\top}} \\frac{d \mathbf{M^f_{\sigma^{-1}}}\mathbf{v} }{d \mathbf{m}}
:param float freq: frequency
:param numpy.ndarray u: solution vector (nF,)
:param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
"""
MeMuI = self.MeMuI
@@ -442,12 +537,15 @@ class Problem_j(BaseFDEMProblem):
def getRHS(self, freq):
"""
.. math ::
Right hand side for the system
\mathbf{RHS} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1}\mathbf{s_m} -i\omega \mathbf{s_e}
:param float freq: Frequency
:rtype: numpy.ndarray (nE, nSrc)
:return: RHS
.. math ::
\mathbf{RHS} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1}\mathbf{s_m} -i\omega \mathbf{s_e}
:param float freq: Frequency
:rtype: numpy.ndarray (nE, nSrc)
:return: RHS
"""
S_m, S_e = self.getSourceTerm(freq)
@@ -462,9 +560,20 @@ class Problem_j(BaseFDEMProblem):
return RHS
def getRHSDeriv_m(self, freq, src, v, adjoint=False):
"""
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.Src src: FDEM source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of rhs deriv with a vector
"""
C = self.mesh.edgeCurl
MeMuI = self.MeMuI
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
if adjoint:
if self._makeASymmetric:
@@ -485,18 +594,19 @@ class Problem_j(BaseFDEMProblem):
class Problem_h(BaseFDEMProblem):
"""
We eliminate \\\(\\\mathbf{j}\\\) using
We eliminate \\\(\\\mathbf{j}\\\) using
.. math ::
.. math ::
\mathbf{j} = \mathbf{C} \mathbf{h} - \mathbf{s_e}
\mathbf{j} = \mathbf{C} \mathbf{h} - \mathbf{s_e}
and solve for \\\(\\\mathbf{h}\\\) using
and solve for \\\(\\\mathbf{h}\\\) using
.. math ::
.. math ::
\\left(\mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}\\right) \mathbf{h} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{s_e}
\\left(\mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}\\right) \mathbf{h} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{s_e}
:param SimPEG.Mesh mesh: mesh
"""
_fieldType = 'h'
@@ -508,13 +618,14 @@ class Problem_h(BaseFDEMProblem):
def getA(self, freq):
"""
.. math ::
System matrix
\mathbf{A} = \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}
.. math::
\mathbf{A} = \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
:return: A
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
:return: A
"""
MeMu = self.MeMu
@@ -524,6 +635,19 @@ class Problem_h(BaseFDEMProblem):
return C.T * (MfRho * C) + 1j*omega(freq)*MeMu
def getADeriv_m(self, freq, u, v, adjoint=False):
"""
Product of the derivative of our system matrix with respect to the model and a vector
.. math::
\\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = \mathbf{C}^{\\top}\\frac{d \mathbf{M^f_{\\rho}}\mathbf{v} }{d\mathbf{m}}
:param float freq: frequency
:param numpy.ndarray u: solution vector (nE,)
:param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
"""
MeMu = self.MeMu
C = self.mesh.edgeCurl
@@ -535,24 +659,35 @@ class Problem_h(BaseFDEMProblem):
def getRHS(self, freq):
"""
.. math ::
Right hand side for the system
\mathbf{RHS} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{s_e}
.. math ::
:param float freq: Frequency
:rtype: numpy.ndarray (nE, nSrc)
:return: RHS
\mathbf{RHS} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{s_e}
:param float freq: Frequency
:rtype: numpy.ndarray
:return: RHS (nE, nSrc)
"""
S_m, S_e = self.getSourceTerm(freq)
C = self.mesh.edgeCurl
MfRho = self.MfRho
RHS = S_m + C.T * ( MfRho * S_e )
return RHS
return S_m + C.T * ( MfRho * S_e )
def getRHSDeriv_m(self, freq, src, v, adjoint=False):
"""
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.Src src: FDEM source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of rhs deriv with a vector
"""
_, S_e = src.eval(self)
C = self.mesh.edgeCurl
MfRho = self.MfRho
@@ -563,7 +698,7 @@ class Problem_h(BaseFDEMProblem):
elif adjoint:
RHSDeriv = MfRhoDeriv.T * (C * v)
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
return RHSDeriv + S_mDeriv(v) + C.T * (MfRho * S_eDeriv(v))
SimPEG/EM/FDEM/FieldsFDEM.py
+512 -9
View File
@@ -7,11 +7,39 @@ from SimPEG.Utils import Zero, Identity
class Fields(SimPEG.Problem.Fields):
"""Fancy Field Storage for a FDEM survey."""
"""
Fancy Field Storage for a FDEM survey. Only one field type is stored for
each problem, the rest are computed. The fields obejct acts like an array and is indexed by
.. code-block:: python
f = problem.fields(m)
e = f[srcList,'e']
b = f[srcList,'b']
If accessing all sources for a given field, use the :code:`:`
.. code-block:: python
f = problem.fields(m)
e = f[:,'e']
b = f[:,'b']
The array returned will be size (nE or nF, nSrcs :math:`\\times` nFrequencies)
"""
knownFields = {}
dtype = complex
class Fields_e(Fields):
"""
Fields object for Problem_e.
:param Mesh mesh: mesh
:param Survey survey: survey
"""
knownFields = {'eSolution':'E'}
aliasFields = {
'e' : ['eSolution','E','_e'],
@@ -30,6 +58,15 @@ class Fields_e(Fields):
self._edgeCurl = self.survey.prob.mesh.edgeCurl
def _ePrimary(self, eSolution, srcList):
"""
Primary electric field from source
:param numpy.ndarray eSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: primary electric field as defined by the sources
"""
ePrimary = np.zeros_like(eSolution)
for i, src in enumerate(srcList):
ep = src.ePrimary(self.prob)
@@ -37,19 +74,67 @@ class Fields_e(Fields):
return ePrimary
def _eSecondary(self, eSolution, srcList):
"""
Secondary electric field is the thing we solved for
:param numpy.ndarray eSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: secondary electric field
"""
return eSolution
def _e(self, eSolution, srcList):
"""
Total electric field is sum of primary and secondary
:param numpy.ndarray eSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: total electric field
"""
return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList)
def _eDeriv_u(self, src, v, adjoint = False):
"""
Derivative of the total electric field with respect to the thing we
solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the derivative of the electric field with respect to the field we solved for with a vector
"""
return Identity()*v
def _eDeriv_m(self, src, v, adjoint = False):
"""
Derivative of the total electric field with respect to the inversion model. Here, we assume that the primary does not depend on the model.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: SimPEG.Utils.Zero
:return: product of the electric field derivative with respect to the inversion model with a vector
"""
# assuming primary does not depend on the model
return Zero()
def _bPrimary(self, eSolution, srcList):
"""
Primary magnetic flux density from source
:param numpy.ndarray eSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: primary magnetic flux density as defined by the sources
"""
bPrimary = np.zeros([self._edgeCurl.shape[0],eSolution.shape[1]],dtype = complex)
for i, src in enumerate(srcList):
bp = src.bPrimary(self.prob)
@@ -57,6 +142,15 @@ class Fields_e(Fields):
return bPrimary
def _bSecondary(self, eSolution, srcList):
"""
Secondary magnetic flux density from eSolution
:param numpy.ndarray eSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: secondary magnetic flux density
"""
C = self._edgeCurl
b = (C * eSolution)
for i, src in enumerate(srcList):
@@ -66,29 +160,84 @@ class Fields_e(Fields):
return b
def _bSecondaryDeriv_u(self, src, v, adjoint = False):
"""
Derivative of the secondary magnetic flux density with respect to the thing we solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the derivative of the secondary magnetic flux density with respect to the field we solved for with a vector
"""
C = self._edgeCurl
if adjoint:
return - 1./(1j*omega(src.freq)) * (C.T * v)
return - 1./(1j*omega(src.freq)) * (C * v)
def _bSecondaryDeriv_m(self, src, v, adjoint = False):
S_mDeriv, _ = src.evalDeriv(self.prob, adjoint)
S_mDeriv = S_mDeriv(v)
"""
Derivative of the secondary magnetic flux density with respect to the inversion model.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the secondary magnetic flux density derivative with respect to the inversion model with a vector
"""
S_mDeriv, _ = src.evalDeriv(self.prob, v, adjoint)
return 1./(1j * omega(src.freq)) * S_mDeriv
def _b(self, eSolution, srcList):
"""
Total magnetic flux density is sum of primary and secondary
:param numpy.ndarray eSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: total magnetic flux density
"""
return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList)
def _bDeriv_u(self, src, v, adjoint=False):
"""
Derivative of the total magnetic flux density with respect to the thing we solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
"""
# Primary does not depend on u
return self._bSecondaryDeriv_u(src, v, adjoint)
def _bDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the total magnetic flux density with respect to the inversion model.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: SimPEG.Utils.Zero
:return: product of the magnetic flux density derivative with respect to the inversion model with a vector
"""
# Assuming the primary does not depend on the model
return self._bSecondaryDeriv_m(src, v, adjoint)
class Fields_b(Fields):
"""
Fields object for Problem_b.
:param Mesh mesh: mesh
:param Survey survey: survey
"""
knownFields = {'bSolution':'F'}
aliasFields = {
'b' : ['bSolution','F','_b'],
@@ -111,6 +260,15 @@ class Fields_b(Fields):
self._Me = self.survey.prob.Me
def _bPrimary(self, bSolution, srcList):
"""
Primary magnetic flux density from source
:param numpy.ndarray bSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: primary electric field as defined by the sources
"""
bPrimary = np.zeros_like(bSolution)
for i, src in enumerate(srcList):
bp = src.bPrimary(self.prob)
@@ -118,19 +276,66 @@ class Fields_b(Fields):
return bPrimary
def _bSecondary(self, bSolution, srcList):
"""
Secondary magnetic flux density is the thing we solved for
:param numpy.ndarray bSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: secondary magnetic flux density
"""
return bSolution
def _b(self, bSolution, srcList):
"""
Total magnetic flux density is sum of primary and secondary
:param numpy.ndarray bSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: total magnetic flux density
"""
return self._bPrimary(bSolution, srcList) + self._bSecondary(bSolution, srcList)
def _bDeriv_u(self, src, v, adjoint=False):
"""
Derivative of the total magnetic flux density with respect to the thing we
solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
"""
return Identity()*v
def _bDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the total magnetic flux density with respect to the inversion model. Here, we assume that the primary does not depend on the model.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: SimPEG.Utils.Zero
:return: product of the magnetic flux density derivative with respect to the inversion model with a vector
"""
# assuming primary does not depend on the model
return Zero()
def _ePrimary(self, bSolution, srcList):
"""
Primary electric field from source
:param numpy.ndarray bSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: primary electric field as defined by the sources
"""
ePrimary = np.zeros([self._edgeCurl.shape[1],bSolution.shape[1]],dtype = complex)
for i,src in enumerate(srcList):
ep = src.ePrimary(self.prob)
@@ -138,6 +343,15 @@ class Fields_b(Fields):
return ePrimary
def _eSecondary(self, bSolution, srcList):
"""
Secondary electric field from bSolution
:param numpy.ndarray bSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: secondary electric field
"""
e = self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * bSolution))
for i,src in enumerate(srcList):
_,S_e = src.eval(self.prob)
@@ -145,12 +359,32 @@ class Fields_b(Fields):
return e
def _eSecondaryDeriv_u(self, src, v, adjoint=False):
"""
Derivative of the secondary electric field with respect to the thing we solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the derivative of the secondary electric field with respect to the field we solved for with a vector
"""
if not adjoint:
return self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * v) )
else:
return self._MfMui.T * (self._edgeCurl * (self._MeSigmaI.T * v))
def _eSecondaryDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the secondary electric field with respect to the inversion model
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the derivative of the secondary electric field with respect to the model with a vector
"""
bSolution = self[[src],'bSolution']
_,S_e = src.eval(self.prob)
Me = self._Me
@@ -166,25 +400,60 @@ class Fields_b(Fields):
elif adjoint:
de_dm = self._MeSigmaIDeriv(w).T * v
_, S_eDeriv = src.evalDeriv(self.prob, adjoint)
Se_Deriv = S_eDeriv(v)
_, S_eDeriv = src.evalDeriv(self.prob, v, adjoint)
de_dm = de_dm - self._MeSigmaI * Se_Deriv
de_dm = de_dm - self._MeSigmaI * S_eDeriv
return de_dm
def _e(self, bSolution, srcList):
"""
Total electric field is sum of primary and secondary
:param numpy.ndarray eSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: total electric field
"""
return self._ePrimary(bSolution, srcList) + self._eSecondary(bSolution, srcList)
def _eDeriv_u(self, src, v, adjoint=False):
"""
Derivative of the total electric field with respect to the thing we solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the derivative of the electric field with respect to the field we solved for with a vector
"""
return self._eSecondaryDeriv_u(src, v, adjoint)
def _eDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the total electric field density with respect to the inversion model.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the electric field derivative with respect to the inversion model with a vector
"""
# assuming primary doesn't depend on model
return self._eSecondaryDeriv_m(src, v, adjoint)
class Fields_j(Fields):
"""
Fields object for Problem_j.
:param Mesh mesh: mesh
:param Survey survey: survey
"""
knownFields = {'jSolution':'F'}
aliasFields = {
'j' : ['jSolution','F','_j'],
@@ -207,6 +476,15 @@ class Fields_j(Fields):
self._Me = self.survey.prob.Me
def _jPrimary(self, jSolution, srcList):
"""
Primary current density from source
:param numpy.ndarray jSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: primary current density as defined by the sources
"""
jPrimary = np.zeros_like(jSolution,dtype = complex)
for i, src in enumerate(srcList):
jp = src.jPrimary(self.prob)
@@ -214,19 +492,66 @@ class Fields_j(Fields):
return jPrimary
def _jSecondary(self, jSolution, srcList):
"""
Secondary current density is the thing we solved for
:param numpy.ndarray jSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: secondary current density
"""
return jSolution
def _j(self, jSolution, srcList):
"""
Total current density is sum of primary and secondary
:param numpy.ndarray jSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: total current density
"""
return self._jPrimary(jSolution, srcList) + self._jSecondary(jSolution, srcList)
def _jDeriv_u(self, src, v, adjoint=False):
"""
Derivative of the total current density with respect to the thing we
solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the derivative of the current density with respect to the field we solved for with a vector
"""
return Identity()*v
def _jDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the total current density with respect to the inversion model. Here, we assume that the primary does not depend on the model.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: SimPEG.Utils.Zero
:return: product of the current density derivative with respect to the inversion model with a vector
"""
# assuming primary does not depend on the model
return Zero()
def _hPrimary(self, jSolution, srcList):
"""
Primary magnetic field from source
:param numpy.ndarray hSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: primary magnetic field as defined by the sources
"""
hPrimary = np.zeros([self._edgeCurl.shape[1],jSolution.shape[1]],dtype = complex)
for i, src in enumerate(srcList):
hp = src.hPrimary(self.prob)
@@ -234,6 +559,15 @@ class Fields_j(Fields):
return hPrimary
def _hSecondary(self, jSolution, srcList):
"""
Secondary magnetic field from bSolution
:param numpy.ndarray jSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: secondary magnetic field
"""
h = self._MeMuI * (self._edgeCurl.T * (self._MfRho * jSolution) )
for i, src in enumerate(srcList):
h[:,i] *= -1./(1j*omega(src.freq))
@@ -242,12 +576,32 @@ class Fields_j(Fields):
return h
def _hSecondaryDeriv_u(self, src, v, adjoint=False):
"""
Derivative of the secondary magnetic field with respect to the thing we solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the derivative of the secondary magnetic field with respect to the field we solved for with a vector
"""
if not adjoint:
return -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * v) )
elif adjoint:
return -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * v))
def _hSecondaryDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the secondary magnetic field with respect to the inversion model
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the derivative of the secondary magnetic field with respect to the model with a vector
"""
jSolution = self[[src],'jSolution']
MeMuI = self._MeMuI
C = self._edgeCurl
@@ -260,7 +614,7 @@ class Fields_j(Fields):
elif adjoint:
hDeriv_m = -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) )
S_mDeriv,_ = src.evalDeriv(self.prob, adjoint)
S_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
if not adjoint:
S_mDeriv = S_mDeriv(v)
@@ -272,17 +626,53 @@ class Fields_j(Fields):
def _h(self, jSolution, srcList):
"""
Total magnetic field is sum of primary and secondary
:param numpy.ndarray eSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: total magnetic field
"""
return self._hPrimary(jSolution, srcList) + self._hSecondary(jSolution, srcList)
def _hDeriv_u(self, src, v, adjoint=False):
"""
Derivative of the total magnetic field with respect to the thing we solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
"""
return self._hSecondaryDeriv_u(src, v, adjoint)
def _hDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the total magnetic field density with respect to the inversion model.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the magnetic field derivative with respect to the inversion model with a vector
"""
# assuming the primary doesn't depend on the model
return self._hSecondaryDeriv_m(src, v, adjoint)
class Fields_h(Fields):
"""
Fields object for Problem_h.
:param Mesh mesh: mesh
:param Survey survey: survey
"""
knownFields = {'hSolution':'E'}
aliasFields = {
'h' : ['hSolution','E','_h'],
@@ -303,6 +693,15 @@ class Fields_h(Fields):
self._MfRho = self.survey.prob.MfRho
def _hPrimary(self, hSolution, srcList):
"""
Primary magnetic field from source
:param numpy.ndarray eSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: primary magnetic field as defined by the sources
"""
hPrimary = np.zeros_like(hSolution,dtype = complex)
for i, src in enumerate(srcList):
hp = src.hPrimary(self.prob)
@@ -310,19 +709,67 @@ class Fields_h(Fields):
return hPrimary
def _hSecondary(self, hSolution, srcList):
"""
Secondary magnetic field is the thing we solved for
:param numpy.ndarray hSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: secondary magnetic field
"""
return hSolution
def _h(self, hSolution, srcList):
"""
Total magnetic field is sum of primary and secondary
:param numpy.ndarray hSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: total magnetic field
"""
return self._hPrimary(hSolution, srcList) + self._hSecondary(hSolution, srcList)
def _hDeriv_u(self, src, v, adjoint=False):
"""
Derivative of the total magnetic field with respect to the thing we
solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
"""
return Identity()*v
def _hDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the total magnetic field with respect to the inversion model. Here, we assume that the primary does not depend on the model.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: SimPEG.Utils.Zero
:return: product of the magnetic field derivative with respect to the inversion model with a vector
"""
# assuming primary does not depend on the model
return Zero()
def _jPrimary(self, hSolution, srcList):
"""
Primary current density from source
:param numpy.ndarray hSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: primary current density as defined by the sources
"""
jPrimary = np.zeros([self._edgeCurl.shape[0], hSolution.shape[1]], dtype = complex)
for i, src in enumerate(srcList):
jp = src.jPrimary(self.prob)
@@ -330,6 +777,15 @@ class Fields_h(Fields):
return jPrimary
def _jSecondary(self, hSolution, srcList):
"""
Secondary current density from eSolution
:param numpy.ndarray hSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: secondary current density
"""
j = self._edgeCurl*hSolution
for i, src in enumerate(srcList):
_,S_e = src.eval(self.prob)
@@ -337,22 +793,69 @@ class Fields_h(Fields):
return j
def _jSecondaryDeriv_u(self, src, v, adjoint=False):
"""
Derivative of the secondary current density with respect to the thing we solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the derivative of the secondary current density with respect to the field we solved for with a vector
"""
if not adjoint:
return self._edgeCurl*v
elif adjoint:
return self._edgeCurl.T*v
def _jSecondaryDeriv_m(self, src, v, adjoint=False):
_,S_eDeriv = src.evalDeriv(self.prob, adjoint)
S_eDeriv = S_eDeriv(v)
"""
Derivative of the secondary current density with respect to the inversion model.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the secondary current density derivative with respect to the inversion model with a vector
"""
_,S_eDeriv = src.evalDeriv(self.prob, v, adjoint)
return -S_eDeriv
def _j(self, hSolution, srcList):
"""
Total current density is sum of primary and secondary
:param numpy.ndarray eSolution: field we solved for
:param list srcList: list of sources
:rtype: numpy.ndarray
:return: total current density
"""
return self._jPrimary(hSolution, srcList) + self._jSecondary(hSolution, srcList)
def _jDeriv_u(self, src, v, adjoint=False):
"""
Derivative of the total current density with respect to the thing we solved for
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of the derivative of the current density with respect to the field we solved for with a vector
"""
return self._jSecondaryDeriv_u(src,v,adjoint)
def _jDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the total current density with respect to the inversion model.
:param SimPEG.EM.FDEM.Src src: source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: SimPEG.Utils.Zero
:return: product of the current density with respect to the inversion model with a vector
"""
# assuming the primary does not depend on the model
return self._jSecondaryDeriv_m(src,v,adjoint)
SimPEG/EM/FDEM/SrcFDEM.py
+274 -19
View File
@@ -1,55 +1,141 @@
from SimPEG import Survey, Problem, Utils, np, sp
from scipy.constants import mu_0
from SimPEG.EM.Utils import *
from SimPEG.Utils import Zero
# from SurveyFDEM import Rx
from SimPEG.Utils import Zero
class BaseSrc(Survey.BaseSrc):
"""
Base source class for FDEM Survey
"""
freq = None
# rxPair = Rx
# rxPair = RxFDEM
integrate = True
def eval(self, prob):
"""
Evaluate the source terms.
- :math:`S_m` : magnetic source term
- :math:`S_e` : electric source term
:param Problem prob: FDEM Problem
:rtype: (numpy.ndarray, numpy.ndarray)
:return: tuple with magnetic source term and electric source term
"""
S_m = self.S_m(prob)
S_e = self.S_e(prob)
return S_m, S_e
def evalDeriv(self, prob, adjoint=False):
return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint)
def evalDeriv(self, prob, v=None, adjoint=False):
"""
Derivatives of the source terms with respect to the inversion model
- :code:`S_mDeriv` : derivative of the magnetic source term
- :code:`S_eDeriv` : derivative of the electric source term
:param Problem prob: FDEM Problem
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: (numpy.ndarray, numpy.ndarray)
:return: tuple with magnetic source term and electric source term derivatives times a vector
"""
if v is not None:
return self.S_mDeriv(prob,v,adjoint), self.S_eDeriv(prob,v,adjoint)
else:
return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint)
def bPrimary(self, prob):
"""
Primary magnetic flux density
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: primary magnetic flux density
"""
return Zero()
def hPrimary(self, prob):
"""
Primary magnetic field
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
return Zero()
def ePrimary(self, prob):
"""
Primary electric field
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: primary electric field
"""
return Zero()
def jPrimary(self, prob):
"""
Primary current density
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: primary current density
"""
return Zero()
def S_m(self, prob):
"""
Magnetic source term
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: magnetic source term on mesh
"""
return Zero()
def S_e(self, prob):
"""
Electric source term
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: electric source term on mesh
"""
return Zero()
def S_mDeriv(self, prob, v, adjoint = False):
"""
Derivative of magnetic source term with respect to the inversion model
:param Problem prob: FDEM Problem
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of magnetic source term derivative with a vector
"""
return Zero()
def S_eDeriv(self, prob, v, adjoint = False):
"""
Derivative of electric source term with respect to the inversion model
:param Problem prob: FDEM Problem
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: product of electric source term derivative with a vector
"""
return Zero()
class RawVec_e(BaseSrc):
"""
RawVec electric source. It is defined by the user provided vector S_e
RawVec electric source. It is defined by the user provided vector S_e
:param numpy.array S_e: electric source term
:param float freq: frequency
:param rxList: receiver list
:param list rxList: receiver list
:param float freq: frequency
:param numpy.array S_e: electric source term
"""
def __init__(self, rxList, freq, S_e): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
@@ -58,16 +144,17 @@ class RawVec_e(BaseSrc):
BaseSrc.__init__(self, rxList)
def S_e(self, prob):
return self._S_e
class RawVec_m(BaseSrc):
"""
RawVec magnetic source. It is defined by the user provided vector S_m
RawVec magnetic source. It is defined by the user provided vector S_m
:param numpy.array S_m: magnetic source term
:param float freq: frequency
:param rxList: receiver list
:param float freq: frequency
:param rxList: receiver list
:param numpy.array S_m: magnetic source term
"""
def __init__(self, rxList, freq, S_m, integrate = True): #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()):
@@ -78,17 +165,24 @@ class RawVec_m(BaseSrc):
BaseSrc.__init__(self, rxList)
def S_m(self, prob):
"""
Magnetic source term
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: magnetic source term on mesh
"""
return self._S_m
class RawVec(BaseSrc):
"""
RawVec source. It is defined by the user provided vectors S_m, S_e
RawVec source. It is defined by the user provided vectors S_m, S_e
:param numpy.array S_m: magnetic source term
:param numpy.array S_e: electric source term
:param float freq: frequency
:param rxList: receiver list
:param rxList: receiver list
:param float freq: frequency
:param numpy.array S_m: magnetic source term
:param numpy.array S_e: electric source term
"""
def __init__(self, rxList, freq, S_m, S_e, integrate = True):
self._S_m = np.array(S_m,dtype=complex)
@@ -109,6 +203,51 @@ class RawVec(BaseSrc):
class MagDipole(BaseSrc):
"""
Point magnetic dipole source calculated by taking the curl of a magnetic
vector potential. By taking the discrete curl, we ensure that the magnetic
flux density is divergence free (no magnetic monopoles!).
This approach uses a primary-secondary in frequency. Here we show the
derivation for E-B formulation noting that similar steps are followed for
the H-J formulation.
.. math::
\mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\
{\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}}
We split up the fields and :math:`\mu^{-1}` into primary (:math:`\mathbf{P}`) and secondary (:math:`\mathbf{S}`) components
- :math:`\mathbf{e} = \mathbf{e^P} + \mathbf{e^S}`
- :math:`\mathbf{b} = \mathbf{b^P} + \mathbf{b^S}`
- :math:`\\boldsymbol{\mu}^{\mathbf{-1}} = \\boldsymbol{\mu}^{\mathbf{-1}^\mathbf{P}} + \\boldsymbol{\mu}^{\mathbf{-1}^\mathbf{S}}`
and define a zero-frequency primary problem, noting that the source is
generated by a divergence free electric current
.. math::
\mathbf{C} \mathbf{e^P} = \mathbf{s_m^P} = 0 \\\\
{\mathbf{C}^T \mathbf{{M_{\mu^{-1}}^f}^P} \mathbf{b^P} - \mathbf{M_{\sigma}^e} \mathbf{e^P} = \mathbf{M^e} \mathbf{s_e^P}}
Since :math:`\mathbf{e^P}` is curl-free, divergence-free, we assume that there is no constant field background, the :math:`\mathbf{e^P} = 0`, so our primary problem is
.. math::
\mathbf{e^P} = 0 \\\\
{\mathbf{C}^T \mathbf{{M_{\mu^{-1}}^f}^P} \mathbf{b^P} = \mathbf{s_e^P}}
Our secondary problem is then
.. math::
\mathbf{C} \mathbf{e^S} + i \omega \mathbf{b^S} = - i \omega \mathbf{b^P} \\\\
{\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b^S} - \mathbf{M_{\sigma}^e} \mathbf{e^S} = -\mathbf{C}^T \mathbf{{M_{\mu^{-1}}^f}^S} \mathbf{b^P}}
:param list rxList: receiver list
:param float freq: frequency
:param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`)
:param string orientation: 'X', 'Y', 'Z'
:param float moment: magnetic dipole moment
:param float mu: background magnetic permeability
"""
#TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0):
@@ -121,6 +260,13 @@ class MagDipole(BaseSrc):
BaseSrc.__init__(self, rxList)
def bPrimary(self, prob):
"""
The primary magnetic flux density from a magnetic vector potential
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
eqLocs = prob._eqLocs
if eqLocs is 'FE':
@@ -152,14 +298,37 @@ class MagDipole(BaseSrc):
return C*a
def hPrimary(self, prob):
"""
The primary magnetic field from a magnetic vector potential
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
b = self.bPrimary(prob)
return h_from_b(prob,b)
def S_m(self, prob):
"""
The magnetic source term
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
b_p = self.bPrimary(prob)
return -1j*omega(self.freq)*b_p
def S_e(self, prob):
"""
The electric source term
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
return Zero()
else:
@@ -179,6 +348,21 @@ class MagDipole(BaseSrc):
class MagDipole_Bfield(BaseSrc):
"""
Point magnetic dipole source calculated with the analytic solution for the
fields from a magnetic dipole. No discrete curl is taken, so the magnetic
flux density may not be strictly divergence free.
This approach uses a primary-secondary in frequency in the same fashion as the MagDipole.
:param list rxList: receiver list
:param float freq: frequency
:param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`)
:param string orientation: 'X', 'Y', 'Z'
:param float moment: magnetic dipole moment
:param float mu: background magnetic permeability
"""
#TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
#TODO: neither does moment
def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0):
@@ -190,6 +374,14 @@ class MagDipole_Bfield(BaseSrc):
BaseSrc.__init__(self, rxList)
def bPrimary(self, prob):
"""
The primary magnetic flux density from the analytic solution for magnetic fields from a dipole
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
eqLocs = prob._eqLocs
if eqLocs is 'FE':
@@ -221,14 +413,35 @@ class MagDipole_Bfield(BaseSrc):
return b
def hPrimary(self, prob):
"""
The primary magnetic field from a magnetic vector potential
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
b = self.bPrimary(prob)
return h_from_b(prob, b)
def S_m(self, prob):
"""
The magnetic source term
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
b = self.bPrimary(prob)
return -1j*omega(self.freq)*b
def S_e(self, prob):
"""
The electric source term
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
return Zero()
else:
@@ -247,6 +460,20 @@ class MagDipole_Bfield(BaseSrc):
class CircularLoop(BaseSrc):
"""
Circular loop magnetic source calculated by taking the curl of a magnetic
vector potential. By taking the discrete curl, we ensure that the magnetic
flux density is divergence free (no magnetic monopoles!).
This approach uses a primary-secondary in frequency in the same fashion as the MagDipole.
:param list rxList: receiver list
:param float freq: frequency
:param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`)
:param string orientation: 'X', 'Y', 'Z'
:param float moment: magnetic dipole moment
:param float mu: background magnetic permeability
"""
#TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
def __init__(self, rxList, freq, loc, orientation='Z', radius = 1., mu=mu_0):
@@ -259,6 +486,13 @@ class CircularLoop(BaseSrc):
BaseSrc.__init__(self, rxList)
def bPrimary(self, prob):
"""
The primary magnetic flux density from a magnetic vector potential
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
eqLocs = prob._eqLocs
if eqLocs is 'FE':
@@ -289,14 +523,35 @@ class CircularLoop(BaseSrc):
return C*a
def hPrimary(self, prob):
"""
The primary magnetic field from a magnetic vector potential
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
b = self.bPrimary(prob)
return 1./self.mu*b
def S_m(self, prob):
"""
The magnetic source term
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
b = self.bPrimary(prob)
return -1j*omega(self.freq)*b
def S_e(self, prob):
"""
The electric source term
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
return Zero()
else:
SimPEG/EM/FDEM/SurveyFDEM.py
+42 -2
View File
@@ -10,6 +10,12 @@ import SrcFDEM as Src
####################################################
class Rx(SimPEG.Survey.BaseRx):
"""
Frequency domain receivers
:param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`)
:param string rxType: reciever type from knownRxTypes
"""
knownRxTypes = {
'exr':['e', 'Ex', 'real'],
@@ -61,6 +67,15 @@ class Rx(SimPEG.Survey.BaseRx):
return self.knownRxTypes[self.rxType][2]
def projectFields(self, src, mesh, u):
"""
Project fields to recievers to get data.
:param Source src: FDEM source
:param Mesh mesh: mesh used
:param Fields u: fields object
:rtype: numpy.ndarray
:return: fields projected to recievers
"""
P = self.getP(mesh)
u_part_complex = u[src, self.projField]
# get the real or imag component
@@ -69,6 +84,16 @@ class Rx(SimPEG.Survey.BaseRx):
return P*u_part
def projectFieldsDeriv(self, src, mesh, u, v, adjoint=False):
"""
Derivative of projected fields with respect to the inversion model times a vector.
:param Source src: FDEM source
:param Mesh mesh: mesh used
:param Fields u: fields object
:param numpy.ndarray v: vector to multiply
:rtype: numpy.ndarray
:return: fields projected to recievers
"""
P = self.getP(mesh)
if not adjoint:
@@ -95,10 +120,13 @@ class Rx(SimPEG.Survey.BaseRx):
class Survey(SimPEG.Survey.BaseSurvey):
"""
docstring for SurveyFDEM
Frequency domain electromagnetic survey
:param list srcList: list of FDEM sources used in the survey
"""
srcPair = Src.BaseSrc
rxPaair = Rx
def __init__(self, srcList, **kwargs):
# Sort these by frequency
@@ -126,6 +154,7 @@ class Survey(SimPEG.Survey.BaseSurvey):
@property
def nSrcByFreq(self):
"""Number of sources at each frequency"""
if getattr(self, '_nSrcByFreq', None) is None:
self._nSrcByFreq = {}
for freq in self.freqs:
@@ -133,11 +162,22 @@ class Survey(SimPEG.Survey.BaseSurvey):
return self._nSrcByFreq
def getSrcByFreq(self, freq):
"""Returns the sources associated with a specific frequency."""
"""
Returns the sources associated with a specific frequency.
:param float freq: frequency for which we look up sources
:rtype: dictionary
:return: sources at the sepcified frequency
"""
assert freq in self._freqDict, "The requested frequency is not in this survey."
return self._freqDict[freq]
def projectFields(self, u):
"""
Project fields to receiver locations
:param Fields u: fields object
:rtype: numpy.ndarray
:return: data
"""
data = SimPEG.Survey.Data(self)
for src in self.srcList:
for rx in src.rxList:
SimPEG/EM/TDEM/BaseTDEM.py
+9 -1
View File
@@ -31,12 +31,20 @@ class FieldsTDEM(Problem.TimeFields):
class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
"""docstring for ProblemTDEM1D"""
"""docstring for BaseTDEMProblem"""
def __init__(self, mesh, mapping=None, **kwargs):
BaseTimeProblem.__init__(self, mesh, mapping=mapping, **kwargs)
_FieldsForward_pair = FieldsTDEM #: used for the forward calculation only
waveformType = "STEPOFF"
current = None
def currentwaveform(self, wave):
self._timeSteps = np.diff(wave[:,0])
self.current = wave[:,1]
self.waveformType = "GENERAL"
def fields(self, m):
if self.verbose: print '%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50)
self.curModel = m
SimPEG/EM/TDEM/SurveyTDEM.py
+24 -6
View File
@@ -99,14 +99,33 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
class SrcTDEM_CircularLoop_MVP(SrcTDEM):
def __init__(self,rxList,loc,radius):
def __init__(self,rxList,loc,radius,waveformType):
self.loc = loc
self.radius = radius
SrcTDEM.__init__(self,rxList)
self.waveformType = waveformType
SrcTDEM.__init__(self,rxList)
def getInitialFields(self, mesh):
"""Circular Loop, magnetic vector potential"""
if self.waveformType == "STEPOFF":
print ">> Step waveform: Non-zero initial condition"
if mesh._meshType is 'CYL':
if mesh.isSymmetric:
MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
else:
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
elif mesh._meshType is 'TENSOR':
MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
else:
raise Exception('Unknown mesh for CircularLoop')
return {"b": mesh.edgeCurl*MVP}
elif self.waveformType == "GENERAL":
print ">> General waveform: Zero initial condition"
return {"b": np.zeros(mesh.nF)}
else:
raise NotImplementedError("Only use STEPOFF or GENERAL")
def getMeS(self, mesh, MfMui):
if mesh._meshType is 'CYL':
if mesh.isSymmetric:
MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
@@ -115,9 +134,8 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM):
elif mesh._meshType is 'TENSOR':
MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
else:
raise Exception('Unknown mesh for CircularLoop')
return {"b": mesh.edgeCurl*MVP}
raise Exception('Unknown mesh for CircularLoop')
return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
class SurveyTDEM(Survey.BaseSurvey):
SimPEG/EM/__init__.py
+1 -1
View File
@@ -1,6 +1,6 @@
# from EM import *
import TDEM
import FDEM
import Base
import Analytics
import Utils
from scipy.constants import mu_0, epsilon_0
SimPEG/Examples/EM_FDEM_1D_Inversion.py
+116
View File
@@ -0,0 +1,116 @@
from SimPEG import *
import SimPEG.EM as EM
from SimPEG.EM import mu_0
def run(plotIt=True):
"""
EM: FDEM: 1D: Inversion
=======================
Here we will create and run a FDEM 1D inversion.
"""
cs, ncx, ncz, npad = 5., 25, 15, 15
hx = [(cs,ncx), (cs,npad,1.3)]
hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
mesh = Mesh.CylMesh([hx,1,hz], '00C')
layerz = -100.
active = mesh.vectorCCz<0.
layer = (mesh.vectorCCz<0.) & (mesh.vectorCCz>=layerz)
actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap
sig_half = 2e-2
sig_air = 1e-8
sig_layer = 1e-2
sigma = np.ones(mesh.nCz)*sig_air
sigma[active] = sig_half
sigma[layer] = sig_layer
mtrue = np.log(sigma[active])
if plotIt:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1, figsize = (3, 6))
plt.semilogx(sigma[active], mesh.vectorCCz[active])
ax.set_ylim(-500, 0)
ax.set_xlim(1e-3, 1e-1)
ax.set_xlabel('Conductivity (S/m)', fontsize = 14)
ax.set_ylabel('Depth (m)', fontsize = 14)
ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
rxOffset=10.
bzi = EM.FDEM.Rx(np.array([[rxOffset, 0., 1e-3]]), 'bzi')
freqs = np.logspace(1,3,10)
srcLoc = np.array([0., 0., 10.])
srcList = []
[srcList.append(EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z')) for freq in freqs]
survey = EM.FDEM.Survey(srcList)
prb = EM.FDEM.Problem_b(mesh, mapping=mapping)
try:
from pymatsolver import MumpsSolver
prb.Solver = MumpsSolver
except ImportError, e:
prb.Solver = SolverLU
prb.pair(survey)
std = 0.05
survey.makeSyntheticData(mtrue, std)
survey.std = std
survey.eps = np.linalg.norm(survey.dtrue)*1e-5
if plotIt:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1, figsize = (6, 6))
ax.semilogx(freqs,survey.dtrue[:freqs.size], 'b.-')
ax.semilogx(freqs,survey.dobs[:freqs.size], 'r.-')
ax.legend(('Noisefree', '$d^{obs}$'), fontsize = 16)
ax.set_xlabel('Time (s)', fontsize = 14)
ax.set_ylabel('$B_z$ (T)', fontsize = 16)
ax.set_xlabel('Time (s)', fontsize = 14)
ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
dmisfit = DataMisfit.l2_DataMisfit(survey)
regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
reg = Regularization.Tikhonov(regMesh)
opt = Optimization.InexactGaussNewton(maxIter = 6)
invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
# Create an inversion object
beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=2)
betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e0)
inv = Inversion.BaseInversion(invProb, directiveList=[beta,betaest])
m0 = np.log(np.ones(mtrue.size)*sig_half)
reg.alpha_s = 1e-3
reg.alpha_x = 1.
prb.counter = opt.counter = Utils.Counter()
opt.LSshorten = 0.5
opt.remember('xc')
mopt = inv.run(m0)
if plotIt:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1, figsize = (3, 6))
plt.semilogx(sigma[active], mesh.vectorCCz[active])
plt.semilogx(np.exp(mopt), mesh.vectorCCz[active])
ax.set_ylim(-500, 0)
ax.set_xlim(1e-3, 1e-1)
ax.set_xlabel('Conductivity (S/m)', fontsize = 14)
ax.set_ylabel('Depth (m)', fontsize = 14)
ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
plt.legend(['$\sigma_{true}$', '$\sigma_{pred}$'],loc='best')
plt.show()
if __name__ == '__main__':
run()
SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py
+43
View File
@@ -0,0 +1,43 @@
from SimPEG import *
import SimPEG.EM as EM
def run(XYZ=None, loc=np.r_[0.,0.,0.], sig=1.0, freq=1.0, orientation='Z', plotIt=True):
"""
EM: Magnetic Dipole in a Whole-Space
====================================
Here we plot the magnetic flux density from a harmonic dipole in a wholespace.
"""
if XYZ is None:
x = np.arange(-100.5,100.5,step = 1.) #(avoid putting measurement points where source is located)
y = np.r_[0]
z = x
XYZ = Utils.ndgrid(x,y,z)
Bx, By, Bz = EM.Analytics.FDEM.MagneticDipoleWholeSpace(XYZ, loc, sig, freq, orientation=orientation)
absB = np.sqrt(Bx*Bx.conj()+By*By.conj()+Bz*Bz.conj()).real
if plotIt:
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
fig, ax = plt.subplots(1,1,figsize=(6,5))
bxplt = Bx.reshape(x.size,z.size)
bzplt = Bz.reshape(x.size,z.size)
pc = ax.pcolor(x,z,absB.reshape(x.size,z.size),norm=LogNorm())
ax.streamplot(x,z,bxplt.real,bzplt.real,color='k',density=1)
ax.set_xlim([x.min(),x.max()])
ax.set_ylim([z.min(),z.max()])
ax.set_xlabel('x')
ax.set_ylabel('z')
cb = plt.colorbar(pc,ax = ax)
cb.set_label('|B| (T)')
plt.show()
return fig, ax
if __name__ == '__main__':
run()
SimPEG/EM/Examples/CylInversion.py → SimPEG/Examples/EM_TDEM_1D_Inversion.py
+19 -10
View File
@@ -1,9 +1,16 @@
from SimPEG import *
import SimPEG.EM as EM
from scipy.constants import mu_0
import matplotlib.pyplot as plt
from SimPEG.EM import mu_0
def run(plotIt=True):
"""
EM: TDEM: 1D: Inversion
=======================
Here we will create and run a TDEM 1D inversion.
"""
cs, ncx, ncz, npad = 5., 25, 15, 15
hx = [(cs,ncx), (cs,npad,1.3)]
@@ -24,6 +31,7 @@ def run(plotIt=True):
if plotIt:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1, figsize = (3, 6))
plt.semilogx(sigma[active], mesh.vectorCCz[active])
ax.set_ylim(-600, 0)
@@ -42,19 +50,18 @@ def run(plotIt=True):
prb.Solver = SolverLU
prb.timeSteps = [(1e-06, 20),(1e-05, 20), (0.0001, 20)]
prb.pair(survey)
dtrue = survey.dpred(mtrue)
survey.dtrue = dtrue
# create observed data
std = 0.05
noise = std*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
survey.dobs = survey.dtrue+noise
survey.std = survey.dobs*0 + std
survey.Wd = 1/(abs(survey.dobs)*std)
survey.dobs = survey.makeSyntheticData(mtrue,std)
survey.std = std
survey.eps = 1e-5*np.linalg.norm(survey.dobs)
if plotIt:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1, figsize = (10, 6))
ax.loglog(rx.times, dtrue, 'b.-')
ax.loglog(rx.times, survey.dtrue, 'b.-')
ax.loglog(rx.times, survey.dobs, 'r.-')
ax.legend(('Noisefree', '$d^{obs}$'), fontsize = 16)
ax.set_xlabel('Time (s)', fontsize = 14)
@@ -67,6 +74,7 @@ def run(plotIt=True):
reg = Regularization.Tikhonov(regMesh)
opt = Optimization.InexactGaussNewton(maxIter = 5)
invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
# Create an inversion object
beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=2)
betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e0)
@@ -81,6 +89,7 @@ def run(plotIt=True):
mopt = inv.run(m0)
if plotIt:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1, figsize = (3, 6))
plt.semilogx(sigma[active], mesh.vectorCCz[active])
plt.semilogx(np.exp(mopt), mesh.vectorCCz[active])
SimPEG/FLOW/Examples/Celia1990.py → SimPEG/Examples/FLOW_Richards_1D_Celia1990.py
+35 -1
View File
@@ -1,8 +1,40 @@
from SimPEG import *
from SimPEG.FLOW import Richards
import matplotlib.pyplot as plt
def run(plotIt=True):
"""
FLOW: Richards: 1D: Celia1990
=============================
There are two different forms of Richards equation that differ
on how they deal with the non-linearity in the time-stepping term.
The most fundamental form, referred to as the
'mixed'-form of Richards Equation Celia1990_
.. math::
\\frac{\partial \\theta(\psi)}{\partial t} - \\nabla \cdot k(\psi) \\nabla \psi - \\frac{\partial k(\psi)}{\partial z} = 0
\quad \psi \in \Omega
where \\\\(\\\\theta\\\\) is water content, and \\\\(\\\\psi\\\\) is pressure head.
This formulation of Richards equation is called the
'mixed'-form because the equation is parameterized in \\\\(\\\\psi\\\\)
but the time-stepping is in terms of \\\\(\\\\theta\\\\).
As noted in Celia1990_ the 'head'-based form of Richards
equation can be written in the continuous form as:
.. math::
\\frac{\partial \\theta}{\partial \psi}\\frac{\partial \psi}{\partial t} - \\nabla \cdot k(\psi) \\nabla \psi - \\frac{\partial k(\psi)}{\partial z} = 0 \quad \psi \in \Omega
However, it can be shown that this does not conserve mass in the discrete formulation.
Here we reproduce the results from Celia1990_ demonstrating the head-based formulation and the mixed-formulation.
.. _Celia1990: http://www.webpages.uidaho.edu/ch/papers/Celia.pdf
"""
M = Mesh.TensorMesh([np.ones(40)])
M.setCellGradBC('dirichlet')
params = Richards.Empirical.HaverkampParams().celia1990
@@ -28,6 +60,7 @@ def run(plotIt=True):
Hs_H120= getFields(120.,'head')
if not plotIt:return
import matplotlib.pyplot as plt
plt.figure(figsize=(13,5))
plt.subplot(121)
plt.plot(40-M.gridCC, Hs_M10[-1],'b-')
@@ -47,6 +80,7 @@ def run(plotIt=True):
plt.xlabel('Depth, cm')
plt.ylabel('Pressure Head, cm')
plt.legend(('$\Delta t$ = 10 sec','$\Delta t$ = 30 sec','$\Delta t$ = 120 sec'))
plt.show()
if __name__ == '__main__':
run()
SimPEG/Examples/DCfwd.py → SimPEG/Examples/Forward_BasicDirectCurrent.py
+6 -4
View File
@@ -1,7 +1,4 @@
from SimPEG import Mesh, Utils, np, SolverLU
import matplotlib.pyplot as plt
import matplotlib
from matplotlib.mlab import griddata
## 2D DC forward modeling example with Tensor and Curvilinear Meshes
@@ -12,7 +9,6 @@ def run(plotIt=True):
tM = Mesh.TensorMesh(sz)
# Curvilinear Mesh
rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
# Step2: Direct Current (DC) operator
def DCfun(mesh, pts):
D = mesh.faceDiv
@@ -39,6 +35,11 @@ def run(plotIt=True):
phirM = AinvrM*rhsrM
if not plotIt: return
import matplotlib.pyplot as plt
import matplotlib
from matplotlib.mlab import griddata
#Step4: Making Figure
fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
label = ["(a)", "(b)"]
@@ -69,6 +70,7 @@ def run(plotIt=True):
else:
axes[i].set_ylabel(" ")
axes[i].set_xlabel("x")
plt.show()
if __name__ == '__main__':
SimPEG/Examples/Linear.py → SimPEG/Examples/Inversion_Linear.py
+36 -24
View File
@@ -1,29 +1,39 @@
from SimPEG import *
class LinearSurvey(Survey.BaseSurvey):
def projectFields(self, u):
return u
class LinearProblem(Problem.BaseProblem):
"""docstring for LinearProblem"""
def run(N=100, plotIt=True):
"""
Inversion: Linear Problem
=========================
surveyPair = LinearSurvey
Here we go over the basics of creating a linear problem and inversion.
def __init__(self, mesh, G, **kwargs):
Problem.BaseProblem.__init__(self, mesh, **kwargs)
self.G = G
"""
def fields(self, m, u=None):
return self.G.dot(m)
class LinearSurvey(Survey.BaseSurvey):
def projectFields(self, u):
return u
def Jvec(self, m, v, u=None):
return self.G.dot(v)
class LinearProblem(Problem.BaseProblem):
def Jtvec(self, m, v, u=None):
return self.G.T.dot(v)
surveyPair = LinearSurvey
def __init__(self, mesh, G, **kwargs):
Problem.BaseProblem.__init__(self, mesh, **kwargs)
self.G = G
def fields(self, m, u=None):
return self.G.dot(m)
def Jvec(self, m, v, u=None):
return self.G.dot(v)
def Jtvec(self, m, v, u=None):
return self.G.T.dot(v)
def run(N, plotIt=True):
np.random.seed(1)
mesh = Mesh.TensorMesh([N])
nk = 20
@@ -52,7 +62,7 @@ def run(N, plotIt=True):
reg = Regularization.Tikhonov(mesh)
dmis = DataMisfit.l2_DataMisfit(survey)
opt = Optimization.InexactGaussNewton(maxIter=20)
opt = Optimization.InexactGaussNewton(maxIter=35)
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
beta = Directives.BetaSchedule()
betaest = Directives.BetaEstimate_ByEig()
@@ -63,16 +73,18 @@ def run(N, plotIt=True):
if plotIt:
import matplotlib.pyplot as plt
plt.figure(1)
for i in range(prob.G.shape[0]):
plt.plot(prob.G[i,:])
plt.figure(2)
plt.plot(M.vectorCCx, survey.mtrue, 'b-')
plt.plot(M.vectorCCx, mrec, 'r-')
fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
for i in range(prob.G.shape[0]):
axes[0].plot(prob.G[i,:])
axes[0].set_title('Columns of matrix G')
axes[1].plot(M.vectorCCx, survey.mtrue, 'b-')
axes[1].plot(M.vectorCCx, mrec, 'r-')
axes[1].legend(('True Model', 'Recovered Model'))
plt.show()
return prob, survey, mesh, mrec
if __name__ == '__main__':
run(100)
run()
SimPEG/Examples/Mesh_Basic_PlotImage.py
+46
View File
@@ -0,0 +1,46 @@
from SimPEG import *
def run(plotIt=True):
"""
Mesh: Basic: PlotImage
======================
You can use M.PlotImage to plot images on all of the Meshes.
"""
M = Mesh.TensorMesh([32,32])
v = Utils.ModelBuilder.randomModel(M.vnC, seed=789)
v = Utils.mkvc(v)
O = Mesh.TreeMesh([32,32])
O.refine(1)
def function(cell):
if (cell.center[0] < 0.75 and cell.center[0] > 0.25 and
cell.center[1] < 0.75 and cell.center[1] > 0.25):return 5
if (cell.center[0] < 0.9 and cell.center[0] > 0.1 and
cell.center[1] < 0.9 and cell.center[1] > 0.1):return 4
return 3
O.refine(function)
P = M.getInterpolationMat(O.gridCC, 'CC')
ov = P * v
if plotIt:
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1,2,figsize=(10,5))
out = M.plotImage(v, grid=True, ax=axes[0])
cb = plt.colorbar(out[0], ax=axes[0]); cb.set_label("Random Field")
axes[0].set_title('TensorMesh')
out = O.plotImage(ov, grid=True, ax=axes[1], clim=[0,1])
cb = plt.colorbar(out[0], ax=axes[1]); cb.set_label("Random Field")
axes[1].set_title('TreeMesh')
plt.show()
if __name__ == '__main__':
run()
SimPEG/Examples/Mesh_Basic_Types.py
+30
View File
@@ -0,0 +1,30 @@
from SimPEG import *
def run(plotIt=True):
"""
Mesh: Basic: Types
==================
Here we show SimPEG used to create three different types of meshes.
"""
sz = [16,16]
tM = Mesh.TensorMesh(sz)
qM = Mesh.TreeMesh(sz)
qM.refine(lambda cell: 4 if np.sqrt(((np.r_[cell.center]-0.5)**2).sum()) < 0.4 else 3)
rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
if plotIt:
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1,3,figsize=(14,5))
opts = {}
tM.plotGrid(ax=axes[0], **opts)
axes[0].set_title('TensorMesh')
qM.plotGrid(ax=axes[1], **opts)
axes[1].set_title('TreeMesh')
rM.plotGrid(ax=axes[2], **opts)
axes[2].set_title('CurvilinearMesh')
plt.show()
if __name__ == '__main__':
run()
SimPEG/Examples/Mesh_Operators_CahnHilliard.py
+105
View File
@@ -0,0 +1,105 @@
from SimPEG import *
def run(plotIt=True, n=60):
"""
Mesh: Operators: Cahn Hilliard
==============================
This example is based on the example in the FiPy_ library.
Please see their documentation for more information about the Cahn-Hilliard equation.
The "Cahn-Hilliard" equation separates a field \\\\( \\\\phi \\\\) into 0 and 1 with smooth transitions.
.. math::
\\frac{\partial \phi}{\partial t} = \\nabla \cdot D \\nabla \left( \\frac{\partial f}{\partial \phi} - \epsilon^2 \\nabla^2 \phi \\right)
Where \\\\( f \\\\) is the energy function \\\\( f = ( a^2 / 2 )\\\\phi^2(1 - \\\\phi)^2 \\\\)
which drives \\\\( \\\\phi \\\\) towards either 0 or 1, this competes with the term
\\\\(\\\\epsilon^2 \\\\nabla^2 \\\\phi \\\\) which is a diffusion term that creates smooth changes in \\\\( \\\\phi \\\\).
The equation can be factored:
.. math::
\\frac{\partial \phi}{\partial t} = \\nabla \cdot D \\nabla \psi \\\\
\psi = \\frac{\partial^2 f}{\partial \phi^2} (\phi - \phi^{\\text{old}}) + \\frac{\partial f}{\partial \phi} - \epsilon^2 \\nabla^2 \phi
Here we will need the derivatives of \\\\( f \\\\):
.. math::
\\frac{\partial f}{\partial \phi} = (a^2/2)2\phi(1-\phi)(1-2\phi)
\\frac{\partial^2 f}{\partial \phi^2} = (a^2/2)2[1-6\phi(1-\phi)]
The implementation below uses backwards Euler in time with an exponentially increasing time step.
The initial \\\\( \\\\phi \\\\) is a normally distributed field with a standard deviation of 0.1 and mean of 0.5.
The grid is 60x60 and takes a few seconds to solve ~130 times. The results are seen below, and you can see the
field separating as the time increases.
.. _FiPy: http://www.ctcms.nist.gov/fipy/examples/cahnHilliard/generated/examples.cahnHilliard.mesh2DCoupled.html
"""
np.random.seed(5)
# Here we are going to rearrange the equations:
# (phi_ - phi)/dt = A*(d2fdphi2*(phi_ - phi) + dfdphi - L*phi_)
# (phi_ - phi)/dt = A*(d2fdphi2*phi_ - d2fdphi2*phi + dfdphi - L*phi_)
# (phi_ - phi)/dt = A*d2fdphi2*phi_ + A*( - d2fdphi2*phi + dfdphi - L*phi_)
# phi_ - phi = dt*A*d2fdphi2*phi_ + dt*A*(- d2fdphi2*phi + dfdphi - L*phi_)
# phi_ - dt*A*d2fdphi2 * phi_ = dt*A*(- d2fdphi2*phi + dfdphi - L*phi_) + phi
# (I - dt*A*d2fdphi2) * phi_ = dt*A*(- d2fdphi2*phi + dfdphi - L*phi_) + phi
# (I - dt*A*d2fdphi2) * phi_ = dt*A*dfdphi - dt*A*d2fdphi2*phi - dt*A*L*phi_ + phi
# (dt*A*d2fdphi2 - I) * phi_ = dt*A*d2fdphi2*phi + dt*A*L*phi_ - phi - dt*A*dfdphi
# (dt*A*d2fdphi2 - I - dt*A*L) * phi_ = (dt*A*d2fdphi2 - I)*phi - dt*A*dfdphi
h = [(0.25,n)]
M = Mesh.TensorMesh([h,h])
# Constants
D = a = epsilon = 1.
I = Utils.speye(M.nC)
# Operators
A = D * M.faceDiv * M.cellGrad
L = epsilon**2 * M.faceDiv * M.cellGrad
duration = 75
elapsed = 0.
dexp = -5
phi = np.random.normal(loc=0.5,scale=0.01,size=M.nC)
ii, jj = 0, 0
PHIS = []
capture = np.logspace(-1,np.log10(duration),8)
while elapsed < duration:
dt = min(100, np.exp(dexp))
elapsed += dt
dexp += 0.05
dfdphi = a**2 * 2 * phi * (1 - phi) * (1 - 2 * phi)
d2fdphi2 = Utils.sdiag(a**2 * 2 * (1 - 6 * phi * (1 - phi)))
MAT = (dt*A*d2fdphi2 - I - dt*A*L)
rhs = (dt*A*d2fdphi2 - I)*phi - dt*A*dfdphi
phi = Solver(MAT)*rhs
if elapsed > capture[jj]:
PHIS += [(elapsed, phi.copy())]
jj += 1
if ii % 10 == 0: print ii, elapsed
ii += 1
if plotIt:
import matplotlib.pyplot as plt
fig, axes = plt.subplots(2,4,figsize=(14,6))
axes = np.array(axes).flatten().tolist()
for ii, ax in zip(np.linspace(0,len(PHIS)-1,len(axes)),axes):
ii = int(ii)
out = M.plotImage(PHIS[ii][1],ax=ax)
ax.axis('off')
ax.set_title('Elapsed Time: %4.1f'%PHIS[ii][0])
plt.show()
if __name__ == '__main__':
run()
SimPEG/Examples/Mesh_QuadTree_Creation.py
+28
View File
@@ -0,0 +1,28 @@
from SimPEG import *
def run(plotIt=True):
"""
Mesh: QuadTree: Creation
========================
You can give the refine method a function, which is evaluated on every cell
of the TreeMesh.
Occasionally it is useful to initially refine to a constant level
(e.g. 3 in this 32x32 mesh). This means the function is first evaluated
on an 8x8 mesh (2^3).
"""
M = Mesh.TreeMesh([32,32])
M.refine(3)
def function(cell):
xyz = cell.center
for i in range(3):
if np.abs(np.sin(xyz[0]*np.pi*2)*0.5 + 0.5 - xyz[1]) < 0.2*i:
return 6-i
return 0
M.refine(function);
if plotIt: M.plotGrid(showIt=True)
if __name__ == '__main__':
run()
SimPEG/Examples/Mesh_QuadTree_FaceDiv.py
+49
View File
@@ -0,0 +1,49 @@
from SimPEG import *
def run(plotIt=True, n=60):
"""
Mesh: QuadTree: FaceDiv
=======================
"""
M = Mesh.TreeMesh([[(1,16)],[(1,16)]], levels=4)
M._refineCell([0,0,0])
M._refineCell([0,0,1])
M._refineCell([4,4,2])
M.__dirty__ = True
M.number()
if plotIt:
import matplotlib.pyplot as plt
fig, axes = plt.subplots(2,1,figsize=(10,10))
M.plotGrid(cells=True, nodes=False, ax=axes[0])
axes[0].axis('off')
axes[0].set_title('Simple QuadTree Mesh')
axes[0].set_xlim([-1,17])
axes[0].set_ylim([-1,17])
for ii, loc in zip(range(M.nC),M.gridCC):
axes[0].text(loc[0]+0.2,loc[1],'%d'%ii, color='r')
axes[0].plot(M.gridFx[:,0],M.gridFx[:,1], 'g>')
for ii, loc in zip(range(M.nFx),M.gridFx):
axes[0].text(loc[0]+0.2,loc[1],'%d'%ii, color='g')
axes[0].plot(M.gridFy[:,0],M.gridFy[:,1], 'm^')
for ii, loc in zip(range(M.nFy),M.gridFy):
axes[0].text(loc[0]+0.2,loc[1]+0.2,'%d'%(ii+M.nFx), color='m')
axes[1].spy(M.faceDiv)
axes[1].set_title('Face Divergence')
axes[1].set_ylabel('Cell Number')
axes[1].set_xlabel('Face Number')
plt.show()
if __name__ == '__main__':
run()
SimPEG/Examples/Mesh_QuadTree_HangingNodes.py
+32
View File
@@ -0,0 +1,32 @@
from SimPEG import *
def run(plotIt=True):
"""
Mesh: QuadTree: Hanging Nodes
=============================
You can give the refine method a function, which is evaluated on every cell
of the TreeMesh.
Occasionally it is useful to initially refine to a constant level
(e.g. 3 in this 32x32 mesh). This means the function is first evaluated
on an 8x8 mesh (2^3).
"""
M = Mesh.TreeMesh([8,8])
def function(cell):
xyz = cell.center
dist = ((xyz - [0.25,0.25])**2).sum()**0.5
if dist < 0.25:
return 3
return 2
M.refine(function);
M.number()
if plotIt:
import matplotlib.pyplot as plt
M.plotGrid(nodes=True, cells=True, facesX=True)
plt.legend(('Grid', 'Cell Centers', 'Nodes', 'Hanging Nodes', 'X faces', 'Hanging X faces'))
plt.show()
if __name__ == '__main__':
run()
SimPEG/Examples/Mesh_Tensor_Creation.py
+35
View File
@@ -0,0 +1,35 @@
from SimPEG import *
def run(plotIt=True):
"""
Mesh: Tensor: Creation
======================
For tensor meshes, there are some functions that can come
in handy. For example, creating mesh tensors can be a bit time
consuming, these can be created speedily by just giving numbers
and sizes of padding. See the example below, that follows this
notation::
h1 = (
(cellSize, numPad, [, increaseFactor]),
(cellSize, numCore),
(cellSize, numPad, [, increaseFactor])
)
.. note::
You can center your mesh by passing a 'C' for the x0[i] position.
A 'N' will make the entire mesh negative, and a '0' (or a 0) will
make the mesh start at zero.
"""
h1 = [(10, 5, -1.3), (5, 20), (10, 3, 1.3)]
M = Mesh.TensorMesh([h1, h1], x0='CN')
if plotIt:
M.plotGrid(showIt=True)
if __name__ == '__main__':
run()
SimPEG/Examples/__init__.py
+106 -1
View File
@@ -1 +1,106 @@
import Linear, DCfwd
# Run this file to add imports.
##### AUTOIMPORTS #####
import EM_FDEM_1D_Inversion
import EM_FDEM_Analytic_MagDipoleWholespace
import EM_TDEM_1D_Inversion
import FLOW_Richards_1D_Celia1990
import Forward_BasicDirectCurrent
import Inversion_Linear
import Mesh_Basic_PlotImage
import Mesh_Basic_Types
import Mesh_Operators_CahnHilliard
import Mesh_QuadTree_Creation
import Mesh_QuadTree_FaceDiv
import Mesh_QuadTree_HangingNodes
import Mesh_Tensor_Creation
__examples__ = ["EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation"]
##### AUTOIMPORTS #####
if __name__ == '__main__':
"""
Run the following to create the examples documentation and add to the imports at the top.
"""
import shutil, os
from SimPEG import Examples
# Create the examples dir in the docs folder.
fName = os.path.realpath(__file__)
docExamplesDir = os.path.sep.join(fName.split(os.path.sep)[:-3] + ['docs', 'examples'])
shutil.rmtree(docExamplesDir)
os.makedirs(docExamplesDir)
# Get all the python examples in this folder
thispath = os.path.sep.join(fName.split(os.path.sep)[:-1])
exfiles = [f[:-3] for f in os.listdir(thispath) if os.path.isfile(os.path.join(thispath, f)) and f.endswith('.py') and not f.startswith('_')]
# Add the imports to the top in the AUTOIMPORTS section
f = file(fName, 'r')
inimports = False
out = ''
for line in f:
if not inimports:
out += line
if line == "##### AUTOIMPORTS #####\n":
inimports = not inimports
if inimports:
out += '\n'.join(["import %s"%_ for _ in exfiles])
out += '\n\n__examples__ = ["' + '", "'.join(exfiles)+ '"]\n'
out += '\n##### AUTOIMPORTS #####\n'
f.close()
f = file(fName, 'w')
f.write(out)
f.close()
def _makeExample(filePath, runFunction):
"""Makes the example given a path of the file and the run function."""
filePath = os.path.realpath(filePath)
name = filePath.split(os.path.sep)[-1].rstrip('.pyc').rstrip('.py')
docstr = runFunction.__doc__
if docstr is None:
doc = '%s\n%s'%(name.replace('_',' '),'='*len(name))
else:
doc = '\n'.join([_[8:].rstrip() for _ in docstr.split('\n')])
out = """.. _examples_%s:
.. --------------------------------- ..
.. ..
.. THIS FILE IS AUTO GENEREATED ..
.. ..
.. SimPEG/Examples/__init__.py ..
.. ..
.. --------------------------------- ..
%s
.. plot::
from SimPEG import Examples
Examples.%s.run()
.. literalinclude:: ../../SimPEG/Examples/%s.py
:language: python
:linenos:
"""%(name,doc,name,name)
rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'examples', name + '.rst']))
print 'Creating: %s.rst'%name
f = open(rst, 'w')
f.write(out)
f.close()
for ex in dir(Examples):
if ex.startswith('_'): continue
E = getattr(Examples,ex)
_makeExample(E.__file__, E.run)
SimPEG/FLOW/Examples/__init__.py
-1
View File
@@ -1 +0,0 @@
import Celia1990
SimPEG/InvProblem.py
+2 -2
View File
@@ -66,8 +66,8 @@ class BaseInvProblem(object):
self.curModel = m0
print """SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
***Done using same solver as the problem***"""
self.opt.bfgsH0 = self.prob.Solver(self.reg.eval2Deriv(self.curModel))
***Done using same Solver and solverOpts as the problem***"""
self.opt.bfgsH0 = self.prob.Solver(self.reg.eval2Deriv(self.curModel), **self.prob.solverOpts)
@property
def warmstart(self):
SimPEG/Maps.py
+54 -45
View File
@@ -10,21 +10,25 @@ class IdentityMap(object):
SimPEG Map
"""
__metaclass__ = Utils.SimPEGMetaClass
mesh = None #: A SimPEG Mesh
def __init__(self, mesh, **kwargs):
def __init__(self, mesh=None, nP=None, **kwargs):
Utils.setKwargs(self, **kwargs)
if nP is not None:
assert type(nP) in [int, long], ' Number of parameters must be an integer.'
self.mesh = mesh
self._nP = nP
@property
def nP(self):
"""
:rtype: int
:return: number of parameters in the model
:return: number of parameters that the mapping accepts
"""
if self._nP is not None:
return self._nP
if self.mesh is None:
return '*'
return self.mesh.nC
@@ -32,11 +36,15 @@ class IdentityMap(object):
@property
def shape(self):
"""
The default shape is (mesh.nC, nP).
The default shape is (mesh.nC, nP) if the mesh is defined.
If this is a meshless mapping (i.e. nP is defined independently)
the shape will be the the shape (nP,nP).
:rtype: (int,int)
:return: shape of the operator as a tuple
"""
if self._nP is not None:
return (self.nP, self.nP)
if self.mesh is None:
return ('*', self.nP)
return (self.mesh.nC, self.nP)
@@ -118,6 +126,7 @@ class IdentityMap(object):
def __str__(self):
return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1])
class ComboMap(IdentityMap):
"""Combination of various maps."""
@@ -475,7 +484,7 @@ class ActiveCells(IdentityMap):
else:
self.valInactive = valInactive.copy()
self.valInactive[self.indActive] = 0
inds = np.nonzero(self.indActive)[0]
self.P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), shape=(self.nC, self.nP))
@@ -708,7 +717,7 @@ class PolyMap(IdentityMap):
Parameterize the model space using a polynomials in a wholespace.
..math::
y = \mathbf{V} c
Define the model as:
@@ -752,10 +761,10 @@ class PolyMap(IdentityMap):
else:
raise(Exception("Input for normal = X or Y or Z"))
#3D
elif self.mesh.dim == 3:
elif self.mesh.dim == 3:
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
Z = self.mesh.gridCC[:,2]
Y = self.mesh.gridCC[:,1]
Z = self.mesh.gridCC[:,2]
if self.normal =='X':
f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X
elif self.normal =='Y':
@@ -766,43 +775,43 @@ class PolyMap(IdentityMap):
raise(Exception("Input for normal = X or Y or Z"))
else:
raise(Exception("Only supports 2D"))
return sig1+(sig2-sig1)*(np.arctan(alpha*f)/np.pi+0.5)
def deriv(self, m):
alpha = self.slope
sig1,sig2, c = m[0],m[1],m[2:]
if self.logSigma:
sig1, sig2 = np.exp(sig1), np.exp(sig2)
#2D
if self.mesh.dim == 2:
if self.mesh.dim == 2:
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
if self.normal =='X':
f = polynomial.polyval(Y, c) - X
V = polynomial.polyvander(Y, len(c)-1)
V = polynomial.polyvander(Y, len(c)-1)
elif self.normal =='Y':
f = polynomial.polyval(X, c) - Y
V = polynomial.polyvander(X, len(c)-1)
V = polynomial.polyvander(X, len(c)-1)
else:
raise(Exception("Input for normal = X or Y or Z"))
raise(Exception("Input for normal = X or Y or Z"))
#3D
elif self.mesh.dim == 3:
elif self.mesh.dim == 3:
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
Z = self.mesh.gridCC[:,2]
if self.normal =='X':
f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X
V = polynomial.polyvander2d(Y, Z, self.order)
V = polynomial.polyvander2d(Y, Z, self.order)
elif self.normal =='Y':
f = polynomial.polyval2d(X, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - Y
V = polynomial.polyvander2d(X, Z, self.order)
V = polynomial.polyvander2d(X, Z, self.order)
elif self.normal =='Z':
f = polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,self.order[1]+1))) - Z
V = polynomial.polyvander2d(X, Y, self.order)
V = polynomial.polyvander2d(X, Y, self.order)
else:
raise(Exception("Input for normal = X or Y or Z"))
@@ -815,16 +824,16 @@ class PolyMap(IdentityMap):
g3 = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*V
return sp.csr_matrix(np.c_[g1,g2,g3])
return sp.csr_matrix(np.c_[g1,g2,g3])
class SplineMap(IdentityMap):
"""SplineMap
Parameterize the boundary of two geological units using a spline interpolation
Parameterize the boundary of two geological units using a spline interpolation
..math::
g = f(x)-y
Define the model as:
@@ -849,7 +858,7 @@ class SplineMap(IdentityMap):
def nP(self):
if self.mesh.dim == 2:
return np.size(self.pts)+2
elif self.mesh.dim == 3:
elif self.mesh.dim == 3:
return np.size(self.pts)*2+2
else:
raise(Exception("Only supports 2D and 3D"))
@@ -866,28 +875,28 @@ class SplineMap(IdentityMap):
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
self.spl = UnivariateSpline(self.pts, c, k=self.order, s=0)
if self.normal =='X':
if self.normal =='X':
f = self.spl(Y) - X
elif self.normal =='Y':
f = self.spl(X) - Y
else:
raise(Exception("Input for normal = X or Y or Z"))
# 3D:
# Comments:
# 3D:
# Comments:
# Make two spline functions and link them using linear interpolation.
# This is not quite direct extension of 2D to 3D case
# Using 2D interpolation is possible
elif self.mesh.dim == 3:
elif self.mesh.dim == 3:
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
Y = self.mesh.gridCC[:,1]
Z = self.mesh.gridCC[:,2]
npts = np.size(self.pts)
npts = np.size(self.pts)
if np.mod(c.size, 2):
raise(Exception("Put even points!"))
self.spl = {"splb":UnivariateSpline(self.pts, c[:npts], k=self.order, s=0),
"splt":UnivariateSpline(self.pts, c[npts:], k=self.order, s=0)}
@@ -902,7 +911,7 @@ class SplineMap(IdentityMap):
raise(Exception("Input for normal = X or Y or Z"))
else:
raise(Exception("Only supports 2D and 3D"))
return sig1+(sig2-sig1)*(np.arctan(alpha*f)/np.pi+0.5)
@@ -912,7 +921,7 @@ class SplineMap(IdentityMap):
if self.logSigma:
sig1, sig2 = np.exp(sig1), np.exp(sig2)
#2D
if self.mesh.dim == 2:
if self.mesh.dim == 2:
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
@@ -921,9 +930,9 @@ class SplineMap(IdentityMap):
elif self.normal =='Y':
f = self.spl(X) - Y
else:
raise(Exception("Input for normal = X or Y or Z"))
raise(Exception("Input for normal = X or Y or Z"))
#3D
elif self.mesh.dim == 3:
elif self.mesh.dim == 3:
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
Z = self.mesh.gridCC[:,2]
@@ -931,7 +940,7 @@ class SplineMap(IdentityMap):
zb = self.ptsv[0]
zt = self.ptsv[1]
flines = (self.spl["splt"](Y)-self.spl["splb"](Y))*(Z-zb)/(zt-zb) + self.spl["splb"](Y)
f = flines - X
f = flines - X
# elif self.normal =='Y':
# elif self.normal =='Z':
else:
@@ -944,7 +953,7 @@ class SplineMap(IdentityMap):
g1 = -(np.arctan(alpha*f)/np.pi + 0.5) + 1.0
g2 = (np.arctan(alpha*f)/np.pi + 0.5)
if self.mesh.dim ==2:
g3 = np.zeros((self.mesh.nC, self.npts))
if self.normal =='Y':
@@ -958,7 +967,7 @@ class SplineMap(IdentityMap):
cb = c.copy()
dy = self.mesh.hy[ind]*1.5
ca[i] = ctemp+dy
cb[i] = ctemp-dy
cb[i] = ctemp-dy
spla = UnivariateSpline(self.pts, ca, k=self.order, s=0)
splb = UnivariateSpline(self.pts, cb, k=self.order, s=0)
fderiv = (spla(X)-splb(X))/(2*dy)
@@ -968,7 +977,7 @@ class SplineMap(IdentityMap):
g3 = np.zeros((self.mesh.nC, self.npts*2))
if self.normal =='X':
# Here we use perturbation to compute sensitivity
for i in range(self.npts*2):
for i in range(self.npts*2):
ctemp = c[i]
ind = np.argmin(abs(self.mesh.vectorCCy-ctemp))
ca = c.copy()
@@ -982,20 +991,20 @@ class SplineMap(IdentityMap):
splbb = UnivariateSpline(self.pts, cb[:self.npts], k=self.order, s=0)
flinesa = (self.spl["splt"](Y)-splba(Y))*(Z-zb)/(zt-zb) + splba(Y) - X
flinesb = (self.spl["splt"](Y)-splbb(Y))*(Z-zb)/(zt-zb) + splbb(Y) - X
#treat top boundary
#treat top boundary
else:
splta = UnivariateSpline(self.pts, ca[self.npts:], k=self.order, s=0)
spltb = UnivariateSpline(self.pts, ca[self.npts:], k=self.order, s=0)
flinesa = (self.spl["splt"](Y)-splta(Y))*(Z-zb)/(zt-zb) + splta(Y) - X
flinesb = (self.spl["splt"](Y)-spltb(Y))*(Z-zb)/(zt-zb) + spltb(Y) - X
fderiv = (flinesa-flinesb)/(2*dy)
flinesb = (self.spl["splt"](Y)-spltb(Y))*(Z-zb)/(zt-zb) + spltb(Y) - X
fderiv = (flinesa-flinesb)/(2*dy)
g3[:,i] = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*fderiv
else :
raise(Exception("Not Implemented for Y and Z, your turn :)"))
return sp.csr_matrix(np.c_[g1,g2,g3])
return sp.csr_matrix(np.c_[g1,g2,g3])
SimPEG/Mesh/MeshIO.py
+416
View File
@@ -0,0 +1,416 @@
import numpy as np, os
from SimPEG import Utils
class TensorMeshIO(object):
@classmethod
def readUBC(TensorMesh, fileName):
"""
Read UBC GIF 3DTensor mesh and generate 3D Tensor mesh in simpegTD
Input:
:param fileName, path to the UBC GIF mesh file
Output:
:param SimPEG TensorMesh object
"""
# Interal function to read cell size lines for the UBC mesh files.
def readCellLine(line):
for seg in line.split():
if '*' in seg:
st = seg
sp = seg.split('*')
re = np.array(sp[0],dtype=int)*(' ' + sp[1])
line = line.replace(st,re.strip())
return np.array(line.split(),dtype=float)
# Read the file as line strings, remove lines with comment = !
msh = np.genfromtxt(fileName,delimiter='\n',dtype=np.str,comments='!')
# Fist line is the size of the model
sizeM = np.array(msh[0].split(),dtype=float)
# Second line is the South-West-Top corner coordinates.
x0 = np.array(msh[1].split(),dtype=float)
# Read the cell sizes
h1 = readCellLine(msh[2])
h2 = readCellLine(msh[3])
h3temp = readCellLine(msh[4])
h3 = h3temp[::-1] # Invert the indexing of the vector to start from the bottom.
# Adjust the reference point to the bottom south west corner
x0[2] = x0[2] - np.sum(h3)
# Make the mesh
tensMsh = TensorMesh([h1,h2,h3],x0)
return tensMsh
@classmethod
def readVTK(TensorMesh, fileName):
"""
Read VTK Rectilinear (vtr xml file) and return SimPEG Tensor mesh and model
Input:
:param vtrFileName, path to the vtr model file to write to
Output:
:return SimPEG TensorMesh object
:return SimPEG model dictionary
"""
# Import
from vtk import vtkXMLRectilinearGridReader as vtrFileReader
from vtk.util.numpy_support import vtk_to_numpy
# Read the file
vtrReader = vtrFileReader()
vtrReader.SetFileName(fileName)
vtrReader.Update()
vtrGrid = vtrReader.GetOutput()
# Sort information
hx = np.abs(np.diff(vtk_to_numpy(vtrGrid.GetXCoordinates())))
xR = vtk_to_numpy(vtrGrid.GetXCoordinates())[0]
hy = np.abs(np.diff(vtk_to_numpy(vtrGrid.GetYCoordinates())))
yR = vtk_to_numpy(vtrGrid.GetYCoordinates())[0]
zD = np.diff(vtk_to_numpy(vtrGrid.GetZCoordinates()))
# Check the direction of hz
if np.all(zD < 0):
hz = np.abs(zD[::-1])
zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[-1]
else:
hz = np.abs(zD)
zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[0]
x0 = np.array([xR,yR,zR])
# Make the SimPEG object
tensMsh = TensorMesh([hx,hy,hz],x0)
# Grap the models
models = {}
for i in np.arange(vtrGrid.GetCellData().GetNumberOfArrays()):
modelName = vtrGrid.GetCellData().GetArrayName(i)
if np.all(zD < 0):
modFlip = vtk_to_numpy(vtrGrid.GetCellData().GetArray(i))
tM = tensMsh.r(modFlip,'CC','CC','M')
modArr = tensMsh.r(tM[:,:,::-1],'CC','CC','V')
else:
modArr = vtk_to_numpy(vtrGrid.GetCellData().GetArray(i))
models[modelName] = modArr
# Return the data
return tensMsh, models
def writeVTK(mesh, fileName, models=None):
"""
Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model.
Input:
:param str, path to the output vtk file
:param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK
:param models, dictionary of numpy.array - Name('s) and array('s). Match number of cells
"""
# Import
from vtk import vtkRectilinearGrid as rectGrid, vtkXMLRectilinearGridWriter as rectWriter, VTK_VERSION
from vtk.util.numpy_support import numpy_to_vtk
# Deal with dimensionalities
if mesh.dim >= 1:
vX = mesh.vectorNx
xD = mesh.nNx
yD,zD = 1,1
vY, vZ = np.array([0,0])
if mesh.dim >= 2:
vY = mesh.vectorNy
yD = mesh.nNy
if mesh.dim == 3:
vZ = mesh.vectorNz
zD = mesh.nNz
# Use rectilinear VTK grid.
# Assign the spatial information.
vtkObj = rectGrid()
vtkObj.SetDimensions(xD,yD,zD)
vtkObj.SetXCoordinates(numpy_to_vtk(vX,deep=1))
vtkObj.SetYCoordinates(numpy_to_vtk(vY,deep=1))
vtkObj.SetZCoordinates(numpy_to_vtk(vZ,deep=1))
# Assign the model('s) to the object
if models is not None:
for item in models.iteritems():
# Convert numpy array
vtkDoubleArr = numpy_to_vtk(item[1],deep=1)
vtkDoubleArr.SetName(item[0])
vtkObj.GetCellData().AddArray(vtkDoubleArr)
# Set the active scalar
vtkObj.GetCellData().SetActiveScalars(models.keys()[0])
# vtkObj.Update()
# Check the extension of the fileName
ext = os.path.splitext(fileName)[1]
if ext is '':
fileName = fileName + '.vtr'
elif ext not in '.vtr':
raise IOError('{:s} is an incorrect extension, has to be .vtr')
# Write the file.
vtrWriteFilter = rectWriter()
if float(VTK_VERSION.split('.')[0]) >=6:
vtrWriteFilter.SetInputData(vtkObj)
else:
vtuWriteFilter.SetInput(vtuObj)
vtrWriteFilter.SetFileName(fileName)
vtrWriteFilter.Update()
def readModelUBC(mesh, fileName):
"""
Read UBC 3DTensor mesh model and generate 3D Tensor mesh model in simpeg
Input:
:param fileName, path to the UBC GIF mesh file to read
:param mesh, TensorMesh object, mesh that coresponds to the model
Output:
:return numpy array, model with TensorMesh ordered
"""
f = open(fileName, 'r')
model = np.array(map(float, f.readlines()))
f.close()
model = np.reshape(model, (mesh.nCz, mesh.nCx, mesh.nCy), order = 'F')
model = model[::-1,:,:]
model = np.transpose(model, (1, 2, 0))
model = Utils.mkvc(model)
return model
def writeModelUBC(mesh, fileName, model):
"""
Writes a model associated with a SimPEG TensorMesh
to a UBC-GIF format model file.
:param str fileName: File to write to
:param simpeg.Mesh.TensorMesh mesh: The mesh
:param numpy.ndarray model: The model
"""
# Reshape model to a matrix
modelMat = mesh.r(model,'CC','CC','M')
# Transpose the axes
modelMatT = modelMat.transpose((2,0,1))
# Flip z to positive down
modelMatTR = Utils.mkvc(modelMatT[::-1,:,:])
np.savetxt(fileName, modelMatTR.ravel())
def writeUBC(mesh, fileName, models=None):
"""
Writes a SimPEG TensorMesh to a UBC-GIF format mesh file.
:param str fileName: File to write to
:param simpeg.Mesh.TensorMesh mesh: The mesh
"""
assert mesh.dim == 3
s = ''
s += '%i %i %i\n' %tuple(mesh.vnC)
origin = mesh.x0 + np.array([0,0,mesh.hz.sum()]) # Have to it in the same operation or use mesh.x0.copy(), otherwise the mesh.x0 is updated.
origin.dtype = float
s += '%.2f %.2f %.2f\n' %tuple(origin)
s += ('%.2f '*mesh.nCx+'\n')%tuple(mesh.hx)
s += ('%.2f '*mesh.nCy+'\n')%tuple(mesh.hy)
s += ('%.2f '*mesh.nCz+'\n')%tuple(mesh.hz[::-1])
f = open(fileName, 'w')
f.write(s)
f.close()
if models is None: return
assert type(models) is dict, 'models must be a dict'
for key in models:
assert type(key) is str, 'The dict key is a file name'
mesh.writeModelUBC(key, models[key])
class TreeMeshIO(object):
def writeUBC(mesh, fileName, models=None):
"""
Write UBC ocTree mesh and model files from a simpeg ocTree mesh and model.
:param str fileName: File to write to
:param simpeg.Mesh.TreeMesh mesh: The mesh
:param dictionary models: The models in a dictionary, where the keys is the name of the of the model file
"""
# Calculate information to write in the file.
# Number of cells in the underlying mesh
nCunderMesh = np.array([h.size for h in mesh.h],dtype=np.int64)
# The top-south-west most corner of the mesh
tswCorn = mesh.x0 + np.array([0,0,np.sum(mesh.h[2])])
# Smallest cell size
smallCell = np.array([h.min() for h in mesh.h])
# Number of cells
nrCells = mesh.nC
## Extract iformation about the cells.
# cell pointers
cellPointers = np.array([c._pointer for c in mesh])
# cell with
cellW = np.array([ mesh._levelWidth(i) for i in cellPointers[:,-1] ])
# Need to shift the pointers to work with UBC indexing
# UBC Octree indexes always the top-left-close (top-south-west) corner first and orders the cells in z(top-down),x,y vs x,y,z(bottom-up).
# Shift index up by 1
ubcCellPt = cellPointers[:,0:-1].copy() + np.array([1.,1.,1.])
# Need reindex the z index to be from the top-left-close corner and to be from the global top.
ubcCellPt[:,2] = ( nCunderMesh[-1] + 2) - (ubcCellPt[:,2] + cellW)
# Reorder the ubcCellPt
ubcReorder = np.argsort(ubcCellPt.view(','.join(3*['float'])),axis=0,order=['f2','f1','f0'])[:,0]
# Make a array with the pointers and the withs, that are order in the ubc ordering
indArr = np.concatenate((ubcCellPt[ubcReorder,:],cellW[ubcReorder].reshape((-1,1)) ),axis=1)
## Write the UBC octree mesh file
with open(fileName,'w') as mshOut:
mshOut.write('{:.0f} {:.0f} {:.0f}\n'.format(nCunderMesh[0],nCunderMesh[1],nCunderMesh[2]))
mshOut.write('{:.4f} {:.4f} {:.4f}\n'.format(tswCorn[0],tswCorn[1],tswCorn[2]))
mshOut.write('{:.3f} {:.3f} {:.3f}\n'.format(smallCell[0],smallCell[1],smallCell[2]))
mshOut.write('{:.0f} \n'.format(nrCells))
np.savetxt(mshOut,indArr,fmt='%i')
## Print the models
# Assign the model('s) to the object
if models is not None:
# indUBCvector = np.argsort(cX0[np.argsort(np.concatenate((cX0[:,0:2],cX0[:,2:3].max() - cX0[:,2:3]),axis=1).view(','.join(3*['float'])),axis=0,order=('f2','f1','f0'))[:,0]].view(','.join(3*['float'])),axis=0,order=('f2','f1','f0'))[:,0]
for item in models.iteritems():
# Save the data
np.savetxt(item[0],item[1][ubcReorder],fmt='%3.5e')
@classmethod
def readUBC(TreeMesh, meshFile):
"""
Read UBC 3D OcTree mesh and/or modelFiles
Input:
:param str meshFile: path to the UBC GIF OcTree mesh file to read
Output:
:return SimPEG.Mesh.TreeMesh mesh: The octree mesh
:return list of ndarray's: models as a list of numpy array's
"""
## Read the file lines
fileLines = np.genfromtxt(meshFile,dtype=str,delimiter='\n')
# Extract the data
nCunderMesh = np.array(fileLines[0].split(),dtype=float)
# I think this is the case?
if np.unique(nCunderMesh).size >1:
raise Exception('SimPEG TreeMeshes have the same number of cell in all directions')
tswCorn = np.array(fileLines[1].split(),dtype=float)
smallCell = np.array(fileLines[2].split(),dtype=float)
nrCells = np.array(fileLines[3].split(),dtype=float)
# Read the index array
indArr = np.genfromtxt(fileLines[4::],dtype=np.int)
## Calculate simpeg parameters
h1,h2,h3 = [np.ones(nr)*sz for nr,sz in zip(nCunderMesh,smallCell)]
x0 = tswCorn - np.array([0,0,np.sum(h3)])
# Need to convert the index array to a points list that complies with SimPEG TreeMesh.
# Shift to start at 0
simpegCellPt = indArr[:,0:-1].copy()
simpegCellPt[:,2] = ( nCunderMesh[-1] + 2) - (simpegCellPt[:,2] + indArr[:,3])
# Need reindex the z index to be from the bottom-left-close corner and to be from the global bottom.
simpegCellPt = simpegCellPt - np.array([1.,1.,1.])
# Calculate the cell level
simpegLevel = np.log2(np.min(nCunderMesh)) - np.log2(indArr[:,3])
# Make a pointer matrix
simpegPointers = np.concatenate((simpegCellPt,simpegLevel.reshape((-1,1))),axis=1)
## Make the tree mesh
mesh = TreeMesh([h1,h2,h3],x0)
mesh._cells = set([mesh._index(p) for p in simpegPointers.tolist()])
# Figure out the reordering
mesh._simpegReorderUBC = np.argsort(np.array([mesh._index(i) for i in simpegPointers.tolist()]))
# mesh._simpegReorderUBC = np.argsort((np.array([[1,1,1,-1]])*simpegPointers).view(','.join(4*['float'])),axis=0,order=['f3','f2','f1','f0'])[:,0]
return mesh
def readModelUBC(mesh, fileName):
"""
Read UBC OcTree model and get vector
Input:
:param fileName, path to the UBC GIF model file to read
Output:
:return numpy array, OcTree model
"""
if type(fileName) is list:
out = {}
for f in fileName:
out[f] = mesh.readModelUBC(f)
return out
assert hasattr(mesh, '_simpegReorderUBC'), 'The file must have been loaded from a UBC format.'
assert mesh.dim == 3
modList = []
modArr = np.loadtxt(fileName)
if len(modArr.shape) == 1:
modList.append(modArr[mesh._simpegReorderUBC])
else:
modList.append(modArr[mesh._simpegReorderUBC,:])
return modList
def writeVTK(mesh, fileName, models=None):
"""
Function to write a VTU file from a SimPEG TreeMesh and model.
"""
import vtk
from vtk import vtkXMLUnstructuredGridWriter as Writer, VTK_VERSION
from vtk.util.numpy_support import numpy_to_vtk, numpy_to_vtkIdTypeArray
if str(type(mesh)).split()[-1][1:-2] not in 'SimPEG.Mesh.TreeMesh.TreeMesh':
raise IOError('mesh is not a SimPEG TreeMesh.')
# Make the data parts for the vtu object
# Points
mesh.number()
ptsMat = mesh._gridN + mesh.x0
vtkPts = vtk.vtkPoints()
vtkPts.SetData(numpy_to_vtk(ptsMat,deep=True))
# Cells
cellConn = np.array([c.nodes for c in mesh],dtype=np.int64)
cellsMat = np.concatenate((np.ones((cellConn.shape[0],1),dtype=np.int64)*cellConn.shape[1],cellConn),axis=1).ravel()
cellsArr = vtk.vtkCellArray()
cellsArr.SetNumberOfCells(cellConn.shape[0])
cellsArr.SetCells(cellConn.shape[0],numpy_to_vtkIdTypeArray(cellsMat,deep=True))
# Make the object
vtuObj = vtk.vtkUnstructuredGrid()
vtuObj.SetPoints(vtkPts)
vtuObj.SetCells(vtk.VTK_VOXEL,cellsArr)
# Add the level of refinement as a cell array
cellSides = np.array([np.array(vtuObj.GetCell(i).GetBounds()).reshape((3,2)).dot(np.array([-1, 1])) for i in np.arange(vtuObj.GetNumberOfCells())])
uniqueLevel, indLevel = np.unique(np.prod(cellSides,axis=1),return_inverse=True)
refineLevelArr = numpy_to_vtk(indLevel.max() - indLevel,deep=1)
refineLevelArr.SetName('octreeLevel')
vtuObj.GetCellData().AddArray(refineLevelArr)
# Assign the model('s) to the object
if models is not None:
for item in models.iteritems():
# Convert numpy array
vtkDoubleArr = numpy_to_vtk(item[1],deep=1)
vtkDoubleArr.SetName(item[0])
vtuObj.GetCellData().AddArray(vtkDoubleArr)
# Make the writer
vtuWriteFilter = Writer()
if float(VTK_VERSION.split('.')[0]) >=6:
vtuWriteFilter.SetInputData(vtuObj)
else:
vtuWriteFilter.SetInput(vtuObj)
vtuWriteFilter.SetFileName(fileName)
# Write the file
vtuWriteFilter.Update()
SimPEG/Mesh/TensorMesh.py
+559 -558
View File
File diff suppressed because it is too large Load Diff
SimPEG/Mesh/TreeMesh.py
+121 -141
View File
@@ -90,24 +90,28 @@
#
from SimPEG import np, sp, Utils, Solver
import matplotlib.pyplot as plt
import matplotlib
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.colors as colors
import matplotlib.cm as cmx
import TreeUtils
try:
import TreeUtils
_IMPORT_TREEUTILS = True
except Exception, e:
_IMPORT_TREEUTILS = False
from InnerProducts import InnerProducts
from TensorMesh import TensorMesh, BaseTensorMesh
from MeshIO import TreeMeshIO
import time
MAX_BITS = 20
class TreeMesh(BaseTensorMesh, InnerProducts):
class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
_meshType = 'TREE'
def __init__(self, h, x0=None, levels=None):
if not _IMPORT_TREEUTILS:
raise Exception('Could not import the Cython code to run the TreeMesh Try:.\n\npython setup.py build_ext --inplace')
assert type(h) is list, 'h must be a list'
assert len(h) in [2,3], "There is only support for TreeMesh in 2D or 3D."
@@ -561,15 +565,18 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
return [p - (p % mod) for p in pointer[:-1]] + [pointer[-1]-1]
def _cellN(self, p):
"""Node location [x,y(,z)] of a single cell, closest to origin, given a pointer."""
p = self._asPointer(p)
return [hi[:p[ii]].sum() for ii, hi in enumerate(self.h)]
def _cellH(self, p):
"""Widths of a single cell given a pointer."""
p = self._asPointer(p)
w = self._levelWidth(p[-1])
return [hi[p[ii]:p[ii]+w].sum() for ii, hi in enumerate(self.h)]
def _cellC(self, p):
"""Cell center of a single cell (without origin correction), given a pointer."""
return (np.array(self._cellH(p))/2.0 + self._cellN(p)).tolist()
def _levelWidth(self, level):
@@ -824,8 +831,10 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
def _numberCells(self, force=False):
if not self.__dirtyCells__ and not force: return
self._cc2i = dict()
self._i2cc = dict()
for ii, c in enumerate(sorted(self._cells)):
self._cc2i[c] = ii
self._i2cc[ii] = c
self.__dirtyCells__ = False
def _numberNodes(self, force=False):
@@ -1701,9 +1710,9 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
"Construct the averaging operator on cell faces to cell centers."
if getattr(self, '_aveF2CC', None) is None:
if self.dim == 2:
self._aveF2CC = 1./self.dim*sp.hstack([self.aveFx2CC, self.aveFy2CC])
self._aveF2CC = 1./self.dim*sp.hstack([self.aveFx2CC, self.aveFy2CC]).tocsr()
elif self.dim == 3:
self._aveF2CC = 1./self.dim*sp.hstack([self.aveFx2CC, self.aveFy2CC, self.aveFz2CC])
self._aveF2CC = 1./self.dim*sp.hstack([self.aveFx2CC, self.aveFy2CC, self.aveFz2CC]).tocsr()
return self._aveF2CC
@property
@@ -1711,9 +1720,9 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
"Construct the averaging operator on cell faces to cell centers."
if getattr(self, '_aveF2CCV', None) is None:
if self.dim == 2:
self._aveF2CCV = sp.block_diag([self.aveFx2CC, self.aveFy2CC])
self._aveF2CCV = sp.block_diag([self.aveFx2CC, self.aveFy2CC]).tocsr()
elif self.dim == 3:
self._aveF2CCV = sp.block_diag([self.aveFx2CC, self.aveFy2CC, self.aveFz2CC])
self._aveF2CCV = sp.block_diag([self.aveFx2CC, self.aveFy2CC, self.aveFz2CC]).tocsr()
return self._aveF2CCV
@property
@@ -1960,11 +1969,18 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
def plotGrid(self, ax=None, showIt=False,
grid=True,
cells=True, cellLine=False,
cells=False, cellLine=False,
nodes=False,
facesX=False, facesY=False, facesZ=False,
edgesX=False, edgesY=False, edgesZ=False):
import matplotlib.pyplot as plt
import matplotlib
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.colors as colors
import matplotlib.cm as cmx
# self.number()
axOpts = {'projection':'3d'} if self.dim == 3 else {}
@@ -1975,24 +1991,28 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
fig = ax.figure
if grid:
X, Y, Z = [], [], []
for ind in self._sortedCells:
p = self._asPointer(ind)
n = self._cellN(p)
h = self._cellH(p)
x = [n[0] , n[0] + h[0], n[0] + h[0], n[0] , n[0]]
y = [n[1] , n[1] , n[1] + h[1], n[1] + h[1], n[1]]
if self.dim == 2:
ax.plot(x,y, 'b-')
X += [n[0] , n[0] + h[0], n[0] + h[0], n[0] , n[0], np.nan]
Y += [n[1] , n[1] , n[1] + h[1], n[1] + h[1], n[1], np.nan]
elif self.dim == 3:
ax.plot(x,y, 'b-', zs=[n[2]]*5)
z = [n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2]]
ax.plot(x,y, 'b-', zs=z)
X += [n[0] , n[0] + h[0], n[0] + h[0], n[0] , n[0], np.nan]*2
Y += [n[1] , n[1] , n[1] + h[1], n[1] + h[1], n[1], np.nan]*2
Z += [n[2]]*5+[np.nan]
Z += [n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2], np.nan]
sides = [0,0], [h[0],0], [0,h[1]], [h[0],h[1]]
for s in sides:
x = [n[0] + s[0], n[0] + s[0]]
y = [n[1] + s[1], n[1] + s[1]]
z = [n[2] , n[2] + h[2]]
ax.plot(x,y, 'b-', zs=z)
X += [n[0] + s[0], n[0] + s[0]]
Y += [n[1] + s[1], n[1] + s[1]]
Z += [n[2] , n[2] + h[2]]
if self.dim == 2:
ax.plot(X,Y, 'b-')
elif self.dim == 3:
ax.plot(X,Y, 'b-', zs=Z)
if self.dim == 2:
if cells:
@@ -2004,11 +2024,13 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
ax.plot(self._gridN[:,0], self._gridN[:,1], 'ms')
ax.plot(self._gridN[self._hangingN.keys(),0], self._gridN[self._hangingN.keys(),1], 'ms', ms=10, mfc='none', mec='m')
if facesX:
ax.plot(self._gridFx[self._hangingFx.keys(),0], self._gridFx[self._hangingFx.keys(),1], 'gs', ms=10, mfc='none', mec='g')
ax.plot(self._gridFx[:,0], self._gridFx[:,1], 'g>')
ax.plot(self._gridFx[self._hangingFx.keys(),0], self._gridFx[self._hangingFx.keys(),1], 'gs', ms=10, mfc='none', mec='g')
if facesY:
ax.plot(self._gridFy[self._hangingFy.keys(),0], self._gridFy[self._hangingFy.keys(),1], 'gs', ms=10, mfc='none', mec='g')
ax.plot(self._gridFy[:,0], self._gridFy[:,1], 'g^')
ax.plot(self._gridFy[self._hangingFy.keys(),0], self._gridFy[self._hangingFy.keys(),1], 'gs', ms=10, mfc='none', mec='g')
ax.set_xlabel('x1')
ax.set_ylabel('x2')
elif self.dim == 3:
if cells:
ax.plot(self.gridCC[:,0], self.gridCC[:,1], 'r.', zs=self.gridCC[:,2])
@@ -2056,7 +2078,6 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
ind = [key, hf[0]]
ax.plot(self._gridEx[ind,0], self._gridEx[ind,1], 'k:', zs=self._gridEx[ind,2])
if edgesY:
ax.plot(self._gridEy[:,0], self._gridEy[:,1], 'k<', zs=self._gridEy[:,2])
ax.plot(self._gridEy[self._hangingEy.keys(),0], self._gridEy[self._hangingEy.keys(),1], 'ks', ms=10, mfc='none', mec='k', zs=self._gridEy[self._hangingEy.keys(),2])
@@ -2072,15 +2093,28 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
for hf in self._hangingEz[key]:
ind = [key, hf[0]]
ax.plot(self._gridEz[ind,0], self._gridEz[ind,1], 'k:', zs=self._gridEz[ind,2])
ax.set_xlabel('x1')
ax.set_ylabel('x2')
ax.set_zlabel('x3')
ax.grid(True)
if showIt:plt.show()
def plotImage(self, I, ax=None, showIt=True, grid=False):
def plotImage(self, I, ax=None, showIt=False, grid=False, clim=None):
if self.dim == 3: raise Exception('Use plot slice?')
import matplotlib.pyplot as plt
import matplotlib
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.colors as colors
import matplotlib.cm as cmx
if ax is None: ax = plt.subplot(111)
jet = cm = plt.get_cmap('jet')
cNorm = colors.Normalize(vmin=I.min(), vmax=I.max())
cNorm = colors.Normalize(
vmin=I.min() if clim is None else clim[0],
vmax=I.max() if clim is None else clim[1])
scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet)
ax.set_xlim((self.x0[0], self.h[0].sum()))
ax.set_ylim((self.x0[1], self.h[1].sum()))
@@ -2089,8 +2123,10 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
ax.add_patch(plt.Rectangle((x0[0], x0[1]), sz[0], sz[1], facecolor=scalarMap.to_rgba(I[ii]), edgecolor='k' if grid else 'none'))
# if text: ax.text(self.center[0],self.center[1],self.num)
scalarMap._A = [] # http://stackoverflow.com/questions/8342549/matplotlib-add-colorbar-to-a-sequence-of-line-plots
plt.colorbar(scalarMap)
ax.set_xlabel('x')
ax.set_ylabel('y')
if showIt: plt.show()
return [scalarMap]
def plotSlice(self, v, vType='CC',
normal='Z', ind=None, grid=True, view='real',
@@ -2102,6 +2138,13 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
assert vType in ['CC','F','E']
assert self.dim == 3
import matplotlib.pyplot as plt
import matplotlib
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.colors as colors
import matplotlib.cm as cmx
szSliceDim = len(getattr(self, 'h'+normal.lower())) #: Size of the sliced dimension
if ind is None: ind = int(szSliceDim/2)
assert type(ind) in [int, long], 'ind must be an integer'
@@ -2181,6 +2224,25 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
if showIt: plt.show()
return tuple(out)
def __len__(self): return self.nC
def __getitem__(self, key):
if isinstance( key, slice ) :
#Get the start, stop, and step from the slice
return [self[ii] for ii in xrange(*key.indices(len(self)))]
elif isinstance( key, int ) :
if key < 0 : #Handle negative indices
key += len( self )
if key >= len( self ) :
raise IndexError, "The index (%d) is out of range."%key
self._numberCells() # no-op if numbered
index = self._i2cc[key]
pointer = self._asPointer(index)
return Cell(self, index, pointer)
else:
raise TypeError, "Invalid argument type."
class Cell(object):
def __init__(self, mesh, index, pointer):
@@ -2188,10 +2250,39 @@ class Cell(object):
self._index = index
self._pointer = pointer
@property
def nodes(self):
"""The node index in _gridN (this may include hanging nodes)."""
M = self.mesh
M._numberNodes()
p = self._pointer
i = self._index
w = M._levelWidth(p[-1])
if M.dim == 2:
n = [
i,
M._index([ p[0] + w, p[1] , p[2]]),
M._index([ p[0] , p[1]+ w, p[2]]),
M._index([ p[0] + w, p[1]+ w, p[2]]),
]
elif self.dim == 3:
n = [
i,
M._index([ p[0] + w, p[1] , p[2] ,p[3]]),
M._index([ p[0] , p[1] + w, p[2] ,p[3]]),
M._index([ p[0] + w, p[1] + w, p[2] ,p[3]]),
M._index([ p[0] , p[1] , p[2] + w,p[3]]),
M._index([ p[0] + w, p[1] , p[2] + w,p[3]]),
M._index([ p[0] , p[1] + w, p[2] + w,p[3]]),
M._index([ p[0] + w, p[1] + w, p[2] + w,p[3]]),
]
return [M._n2i[_] for _ in n]
@property
def center(self):
if getattr(self, '_center', None) is None:
self._center = self.mesh._cellC(self._pointer)
self._center = np.array(self.mesh._cellC(self._pointer))
return self._center
@property
def h(self): return self.mesh._cellH(self._pointer)
@@ -2245,114 +2336,3 @@ class NotBalancedException(TreeException):
pass
class CellLookUpException(TreeException):
pass
if __name__ == '__main__':
def topo(x):
return np.sin(x*(2.*np.pi))*0.3 + 0.5
def function(cell):
r = cell.center - np.array([0.5]*len(cell.center))
dist = np.sqrt(r.dot(r))
# dist2 = np.abs(cell.center[-1] - topo(cell.center[0]))
# dist = min([dist1,dist2])
# if dist < 0.05:
# return 5
if dist < 0.1:
return 5
if dist < 0.2:
return 4
if dist < 0.4:
return 3
return 2
# T = TreeMesh([[(1,128)],[(1,128)],[(1,128)]],levels=7)
# T = TreeMesh([128,128,128])
# T = TreeMesh([64,64],levels=6)
T = TreeMesh([4,4,4])
# T = TreeMesh([[(1,128)],[(1,128)]],levels=7)
# T.refine(lambda xc:2, balance=False)
# T._index([0,0,0])
# T._pointer(0)
# tic = time.time()
T.refine(function)#, balance=False)
# print time.time() - tic
# print T.nC
T.plotSlice(np.log(T.vol))#np.random.rand(T.nC))
plt.show()
blah
# T.plotImage(np.arange(len(T.vol)),showIt=True)
# print T.getFaceInnerProduct()
# print T.gridFz
# T._refineCell([8,0,1])
# T._refineCell([8,0,2])
# T._refineCell([12,0,2])
# T._refineCell([8,4,2])
# T._refineCell([6,0,3])
# T._refineCell([8,8,1])
# T._refineCell([0,0,0,1])
# T.__dirty__ = True
# print T.gridFx.shape[0], T.nFx
ax = plt.subplot(211)
ax.spy(T.edgeCurl)
# print Mesh.TensorMesh([2,2,2]).edgeCurl.todense()
# print T.edgeCurl.todense()
# print Mesh.TensorMesh([2,2,2]).edgeCurl.todense() - T.edgeCurl.todense()
# print T.gridEy - Mesh.TensorMesh([2,2,2]).gridEy
# print T.edge
# T.plotGrid(ax=ax)
# R = deflationMatrix(T._facesX, T._hangingFx, T._fx2i)
# print R
ax = plt.subplot(212)#, projection='3d')
ax.spy(Mesh.TensorMesh([2,2,2]).edgeCurl)
# ax = plt.subplot(313)
# ax.spy(T.faceDiv[:,:T.nFx] * R)
# T.balance()
# T.plotGrid(ax=ax)
# cx = T._getNextCell([0,0,1],direction=0,positive=True)
# print cx
# # print [T._asPointer(_) for _ in cx]
# cx = T._getNextCell([8,0,3],direction=0,positive=False)
# print T._asPointer(cx)
# cx = T._getNextCell([8,8,1],direction=1,positive=False)
# print cx, #[T._asPointer(_) for _ in cx]
# cm = T._getNextCell([64,80,4],direction=0,positive=False)
# cy = T._getNextCell([64,80,4],direction=1,positive=True)
# cp = T._getNextCell([64,80,4],direction=1,positive=False)
# ax.plot( T._cellN([4,0,1])[0],T._cellN([4,0,1])[1], 'yd')
# ax.plot( T._cellN(cx)[0],T._cellN(cx)[1], 'ys')
# ax.plot( T._cellN(cm)[0],T._cellN(cm)[1], 'ys')
# ax.plot( T._cellN(cy)[0],T._cellN(cy)[1], 'ys')
# ax.plot( T._cellN(cp[0])[0],T._cellN(cp[0])[1], 'ys')
# ax.plot( T._cellN(cp[1])[0],T._cellN(cp[1])[1], 'ys')
# print T.nN
plt.show()
SimPEG/Mesh/View.py
+7 -61
View File
@@ -1,8 +1,11 @@
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from mpl_toolkits.mplot3d import Axes3D
from SimPEG.Utils import mkvc, animate
from SimPEG.Utils import mkvc
try:
import matplotlib.pyplot as plt
import matplotlib
from mpl_toolkits.mplot3d import Axes3D
except ImportError, e:
print 'Trouble importing matplotlib.'
class TensorView(object):
@@ -479,63 +482,6 @@ class TensorView(object):
ax.grid(True)
if showIt: plt.show()
def slicer(mesh, var, imageType='CC', normal='z', index=0, ax=None, clim=None):
assert normal in 'xyz', 'normal must be x, y, or z'
if ax is None: ax = plt.subplot(111)
I = mesh.r(var,'CC','CC','M')
axes = [p for p in 'xyz' if p not in normal.lower()]
if normal is 'x': I = I[index,:,:]
if normal is 'y': I = I[:,index,:]
if normal is 'z': I = I[:,:,index]
if clim is None: clim = [I.min(),I.max()]
p = ax.pcolormesh(getattr(mesh,'vectorN'+axes[0]),getattr(mesh,'vectorN'+axes[1]),I.T,vmin=clim[0],vmax=clim[1])
ax.axis('tight')
ax.set_xlabel(axes[0])
ax.set_ylabel(axes[1])
return p
def videoSlicer(mesh,var,imageType='CC',normal='z',figsize=(10,8)):
assert mesh.dim > 2, 'This is for 3D meshes only.'
# First set up the figure, the axis, and the plot element we want to animate
fig = plt.figure(figsize=figsize)
ax = plt.axes()
clim = [var.min(),var.max()]
plt.colorbar(mesh.slicer(var, imageType=imageType, normal=normal, index=0, ax=ax, clim=clim))
tlt = plt.title(normal)
def animateFrame(i):
mesh.slicer(var, imageType=imageType, normal=normal, index=i, ax=ax, clim=clim)
tlt.set_text(normal.upper()+('-Slice: %d, %4.4f' % (i,getattr(mesh,'vectorCC'+normal)[i])))
return animate(fig, animateFrame, frames=mesh.vnC['xyz'.index(normal)])
def video(mesh, var, function, figsize=(10, 8), colorbar=True, skip=1):
"""
Call a function for a list of models to create a video.
::
def function(var, ax, clim, tlt, i):
tlt.set_text('%d'%i)
return mesh.plotImage(var, imageType='CC', ax=ax, clim=clim)
mesh.video([model1, model2, ..., modeln],function)
"""
# First set up the figure, the axis, and the plot element we want to animate
fig = plt.figure(figsize=figsize)
ax = plt.axes()
VAR = np.concatenate(var)
clim = [VAR.min(),VAR.max()]
tlt = plt.title('')
if colorbar:
plt.colorbar(function(var[0],ax,clim,tlt,0))
frames = np.arange(0,len(var),skip)
def animateFrame(j):
i = frames[j]
function(var[i],ax,clim,tlt,i)
return animate(fig, animateFrame, frames=len(frames))
class CylView(object):
SimPEG/Problem.py
+2 -2
View File
@@ -32,8 +32,8 @@ class BaseProblem(object):
val._assertMatchesPair(self.mapPair)
self._mapping = val
else:
self._mapping = self.PropMap(val)
self._mapping = self.PropMap(val)
def __init__(self, mesh, mapping=None, **kwargs):
Utils.setKwargs(self, **kwargs)
assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
SimPEG/Regularization.py
+47 -87
View File
@@ -20,12 +20,13 @@ class BaseRegularization(object):
mesh = None #: A SimPEG.Mesh instance.
mref = None #: Reference model.
def __init__(self, mesh, mapping=None, **kwargs):
def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
Utils.setKwargs(self, **kwargs)
self.mesh = mesh
assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
self.mapping = mapping or Maps.IdentityMap(mesh)
self.mapping = mapping or self.mapPair(mesh)
self.mapping._assertMatchesPair(self.mapPair)
self.indActive = indActive
@property
def parent(self):
@@ -112,89 +113,8 @@ class BaseRegularization(object):
return mD.T * ( self.W.T * ( self.W * ( mD * v) ) )
class Tikhonov(BaseRegularization):
"""**Tikhonov Regularization**
Here we will define regularization of a model, m, in general however, this should be thought of as (m-m_ref) but otherwise it is exactly the same:
.. math::
R(m) = \int_\Omega \\frac{\\alpha_x}{2}\left(\\frac{\partial m}{\partial x}\\right)^2 + \\frac{\\alpha_y}{2}\left(\\frac{\partial m}{\partial y}\\right)^2 \partial v
Our discrete gradient operator works on cell centers and gives the derivative on the cell faces, which is not where we want to be evaluating this integral. We need to average the values back to the cell-centers before we integrate. To avoid null spaces, we square first and then average. In 2D with ij notation it looks like this:
.. math::
R(m) \\approx \sum_{ij} \left[\\frac{\\alpha_x}{2}\left[\left(\\frac{m_{i+1,j} - m_{i,j}}{h}\\right)^2 + \left(\\frac{m_{i,j} - m_{i-1,j}}{h}\\right)^2\\right]
+ \\frac{\\alpha_y}{2}\left[\left(\\frac{m_{i,j+1} - m_{i,j}}{h}\\right)^2 + \left(\\frac{m_{i,j} - m_{i,j-1}}{h}\\right)^2\\right]
\\right]h^2
If we let D_1 be the derivative matrix in the x direction
.. math::
\mathbf{D}_1 = \mathbf{I}_2\otimes\mathbf{d}_1
.. math::
\mathbf{D}_2 = \mathbf{d}_2\otimes\mathbf{I}_1
Where d_1 is the one dimensional derivative:
.. math::
\mathbf{d}_1 = \\frac{1}{h} \left[ \\begin{array}{cccc}
-1 & 1 & & \\\\
& \ddots & \ddots&\\\\
& & -1 & 1\end{array} \\right]
.. math::
R(m) \\approx \mathbf{v}^\\top \left[\\frac{\\alpha_x}{2}\mathbf{A}_1 (\mathbf{D}_1 m) \odot (\mathbf{D}_1 m) + \\frac{\\alpha_y}{2}\mathbf{A}_2 (\mathbf{D}_2 m) \odot (\mathbf{D}_2 m) \\right]
Recall that this is really a just point wise multiplication, or a diagonal matrix times a vector. When we multiply by something in a diagonal we can interchange and it gives the same results (i.e. it is point wise)
.. math::
\mathbf{a\odot b} = \\text{diag}(\mathbf{a})\mathbf{b} = \\text{diag}(\mathbf{b})\mathbf{a} = \mathbf{b\odot a}
and the transpose also is true (but the sizes have to make sense...):
.. math::
\mathbf{a}^\\top\\text{diag}(\mathbf{b}) = \mathbf{b}^\\top\\text{diag}(\mathbf{a})
So R(m) can simplify to:
.. math::
R(m) \\approx \mathbf{m}^\\top \left[\\frac{\\alpha_x}{2}\mathbf{D}_1^\\top \\text{diag}(\mathbf{A}_1^\\top\mathbf{v}) \mathbf{D}_1 + \\frac{\\alpha_y}{2}\mathbf{D}_2^\\top \\text{diag}(\mathbf{A}_2^\\top \mathbf{v}) \mathbf{D}_2 \\right] \mathbf{m}
We will define W_x as:
.. math::
\mathbf{W}_x = \sqrt{\\alpha_x}\\text{diag}\left(\sqrt{\mathbf{A}_1^\\top\mathbf{v}}\\right) \mathbf{D}_1
And then W as a tall matrix of all of the different regularization terms:
.. math::
\mathbf{W} = \left[ \\begin{array}{c}
\mathbf{W}_s\\\\
\mathbf{W}_x\\\\
\mathbf{W}_y\end{array} \\right]
Then we can write
.. math::
R(m) \\approx \\frac{1}{2}\mathbf{m^\\top W^\\top W m}
"""
"""
smoothModel = True #: SMOOTH and SMOOTH_MOD_DIF options
alpha_s = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Ws'], "Smallness weight")
@@ -205,14 +125,18 @@ class Tikhonov(BaseRegularization):
alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction")
alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction")
def __init__(self, mesh, mapping=None, **kwargs):
def __init__(self, mesh, mapping=None, indActive = None, **kwargs):
BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs)
self.indActive = indActive
@property
def Ws(self):
"""Regularization matrix Ws"""
if getattr(self,'_Ws', None) is None:
self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
if self.indActive is not None:
Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
self._Ws = Pac.T * self._Ws * Pac
return self._Ws
@property
@@ -221,6 +145,13 @@ class Tikhonov(BaseRegularization):
if getattr(self, '_Wx', None) is None:
Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFx].T*self.mesh.vol
self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx
if self.indActive is not None:
indActive_Fx = (self.mesh.aveFx2CC.T * self.indActive) == 1
Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
Pafx = Utils.speye(self.mesh.nFx)[:,indActive_Fx]
self._Wx = Pafx.T*self._Wx*Pac
return self._Wx
@property
@@ -229,6 +160,13 @@ class Tikhonov(BaseRegularization):
if getattr(self, '_Wy', None) is None:
Ave_y_vol = self.mesh.aveF2CC[:,self.mesh.nFx:np.sum(self.mesh.vnF[:2])].T*self.mesh.vol
self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady
if self.indActive is not None:
indActive_Fy = (self.mesh.aveFy2CC.T * self.indActive) == 1
Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
Pafy = Utils.speye(self.mesh.nFy)[:,indActive_Fy]
self._Wy = Pafy.T*self._Wy*Pac
return self._Wy
@property
@@ -237,6 +175,13 @@ class Tikhonov(BaseRegularization):
if getattr(self, '_Wz', None) is None:
Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.vnF[:2]):].T*self.mesh.vol
self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz
if self.indActive is not None:
indActive_Fz = (self.mesh.aveFz2CC.T * self.indActive) == 1
Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
Pafz = Utils.speye(self.mesh.nFz)[:,indActive_Fz]
self._Wz = Pafz.T*self._Wz*Pac
return self._Wz
@property
@@ -244,6 +189,11 @@ class Tikhonov(BaseRegularization):
"""Regularization matrix Wxx"""
if getattr(self, '_Wxx', None) is None:
self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx
if self.indActive is not None:
Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
self._Wxx = Pac.T*self._Wxx*Pac
return self._Wxx
@property
@@ -251,6 +201,11 @@ class Tikhonov(BaseRegularization):
"""Regularization matrix Wyy"""
if getattr(self, '_Wyy', None) is None:
self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady
if self.indActive is not None:
Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
self._Wyy = Pac.T*self._Wyy*Pac
return self._Wyy
@property
@@ -258,6 +213,11 @@ class Tikhonov(BaseRegularization):
"""Regularization matrix Wzz"""
if getattr(self, '_Wzz', None) is None:
self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz
if self.indActive is not None:
Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
self._Wzz = Pac.T*self._Wzz*Pac
return self._Wzz
@property
@@ -311,7 +271,7 @@ class Tikhonov(BaseRegularization):
if self.smoothModel == True:
mD1 = self.mapping.deriv(m)
mD2 = self.mapping.deriv(m - self.mref)
r1 = self.Wsmooth * ( self.mapping * (m))
r1 = self.Wsmooth * ( self.mapping * (m))
r2 = self.Ws * ( self.mapping * (m - self.mref) )
out1 = mD1.T * ( self.Wsmooth.T * r1 )
out2 = mD2.T * ( self.Ws.T * r2 )
SimPEG/Survey.py
+1
View File
@@ -205,6 +205,7 @@ class BaseSurvey(object):
__metaclass__ = Utils.SimPEGMetaClass
std = None #: Estimated Standard Deviations
eps = None #: Estimated Noise Floor
dobs = None #: Observed data
dtrue = None #: True data, if data is synthetic
mtrue = None #: True model, if data is synthetic
SimPEG/Tests.py
+1 -1
View File
@@ -1,5 +1,4 @@
import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import norm
from SimPEG.Utils import mkvc, sdiag, diagEst
from SimPEG import Utils
@@ -311,6 +310,7 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
if plotIt:
import matplotlib.pyplot as plt
ax = ax or plt.subplot(111)
ax.loglog(h, E0, 'b')
ax.loglog(h, E1, 'g--')
SimPEG/Utils/SolverUtils.py
+17 -4
View File
@@ -26,7 +26,14 @@ def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6):
def __init__(self, A, **kwargs):
self.A = A.tocsc()
self.checkAccuracy = kwargs.get("checkAccuracy", checkAccuracy)
if kwargs.has_key("checkAccuracy"): del kwargs["checkAccuracy"]
self.accuracyTol = kwargs.get("accuracyTol", accuracyTol)
if kwargs.has_key("accuracyTol"): del kwargs["accuracyTol"]
self.kwargs = kwargs
if factorize:
self.solver = fun(self.A, **kwargs)
@@ -57,8 +64,8 @@ def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6):
else:
X[:,i] = fun(self.A, b[:,i], **self.kwargs)
if checkAccuracy:
_checkAccuracy(self.A, b, X, accuracyTol)
if self.checkAccuracy:
_checkAccuracy(self.A, b, X, self.accuracyTol)
return X
def clean(self):
@@ -81,6 +88,12 @@ def SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-5):
def __init__(self, A, **kwargs):
self.A = A
self.checkAccuracy = kwargs.get("checkAccuracy", checkAccuracy)
if kwargs.has_key("checkAccuracy"): del kwargs["checkAccuracy"]
self.accuracyTol = kwargs.get("accuracyTol", accuracyTol)
if kwargs.has_key("accuracyTol"): del kwargs["accuracyTol"]
self.kwargs = kwargs
def __mul__(self, b):
@@ -108,8 +121,8 @@ def SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-5):
else:
X[:,i] = out
if checkAccuracy:
_checkAccuracy(self.A, b, X, accuracyTol)
if self.checkAccuracy:
_checkAccuracy(self.A, b, X, self.accuracyTol)
return X
def clean(self):
SimPEG/Utils/__init__.py
+1 -2
View File
@@ -1,9 +1,8 @@
from matutils import *
from codeutils import *
from meshutils import exampleLrmGrid, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCTensorModel, readVTRFile, writeVTRFile
from meshutils import *
from curvutils import volTetra, faceInfo, indexCube
from interputils import interpmat
from ipythonutils import easyAnimate as animate
from CounterUtils import *
import ModelBuilder
import SolverUtils
SimPEG/Utils/codeutils.py
+1 -1
View File
@@ -17,7 +17,7 @@ def memProfileWrapper(towrap, *funNames):
For example::
foo_mem = memProfile(foo,'my_func')
foo_mem = memProfileWrapper(foo,['my_func'])
fooi = foo_mem()
for i in range(5):
fooi.my_func()
SimPEG/Utils/ipythonutils.py
-28
View File
@@ -1,28 +0,0 @@
from tempfile import NamedTemporaryFile
import matplotlib.pyplot as plt
from matplotlib import animation
# http://jakevdp.github.io/blog/2013/05/12/embedding-matplotlib-animations/
# http://www.renevolution.com/how-to-install-ffmpeg-on-mac-os-x/
VIDEO_TAG = """<video controls loop>
<source src="data:video/x-m4v;base64,{0}" type="video/mp4">
Your browser does not support the video tag.
</video>"""
def anim_to_html(anim):
if not hasattr(anim, '_encoded_video'):
with NamedTemporaryFile(suffix='.mp4') as f:
anim.save(f.name, fps=20, extra_args=['-vcodec', 'libx264', '-pix_fmt', 'yuv420p'])
video = open(f.name, "rb").read()
anim._encoded_video = video.encode("base64")
return VIDEO_TAG.format(anim._encoded_video)
def display_animation(anim):
plt.close(anim._fig)
return anim_to_html(anim)
animation.Animation._repr_html_ = display_animation
easyAnimate = animation.FuncAnimation
SimPEG/Utils/matutils.py
+12 -1
View File
@@ -2,7 +2,6 @@ import numpy as np
import scipy.sparse as sp
from codeutils import isScalar
def mkvc(x, numDims=1):
"""Creates a vector with the number of dimension specified
@@ -26,6 +25,9 @@ def mkvc(x, numDims=1):
if hasattr(x, 'tovec'):
x = x.tovec()
if isinstance(x, Zero):
return x
assert isinstance(x, np.ndarray), "Vector must be a numpy array"
if numDims == 1:
@@ -37,6 +39,9 @@ def mkvc(x, numDims=1):
def sdiag(h):
"""Sparse diagonal matrix"""
if isinstance(h, Zero):
return Zero()
return sp.spdiags(mkvc(h), 0, h.size, h.size, format="csr")
def sdInv(M):
@@ -417,6 +422,12 @@ class Zero(object):
def __ge__(self, v):return 0 >= v
def __gt__(self, v):return 0 > v
@property
def transpose(self): return Zero()
@property
def T(self): return Zero()
class Identity(object):
_positive = True
def __init__(self, positive=True):
SimPEG/Utils/meshutils.py
-217
View File
@@ -102,223 +102,6 @@ def closestPoints(mesh, pts, gridLoc='CC'):
return nodeInds
def readUBCTensorMesh(fileName):
"""
Read UBC GIF 3DTensor mesh and generate 3D Tensor mesh in simpegTD
Input:
:param fileName, path to the UBC GIF mesh file
Output:
:param SimPEG TensorMesh object
:return
"""
# Interal function to read cell size lines for the UBC mesh files.
def readCellLine(line):
for seg in line.split():
if '*' in seg:
st = seg
sp = seg.split('*')
re = np.array(sp[0],dtype=int)*(' ' + sp[1])
line = line.replace(st,re.strip())
return np.array(line.split(),dtype=float)
# Read the file as line strings, remove lines with comment = !
msh = np.genfromtxt(fileName,delimiter='\n',dtype=np.str,comments='!')
# Fist line is the size of the model
sizeM = np.array(msh[0].split(),dtype=float)
# Second line is the South-West-Top corner coordinates.
x0 = np.array(msh[1].split(),dtype=float)
# Read the cell sizes
h1 = readCellLine(msh[2])
h2 = readCellLine(msh[3])
h3temp = readCellLine(msh[4])
h3 = h3temp[::-1] # Invert the indexing of the vector to start from the bottom.
# Adjust the reference point to the bottom south west corner
x0[2] = x0[2] - np.sum(h3)
# Make the mesh
from SimPEG import Mesh
tensMsh = Mesh.TensorMesh([h1,h2,h3],x0)
return tensMsh
def readUBCTensorModel(fileName, mesh):
"""
Read UBC 3DTensor mesh model and generate 3D Tensor mesh model in simpeg
Input:
:param fileName, path to the UBC GIF mesh file to read
:param mesh, TensorMesh object, mesh that coresponds to the model
Output:
:return numpy array, model with TensorMesh ordered
"""
f = open(fileName, 'r')
model = np.array(map(float, f.readlines()))
f.close()
model = np.reshape(model, (mesh.nCz, mesh.nCx, mesh.nCy), order = 'F')
model = model[::-1,:,:]
model = np.transpose(model, (1, 2, 0))
model = mkvc(model)
return model
def writeUBCTensorMesh(fileName, mesh):
"""
Writes a SimPEG TensorMesh to a UBC-GIF format mesh file.
:param str fileName: File to write to
:param simpeg.Mesh.TensorMesh mesh: The mesh
"""
assert mesh.dim == 3
s = ''
s += '%i %i %i\n' %tuple(mesh.vnC)
origin = mesh.x0 + np.array([0,0,mesh.hz.sum()]) # Have to it in the same operation or use mesh.x0.copy(), otherwise the mesh.x0 is updated.
origin.dtype = float
s += '%.2f %.2f %.2f\n' %tuple(origin)
s += ('%.2f '*mesh.nCx+'\n')%tuple(mesh.hx)
s += ('%.2f '*mesh.nCy+'\n')%tuple(mesh.hy)
s += ('%.2f '*mesh.nCz+'\n')%tuple(mesh.hz[::-1])
f = open(fileName, 'w')
f.write(s)
f.close()
def writeUBCTensorModel(fileName, mesh, model):
"""
Writes a model associated with a SimPEG TensorMesh
to a UBC-GIF format model file.
:param str fileName: File to write to
:param simpeg.Mesh.TensorMesh mesh: The mesh
:param numpy.ndarray model: The model
"""
# Reshape model to a matrix
modelMat = mesh.r(model,'CC','CC','M')
# Transpose the axes
modelMatT = modelMat.transpose((2,0,1))
# Flip z to positive down
modelMatTR = mkvc(modelMatT[::-1,:,:])
np.savetxt(fileName, modelMatTR.ravel())
def readVTRFile(fileName):
"""
Read VTK Rectilinear (vtr xml file) and return SimPEG Tensor mesh and model
Input:
:param vtrFileName, path to the vtr model file to write to
Output:
:return SimPEG TensorMesh object
:return SimPEG model dictionary
"""
# Import
from vtk import vtkXMLRectilinearGridReader as vtrFileReader
from vtk.util.numpy_support import vtk_to_numpy
# Read the file
vtrReader = vtrFileReader()
vtrReader.SetFileName(fileName)
vtrReader.Update()
vtrGrid = vtrReader.GetOutput()
# Sort information
hx = np.abs(np.diff(vtk_to_numpy(vtrGrid.GetXCoordinates())))
xR = vtk_to_numpy(vtrGrid.GetXCoordinates())[0]
hy = np.abs(np.diff(vtk_to_numpy(vtrGrid.GetYCoordinates())))
yR = vtk_to_numpy(vtrGrid.GetYCoordinates())[0]
zD = np.diff(vtk_to_numpy(vtrGrid.GetZCoordinates()))
# Check the direction of hz
if np.all(zD < 0):
hz = np.abs(zD[::-1])
zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[-1]
else:
hz = np.abs(zD)
zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[0]
x0 = np.array([xR,yR,zR])
# Make the SimPEG object
from SimPEG import Mesh
tensMsh = Mesh.TensorMesh([hx,hy,hz],x0)
# Grap the models
modelDict = {}
for i in np.arange(vtrGrid.GetCellData().GetNumberOfArrays()):
modelName = vtrGrid.GetCellData().GetArrayName(i)
if np.all(zD < 0):
modFlip = vtk_to_numpy(vtrGrid.GetCellData().GetArray(i))
tM = tensMsh.r(modFlip,'CC','CC','M')
modArr = tensMsh.r(tM[:,:,::-1],'CC','CC','V')
else:
modArr = vtk_to_numpy(vtrGrid.GetCellData().GetArray(i))
modelDict[modelName] = modArr
# Return the data
return tensMsh, modelDict
def writeVTRFile(fileName,mesh,model=None):
"""
Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model.
Input:
:param str, path to the output vtk file
:param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK
:param model, dictionary of numpy.array - Name('s) and array('s). Match number of cells
"""
# Import
from vtk import vtkRectilinearGrid as rectGrid, vtkXMLRectilinearGridWriter as rectWriter
from vtk.util.numpy_support import numpy_to_vtk
# Deal with dimensionalities
if mesh.dim >= 1:
vX = mesh.vectorNx
xD = mesh.nNx
yD,zD = 1,1
vY, vZ = np.array([0,0])
if mesh.dim >= 2:
vY = mesh.vectorNy
yD = mesh.nNy
if mesh.dim == 3:
vZ = mesh.vectorNz
zD = mesh.nNz
# Use rectilinear VTK grid.
# Assign the spatial information.
vtkObj = rectGrid()
vtkObj.SetDimensions(xD,yD,zD)
vtkObj.SetXCoordinates(numpy_to_vtk(vX,deep=1))
vtkObj.SetYCoordinates(numpy_to_vtk(vY,deep=1))
vtkObj.SetZCoordinates(numpy_to_vtk(vZ,deep=1))
# Assign the model('s) to the object
for item in model.iteritems():
# Convert numpy array
vtkDoubleArr = numpy_to_vtk(item[1],deep=1)
vtkDoubleArr.SetName(item[0])
vtkObj.GetCellData().AddArray(vtkDoubleArr)
# Set the active scalar
vtkObj.GetCellData().SetActiveScalars(model.keys()[0])
vtkObj.Update()
# Check the extension of the fileName
ext = os.path.splitext(fileName)[1]
if ext is '':
fileName = fileName + '.vtr'
elif ext not in '.vtr':
raise IOError('{:s} is an incorrect extension, has to be .vtr')
# Write the file.
vtrWriteFilter = rectWriter()
vtrWriteFilter.SetInput(vtkObj)
vtrWriteFilter.SetFileName(fileName)
vtrWriteFilter.Update()
def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
"""
Extracts Core Mesh from Global mesh
SimPEG/__init__.py
+1 -1
View File
@@ -15,7 +15,7 @@ import Directives
import Inversion
import Tests
__version__ = '0.1.3'
__version__ = '0.1.9'
__author__ = 'Rowan Cockett'
__license__ = 'MIT'
__copyright__ = 'Copyright 2014 Rowan Cockett'