mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-05 07:05:16 +08:00
docs for regmesh, cellGrad--> cellDiff, faceDiv--> faceDiff for regmesh
This commit is contained in:
Lindsey Heagy
+174
-64
@@ -1,6 +1,16 @@
|
||||
import Utils, Maps, Mesh, numpy as np, scipy.sparse as sp
|
||||
|
||||
class RegularizationMesh(object):
|
||||
"""
|
||||
**Regularization Mesh**
|
||||
|
||||
This contains the operators used in the regularization. Note that these
|
||||
are not necessarily true differential operators, but are constructed from
|
||||
a SimPEG Mesh.
|
||||
|
||||
:param Mesh mesh: problem mesh
|
||||
:param numpy.array indActive: bool array, size nC, that is True where we have active cells. Used to reduce the operators so we regularize only on active cells
|
||||
"""
|
||||
|
||||
def __init__(self, mesh, indActive=None):
|
||||
self.mesh = mesh
|
||||
@@ -9,12 +19,22 @@ class RegularizationMesh(object):
|
||||
|
||||
@property
|
||||
def vol(self):
|
||||
"""
|
||||
reduced volume vector
|
||||
:rtype: numpy.array
|
||||
:return: reduced cell volume
|
||||
"""
|
||||
if getattr(self, '_vol', None) is None:
|
||||
self._vol = self._Pac.T * self.mesh.vol
|
||||
return self._vol
|
||||
|
||||
@property
|
||||
def nC(self):
|
||||
"""
|
||||
reduced number of cells
|
||||
:rtype: int
|
||||
:return: number of cells being regularized
|
||||
"""
|
||||
if getattr(self, '_nC', None) is None:
|
||||
if self.indActive is None:
|
||||
self._nC = self.mesh.nC
|
||||
@@ -24,6 +44,11 @@ class RegularizationMesh(object):
|
||||
|
||||
@property
|
||||
def dim(self):
|
||||
"""
|
||||
dimension of regularization mesh (1D, 2D, 3D)
|
||||
:rtype: int
|
||||
:return: dimension
|
||||
"""
|
||||
if getattr(self, '_dim', None) is None:
|
||||
self._dim = self.mesh.dim
|
||||
return self._dim
|
||||
@@ -31,6 +56,11 @@ class RegularizationMesh(object):
|
||||
|
||||
@property
|
||||
def _Pac(self):
|
||||
"""
|
||||
projection matrix that takes from the reduced space of active cells to full modelling space (ie. nC x nindActive)
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: active cell projection matrix
|
||||
"""
|
||||
if getattr(self, '__Pac', None) is None:
|
||||
if self.indActive is None:
|
||||
self.__Pac = Utils.speye(self.mesh.nC)
|
||||
@@ -40,6 +70,11 @@ class RegularizationMesh(object):
|
||||
|
||||
@property
|
||||
def _Pafx(self):
|
||||
"""
|
||||
projection matrix that takes from the reduced space of active x-faces to full modelling space (ie. nFx x nindActive_Fx )
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: active face-x projection matrix
|
||||
"""
|
||||
if getattr(self, '__Pafx', None) is None:
|
||||
if self.indActive is None:
|
||||
self.__Pafx = Utils.speye(self.mesh.nFx)
|
||||
@@ -50,6 +85,11 @@ class RegularizationMesh(object):
|
||||
|
||||
@property
|
||||
def _Pafy(self):
|
||||
"""
|
||||
projection matrix that takes from the reduced space of active y-faces to full modelling space (ie. nFy x nindActive_Fy )
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: active face-y projection matrix
|
||||
"""
|
||||
if getattr(self, '__Pafy', None) is None:
|
||||
if self.indActive is None:
|
||||
self.__Pafy = Utils.speye(self.mesh.nFy)
|
||||
@@ -60,6 +100,11 @@ class RegularizationMesh(object):
|
||||
|
||||
@property
|
||||
def _Pafz(self):
|
||||
"""
|
||||
projection matrix that takes from the reduced space of active z-faces to full modelling space (ie. nFz x nindActive_Fz )
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: active face-z projection matrix
|
||||
"""
|
||||
if getattr(self, '__Pafz', None) is None:
|
||||
if self.indActive is None:
|
||||
self.__Pafz = Utils.speye(self.mesh.nFz)
|
||||
@@ -70,108 +115,173 @@ class RegularizationMesh(object):
|
||||
|
||||
@property
|
||||
def aveFx2CC(self):
|
||||
"""
|
||||
averaging from active cell centers to active x-faces
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: averaging from active cell centers to active x-faces
|
||||
"""
|
||||
if getattr(self, '_aveFx2CC', None) is None:
|
||||
self._aveFx2CC = self._Pac.T * self.mesh.aveFx2CC * self._Pafx
|
||||
return self._aveFx2CC
|
||||
|
||||
@property
|
||||
def aveCC2Fx(self):
|
||||
"""
|
||||
averaging from active x-faces to active cell centers
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: averaging matrix from active x-faces to active cell centers
|
||||
"""
|
||||
if getattr(self, '_aveCC2Fx', None) is None:
|
||||
self._aveCC2Fx = Utils.sdiag(1./(self.aveFx2CC.T).sum(1)) * self.aveFx2CC.T
|
||||
return self._aveCC2Fx
|
||||
|
||||
@property
|
||||
def aveFy2CC(self):
|
||||
"""
|
||||
averaging from active cell centers to active y-faces
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: averaging from active cell centers to active y-faces
|
||||
"""
|
||||
if getattr(self, '_aveFy2CC', None) is None:
|
||||
self._aveFy2CC = self._Pac.T * self.mesh.aveFy2CC * self._Pafy
|
||||
return self._aveFy2CC
|
||||
|
||||
@property
|
||||
def aveCC2Fy(self):
|
||||
"""
|
||||
averaging from active y-faces to active cell centers
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: averaging matrix from active y-faces to active cell centers
|
||||
"""
|
||||
if getattr(self, '_aveCC2Fy', None) is None:
|
||||
self._aveCC2Fy = Utils.sdiag(1./(self.aveFy2CC.T).sum(1)) * self.aveFy2CC.T
|
||||
return self._aveCC2Fy
|
||||
|
||||
@property
|
||||
def aveFz2CC(self):
|
||||
"""
|
||||
averaging from active cell centers to active z-faces
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: averaging from active cell centers to active z-faces
|
||||
"""
|
||||
if getattr(self, '_aveFz2CC', None) is None:
|
||||
self._aveFz2CC = self._Pac.T * self.mesh.aveFz2CC * self._Pafz
|
||||
return self._aveFz2CC
|
||||
|
||||
@property
|
||||
def aveCC2Fz(self):
|
||||
"""
|
||||
averaging from active z-faces to active cell centers
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: averaging matrix from active z-faces to active cell centers
|
||||
"""
|
||||
if getattr(self, '_aveCC2Fz', None) is None:
|
||||
self._aveCC2Fz = Utils.sdiag(1./(self.aveFz2CC.T).sum(1)) * self.aveFz2CC.T
|
||||
return self._aveCC2Fz
|
||||
|
||||
@property
|
||||
def cellGradx(self):
|
||||
if getattr(self, '_cellGradx', None) is None:
|
||||
self._cellGradx = self._Pafx.T * self.mesh.cellGradx * self._Pac
|
||||
return self._cellGradx
|
||||
def cellDiffx(self):
|
||||
"""
|
||||
cell centered difference in the x-direction
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: differencing matrix for active cells in the x-direction
|
||||
"""
|
||||
if getattr(self, '_cellDiffx', None) is None:
|
||||
self._cellDiffx = self._Pafx.T * self.mesh.cellGradx * self._Pac
|
||||
return self._cellDiffx
|
||||
|
||||
@property
|
||||
def cellGrady(self):
|
||||
if getattr(self, '_cellGrady', None) is None:
|
||||
self._cellGrady = self._Pafy.T * self.mesh.cellGrady * self._Pac
|
||||
return self._cellGrady
|
||||
def cellDiffy(self):
|
||||
"""
|
||||
cell centered difference in the y-direction
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: differencing matrix for active cells in the y-direction
|
||||
"""
|
||||
if getattr(self, '_cellDiffy', None) is None:
|
||||
self._cellDiffy = self._Pafy.T * self.mesh.cellGrady * self._Pac
|
||||
return self._cellDiffy
|
||||
|
||||
@property
|
||||
def cellGradz(self):
|
||||
if getattr(self, '_cellGradz', None) is None:
|
||||
self._cellGradz = self._Pafz.T * self.mesh.cellGradz * self._Pac
|
||||
return self._cellGradz
|
||||
def cellDiffz(self):
|
||||
"""
|
||||
cell centered difference in the z-direction
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: differencing matrix for active cells in the z-direction
|
||||
"""
|
||||
if getattr(self, '_cellDiffz', None) is None:
|
||||
self._cellDiffz = self._Pafz.T * self.mesh.cellGradz * self._Pac
|
||||
return self._cellDiffz
|
||||
|
||||
@property
|
||||
def faceDivx(self):
|
||||
if getattr(self, '_faceDivx', None) is None:
|
||||
self._faceDivx = self._Pac.T * self.mesh.faceDivx * self._Pafx
|
||||
return self._faceDivx
|
||||
def faceDiffx(self):
|
||||
"""
|
||||
x-face differences
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: differencing matrix for active faces in the x-direction
|
||||
"""
|
||||
if getattr(self, '_faceDiffx', None) is None:
|
||||
self._faceDiffx = self._Pac.T * self.mesh.faceDivx * self._Pafx
|
||||
return self._faceDiffx
|
||||
|
||||
@property
|
||||
def faceDivy(self):
|
||||
if getattr(self, '_faceDivy', None) is None:
|
||||
self._faceDivy = self._Pac.T * self.mesh.faceDivy * self._Pafy
|
||||
return self._faceDivy
|
||||
def faceDiffy(self):
|
||||
"""
|
||||
y-face differences
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: differencing matrix for active faces in the y-direction
|
||||
"""
|
||||
if getattr(self, '_faceDiffy', None) is None:
|
||||
self._faceDiffy = self._Pac.T * self.mesh.faceDivy * self._Pafy
|
||||
return self._faceDiffy
|
||||
|
||||
@property
|
||||
def faceDivz(self):
|
||||
if getattr(self, '_faceDivz', None) is None:
|
||||
self._faceDivz = self._Pac.T * self.mesh.faceDivz * self._Pafz
|
||||
return self._faceDivz
|
||||
def faceDiffz(self):
|
||||
"""
|
||||
z-face differences
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: differencing matrix for active faces in the z-direction
|
||||
"""
|
||||
if getattr(self, '_faceDiffz', None) is None:
|
||||
self._faceDiffz = self._Pac.T * self.mesh.faceDivz * self._Pafz
|
||||
return self._faceDiffz
|
||||
|
||||
@property
|
||||
def cellGradxStencil(self):
|
||||
"""Cell centered Gradient in the x dimension used for
|
||||
regularization. The gradient operator is square (nC-by-nC)"""
|
||||
def cellDiffxStencil(self):
|
||||
"""
|
||||
cell centered difference stencil (no cell lengths include) in the x-direction
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: differencing matrix for active cells in the x-direction
|
||||
"""
|
||||
if getattr(self, '_cellDiffxStencil', None) is None:
|
||||
|
||||
# if self.dim < 3: return None
|
||||
if getattr(self, '_cellGradxStencil', None) is None:
|
||||
|
||||
self._cellGradxStencil = self._Pafx.T * self.mesh._cellGradxStencil() * self._Pac
|
||||
return self._cellGradxStencil
|
||||
self._cellDiffxStencil = self._Pafx.T * self.mesh._cellGradxStencil() * self._Pac
|
||||
return self._cellDiffxStencil
|
||||
|
||||
@property
|
||||
def cellGradyStencil(self):
|
||||
"""Cell centered Gradient in the x dimension used for
|
||||
regularization. The gradient operator is square (nC-by-nC)"""
|
||||
|
||||
def cellDiffyStencil(self):
|
||||
"""
|
||||
cell centered difference stencil (no cell lengths include) in the y-direction
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: differencing matrix for active cells in the y-direction
|
||||
"""
|
||||
if self.dim < 2: return None
|
||||
if getattr(self, '_cellGradyStencil', None) is None:
|
||||
if getattr(self, '_cellDiffyStencil', None) is None:
|
||||
|
||||
self._cellGradyStencil = self._Pafy.T * self.mesh._cellGradyStencil() * self._Pac
|
||||
return self._cellGradyStencil
|
||||
self._cellDiffyStencil = self._Pafy.T * self.mesh._cellGradyStencil() * self._Pac
|
||||
return self._cellDiffyStencil
|
||||
|
||||
@property
|
||||
def cellGradzStencil(self):
|
||||
"""Cell centered Gradient in the x dimension used for
|
||||
regularization. The gradient operator is square (nC-by-nC)"""
|
||||
|
||||
def cellDiffzStencil(self):
|
||||
"""
|
||||
cell centered difference stencil (no cell lengths include) in the y-direction
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: differencing matrix for active cells in the y-direction
|
||||
"""
|
||||
if self.dim < 3: return None
|
||||
if getattr(self, '_cellGradzStencil', None) is None:
|
||||
if getattr(self, '_cellDiffzStencil', None) is None:
|
||||
|
||||
self._cellGradzStencil = self._Pafz.T * self.mesh._cellGradzStencil() * self._Pac
|
||||
return self._cellGradzStencil
|
||||
self._cellDiffzStencil = self._Pafz.T * self.mesh._cellGradzStencil() * self._Pac
|
||||
return self._cellDiffzStencil
|
||||
|
||||
|
||||
class BaseRegularization(object):
|
||||
@@ -317,7 +427,7 @@ class Tikhonov(BaseRegularization):
|
||||
"""Regularization matrix Wx"""
|
||||
if getattr(self, '_Wx', None) is None:
|
||||
Ave_x_vol = self.regmesh.aveCC2Fx * self.regmesh.vol
|
||||
self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.regmesh.cellGradx
|
||||
self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.regmesh.cellDiffx
|
||||
return self._Wx
|
||||
|
||||
@property
|
||||
@@ -325,7 +435,7 @@ class Tikhonov(BaseRegularization):
|
||||
"""Regularization matrix Wy"""
|
||||
if getattr(self, '_Wy', None) is None:
|
||||
Ave_y_vol = self.regmesh.aveCC2Fy * self.regmesh.vol
|
||||
self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.regmesh.cellGrady
|
||||
self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.regmesh.cellDiffy
|
||||
return self._Wy
|
||||
|
||||
@property
|
||||
@@ -333,28 +443,28 @@ class Tikhonov(BaseRegularization):
|
||||
"""Regularization matrix Wz"""
|
||||
if getattr(self, '_Wz', None) is None:
|
||||
Ave_z_vol = self.regmesh.aveCC2Fz * self.regmesh.vol
|
||||
self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.regmesh.cellGradz
|
||||
self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.regmesh.cellDiffz
|
||||
return self._Wz
|
||||
|
||||
@property
|
||||
def Wxx(self):
|
||||
"""Regularization matrix Wxx"""
|
||||
if getattr(self, '_Wxx', None) is None:
|
||||
self._Wxx = Utils.sdiag((self.regmesh.vol*self.alpha_xx)**0.5)*self.regmesh.faceDivx*self.regmesh.cellGradx
|
||||
self._Wxx = Utils.sdiag((self.regmesh.vol*self.alpha_xx)**0.5)*self.regmesh.faceDiffx*self.regmesh.cellDiffx
|
||||
return self._Wxx
|
||||
|
||||
@property
|
||||
def Wyy(self):
|
||||
"""Regularization matrix Wyy"""
|
||||
if getattr(self, '_Wyy', None) is None:
|
||||
self._Wyy = Utils.sdiag((self.regmesh.vol*self.alpha_yy)**0.5)*self.regmesh.faceDivy*self.regmesh.cellGrady
|
||||
self._Wyy = Utils.sdiag((self.regmesh.vol*self.alpha_yy)**0.5)*self.regmesh.faceDiffy*self.regmesh.cellDiffy
|
||||
return self._Wyy
|
||||
|
||||
@property
|
||||
def Wzz(self):
|
||||
"""Regularization matrix Wzz"""
|
||||
if getattr(self, '_Wzz', None) is None:
|
||||
self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDivz*self.regmesh.cellGradz
|
||||
self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDiffz*self.regmesh.cellDiffz
|
||||
return self._Wzz
|
||||
|
||||
@property
|
||||
@@ -445,21 +555,21 @@ class Simple(BaseRegularization):
|
||||
def Wx(self):
|
||||
"""Regularization matrix Wx"""
|
||||
if getattr(self, '_Wx', None) is None:
|
||||
self._Wx = Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x)**0.5)*self.regmesh.cellGradxStencil
|
||||
self._Wx = Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x)**0.5)*self.regmesh.cellDiffxStencil
|
||||
return self._Wx
|
||||
|
||||
@property
|
||||
def Wy(self):
|
||||
"""Regularization matrix Wy"""
|
||||
if getattr(self, '_Wy', None) is None:
|
||||
self._Wy = Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol * self.alpha_y)**0.5)*self.regmesh.cellGradyStencil
|
||||
self._Wy = Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol * self.alpha_y)**0.5)*self.regmesh.cellDiffyStencil
|
||||
return self._Wy
|
||||
|
||||
@property
|
||||
def Wz(self):
|
||||
"""Regularization matrix Wz"""
|
||||
if getattr(self, '_Wz', None) is None:
|
||||
self._Wz = Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z)**0.5)*self.regmesh.cellGradzStencil
|
||||
self._Wz = Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z)**0.5)*self.regmesh.cellDiffzStencil
|
||||
return self._Wz
|
||||
|
||||
@property
|
||||
@@ -563,15 +673,15 @@ class SparseRegularization(Simple):
|
||||
"""Regularization matrix Wx"""
|
||||
|
||||
if getattr(self, 'm', None) is None:
|
||||
self.Rx = Utils.speye(self.regmesh.cellGradxStencil.shape[0])
|
||||
self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
|
||||
|
||||
else:
|
||||
f_m = self.regmesh.cellGradxStencil * self.m
|
||||
f_m = self.regmesh.cellDiffxStencil * self.m
|
||||
self.rx = self.R( f_m , self.qx, self.eps)
|
||||
self.Rx = Utils.sdiag( self.rx )
|
||||
|
||||
if getattr(self, '_Wx', None) is None:
|
||||
self._Wx = Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma)**0.5)*self.Rx*self.regmesh.cellGradxStencil
|
||||
self._Wx = Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma)**0.5)*self.Rx*self.regmesh.cellDiffxStencil
|
||||
return self._Wx
|
||||
|
||||
@property
|
||||
@@ -579,15 +689,15 @@ class SparseRegularization(Simple):
|
||||
"""Regularization matrix Wy"""
|
||||
|
||||
if getattr(self, 'm', None) is None:
|
||||
self.Ry = Utils.speye(self.regmesh.cellGradyStencil.shape[0])
|
||||
self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
|
||||
|
||||
else:
|
||||
f_m = self.regmesh.cellGradyStencil * self.m
|
||||
f_m = self.regmesh.cellDiffyStencil * self.m
|
||||
self.ry = self.R( f_m , self.qy, self.eps)
|
||||
self.Ry = Utils.sdiag( self.ry )
|
||||
|
||||
if getattr(self, '_Wy', None) is None:
|
||||
self._Wy = Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma)**0.5)*self.Ry*self.regmesh.cellGradyStencil
|
||||
self._Wy = Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma)**0.5)*self.Ry*self.regmesh.cellDiffyStencil
|
||||
return self._Wy
|
||||
|
||||
@property
|
||||
@@ -595,15 +705,15 @@ class SparseRegularization(Simple):
|
||||
"""Regularization matrix Wz"""
|
||||
|
||||
if getattr(self, 'm', None) is None:
|
||||
self.Rz = Utils.speye(self.regmesh.cellGradzStencil.shape[0])
|
||||
self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
|
||||
|
||||
else:
|
||||
f_m = self.regmesh.cellGradzStencil * self.m
|
||||
f_m = self.regmesh.cellDiffzStencil * self.m
|
||||
self.rz = self.R( f_m , self.qz, self.eps)
|
||||
self.Rz = Utils.sdiag( self.rz )
|
||||
|
||||
if getattr(self, '_Wz', None) is None:
|
||||
self._Wz = Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma)**0.5)*self.Rz*self.regmesh.cellGradzStencil
|
||||
self._Wz = Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma)**0.5)*self.Rz*self.regmesh.cellDiffzStencil
|
||||
return self._Wz
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user