mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-27 11:27:23 +08:00
Merge branch 'master' into em/dev
This commit is contained in:
Lindsey Heagy
@@ -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'))
|
||||
|
||||
@@ -1 +0,0 @@
|
||||
import CylInversion
|
||||
@@ -1,9 +1,16 @@
|
||||
from SimPEG import *
|
||||
import SimPEG.EM as EM
|
||||
from scipy.constants import mu_0
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
|
||||
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)]
|
||||
@@ -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)
|
||||
@@ -53,6 +61,7 @@ def run(plotIt=True):
|
||||
survey.Wd = 1/(abs(survey.dobs)*std)
|
||||
|
||||
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.dobs, 'r.-')
|
||||
@@ -81,6 +90,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])
|
||||
@@ -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()
|
||||
@@ -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__':
|
||||
@@ -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()
|
||||
@@ -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()
|
||||
@@ -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()
|
||||
@@ -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()
|
||||
@@ -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()
|
||||
@@ -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()
|
||||
@@ -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()
|
||||
@@ -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()
|
||||
|
||||
+103
-1
@@ -1 +1,103 @@
|
||||
import Linear, DCfwd
|
||||
# Run this file to add imports.
|
||||
|
||||
##### AUTOIMPORTS #####
|
||||
import EM_FDEM_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", "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.
|
||||
docExamplesDir = os.path.sep.join(os.path.realpath(__file__).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(__file__.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(__file__, '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(__file__, '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)
|
||||
|
||||
@@ -1 +0,0 @@
|
||||
import Celia1990
|
||||
+69
-25
@@ -90,13 +90,14 @@
|
||||
#
|
||||
|
||||
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
|
||||
import time
|
||||
@@ -108,6 +109,8 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
|
||||
_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."
|
||||
|
||||
@@ -1960,11 +1963,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 +1985,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 +2018,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 +2072,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 +2087,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 +2117,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 +2132,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'
|
||||
@@ -2191,7 +2228,7 @@ class Cell(object):
|
||||
@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)
|
||||
@@ -2248,6 +2285,13 @@ class CellLookUpException(TreeException):
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib
|
||||
from mpl_toolkits.mplot3d import Axes3D
|
||||
import matplotlib.colors as colors
|
||||
import matplotlib.cm as cmx
|
||||
|
||||
def topo(x):
|
||||
return np.sin(x*(2.*np.pi))*0.3 + 0.5
|
||||
|
||||
|
||||
+7
-61
@@ -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):
|
||||
|
||||
|
||||
@@ -24,7 +24,7 @@ class BaseRegularization(object):
|
||||
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)
|
||||
|
||||
@property
|
||||
@@ -115,86 +115,7 @@ class BaseRegularization(object):
|
||||
|
||||
|
||||
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")
|
||||
@@ -311,7 +232,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 )
|
||||
|
||||
+1
-1
@@ -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--')
|
||||
|
||||
@@ -3,7 +3,6 @@ from codeutils import *
|
||||
from meshutils import exampleLrmGrid, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCTensorModel, readVTRFile, writeVTRFile
|
||||
from curvutils import volTetra, faceInfo, indexCube
|
||||
from interputils import interpmat
|
||||
from ipythonutils import easyAnimate as animate
|
||||
from CounterUtils import *
|
||||
import ModelBuilder
|
||||
import SolverUtils
|
||||
|
||||
@@ -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
|
||||
+1
-1
@@ -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'
|
||||
|
||||
Reference in New Issue
Block a user