1//  run with MPI:  ff-mpirun -np 4 script.edp
2// NBPROC 4
3
4load "PETSc"
5macro dimension()2// EOM
6include "macro_ddm.idp"
7
8real c = 6.25;
9int N  = 80;
10
11// domain: unit square
12border aa(t=0,1) { x=t;   y=0;   };
13border bb(t=0,1) { x=1;   y=t;   };
14border cc(t=0,1) { x=1-t; y=1;   };
15border dd(t=0,1) { x=0;   y=1-t; };
16
17mesh M = buildmesh(aa(N)+bb(N)+cc(N)+dd(N));
18
19load "Element_P3"
20func Pk = P1;
21Mat A;
22buildMat(M, getARGV("-split", 1), A, Pk, mpiCommWorld)
23fespace Vh(M, Pk);
24Vh u;
25func BC = cos(pi*x)*cos(pi*y);
26varf vInit(w, v) = on(aa, bb, cc, dd, w = BC);
27varf vJ(w, v) = int2d(M)(dx(w)*dx(v) + dy(w)*dy(v) - c*exp(u)*w*v) + on(aa, bb, cc, dd, w = 0);
28varf vRes(w, v) = int2d(M)(dx(u)*dx(v) + dy(u)*dy(v) - c*exp(u)*v) + on(aa, bb, cc, dd, w = u);
29func real[int] funcRes(real[int]& inPETSc) {
30    changeNumbering(A, u[], inPETSc, inverse = true, exchange = true);
31    real[int] out(Vh.ndof);
32    out = vRes(0, Vh, tgv = -2);
33    real[int] outPETSc;
34    changeNumbering(A, out, outPETSc);
35    return outPETSc;
36}
37func int funcJ(real[int]& inPETSc) {
38    changeNumbering(A, u[], inPETSc, inverse = true, exchange = true);
39    A = vJ(Vh, Vh, tgv = -2);
40    return 0;
41}
42if(hasType("PC", "hpddm")) {
43    bool stage = usedARGV("-stage") != -1;
44    for(int i = 0; i < 3; ++i) {
45        string params;
46        if(i == 0) {
47            params = "";
48            if(stage)
49                PetscLogStagePush("Everything");
50        }
51        else if(i == 1) {
52            params = "-pc_hpddm_coarse_ksp_reuse_preconditioner";
53            if(stage)
54                PetscLogStagePush("Fine level");
55        }
56        else if(i == 2) {
57            params = "-pc_hpddm_levels_1_ksp_reuse_preconditioner";
58            if(stage)
59                PetscLogStagePush("Nothing");
60        }
61        set(A, sparams = "-pc_type hpddm -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_levels_1_eps_nev 10 -pc_hpddm_define_subdomains -pc_hpddm_has_neumann -ksp_converged_reason -pc_hpddm_levels_1_pc_type asm " + params);
62        real[int] bPETSc;
63        u[] = vInit(0, Vh, tgv = -2);
64        changeNumbering(A, u[], bPETSc);
65        real[int] xPETSc = bPETSc;
66        SNESSolve(A, funcJ, funcRes, bPETSc, xPETSc, sparams = "-snes_monitor -snes_type newtonls -snes_converged_reason -ksp_converged_reason -snes_view");
67        changeNumbering(A, u[], xPETSc, inverse = true, exchange = true);
68        macro def(u)u//
69        plotMPI(M, u, Pk, def, real, cmm = "Global solution")
70        set(A, sparams = "-pc_type none");
71        if(stage)
72            PetscLogStagePop();
73    }
74}
75