diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py new file mode 100644 index 00000000..c5daa96d --- /dev/null +++ b/SimPEG/DiffOperators.py @@ -0,0 +1,260 @@ +import numpy as np +from scipy import sparse as sp +from sputils import sdiag, speye, kron3, spzeros +from utils import mkvc + + +def ddx(n): + """Define 1D derivatives, inner, this means we go from n+1 to n+1""" + return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr") + + +def checkBC(bc): + """ Checks if boundary condition 'bc' is valid. """ + if(type(bc) is str): + bc = [bc, bc] + assert type(bc) is list, 'bc must be a list' + assert len(bc) == 2, 'bc must have two elements' + + for bc_i in bc: + assert type(bc_i) is str, "each bc must be a string" + assert bc_i in ['dirichlet', 'neumann'], "each bc must be either, 'dirichlet' or 'neumann'" + return bc + + +def ddxCellGrad(n, bc): + """Create 1D derivative operator from cell-centres to nodes this means we go from n to n+1""" + bc = checkBC(bc) + + D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n, format="csr") + # Set the first side + if(bc[0] == 'dirichlet'): + D[0, 0] = 2 + elif(bc[0] == 'neumann'): + D[0, 0] = 0 + # Set the second side + if(bc[1] == 'dirichlet'): + D[-1, -1] = -2 + elif(bc[1] == 'neumann'): + D[-1, -1] = 0 + return D + + +def av(n): + """Define 1D averaging operator from cell-centres to nodes.""" + return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr") + + +class DiffOperators(object): + """ + Class creates the differential operators that you need! + """ + def __init__(self): + raise Exception('DiffOperators is a base class providing differential operators on meshes and cannot run on its own. Inherit to your favorite Mesh class.') + + def faceDiv(): + doc = "Construct divergence operator (face-stg to cell-centres)." + + def fget(self): + if(self._faceDiv is None): + # The number of cell centers in each direction + n = self.n + # Compute faceDivergence operator on faces + if(self.dim == 1): + D = ddx(n[0]) + elif(self.dim == 2): + D1 = sp.kron(speye(n[1]), ddx(n[0])) + D2 = sp.kron(ddx(n[1]), speye(n[0])) + D = sp.hstack((D1, D2), format="csr") + elif(self.dim == 3): + D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0])) + D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0])) + D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0])) + D = sp.hstack((D1, D2, D3), format="csr") + # Compute areas of cell faces & volumes + S = self.area + V = self.vol + self._faceDiv = sdiag(1/V)*D*sdiag(S) + + return self._faceDiv + return locals() + _faceDiv = None + faceDiv = property(**faceDiv()) + + def nodalGrad(): + doc = "Construct gradient operator (nodes to edges)." + + def fget(self): + if(self._nodalGrad is None): + # The number of cell centers in each direction + n = self.n + # Compute divergence operator on faces + if(self.dim == 1): + G = ddx(n[0]) + elif(self.dim == 2): + D1 = sp.kron(speye(n[1]+1), ddx(n[0])) + D2 = sp.kron(ddx(n[1]), speye(n[0]+1)) + G = sp.vstack((D1, D2), format="csr") + elif(self.dim == 3): + D1 = kron3(speye(n[2]+1), speye(n[1]+1), ddx(n[0])) + D2 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0]+1)) + D3 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]+1)) + G = sp.vstack((D1, D2, D3), format="csr") + # Compute lengths of cell edges + L = self.edge + self._nodalGrad = sdiag(1/L)*G + return self._nodalGrad + return locals() + _nodalGrad = None + nodalGrad = property(**nodalGrad()) + + def setCellGradBC(self, BC): + """ + Function that sets the boundary conditions for cell-centred derivative operators. + + Examples: + + BC = 'neumann' # Neumann in all directions + BC = ['neumann', 'dirichlet', 'neumann'] # 3D, Dirichlet in y Neumann else + BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet'] # 3D, Neumann in x on bottom of domain, + # Dirichlet else + + """ + if(type(BC) is str): + BC = [BC for _ in self.n] # Repeat the str self.dim times + elif(type(BC) is list): + assert len(BC) == self.dim, 'BC list must be the size of your mesh' + else: + raise Exception("BC must be a str or a list.") + + for i, bc_i in enumerate(BC): + BC[i] = checkBC(bc_i) + + self._cellGrad = None # ensure we create a new gradient next time we call it + self._cellGradBC = BC + return BC + _cellGradBC = 'neumann' + + def cellGrad(): + doc = "The cell centered Gradient, takes you to cell faces." + + def fget(self): + if(self._cellGrad is None): + BC = self.setCellGradBC(self._cellGradBC) + n = self.n + if(self.dim == 1): + G = ddxCellGrad(n[0], BC[0]) + elif(self.dim == 2): + G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC[0])) + G2 = sp.kron(ddxCellGrad(n[1], BC[1]), speye(n[0])) + G = sp.vstack((G1, G2), format="csr") + elif(self.dim == 3): + G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC[0])) + G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC[1]), speye(n[0])) + G3 = kron3(ddxCellGrad(n[2], BC[2]), speye(n[1]), speye(n[0])) + G = sp.vstack((G1, G2, G3), format="csr") + # Compute areas of cell faces & volumes + S = self.area + V = self.vol + self._cellGrad = sdiag(S)*G*sdiag(1/V) + return self._cellGrad + return locals() + _cellGrad = None + cellGrad = property(**cellGrad()) + + def edgeCurl(): + doc = "Construct the 3D curl operator." + + def fget(self): + if(self._edgeCurl is None): + # The number of cell centers in each direction + n1 = np.size(self.hx) + n2 = np.size(self.hy) + n3 = np.size(self.hz) + + # Compute lengths of cell edges + L = self.edge + + # Compute areas of cell faces + S = self.area + + # Compute divergence operator on faces + d1 = ddx(n1) + d2 = ddx(n2) + d3 = ddx(n3) + + D32 = kron3(d3, speye(n2), speye(n1+1)) + D23 = kron3(speye(n3), d2, speye(n1+1)) + D31 = kron3(d3, speye(n2+1), speye(n1)) + D13 = kron3(speye(n3), speye(n2+1), d1) + D21 = kron3(speye(n3+1), d2, speye(n1)) + D12 = kron3(speye(n3+1), speye(n2), d1) + + O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1]) + O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1]) + O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1]) + + C = sp.vstack((sp.hstack((O1, -D32, D23)), + sp.hstack((D31, O2, -D13)), + sp.hstack((-D21, D12, O3))), format="csr") + + self._edgeCurl = sdiag(1/S)*(C*sdiag(L)) + return self._edgeCurl + return locals() + _edgeCurl = None + edgeCurl = property(**edgeCurl()) + + def faceAve(): + doc = "Construct the averaging operator on cell faces to cell centers." + + def fget(self): + if(self._faceAve is None): + n = self.n + if(self.dim == 1): + self._faceAve = av(n[0]) + elif(self.dim == 2): + self._faceAve = sp.hstack((sp.kron(speye(n[1]), av(n[0])), + sp.kron(av(n[1]), speye(n[0]))), format="csr") + elif(self.dim == 3): + self._faceAve = sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])), + kron3(speye(n[2]), av(n[1]), speye(n[0])), + kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr") + return self._faceAve + return locals() + _faceAve = None + faceAve = property(**faceAve()) + + def edgeAve(): + doc = "Construct the averaging operator on cell edges to cell centers." + + def fget(self): + if(self._edgeAve is None): + # The number of cell centers in each direction + n = self.n + if(self.dim == 1): + raise Exception('Edge Averaging does not make sense in 1D: Use Identity?') + elif(self.dim == 2): + self._edgeAve = sp.hstack((sp.kron(av(n[1]), speye(n[0])), + sp.kron(speye(n[1]), av(n[0]))), format="csr") + elif(self.dim == 3): + self._edgeAve = sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])), + kron3(av(n[2]), speye(n[1]), av(n[0])), + kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr") + return self._edgeAve + return locals() + _edgeAve = None + edgeAve = property(**edgeAve()) + + def getEdgeMass(self, materialProp=None): + """mass matix for products of edge functions w'*M(materialProp)*e""" + if(materialProp is None): + materialProp = np.ones(self.nC) + Av = self.edgeAve + return sdiag(Av.T * (self.vol * mkvc(materialProp))) + + def getFaceMass(self, materialProp=None): + """mass matix for products of edge functions w'*M(materialProp)*e""" + if(materialProp is None): + materialProp = np.ones(self.nC) + Av = self.faceAve + return sdiag(Av.T*(self.vol*mkvc(materialProp))) diff --git a/SimPEG/TensorMesh.py b/SimPEG/TensorMesh.py index 5f5e44b3..8395c54d 100644 --- a/SimPEG/TensorMesh.py +++ b/SimPEG/TensorMesh.py @@ -1,10 +1,11 @@ import numpy as np from BaseMesh import BaseMesh from TensorView import TensorView -from utils import ndgrid +from DiffOperators import DiffOperators +from utils import ndgrid, mkvc -class TensorMesh(BaseMesh, TensorView): +class TensorMesh(BaseMesh, TensorView, DiffOperators): """ TensorMesh is a mesh class that deals with tensor product meshes. @@ -21,15 +22,16 @@ class TensorMesh(BaseMesh, TensorView): """ def __init__(self, h, x0=None): - super(TensorMesh, self).__init__(np.array([len(x) for x in h]), x0) + super(TensorMesh, self).__init__(np.array([x.size for x in h]), x0) assert len(h) == len(self.x0), "Dimension mismatch. x0 != len(h)" for i, h_i in enumerate(h): assert type(h_i) == np.ndarray, ("h[%i] is not a numpy array." % i) + assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i) # Ensure h contains 1D vectors - self._h = [x.ravel() for x in h] + self._h = [mkvc(x) for x in h] def h(): doc = "h is a list containing the cell widths of the tensor mesh in each dimension." @@ -186,6 +188,79 @@ class TensorMesh(BaseMesh, TensorView): def getCellNumbering(self): pass + # --------------- Geometries --------------------- + def vol(): + doc = "Construct cell volumes of the 3D model as 1d array." + + def fget(self): + if(self._vol is None): + vh = self.h + # Compute cell volumes + if(self.dim == 1): + self._vol = mkvc(vh[0]) + elif(self.dim == 2): + # Cell sizes in each direction + self._vol = mkvc(np.outer(vh[0], vh[1])) + elif(self.dim == 3): + # Cell sizes in each direction + self._vol = mkvc(np.outer(mkvc(np.outer(vh[0], vh[1])), vh[2])) + return self._vol + return locals() + _vol = None + vol = property(**vol()) + + def area(): + doc = "Construct face areas of the 3D model as 1d array." + + def fget(self): + if(self._area is None): + # Ensure that we are working with column vectors + vh = self.h + # The number of cell centers in each direction + n = self.n + # Compute areas of cell faces + if(self.dim == 1): + self._area = np.ones(n[0]+1) + elif(self.dim == 2): + area1 = np.outer(np.ones(n[0]+1), vh[1]) + area2 = np.outer(vh[0], np.ones(n[1]+1)) + self._area = np.r_[mkvc(area1), mkvc(area2)] + elif(self.dim == 3): + area1 = np.outer(np.ones(n[0]+1), mkvc(np.outer(vh[1], vh[2]))) + area2 = np.outer(vh[0], mkvc(np.outer(np.ones(n[1]+1), vh[2]))) + area3 = np.outer(vh[0], mkvc(np.outer(vh[1], np.ones(n[2]+1)))) + self._area = np.r_[mkvc(area1), mkvc(area2), mkvc(area3)] + return self._area + return locals() + _area = None + area = property(**area()) + + def edge(): + doc = "Construct edge legnths of the 3D model as 1d array." + + def fget(self): + if(self._edge is None): + # Ensure that we are working with column vectors + vh = self.h + # The number of cell centers in each direction + n = self.n + # Compute edge lengths + if(self.dim == 1): + self._edge = mkvc(vh[0]) + elif(self.dim == 2): + l1 = np.outer(vh[0], np.ones(n[1]+1)) + l2 = np.outer(np.ones(n[0]+1), vh[1]) + self._edge = np.r_[mkvc(l1), mkvc(l2)] + elif(self.dim == 3): + l1 = np.outer(vh[0], mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1)))) + l2 = np.outer(np.ones(n[0]+1), mkvc(np.outer(vh[1], np.ones(n[2]+1)))) + l3 = np.outer(np.ones(n[0]+1), mkvc(np.outer(np.ones(n[1]+1), vh[2]))) + self._edge = np.r_[mkvc(l1), mkvc(l2), mkvc(l3)] + return self._edge + return locals() + _edge = None + edge = property(**edge()) + if __name__ == '__main__': print('Welcome to tensor mesh!') diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 2d1b1e99..bcbba681 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -1 +1,2 @@ from TensorMesh import TensorMesh +import utils diff --git a/SimPEG/getDIV.py b/SimPEG/getDIV.py deleted file mode 100644 index 6c45c57c..00000000 --- a/SimPEG/getDIV.py +++ /dev/null @@ -1,75 +0,0 @@ -import numpy as np -from scipy import sparse -from utils import mkvc -from sputils import ddx, sdiag, speye, kron3 - -def getvol(h): - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # Compute cell volumes - v12 = h1.T*h2 - V = mkvc(v12.reshape(-1,1)*h3) - - return V - -def getarea(h): - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - # Compute areas of cell faces - area1 = np.ones((n1+1,1))*mkvc(h2.T*h3) - area2 = h1.T*mkvc(np.ones((n2+1,1))*h3) - area3 = h1.T*mkvc(h2.T*np.ones(n3+1)) - area = np.hstack((np.hstack((mkvc(area1), mkvc(area2))), mkvc(area3))) - - return area - -def getDivMatrix(h): - """Consturct the 3D divergence operator on Faces.""" - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - - # Compute areas of cell faces - #area1 = np.ones((n1+1,1))*mkvc(h2.T*h3) - #area2 = h1.T*mkvc(np.ones((n2+1,1))*h3) - #area3 = h1.T*mkvc(h2.T*np.ones(n3+1)) - #area = np.hstack((np.hstack((mkvc(area1), mkvc(area2))), mkvc(area3))) - area = getarea(h) - - S = sdiag(area) - - # Compute cell volumes - #v12 = h1.T*h2 - #V = mkvc(v12.reshape(-1,1)*h3) - V = getvol(h) - - # Compute divergence operator on faces - d1 = ddx(n1) - d2 = ddx(n2) - d3 = ddx(n3) - D1 = kron3(speye(n3), speye(n2), d1) - D2 = kron3(speye(n3), d2, speye(n1)) - D3 = kron3(d3, speye(n2), speye(n1)) - - D = sparse.hstack((sparse.hstack((D1, D2)), D3)) - return sdiag(1/V)*D*S - diff --git a/SimPEG/getDiffOpps.py b/SimPEG/getDiffOpps.py deleted file mode 100644 index bcfe2c8a..00000000 --- a/SimPEG/getDiffOpps.py +++ /dev/null @@ -1,207 +0,0 @@ -import numpy as np -from scipy import sparse -from utils import mkvc -from sputils import * -#from sputils import ddx, sdiag, speye, kron3, spzeros, appendBottom3, - -def getVol(h): - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # Compute cell volumes - V = mkvc(np.outer(mkvc(np.outer(h1,h2)),h3)) - - return V - -def getArea(h): - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - # Compute areas of cell faces - area1 = mkvc(np.outer(np.ones(n1+1),np.outer(h2,h3))) - area2 = mkvc(np.outer(h1,mkvc(np.outer(np.ones(n2+1),h3)))) - area3 = mkvc(np.outer(h1,mkvc(np.outer(h2,np.ones(n3+1))))) - area = np.hstack((np.hstack((area1, area2)), area3)) - - return area - -def getLength(h): - - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - - # compute the length of each edge - Length1 = mkvc(np.outer(h1,mkvc(np.outer(np.ones(n2+1),np.ones(n3+1))))) - Length2 = mkvc(np.outer(np.ones(n1+1),mkvc(np.outer(h2,np.ones(n3+1))))) - Length3 = mkvc(np.outer(np.ones(n1+1),mkvc(np.outer(np.ones(n2+1),h3)))) - - Length = np.hstack((np.hstack((Length1, Length2)), Length3)) - - return Length - - -def getDivMatrix(h): - """Consturct the 3D divergence operator on Faces.""" - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - - area = getArea(h) - S = sdiag(area) - - # Compute cell volumes - V = getVol(h) - - # Compute divergence operator on faces - d1 = ddx(n1) - d2 = ddx(n2) - d3 = ddx(n3) - D1 = kron3(speye(n3), speye(n2), d1) - D2 = kron3(speye(n3), d2, speye(n1)) - D3 = kron3(d3, speye(n2), speye(n1)) - - D = sparse.hstack((sparse.hstack((D1, D2)), D3)) - return sdiag(1/V)*D*S - - -def getCurlMatrix(h): - """Edge CURL """ - - # Cell sizes in each direction - h1 = h[0]; h2 = h[1]; h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1); n2 = np.size(h2); n3 = np.size(h3) - - d1 = ddx(n1); d2 = ddx(n2); d3 = ddx(n3) - # derivatives on x-edge variables - D32 = kron3(d3, speye(n2), speye(n1+1)) - D23 = kron3(speye(n3), d2, speye(n1+1)) - D31 = kron3(d3, speye(n2+1), speye(n1)) - D13 = kron3(speye(n3), speye(n2+1), d1) - D21 = kron3(speye(n3+1), d2, speye(n1)) - D12 = kron3(speye(n3+1), speye(n2), d1) - - O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1]) - O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1]) - O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1]) - - CURL = appendBottom3( - appendRight3(O1, -D32, D23), - appendRight3(D31, O2, -D13), - appendRight3(-D21, D12, O3)) - - - area = getArea(h) - S = sdiag(1/area) - - # Compute edge length - lngth = getLength(h) - L = sdiag(lngth) - - return S*(CURL*L) - - -def getNodalGradient(h): - """Nodal Gradients""" - - # Cell sizes in each direction - h1 = h[0]; h2 = h[1]; h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1); n2 = np.size(h2); n3 = np.size(h3) - - - D1 = kron3(speye(n3+1), speye(n2+1), ddx(n1)) - D2 = kron3(speye(n3+1), ddx(n2), speye(n1+1)) - D3 = kron3(ddx(n3), speye(n2+1), speye(n1+1)) - - # topological gradient - GRAD = appendBottom3(D1, D2, D3) - - # scale for non-uniform mesh - # Compute edge length - lngth = getLength(h) - L = sdiag(1/lngth) - - return L*GRAD - - -def getEdgeToCellAverge(h): - - """Average from Edge to Cell center """ - - # Cell sizes in each direction - h1 = h[0]; h2 = h[1]; h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1); n2 = np.size(h2); n3 = np.size(h3) - - a1 = av(n1); a2 = av(n2); a3 = av(n3) - # derivatives on x-edge variables - A1 = kron3(a3, a2, speye(n1)) - A2 = kron3(a3, speye(n2), a1) - A3 = kron3(speye(n3), a2, a1) - - return appendRight3(A1, A2, A3) - - -def getFaceToCellAverge(h): - - """Average from Edge to Cell center """ - - # Cell sizes in each direction - h1 = h[0]; h2 = h[1]; h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1); n2 = np.size(h2); n3 = np.size(h3) - - a1 = av(n1); a2 = av(n2); a3 = av(n3) - # derivatives on x-edge variables - A1 = kron3(speye(n3), speye(n2), a1) - A2 = kron3(speye(n3), a2, speye(n1)) - A3 = kron3(a3, speye(n2), speye(n1)) - - return appendRight3(A1, A2, A3) - - -def getEdgeMassMatrix(h,sigma): - # mass matix for products of edge functions w'*M(sigma)*e - - Av = getEdgeToCellAverge(h) - v = getVol(h) - sigma = mkvc(sigma) - - return sdiag(Av.T*(v*sigma)) - -def getFaceMassMatrix(h,sigma): - # mass matix for products of edge functions w'*M(sigma)*e - - Av = getFaceToCellAverge(h) - v = getVol(h) - sigma = mkvc(sigma) - - return sdiag(Av.T*(v*sigma)) \ No newline at end of file diff --git a/SimPEG/sputils.py b/SimPEG/sputils.py index 7371e46b..b078fff2 100644 --- a/SimPEG/sputils.py +++ b/SimPEG/sputils.py @@ -1,35 +1,30 @@ import numpy as np -from scipy import sparse -from numpy import ones +from scipy import sparse as sp -def ddx(n): - """Define 1D derivatives""" - return sparse.spdiags((np.ones((n+1,1))*[-1,1]).T, [0,1], n, n+1) - def sdiag(h): """Sparse diagonal matrix""" - return sparse.spdiags(h, 0, np.size(h), np.size(h)) + return sp.spdiags(h, 0, np.size(h), np.size(h), format="csr") + def speye(n): """Sparse identity""" - return sparse.identity(n) + return sp.identity(n, format="csr") + def kron3(A, B, C): - """Two kron prods""" - return sparse.kron(sparse.kron(A, B), C) - -def av(n): - """Define 1D average""" - return 0.5*(sparse.spdiags(ones(n+1), 0, n, n+1) + sparse.spdiags(ones(n+1), 1, n, n+1)) - + """Three kron prods""" + return sp.kron(sp.kron(A, B), C, format="csr") + + def spzeros(n1, n2): """spzeros""" - return sparse.coo_matrix((n1, n2)) - + return sp.coo_matrix((n1, n2)).tocsr() + + def appendBottom(A, B): """append on bottom""" - C = sparse.vstack((A, B)) + C = sp.vstack((A, B)) C = C.tocsr() return C @@ -43,7 +38,7 @@ def appendBottom3(A, B, C): def appendRight(A, B): """append on right""" - C = sparse.hstack((A, B)) + C = sp.hstack((A, B)) C = C.tocsr() return C @@ -57,9 +52,9 @@ def appendRight3(A, B, C): def blkDiag(A, B): """blockdigonal""" - O12 = sparse.coo_matrix((np.shape(A)[0], np.shape(B)[1])) - O21 = sparse.coo_matrix((np.shape(B)[0], np.shape(A)[1])) - C = sparse.vstack((sparse.hstack((A, O12)), sparse.hstack((O21, B)))) + O12 = sp.coo_matrix((np.shape(A)[0], np.shape(B)[1])) + O21 = sp.coo_matrix((np.shape(B)[0], np.shape(A)[1])) + C = sp.vstack((sp.hstack((A, O12)), sp.hstack((O21, B)))) C = C.tocsr() return C @@ -69,8 +64,3 @@ def blkDiag3(A, B, C): ABC = blkDiag(blkDiag(A, B), C) ABC = ABC.tocsr() return ABC - - - - - \ No newline at end of file diff --git a/SimPEG/tests/OrderTest.py b/SimPEG/tests/OrderTest.py new file mode 100644 index 00000000..7f9f97ca --- /dev/null +++ b/SimPEG/tests/OrderTest.py @@ -0,0 +1,113 @@ +import sys +sys.path.append('../../') +from SimPEG import TensorMesh +import numpy as np +import unittest + + +class OrderTest(unittest.TestCase): + """ + + OrderTest is a base class for testing convergence orders with respect to mesh + sizes of integral/differential operators. + + Mathematical Problem: + + Given are an operator A and its discretization A[h]. For a given test function f + and h --> 0 we compare: + + error(h) = \| A[h](f) - A(f) \|_{\infty} + + Note that you can provide any norm. + + Test is passed when estimated rate order of convergence is at least 90% of the + estimated rate supplied by the user. + + Minimal example for a curl operator: + + class TestCURL(OrderTest): + name = "Curl" + + def getError(self): + # For given Mesh, generate A[h], f and A(f) and return norm of error. + + + fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x)) + sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y)) + + + Ex = fun(self.M.gridEx[:, 1]) + Ey = fun(self.M.gridEy[:, 2]) + Ez = fun(self.M.gridEz[:, 0]) + f = np.concatenate((Ex, Ey, Ez)) + + Fx = sol(self.M.gridFx[:, 2]) + Fy = sol(self.M.gridFy[:, 0]) + Fz = sol(self.M.gridFz[:, 1]) + Af = np.concatenate((Fx, Fy, Fz)) + + # Generate DIV matrix + Ah = self.M.edgeCurl + + curlE = Ah*E + err = np.linalg.norm((Ah*f -Af), np.inf) + return err + + def test_order(self): + # runs the test + self.orderTest() + + See also: test_operatorOrder.py + + """ + + name = "Order Test" + expectedOrder = 2 + meshSizes = [4, 8, 16, 32] + + def setupMesh(self, nc): + """ + For a given number of cells nc, generate a TensorMesh with uniform cells with edge length h=1/nc. + """ + h1 = np.ones(nc)/nc + h2 = np.ones(nc)/nc + h3 = np.ones(nc)/nc + h = [h1, h2, h3] + self.M = TensorMesh(h) + + def getError(self): + """For given h, generate A[h], f and A(f) and return norm of error.""" + return 1. + + def orderTest(self): + """ + For number of cells specified in meshSizes setup mesh, call getError + and prints mesh size, error, ratio between current and previous error, + and estimated order of convergence. + + + """ + order = [] + err_old = 0. + nc_old = 0. + for ii, nc in enumerate(self.meshSizes): + self.setupMesh(nc) + err = self.getError() + if ii == 0: + print '' + print 'Testing order of: ' + self.name + print '_____________________________________________' + print ' h | error | e(i-1)/e(i) | order' + print '~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~' + print '%4i | %8.2e |' % (nc, err) + else: + order.append(np.log(err/err_old)/np.log(float(nc_old)/float(nc))) + print '%4i | %8.2e | %6.4f | %6.4f' % (nc, err, err_old/err, order[-1]) + err_old = err + nc_old = nc + print '---------------------------------------------' + self.assertTrue(len(np.where(np.array(order) > 0.9*self.expectedOrder)[0]) > np.floor(0.75*len(order))) + + +if __name__ == '__main__': + unittest.main() diff --git a/SimPEG/tests/test_div.py b/SimPEG/tests/test_div.py deleted file mode 100644 index bc0c1a9d..00000000 --- a/SimPEG/tests/test_div.py +++ /dev/null @@ -1,45 +0,0 @@ -import numpy as np - -import sys -sys.path.append('../') -from TensorMesh import TensorMesh -from getDIV import getDivMatrix, getarea, getvol - -# Define the mesh -err=0. -for i in range(4): - icount=i+1; - nc = 2*icount; - h1 = np.pi/nc*np.ones((1,nc)) - h2 = np.pi/nc*np.ones((1,nc)) - h3 = np.pi/nc*np.ones((1,nc)) - h = [h1, h2, h3] - x0 = -np.pi/2*np.ones((3, 1)) - M = TensorMesh(h, x0) - #n = M.plotGrid() - - # Generate DIV matrix - DIV = getDivMatrix(h) - - #Test function - fun = lambda x: np.sin(x) - Fx = fun(M.gridFx[:,0]) - Fy = fun(M.gridFy[:,1]) - Fz = fun(M.gridFz[:,2]) - - F = np.concatenate((Fx,Fy,Fz)) - divF = DIV*F - sol = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) - divF_anal = sol(M.gridCC[:,0], M.gridCC[:,1], M.gridCC[:,2]) - - area = getarea(h) - vol = getvol(h) - err = np.linalg.norm((divF-divF_anal)*np.sqrt(vol), 2) - if icount == 1: - err1 = err - print 'h | 2 norm | error ratio' - print '---------------------------------------' - print '%6.4f | %8.2e |'% (h1[0,0], err) - else: - print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err1/err) - diff --git a/SimPEG/tests/test_tensorMesh.py b/SimPEG/tests/test_tensorMesh.py index e1b74374..c0b194db 100644 --- a/SimPEG/tests/test_tensorMesh.py +++ b/SimPEG/tests/test_tensorMesh.py @@ -3,15 +3,18 @@ import unittest import sys sys.path.append('../') from TensorMesh import TensorMesh +from OrderTest import OrderTest +from scipy.sparse.linalg import dsolve -class TestSequenceFunctions(unittest.TestCase): +class BasicTensorMeshTests(unittest.TestCase): def setUp(self): a = np.array([1, 1, 1]) b = np.array([1, 2]) - x0 = np.array([3, 5]) - self.mesh2 = TensorMesh([a, b], x0) + c = np.array([1, 4]) + self.mesh2 = TensorMesh([a, b], np.array([3, 5])) + self.mesh3 = TensorMesh([a, b, c]) def test_vectorN_2D(self): testNx = np.array([3, 4, 5, 6]) @@ -29,6 +32,140 @@ class TestSequenceFunctions(unittest.TestCase): ytest = np.all(self.mesh2.vectorCCy == testNy) self.assertTrue(xtest and ytest) + def test_area_3D(self): + test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2]) + t1 = np.all(self.mesh3.area == test_area) + self.assertTrue(t1) + def test_vol_3D(self): + test_vol = np.array([1, 1, 1, 2, 2, 2, 4, 4, 4, 8, 8, 8]) + t1 = np.all(self.mesh3.vol == test_vol) + self.assertTrue(t1) + + def test_vol_2D(self): + test_vol = np.array([1, 1, 1, 2, 2, 2]) + t1 = np.all(self.mesh2.vol == test_vol) + self.assertTrue(t1) + + def test_edge_3D(self): + test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]) + t1 = np.all(self.mesh3.edge == test_edge) + self.assertTrue(t1) + + def test_edge_2D(self): + test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2]) + t1 = np.all(self.mesh2.edge == test_edge) + self.assertTrue(t1) + +class TestCurl(OrderTest): + name = "Curl" + + def getError(self): + fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x)) + sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y)) + + Ex = fun(self.M.gridEx[:, 1]) + Ey = fun(self.M.gridEy[:, 2]) + Ez = fun(self.M.gridEz[:, 0]) + E = np.concatenate((Ex, Ey, Ez)) + + Fx = sol(self.M.gridFx[:, 2]) + Fy = sol(self.M.gridFy[:, 0]) + Fz = sol(self.M.gridFz[:, 1]) + curlE_anal = np.concatenate((Fx, Fy, Fz)) + + # Generate DIV matrix + CURL = self.M.edgeCurl + + curlE = CURL*E + err = np.linalg.norm((curlE-curlE_anal), np.inf) + return err + + def test_order(self): + self.orderTest() + + +class TestFaceDiv(OrderTest): + name = "Face Divergence" + + def getError(self): + DIV = self.M.faceDiv + + #Test function + fun = lambda x: np.sin(x) + Fx = fun(self.M.gridFx[:, 0]) + Fy = fun(self.M.gridFy[:, 1]) + Fz = fun(self.M.gridFz[:, 2]) + + F = np.concatenate((Fx, Fy, Fz)) + divF = DIV*F + sol = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) + divF_anal = sol(self.M.gridCC[:, 0], self.M.gridCC[:, 1], self.M.gridCC[:, 2]) + + err = np.linalg.norm((divF-divF_anal), np.inf) + + return err + + def test_order(self): + self.orderTest() + + +class TestNodalGrad(OrderTest): + name = "Nodal Gradient" + + def getError(self): + GRAD = self.M.nodalGrad + #Test function + fun = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) + sol = lambda x: -np.sin(x) # i (sin(x)) + j (sin(y)) + k (sin(z)) + + phi = fun(self.M.gridN[:, 0], self.M.gridN[:, 1], self.M.gridN[:, 2]) + gradE = GRAD*phi + + Ex = sol(self.M.gridEx[:, 0]) + Ey = sol(self.M.gridEy[:, 1]) + Ez = sol(self.M.gridEz[:, 2]) + + gradE_anal = np.concatenate((Ex, Ey, Ez)) + err = np.linalg.norm((gradE-gradE_anal), np.inf) + + return err + + def test_order(self): + self.orderTest() + + +class TestPoissonEqn(OrderTest): + name = "Poisson Equation" + meshSizes = [16, 20, 24] + + def getError(self): + # Create some functions to integrate + fun = lambda x: np.sin(2*np.pi*x[:, 0])*np.sin(2*np.pi*x[:, 1])*np.sin(2*np.pi*x[:, 2]) + sol = lambda x: -3.*((2*np.pi)**2)*fun(x) + + self.M.setCellGradBC('dirichlet') + + D = self.M.faceDiv + G = self.M.cellGrad + if self.forward: + sA = sol(self.M.gridCC) + sN = D*G*fun(self.M.gridCC) + err = np.linalg.norm((sA - sN), np.inf) + else: + fA = fun(self.M.gridCC) + fN = dsolve.spsolve(D*G, sol(self.M.gridCC)) + err = np.linalg.norm((fA - fN), np.inf) + return err + + def test_orderForward(self): + self.name = "Poisson Equation - Forward" + self.forward = True + self.orderTest() + + def test_orderBackward(self): + self.name = "Poisson Equation - Backward" + self.forward = False + self.orderTest() if __name__ == '__main__': unittest.main()