-
Notifications
You must be signed in to change notification settings - Fork 17
/
test_poisson_ball.py
82 lines (62 loc) · 1.99 KB
/
test_poisson_ball.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
# -*- coding: utf-8 -*-
import dolfin
import mshr
import helpers
import pyamg
import pyfvm
from pyfvm.form_language import integrate, n_dot_grad, dS, dV
from sympy import pi, sin, cos
import unittest
import voropy
def exact_sol(x):
return cos(pi/2 * (x[0]**2 + x[1]**2 + x[2]**2))
class Poisson(object):
def apply(self, u):
def rhs(x):
z = pi/2 * (x[0]**2 + x[1]**2 + x[2]**2)
return 2*pi * (1.5 * sin(z) + z * cos(z))
return integrate(lambda x: -n_dot_grad(u(x)), dS) - \
integrate(rhs, dV)
def dirichlet(self, u):
return [
(lambda x: u(x) - exact_sol(x), 'boundary')
]
def get_mesh(k):
h = 0.5**(k+2)
c = mshr.Sphere(dolfin.Point(0., 0., 0.), 1.0, int(2*pi / h))
m = mshr.generate_mesh(c, 2.0 / h)
return voropy.mesh_tetra.MeshTetra(
m.coordinates(),
m.cells(),
mode='geometric'
)
class ConvergencePoisson3dBallTest(unittest.TestCase):
def setUp(self):
return
@staticmethod
def solve(verbose=False):
def solver(mesh):
matrix, rhs = pyfvm.discretize_linear(Poisson(), mesh)
ml = pyamg.smoothed_aggregation_solver(matrix)
u = ml.solve(rhs, tol=1e-10)
return u
return helpers.perform_convergence_tests(
solver,
exact_sol,
get_mesh,
range(3),
verbose=verbose
)
def test(self):
H, error_norm_1, error_norm_inf, order_1, order_inf = self.solve()
expected_order = 2
tol = 2.0e-1
self.assertGreater(order_1[-1], expected_order - tol)
self.assertGreater(order_inf[-1], expected_order - tol)
return
if __name__ == '__main__':
from matplotlib import pyplot as plt
H, error_norm_1, error_norm_inf, order_1, order_inf = \
ConvergencePoisson3dBallTest.solve(verbose=True)
helpers.plot_error_data(H, error_norm_1, error_norm_inf)
plt.show()