which leads to the final discretization
$$ \left(\frac{\partial}{\partial x} k \frac{\partial u} {\partial x}\right)_i \ \approx \frac{k_{i+\frac{1}{2}}\left(u_{i+1} - u_{i}\right) - k_{i-\frac{1}{2}}\left(u_{i} - u_{i-1}\right)}{h^2} + \mathcal{O}(h^2), \quad i = 1, \ldots, n-1, \quad u_0 = u_n = 0. $$In two dimensions, $k(x, y) \in \mathbb{R}^{2 \times 2}$ for every point $(x, y) \in \Omega$ is diffusion tensor.
If $k(x, y)$ is a diagonal matrix $$ k(x,y) = \begin{bmatrix} K_1(x, y) & 0\\ 0 & K_2(x, y) \end{bmatrix} $$ we have
$$ \mathrm{div}(k(x,y) \nabla u) = \frac{\partial}{\partial x} K_1 (x, y) \frac{\partial u}{\partial x} + \frac{\partial}{\partial y} K_2 (x, y) \frac{\partial u}{\partial y}, $$and we discretize each term and get block tridiagonal matrix with tridiagonal blocks.
For the simplest case $K_1 = K_2 = I$ we get 2D Poisson problem, and the matrix can be written
$$\Delta_{2D} = \Delta_{1D} \otimes I + I \otimes \Delta_{1D},$$where $\Delta_{1D} = \frac{1}{h^2}\mathrm{tridiag}(-1, 2, -1)$ is a one-dimensional Laplace operator.
This is a linear system $$ Ac = b $$ with unknown vector $c = [c_1, \dots,c_n]^{\top}$, matrix $$ A = [a_{ij}], \; a_{ij} = (\nabla \psi_i, k\nabla \psi_j) $$ and right-hand side $$ b = [b_i], \; b_i = (\psi_i, f) $$
Is matrix $A$ symmetric and positive-definite?
In [1]:
import numpy as np
import time
n = 4000
a = np.random.randn(n, n)
v = np.random.randn(n)
t = time.time()
np.dot(a, v)
t = time.time() - t
print('Time: {0: 3.1e}, Efficiency: {1: 3.1e} Gflops'.\
format(t, ((2 * n ** 2)/t) / 10 ** 9))
In [2]:
import scipy as sp
import scipy.sparse
n = 4000
r = 100
ex = np.ones(n);
a = sp.sparse.spdiags(np.vstack((ex, -2*ex, ex)), [-1, 0, 1], n, n, 'csr');
rhs = np.random.randn(n, r)
t = time.time()
a.dot(rhs)
t = time.time() - t
print('Time: {0: 3.1e}, Efficiency: {1: 3.1e} Gflops'.\
format(t, (3 * n * r) / t / 10 ** 9))
Next lecture considers methods to solve large sparse linear systems
In [1]:
%matplotlib inline
from __future__ import print_function
import fenics
import matplotlib.pyplot as plt
In [3]:
import dolfin
import mshr
import math
domain_vertices = [dolfin.Point(0.0, 0.0),
dolfin.Point(10.0, 0.0),
dolfin.Point(10.0, 2.0),
dolfin.Point(8.0, 2.0),
dolfin.Point(7.5, 1.0),
dolfin.Point(2.5, 1.0),
dolfin.Point(2.0, 4.0),
dolfin.Point(0.0, 4.0),
dolfin.Point(0.0, 0.0)]
p = mshr.Polygon(domain_vertices);
rect_mesh = mshr.generate_mesh(p, 20)
fenics.plot(rect_mesh)
Out[3]:
In [4]:
V = fenics.FunctionSpace(rect_mesh, 'P', 1)
u_D = fenics.Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
def boundary(x, on_boundary):
return on_boundary
bc = fenics.DirichletBC(V, u_D, boundary)
u = fenics.TrialFunction(V)
v = fenics.TestFunction(V)
f = fenics.Constant(-6.0) # Or f = Expression(’-6’, degree=0)
# Left-hand side
a = fenics.dot(fenics.grad(u), fenics.grad(v))*fenics.dx
# Right-hand side
L = f*v*fenics.dx
u = fenics.Function(V)
fenics.solve(a == L, u, bc)
fenics.plot(u)
Out[4]:
In [5]:
error_L2 = fenics.errornorm(u_D, u, 'L2')
print("Error in L2 norm = {}".format(error_L2))
error_H1 = fenics.errornorm(u_D, u, 'H1')
print("Error in H1 norm = {}".format(error_H1))
In [1]:
%matplotlib inline
import fipy
cellSize = 0.05
radius = 1.
mesh = fipy.Gmsh2D('''
cellSize = %(cellSize)g;
radius = %(radius)g;
Point(1) = {0, 0, 0, cellSize};
Point(2) = {-radius, 0, 0, cellSize};
Point(3) = {0, radius, 0, cellSize};
Point(4) = {radius, 0, 0, cellSize};
Point(5) = {0, -radius, 0, cellSize};
Circle(6) = {2, 1, 3};
Circle(7) = {3, 1, 4};
Circle(8) = {4, 1, 5};
Circle(9) = {5, 1, 2};
Line Loop(10) = {6, 7, 8, 9};
Plane Surface(11) = {10};
''' % locals())
In [2]:
phi = fipy.CellVariable(name = "solution variable",
mesh = mesh,
value = 0.)
In [3]:
# viewer = fipy.Viewer(vars=phi, datamin=-1, datamax=1.)
In [4]:
D = 1.
eq = fipy.TransientTerm() == fipy.DiffusionTerm(coeff=D)
In [5]:
X, Y = mesh.faceCenters
phi.constrain(X, mesh.exteriorFaces)
timeStepDuration = 10 * 0.9 * cellSize**2 / (2 * D)
steps = 10
for step in range(steps):
eq.solve(var=phi, dt=timeStepDuration)
# if viewer is not None:
# viewer
In [6]:
viewer
Out[6]: