From 5e5c7ba0fb408645eb5c2fa85340b3af75f01098 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 1 Mar 2016 17:31:37 -0800 Subject: [PATCH] docs for regmesh, cellGrad--> cellDiff, faceDiv--> faceDiff for regmesh --- SimPEG/Regularization.py | 238 ++++++++++++++++++++++++++++----------- 1 file changed, 174 insertions(+), 64 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index d0fc21b1..acd6c68d 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -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