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
|
'''
This module test the ksp class
'''
import pytest
pytest.importorskip("ngsolve")
from math import sqrt
from ngsolve import Mesh, H1, BilinearForm, LinearForm, grad, Integrate
from ngsolve import x, y, dx, GridFunction
from netgen.geom2d import unit_square
import netgen.meshing as ngm
from mpi4py.MPI import COMM_WORLD
from ngsPETSc import KrylovSolver
def test_ksp_preonly_lu():
'''
Testing the mapping PETSc KSP using a direct solver
'''
if COMM_WORLD.rank == 0:
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1).Distribute(COMM_WORLD))
else:
mesh = Mesh(ngm.Mesh.Receive(COMM_WORLD))
fes = H1(mesh, order=3, dirichlet="left|right|top|bottom")
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx)
solver = KrylovSolver(a, fes.FreeDofs(),
solverParameters={'ksp_type': 'preonly',
'ksp_monitor': '',
'pc_type': 'lu',})
f = LinearForm(fes)
f += 32 * (y*(1-y)+x*(1-x)) * v * dx
f.Assemble()
gfu = GridFunction(fes)
solver.solve(f.vec, gfu.vec)
exact = 16*x*(1-x)*y*(1-y)
assert sqrt(Integrate((gfu-exact)**2, mesh))<1e-4
def test_ksp_cg_gamg():
'''
Testing the mapping PETSc KSP using GAMG
'''
if COMM_WORLD.rank == 0:
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1).Distribute(COMM_WORLD))
else:
mesh = Mesh(ngm.Mesh.Receive(COMM_WORLD))
fes = H1(mesh, order=3, dirichlet="left|right|top|bottom")
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx).Assemble()
solver = KrylovSolver(a,fes.FreeDofs(),
solverParameters={'ksp_type': 'cg',
'ksp_monitor': '',
'pc_type': 'gamg'})
f = LinearForm(fes)
f += 32 * (y*(1-y)+x*(1-x)) * v * dx
f.Assemble()
gfu = GridFunction(fes)
solver.solve(f.vec, gfu.vec)
exact = 16*x*(1-x)*y*(1-y)
assert sqrt(Integrate((gfu-exact)**2, mesh))<1e-4
if __name__ == '__main__':
test_ksp_cg_gamg()
test_ksp_preonly_lu()
|
