1//  run with MPI:  ff-mpirun -np 4 script.edp
2// NBPROC 4
3
4load "hpddm"
5macro dimension()2// EOM            // 2D or 3D
6include "macro_ddm.idp"             // additional DDM functions
7
8macro def(i)i// EOM                 // scalar field definition
9macro init(i)i// EOM                // scalar field initialization
10macro grad(u)[dx(u), dy(u)]// EOM   // two-dimensional gradient
11func Pk = P2;                       // finite element space
12
13real k = getARGV("-waven", 40.0);
14func f = 100 * exp(-10 * k * ((x-0.5)^2 + (y-0.5)^2));
15
16real lambda = 2 * pi / k;
17
18real epsilonA = 0;
19real epsilonE = 0;
20
21int Dirichlet = 1;
22int Robin = 2;
23
24int s = getARGV("-split", 3);
25real nw = 15.0 / s;
26meshN ThCoarse = square(nw * 1/lambda, nw * 1/lambda);
27
28{
29    int[int] chlab = [1, Dirichlet, 2, Dirichlet, 3, Dirichlet, 4, Dirichlet];
30    ThCoarse = change(ThCoarse, refe = chlab);
31}
32
33meshN Th;
34fespace Wh(Th, Pk);
35real[int] part;
36fespace PhCoarse(ThCoarse, P0);
37PhCoarse partCoarse;
38if(mpirank == 0)
39    partitionerSeq(partCoarse[], ThCoarse, mpisize);
40partitionerPar(partCoarse[], ThCoarse, mpiCommWorld, mpisize);
41int[int][int] intersection;   // local-to-neighbors renumbering
42real[int] D;                  // partition of unity
43{
44    Th = trunc(ThCoarse, 1, split = s);
45    fespace Ph(Th, P0);
46    Ph part;
47    part = partCoarse;
48    buildWithPartitioning(Th, part[], 1, intersection, D, Pk, mpiCommWorld);
49}
50
51real epsilon = epsilonA;
52varf vPb(u, v) = intN(Th)(-(k^2 - 1i*epsilon)*u*v + grad(u)'*grad(v))
53               + intN1(Th, Robin)(1i*k*u*v)
54               - intN(Th)(f*v)
55               + on(Dirichlet, u = 0);
56matrix<complex> Mat = vPb(Wh, Wh, sym = 1);
57fespace WhCoarse(ThCoarse, Pk);
58matrix<complex> MatCoarse;
59{
60    meshN ThBackup = Th;
61    Th = ThCoarse;
62    epsilon = epsilonE;
63    MatCoarse = vPb(Wh, Wh, sym = 1);
64    epsilon = epsilonA;
65    Th = ThBackup;
66}
67set(MatCoarse, solver = sparsesolverSym, sym = 1);
68schwarz<complex> ACoarseSeq(MatCoarse);
69set(ACoarseSeq, prefix = "level_2_");
70matrix R = interpolate(Wh, WhCoarse);
71complex[int] rhs = vPb(0, Wh);
72
73schwarz<complex> A(Mat, intersection, D);
74set(A, sparams = "-hpddm_verbosity " + (!mpirank ? "1" : "0") + " -hpddm_schwarz_method ras -hpddm_variant right -hpddm_gmres_restart 200 -hpddm_max_it 200");
75
76Wh<complex> u;
77u[] = A^-1 * rhs;
78plotMPI(Th, real(u), Pk, def, real, cmm = "Global solution with a one-level method")
79
80set(A, sparams = "-hpddm_verbosity " + (!mpirank ? "3" : "0") + " -hpddm_schwarz_coarse_correction deflated -hpddm_variant flexible -hpddm_level_2_verbosity 0");
81func complex[int] correctionExact(complex[int]& in) {
82    complex[int] out(in.n), tmp(WhCoarse.ndof), tmpReduced(WhCoarse.ndof);
83    for[i, di: D] in[i] *= di;
84    tmp = R' * in;
85    mpiAllReduce(tmp, tmpReduced, mpiCommWorld, mpiSUM);
86    tmp = ACoarseSeq^-1 * tmpReduced;
87    out = R * tmp;
88    return out;
89}
90if(mpisize > 1 && isSetOption("schwarz_coarse_correction")) // two-level Schwarz methods
91    attachCoarseOperator(mpiCommWorld, A, correctionExact);
92
93u[] = 0;
94u[] = A^-1 * rhs;
95plotMPI(Th, real(u), Pk, def, real, cmm = "Global solution with a redundant coarse correction")
96
97int[int][int] intersectionCoarse;
98real[int] DCoarse;
99{
100    buildWithPartitioning(ThCoarse, partCoarse[], 1, intersectionCoarse, DCoarse, Pk, mpiCommWorld);
101    meshN ThBackup = Th;
102    Th = ThCoarse;
103    epsilon = epsilonE;
104    MatCoarse = vPb(Wh, Wh, sym = 1);
105    epsilon = epsilonA;
106    Th = ThBackup;
107}
108R = interpolate(Wh, WhCoarse);
109schwarz<complex> ACoarse(MatCoarse, intersectionCoarse, DCoarse);
110set(ACoarse, sparams = "-hpddm_level_2_verbosity " + (!mpirank ? "1" : "0") + " -hpddm_level_2_schwarz_method ras -hpddm_level_2_tol 1e-1 -hpddm_level_2_gmres_restart 200 -hpddm_level_2_max_it 200 -hpddm_level_2_krylov_method gcrodr -hpddm_level_2_recycle 10", prefix = "level_2_");
111macro mplot()wait = 1, cmm = "Transfers on the correction at iteration " + ijk//
112int ijk = 1;
113func complex[int] correctionInexact(complex[int]& in) {
114    complex[int] out(in.n), tmp(WhCoarse.ndof), tmpReduced(WhCoarse.ndof);
115    if(!NoGraphicWindow && ijk > 0 && ijk < 10) {
116        Wh outW;
117        for[i, xi: in] outW[][i] = real(xi);
118        plotMPI(Th, outW, Pk, def, real, mplot)
119    }
120    tmpReduced = R' * in;
121    exchange(ACoarse, tmpReduced, scaled = true);
122    if(!NoGraphicWindow && ijk > 0 && ijk < 10) {
123        WhCoarse outWCoarse;
124        for[i, xi: tmpReduced] outWCoarse[][i] = real(tmpReduced[i]);
125        plotMPI(ThCoarse, outWCoarse, Pk, def, real, mplot)
126        ++ijk;
127    }
128    tmp = ACoarse^-1 * tmpReduced;
129    out = R * tmp;
130    exchange(A, out);
131    return out;
132}
133if(mpisize > 1 && isSetOption("schwarz_coarse_correction")) // two-level Schwarz methods
134    attachCoarseOperator(mpiCommWorld, A, correctionInexact);
135else
136    exit(0);
137
138u[] = 0;
139u[] = A^-1 * rhs;
140plotMPI(Th, real(u), Pk, def, real, cmm = "Global solution with a one-level inner coarse correction (two levels in total)")
141