1//  run with MPI:  ff-mpirun -np 4 script.edp
2// NBPROC 4
3
4load "PETSc"                        // PETSc 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
13int s = getARGV("-split", 1);                                   // refinement factor
14
15meshN Th = square(1, 1);
16fespace Wh(Th, Pk);           // local finite element space
17int[int][int] intersection;   // local-to-neighbors renumbering
18real[int] D;                  // partition of unity
19{
20    int[int] l = [1, 2, 2, 2];
21    Th = square(getARGV("-global", 40), getARGV("-global", 40), label = l);    // global mesh
22    build(Th, s, intersection, D, Pk, mpiCommWorld)
23}
24
25real[int] rhs(Wh.ndof);                     // local right-hand side
26matrix<real> Loc;                           // local operator
27{                                           // local weak form
28    fespace Ph(Th, P0);
29    Ph kappa = x < 0.25 ? 10.0 : 1.0;
30    varf vPb(u, v) = intN(Th)(-1.0 * kappa * grad(u)' * grad(v)) + on(1, u = 0.0);
31    Loc = vPb(Wh, Wh, tgv = -2);
32    rhs = vPb(0, Wh, tgv = -2);
33}
34
35Mat T(Loc, intersection, D);
36func real[int] funcRes(real t, real[int]& in, real[int]& inT) {
37    real[int] out(in.n);
38    T = Loc;
39    MatMult(T, in, out);
40    out = inT - out;
41    return out;
42}
43real shift = 0;
44matrix Id;
45{
46    real[int] D(Loc.n);
47    D = 1.0;
48    Id = D;
49}
50func int funcJ(real t, real[int]& in, real[int]& inT, real a) {
51    matrix B = (-1.0) * Loc + a * Id;
52    T = B;
53    shift = a;
54    return 0;
55}
56Wh w;
57func int funcM(int s, real t, real[int]& u) {
58    changeNumbering(T, w[], u, exchange = true, inverse = true);
59    plotMPI(Th, w, Pk, def, real, cmm = "Global solution at step " + s + " (time " + t + ")")
60}
61w = (0.5 - x)^2 + (0.5 - y)^2 < 0.2 ? 1.0 : 0.0;
62real[int] wPETSc;
63changeNumbering(T, w[], wPETSc);
64TSSolve(T, funcJ, funcRes, wPETSc, sparams = "-ts_type beuler -ts_dt 0.1 -ts_max_time 100 -ts_exact_final_time interpolate -ts_max_snes_failures -1 -ts_view -pc_type lu -ts_adapt_type basic -ts_rtol 1e-3", monitor = funcM);
65