From 4c82ce7dc27ba595321232e5c303dcc41d439f99 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 13 Nov 2013 15:14:33 -0800 Subject: [PATCH] Documentation for cellGradient. Boundary Conditions. Testing in 2D & 3D. --- SimPEG/mesh/DiffOperators.py | 86 +++++++++++++++++++++++- SimPEG/tests/test_operators.py | 117 ++++++++++++++++++++++++++++++++- 2 files changed, 199 insertions(+), 4 deletions(-) diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py index fa3f42e0..27fe3cc2 100644 --- a/SimPEG/mesh/DiffOperators.py +++ b/SimPEG/mesh/DiffOperators.py @@ -4,7 +4,13 @@ from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros, ddx, av def checkBC(bc): - """ Checks if boundary condition 'bc' is valid. """ + """ + + Checks if boundary condition 'bc' is valid. + + Each bc must be either 'dirichlet' or 'neumann' + + """ if(type(bc) is str): bc = [bc, bc] assert type(bc) is list, 'bc must be a list' @@ -17,7 +23,33 @@ def checkBC(bc): def ddxCellGrad(n, bc): - """Create 1D derivative operator from cell-centres to nodes this means we go from n to n+1""" + """ + Create 1D derivative operator from cell-centers to nodes this means we go from n to n+1 + + For Cell-Centered **Dirichlet**, use a ghost point:: + + (u_1 - u_g)/hf = grad + + u_g u_1 u_2 + * | * | * ... + ^ + 0 + + u_g = - u_1 + grad = 2*u1/dx + negitive on the other side. + + For Cell-Centered **Neumann**, use a ghost point:: + + (u_1 - u_g)/hf = 0 + + u_g u_1 u_2 + * | * | * ... + + u_g = u_1 + grad = 0; put a zero in. + + """ bc = checkBC(bc) D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n, format="csr") @@ -33,6 +65,56 @@ def ddxCellGrad(n, bc): D[-1, -1] = 0 return D +def ddxCellGradBC(n, bc): + """ + + Create 1D derivative operator from cell-centers to nodes this means we go from n to n+1 + + For Cell-Centered **Dirichlet**, use a ghost point:: + + (u_1 - u_g)/hf = grad + + u_g u_1 u_2 + * | * | * ... + ^ + u_b + + We know the value at the boundary (u_b):: + + (u_g+u_1)/2 = u_b (the average) + u_g = 2*u_b - u_1 + + So plug in to gradient: + + (u_1 - (2*u_b - u_1))/hf = grad + 2*(u_1-u_b)/hf = grad + + Separate, because BC are known (and can move to RHS later):: + + ( 2/hf )*u_1 + ( -2/hf )*u_b = grad + + ( ^ ) JUST RETURN THIS + + + """ + bc = checkBC(bc) + + ij = (np.array([0, n+1]),np.array([0, 1])) + vals = np.zeros(2) + + # Set the first side + if(bc[0] == 'dirichlet'): + vals[0] = -2 + elif(bc[0] == 'neumann'): + vals[0] = 0 + # Set the second side + if(bc[1] == 'dirichlet'): + vals[1] = 2 + elif(bc[1] == 'neumann'): + vals[1] = 0 + D = sp.csr_matrix((vals, ij), shape=(n+1,2)) + return D + class DiffOperators(object): """ diff --git a/SimPEG/tests/test_operators.py b/SimPEG/tests/test_operators.py index 26233626..d67a382b 100644 --- a/SimPEG/tests/test_operators.py +++ b/SimPEG/tests/test_operators.py @@ -1,6 +1,7 @@ import numpy as np import unittest from TestUtils import OrderTest +import matplotlib.pyplot as plt MESHTYPES = ['uniformTensorMesh', 'uniformLOM', 'rotateLOM'] call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1]) @@ -48,8 +49,120 @@ class TestCurl(OrderTest): self.orderTest() -class TestFaceDiv(OrderTest): - name = "Face Divergence" +class TestCellGrad2D_Dirichlet(OrderTest): + name = "Cell Grad 2D - Dirichlet" + meshTypes = ['uniformTensorMesh'] + meshDimension = 2 + meshSizes = [8, 16, 32, 64] + + def getError(self): + #Test function + fx = lambda x, y: 2*np.pi*np.cos(2*np.pi*x)*np.sin(2*np.pi*y) + fy = lambda x, y: 2*np.pi*np.cos(2*np.pi*y)*np.sin(2*np.pi*x) + sol = lambda x, y: np.sin(2*np.pi*x)*np.sin(2*np.pi*y) + + xc = call2(sol, self.M.gridCC) + + Fc = cartF2(self.M, fx, fy) + gradX_anal = self.M.projectFaceVector(Fc) + + self.M.setCellGradBC('dirichlet') + gradX = self.M.cellGrad.dot(xc) + + err = np.linalg.norm((gradX-gradX_anal), np.inf) + + return err + + def test_order(self): + self.orderTest() + + +class TestCellGrad3D_Dirichlet(OrderTest): + name = "Cell Grad 3D - Dirichlet" + meshTypes = ['uniformTensorMesh'] + meshDimension = 3 + meshSizes = [8, 16, 32, 64] + + def getError(self): + #Test function + fx = lambda x, y, z: 2*np.pi*np.cos(2*np.pi*x)*np.sin(2*np.pi*y)*np.sin(2*np.pi*z) + fy = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*x)*np.cos(2*np.pi*y)*np.sin(2*np.pi*z) + fz = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*x)*np.sin(2*np.pi*y)*np.cos(2*np.pi*z) + sol = lambda x, y, z: np.sin(2*np.pi*x)*np.sin(2*np.pi*y)*np.sin(2*np.pi*z) + + xc = call3(sol, self.M.gridCC) + + Fc = cartF3(self.M, fx, fy, fz) + gradX_anal = self.M.projectFaceVector(Fc) + + self.M.setCellGradBC('dirichlet') + gradX = self.M.cellGrad.dot(xc) + + err = np.linalg.norm((gradX-gradX_anal), np.inf) + + return err + + def test_order(self): + self.orderTest() + +class TestCellGrad2D_Neumann(OrderTest): + name = "Cell Grad 2D - Neumann" + meshTypes = ['uniformTensorMesh'] + meshDimension = 2 + meshSizes = [8, 16, 32, 64] + + def getError(self): + #Test function + fx = lambda x, y: -2*np.pi*np.sin(2*np.pi*x)*np.cos(2*np.pi*y) + fy = lambda x, y: -2*np.pi*np.sin(2*np.pi*y)*np.cos(2*np.pi*x) + sol = lambda x, y: np.cos(2*np.pi*x)*np.cos(2*np.pi*y) + + xc = call2(sol, self.M.gridCC) + + Fc = cartF2(self.M, fx, fy) + gradX_anal = self.M.projectFaceVector(Fc) + + self.M.setCellGradBC('neumann') + gradX = self.M.cellGrad.dot(xc) + + err = np.linalg.norm((gradX-gradX_anal), np.inf) + + return err + + def test_order(self): + self.orderTest() + + +class TestCellGrad3D_Neumann(OrderTest): + name = "Cell Grad 3D - Neumann" + meshTypes = ['uniformTensorMesh'] + meshDimension = 3 + meshSizes = [8, 16, 32, 64] + + def getError(self): + #Test function + fx = lambda x, y, z: -2*np.pi*np.sin(2*np.pi*x)*np.cos(2*np.pi*y)*np.cos(2*np.pi*z) + fy = lambda x, y, z: -2*np.pi*np.cos(2*np.pi*x)*np.sin(2*np.pi*y)*np.cos(2*np.pi*z) + fz = lambda x, y, z: -2*np.pi*np.cos(2*np.pi*x)*np.cos(2*np.pi*y)*np.sin(2*np.pi*z) + sol = lambda x, y, z: np.cos(2*np.pi*x)*np.cos(2*np.pi*y)*np.cos(2*np.pi*z) + + xc = call3(sol, self.M.gridCC) + + Fc = cartF3(self.M, fx, fy, fz) + gradX_anal = self.M.projectFaceVector(Fc) + + self.M.setCellGradBC('neumann') + gradX = self.M.cellGrad.dot(xc) + + err = np.linalg.norm((gradX-gradX_anal), np.inf) + + return err + + def test_order(self): + self.orderTest() + +class TestFaceDiv3D(OrderTest): + name = "Face Divergence 3D" meshTypes = MESHTYPES meshSizes = [8, 16, 32]