1//  run with MPI:  ff-mpirun -np 4 script.edp
2// NBPROC 4
3
4load "hpddm"                        // HPDDM plugin
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 = P1;                       // finite element space
12
13string deflation = getARGV("-deflation", "geneo");              // coarse space construction
14int overlap = getARGV("-overlap", 1);                           // geometric overlap between subdomains
15int fakeInterface = getARGV("-interface", 10);                  // interface between subdomains
16int s = getARGV("-split", 1);                                   // refinement factor
17
18mpiComm comm;
19int p = getARGV("-hpddm_level_2_p", 1);
20bool excluded = splitComm(mpiCommWorld, p, comm, topology = getARGV("-hpddm_level_2_topology", 0), exclude = (usedARGV("-hpddm_level_2_exclude") != -1));
21
22if(verbosity > 0 && mpirank == 0) {
23    cout << " --- " << mpirank << "/" << mpisize;
24    cout << " - diffusion-2d.edp - input parameters: refinement factor = " << s << " - overlap = " << overlap << endl;
25}
26
27meshN Th = square(1, 1);
28fespace Wh(Th, Pk);              // local finite element space
29int[int][int] intersection(0);   // local-to-neighbors renumbering
30real[int] D;                     // partition of unity
31{
32    meshN ThBorder;
33    Th = square(getARGV("-global", 40), getARGV("-global", 40));    // global mesh
34    buildOverlap(Th, ThBorder, fakeInterface, s, overlap, intersection, D, Pk, comm, excluded)
35}
36
37matrix<real> Mat;                           // local operator
38varf vPb(u, v) = intN(Th)(grad(u)' * grad(v)) + intN(Th)(v) + on(1, u = 1.0);
39Mat = vPb(Wh, Wh, tgv = -1);
40real[int] rhs = vPb(0, Wh, tgv = -1);
41
42schwarz A(Mat, intersection, D);
43set(A, sparams = "-hpddm_schwarz_method ras -hpddm_schwarz_coarse_correction deflated -hpddm_geneo_nu 10");
44
45matrix<real> Opt;                           // local operator with optimized boundary conditions
46pair ret;
47{
48    int solver = getOption("schwarz_method");
49    if(solver == 1 || solver == 2 || solver == 4) { // optimized Schwarz methods
50        fespace Ph(Th, P0);
51        real kZero = getARGV("-kZero", 10.0);
52        Ph transmission = kZero;
53        varf vOptimized(u, v) = intN(Th)(grad(u)' * grad(v)) + intN1(Th, fakeInterface)(transmission * (u * v)) + on(1, u = 1.0);
54        Opt = vOptimized(Wh, Wh, tgv = -1);
55    }
56    if(mpisize > 1 && isSetOption("schwarz_coarse_correction")) { // two-level Schwarz methods
57        if(excluded)
58            attachCoarseOperator(mpiCommWorld, A/*, A = noPen, B = overlapRestriction, threshold = 2. * h[].max / diam*/);
59        else {
60            varf vPbNoPen(u, v) = intN(Th)(grad(u)' * grad(v)) + on(1, u = 0.0);
61            matrix<real> noPen = vPbNoPen(Wh, Wh, sym = 1);
62            if(deflation == "geneo") // standard GenEO, no need for RHS -> deduced from LHS (Neumann matrix)
63                attachCoarseOperator(mpiCommWorld, A, A = noPen/*, threshold = 2. * h[].max / diam,*/, ret = ret);
64            else if(deflation == "dtn") {
65                varf vMass(def(u), def(v)) = intN1(Th, fakeInterface)(u * v);
66                matrix<real> massMatrix = vMass(Wh, Wh, sym = 1);
67                attachCoarseOperator(mpiCommWorld, A, A = noPen, B = massMatrix, pattern = Opt/*, threshold = k,*/, ret = ret);
68            }
69            else if(deflation == "geneo-2") // GenEO-2 for optimized Schwarz methods, need for RHS (LHS is still Neumann matrix)
70                attachCoarseOperator(mpiCommWorld, A, A = noPen, B = Opt, pattern = Opt/*, threshold = 2. * h[].max / diam,*/, ret = ret);
71        }
72    }
73}
74
75Wh<real> def(u);    // local solution
76
77if(Opt.n > 0)       // optimized Schwarz methods
78    DDM(A, rhs, u[], excluded = excluded, ret = ret, O = Opt);
79else
80    u[] = A^-1 * rhs;
81
82real[int] err = A * u[]; // global matrix-vector product
83err -= rhs;
84
85plotMPI(Th, u, Pk, def, real, cmm = "Global solution")
86u[] = err;
87plotMPI(Th, u, Pk, def, real, cmm = "Global residual")
88