//  run with MPI:  ff-mpirun -np 4 script.edp
// NBPROC 4

load "PETSc"                        // PETSc plugin
include "macro_ddm.idp"             // additional DDM functions

macro def(i)i// EOM                         // scalar field definition
macro init(i)i// EOM                        // scalar field initialization
macro grad(u)[dx(u), dy(u)]// EOM           // two-dimensional gradient
func Pk = P1;                               // finite element space

macro ThRefinementFactor()getARGV("-split", 1)//
int[int] l = [1, 2, 2, 2];
mesh Th = square(getARGV("-global", 40), getARGV("-global", 40), label = l);
fespace Wh(Th, Pk);           // local finite element space
Mat T;
MatCreate(Th, T, Pk);

real[int] rhs(Wh.ndof);                     // local right-hand side
matrix<real> Loc;                           // local operator
{                                           // local weak form
    fespace Ph(Th, P0);
    Ph kappa = x < 0.25 ? 10.0 : 1.0;
    varf vPb(u, v) = int2d(Th)(-1.0 * kappa * grad(u)' * grad(v)) + on(1, u = 0.0);
    Loc = vPb(Wh, Wh, tgv = -2);
    rhs = vPb(0, Wh, tgv = -2);
}

func real[int] funcRes(real t, real[int]& in, real[int]& inT) {
    real[int] out(in.n);
    T = Loc;
    MatMult(T, in, out);
    out = inT - out;
    return out;
}
real shift = 0;
matrix Id;
{
    real[int] D(Loc.n);
    D = 1.0;
    Id = D;
}
func int funcJ(real t, real[int]& in, real[int]& inT, real a) {
    matrix B = (-1.0) * Loc + a * Id;
    T = B;
    shift = a;
    return 0;
}
Wh w;
func int funcM(int s, real t, real[int]& u) {
    ChangeNumbering(T, w[], u, exchange = true, inverse = true);
    plotMPI(Th, w, Pk, def, real, cmm = "Global solution at step " + s + " (time " + t + ")");
}
w = (0.5 - x)^2 + (0.5 - y)^2 < 0.2 ? 1.0 : 0.0;
real[int] wPETSc;
ChangeNumbering(T, w[], wPETSc);
TSSolve(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);
