/
diffusion-periodic-2d-PETSc.edp
64 lines (54 loc) · 2.33 KB
/
diffusion-periodic-2d-PETSc.edp
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
// run with MPI: ff-mpirun -np 4 script.edp
// NBPROC 4
load "PETSc" // PETSc plugin
macro dimension()2// EOM // 2D or 3D
include "macro_ddm.idp" // additional DDM functions
macro grad(u)[dx(u), dy(u)]// EOM // two-dimensional gradient
int[int] labPeriodic = [2, 4, 1, 3];
macro Pk() P2, periodic=[[labPeriodic[0],x+y], [labPeriodic[1],x+y], [labPeriodic[2],x-y], [labPeriodic[3],x-y]]// EOM
real r = 0.25;
border a(t=0,1) { x = -t+1; y = t; label = 1; };
border b(t=0,1) { x = -t; y = 1-t; label = 2; };
border c(t=0,1) { x = t-1; y = -t; label = 3; };
border d(t=0,1) { x = t; y = -1+t; label = 4; };
border e(t=0,2*pi) { x = r*cos(t); y = -r*sin(t); label = 0; };
mesh Th = buildmesh(a(getARGV("-global", 40)) + b(getARGV("-global", 40)) + c(getARGV("-global", 40)) + d(getARGV("-global", 40)) + e(getARGV("-global", 40)));
Mat A;
int[int] n2o;
macro ThPeriodicity()labPeriodic//
macro ThN2O()n2o//
DmeshCreate(Th);
MatCreate(Th, A, Pk);
fespace Wh(Th, Pk); // local finite element space
func f = (y+x+1) * (y+x-1) * (y-x+1) * (y-x-1);
varf vPb(u, v) = int2d(Th)(grad(u)' * grad(v)) - int2d(Th)((0.39 - f) * v) + on(0, u = 0.0);
real[int] rhs = vPb(0, Wh);
set(A, sparams = "-ksp_view");
Wh<real> u; // local solution
A = vPb(Wh, Wh);
u[] = A^-1 * rhs;
real[int] err = A * u[]; // global matrix-vector product
exchange(A, rhs, scaled = true);
err -= rhs;
macro def(u)u//
plotMPI(Th, u, Pk, def, real, cmm = "Global solution");
u[] = err;
plotMPI(Th, u, Pk, def, real, cmm = "Global residual");
set(A, sparams = "-pc_type gamg -ksp_type gmres -pc_gamg_threshold -1.0 -ksp_max_it 200");
u[] = 0.0;
u[] = A^-1 * rhs;
plotMPI(Th, u, Pk, def, real, cmm = "Global solution");
if(!NoGraphicWindow) {
real[int] sol;
ChangeNumbering(A, u[], sol);
ChangeNumbering(A, u[], sol, inverse = true);
mesh ThPlt = buildmesh(a(getARGV("-global", 40)) + b(getARGV("-global", 40)) + c(getARGV("-global", 40)) + d(getARGV("-global", 40)) + e(getARGV("-global", 40)));
fespace WhPlt(ThPlt, Pk);
WhPlt plt;
WhPlt reduce;
int[int] rest = restrict(Wh, WhPlt, n2o);
plt[](rest) = u[];
mpiReduce(plt[], reduce[], processor(0, mpiCommWorld), mpiSUM);
if(mpirank == 0)
plot(reduce);
}