mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-20 12:40:44 +08:00
added comments
This commit is contained in:
Lars Ruthotto
+13
-9
@@ -10,6 +10,7 @@ def ddx(n):
|
||||
|
||||
|
||||
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'
|
||||
@@ -22,7 +23,7 @@ def checkBC(bc):
|
||||
|
||||
|
||||
def ddxCellGrad(n, bc):
|
||||
"""Define 1D derivatives, outer, this means we go from n to n+1"""
|
||||
"""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")
|
||||
@@ -40,7 +41,7 @@ def ddxCellGrad(n, bc):
|
||||
|
||||
|
||||
def av(n):
|
||||
"""Define 1D averaging operator"""
|
||||
"""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")
|
||||
|
||||
|
||||
@@ -49,10 +50,10 @@ class DiffOperators(object):
|
||||
Class creates the differential operators that you need!
|
||||
"""
|
||||
def __init__(self):
|
||||
raise Exception('You should use a Mesh class.')
|
||||
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 the 3D Divergence operator on Faces."
|
||||
doc = "Construct divergence operator (face-stg to cell-centres)."
|
||||
|
||||
def fget(self):
|
||||
if(self._faceDiv is None):
|
||||
@@ -81,7 +82,7 @@ class DiffOperators(object):
|
||||
faceDiv = property(**faceDiv())
|
||||
|
||||
def nodalGrad():
|
||||
doc = "Construct the 3D nodal gradient operator."
|
||||
doc = "Construct gradient operator (nodes to edges)."
|
||||
|
||||
def fget(self):
|
||||
if(self._nodalGrad is None):
|
||||
@@ -109,11 +110,14 @@ class DiffOperators(object):
|
||||
|
||||
def setCellGradBC(self, BC):
|
||||
"""
|
||||
e.g.
|
||||
Function that sets the boundary conditions for cell-centred derivative operators.
|
||||
|
||||
BC = 'neumann'
|
||||
BC = ['neumann', 'dirichlet', 'neumann']
|
||||
BC = [['neumann', 'dirichlet'], 'dirichlet', 'neumann']
|
||||
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):
|
||||
|
||||
@@ -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))
|
||||
+1
-1
@@ -13,7 +13,7 @@ def speye(n):
|
||||
|
||||
|
||||
def kron3(A, B, C):
|
||||
"""Two kron prods"""
|
||||
"""Three kron prods"""
|
||||
return sp.kron(sp.kron(A, B), C, format="csr")
|
||||
|
||||
|
||||
|
||||
@@ -6,14 +6,69 @@ import unittest
|
||||
|
||||
|
||||
class OrderTest(unittest.TestCase):
|
||||
"""Order test sets up the basics for testing order of decrease for a function on a mesh."""
|
||||
"""
|
||||
|
||||
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):
|
||||
# Define the mesh
|
||||
"""
|
||||
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
|
||||
@@ -21,10 +76,17 @@ class OrderTest(unittest.TestCase):
|
||||
self.M = TensorMesh(h)
|
||||
|
||||
def getError(self):
|
||||
"""Overwrite this function with the guts of the test."""
|
||||
"""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.
|
||||
@@ -34,16 +96,16 @@ class OrderTest(unittest.TestCase):
|
||||
if ii == 0:
|
||||
print ''
|
||||
print 'Testing order of: ' + self.name
|
||||
print '__________________________________________'
|
||||
print ' h | inf norm | ratio | order'
|
||||
print '~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~|~~~~~~~~~~'
|
||||
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])
|
||||
print '%4i | %8.2e | %6.4f | %6.4f' % (nc, err, err_old/err, order[-1])
|
||||
err_old = err
|
||||
nc_old = nc
|
||||
print '------------------------------------------'
|
||||
print '---------------------------------------------'
|
||||
self.assertTrue(len(np.where(np.array(order) > 0.9*self.expectedOrder)[0]) > np.floor(0.75*len(order)))
|
||||
|
||||
|
||||
|
||||
@@ -1,119 +0,0 @@
|
||||
import numpy as np
|
||||
from OrderTest import OrderTest
|
||||
import unittest
|
||||
from scipy.sparse.linalg import dsolve
|
||||
|
||||
|
||||
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()
|
||||
@@ -3,9 +3,11 @@ 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])
|
||||
@@ -55,6 +57,115 @@ class TestSequenceFunctions(unittest.TestCase):
|
||||
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()
|
||||
|
||||
Reference in New Issue
Block a user