wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-27 19:48:52 +08:00
Code Issues Packages Projects Releases Wiki Activity

implement and test harmonic averaged derivatives

Browse Source
This commit is contained in:
rowanc1
2014-07-03 11:00:25 -07:00
parent bb0b9c1de4
commit 42278d35fe
3 changed files with 47 additions and 38 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Mesh/InnerProducts.py
+4 -6
View File
@@ -148,18 +148,16 @@ class InnerProducts(object):
:rtype: scipy.csr_matrix
:return: dMdm, the derivative of the inner product matrix (nE, nC*nA)
"""
tensorType = TensorType(self, prop)
fast = None
if hasattr(self, '_fastInnerProductDeriv') and doFast:
fast = self._fastInnerProductDeriv(projType, tensorType, invProp=invProp, invMat=invMat)
fast = self._fastInnerProductDeriv(projType, prop, invProp=invProp, invMat=invMat)
if fast is not None:
return fast
if invProp or invMat:
raise NotImplementedError('inverting the property or the matrix is not yet implemented for this mesh/tensorType. You should write it!')
if fast is not None:
return fast
tensorType = TensorType(self, prop)
P = self._getInnerProductProjectionMatrices(projType, tensorType=tensorType)
def innerProductDeriv(v):
return self._getInnerProductDerivFunction(tensorType, P, projType, v)
SimPEG/Mesh/TensorMesh.py
+31 -20
View File
@@ -300,46 +300,57 @@ class BaseTensorMesh(BaseRectangularMesh):
else:
return M
def _fastInnerProductDeriv(self, projType, tensorType):
def _fastInnerProductDeriv(self, projType, prop, invProp=False, invMat=False):
"""
:param str projType: 'E' or 'F'
:param TensorType tensorType: type of the tensor
:param bool invProp: inverts the material property
:param bool invMat: inverts the matrix
:rtype: function
:return: dMdmu, the derivative of the inner product matrix
"""
assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
if tensorType == -1:
return None
tensorType = Utils.TensorType(self, prop)
dMdprop = None
if invMat:
MI = self._fastInnerProduct(projType, prop, invProp=invProp, invMat=invMat)
if tensorType == 0:
Av = getattr(self, 'ave'+projType+'2CC')
V = Utils.sdiag(self.vol)
ones = sp.csr_matrix((np.ones(self.nC), (range(self.nC), np.zeros(self.nC))), shape=(self.nC,1))
# if v is None:
# return self.dim * Av.T * V * ones
dim_x_AvT_x_V_x_ones = self.dim * Av.T * V * ones
def scalarInnerProductDeriv(v):
return Utils.sdiag(v) * dim_x_AvT_x_V_x_ones
return scalarInnerProductDeriv
if not invMat and not invProp:
dMdprop = self.dim * Av.T * V * ones
elif invMat and invProp:
dMdprop = self.dim * Utils.sdiag(MI.diagonal()**2) * Av.T * V * ones * Utils.sdiag(1./prop**2)
if tensorType == 1:
Av = getattr(self, 'ave'+projType+'2CC')
V = Utils.sdiag(self.vol)
dim_x_AvT_x_V = self.dim * Av.T * V
def isotropicInnerProductDeriv(v=None):
if v is None:
print 'Depreciation Warning: TensorMesh.isotropicInnerProductDeriv. You should be supplying a vector. Returning: d * Av.T * V'
return dim_x_AvT_x_V
return Utils.sdiag(v) * dim_x_AvT_x_V
return isotropicInnerProductDeriv
if not invMat and not invProp:
dMdprop = self.dim * Av.T * V
elif invMat and invProp:
dMdprop = self.dim * Utils.sdiag(MI.diagonal()**2) * Av.T * V * Utils.sdiag(1./prop**2)
if tensorType == 2: # anisotropic
Av = getattr(self, 'ave'+projType+'2CCV')
V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol))
AvT_x_V = Av.T * V
def anisotropicInnerProductDeriv(v):
return Utils.sdiag(v) * AvT_x_V
return anisotropicInnerProductDeriv
if not invMat and not invProp:
dMdprop = Av.T * V
elif invMat and invProp:
dMdprop = Utils.sdiag(MI.diagonal()**2) * Av.T * V * Utils.sdiag(1./prop**2)
if dMdprop is not None:
def innerProductDeriv(v=None):
if v is None:
print 'Depreciation Warning: TensorMesh.innerProductDeriv. You should be supplying a vector. Use: sdiag(u)*dMdprop'
return dMdprop
return Utils.sdiag(v) * dMdprop
return innerProductDeriv
else:
return None
SimPEG/Tests/test_massMatricesDerivs.py
+12 -12
View File
@@ -20,7 +20,7 @@ class TestInnerProductsDerivs(unittest.TestCase):
M = mesh.getFaceInnerProduct(sig, invProp=invProp, invMat=invMat)
Md = mesh.getFaceInnerProductDeriv(sig, invProp=invProp, invMat=invMat, doFast=fast)
return M*v, Md(v)
print meshType, 'Face', h, rep, fast
print meshType, 'Face', h, rep, fast, ('harmonic' if invProp and invMat else 'standard')
return checkDerivative(fun, sig, num=5, plotIt=False)
def doTestEdge(self, h, rep, fast, meshType, invProp=False, invMat=False):
@@ -37,7 +37,7 @@ class TestInnerProductsDerivs(unittest.TestCase):
M = mesh.getEdgeInnerProduct(sig, invProp=invProp, invMat=invMat)
Md = mesh.getEdgeInnerProductDeriv(sig, invProp=invProp, invMat=invMat, doFast=fast)
return M*v, Md(v)
print meshType, 'Edge', h, rep, fast
print meshType, 'Edge', h, rep, fast, ('harmonic' if invProp and invMat else 'standard')
return checkDerivative(fun, sig, num=5, plotIt=False)
def test_FaceIP_1D_float(self):
@@ -118,22 +118,22 @@ class TestInnerProductsDerivs(unittest.TestCase):
# def test_FaceIP_1D_float_fast_harmonic(self):
# self.assertTrue(self.doTestFace([10],0, True, 'Tensor', invProp=True, invMat=True))
# def test_FaceIP_2D_float_fast_harmonic(self):
# self.assertTrue(self.doTestFace([10, 4],0, True, 'Tensor', invProp=True, invMat=True))
# def test_FaceIP_3D_float_fast_harmonic(self):
# self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Tensor', invProp=True, invMat=True))
def test_FaceIP_1D_float_fast_harmonic(self):
self.assertTrue(self.doTestFace([10],0, True, 'Tensor', invProp=True, invMat=True))
def test_FaceIP_2D_float_fast_harmonic(self):
self.assertTrue(self.doTestFace([10, 4],0, True, 'Tensor', invProp=True, invMat=True))
def test_FaceIP_3D_float_fast_harmonic(self):
self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Tensor', invProp=True, invMat=True))
def test_FaceIP_1D_isotropic_fast_harmonic(self):
self.assertTrue(self.doTestFace([10],1, True, 'Tensor', invProp=True, invMat=True))
def test_FaceIP_2D_isotropic_fast_harmonic(self):
self.assertTrue(self.doTestFace([10, 4],1, True, 'Tensor', invProp=True, invMat=True))
def test_FaceIP_3D_isotropic_fast_harmonic(self):
self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'Tensor', invProp=True, invMat=True))
# def test_FaceIP_2D_anisotropic_fast_harmonic(self):
# self.assertTrue(self.doTestFace([10, 4],2, True, 'Tensor', invProp=True, invMat=True))
# def test_FaceIP_3D_anisotropic_fast_harmonic(self):
# self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Tensor', invProp=True, invMat=True))
def test_FaceIP_2D_anisotropic_fast_harmonic(self):
self.assertTrue(self.doTestFace([10, 4],2, True, 'Tensor', invProp=True, invMat=True))
def test_FaceIP_3D_anisotropic_fast_harmonic(self):
self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Tensor', invProp=True, invMat=True))