wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-14 11:18:18 +08:00
Code Issues Packages Projects Releases Wiki Activity

Documentation for cellGradient. Boundary Conditions. Testing in 2D & 3D.

Browse Source
This commit is contained in:
Rowan Cockett
2013-11-13 15:14:33 -08:00
parent 56ee22d1ce
commit 4c82ce7dc2
2 changed files with 199 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/mesh/DiffOperators.py
+84 -2
View File
@@ -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):
"""
SimPEG/tests/test_operators.py
+115 -2
View File
@@ -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]