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

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

macro grad(u)[dx(u), dy(u)]// EOM   // two-dimensional gradient
func Pk = P1;                       // finite element space

mesh Th = square(getARGV("-global", 40), getARGV("-global", 40)); // global mesh
Mat A;
macro ThRefinementFactor()getARGV("-split", 1)//
MatCreate(Th, A, Pk);

fespace Wh(Th, Pk);                 // local finite element space
varf vPb(u, v) = int2d(Th)(grad(u)' * grad(v));
func f = 1 + x - y;
varf vL(u, v)  = int2d(Th)(f * v);
varf vb(u, v)  = int2d(Th)(v);
real[int] F = vL(0, Wh);
real[int] B = vb(0, Wh);
real[int] pB, pF, pSol;
A = vPb(Wh, Wh);
ChangeNumbering(A, B, pB);
ChangeNumbering(A, F, pF);

Mat N = [[ A  , pB ],
         [ pB', 0  ]];
if(mpirank == 0) {
    pF.resize(pF.n + 1);
    pF = [pF(0:pF.n - 2), 0];
    F.resize(F.n + 1);
    F = [F(0:F.n - 2), 0];
}
{
    set(N, sparams = "-pc_type lu");
    real[int] x = N^-1 * F;
    if(mpirank == 0)
        x.resize(x.n - 1);
    Wh sol;
    sol[] = x;
    macro params()cmm = "Global solution #1", wait = 1, dim = 3// EOM
    macro def(u)u//
    plotMPI(Th, sol, Pk, def, real, params);
}
{
    Mat C;
    MatConvert(N, C);
    set(C, sparams = "-pc_type lu");
    KSPSolve(C, pF, pSol);
    if(mpirank == 0)
        pSol.resize(pSol.n - 1);
    Wh sol;
    ChangeNumbering(A, sol[], pSol, inverse = true, exchange = true);
    macro params()cmm = "Global solution #2", wait = 1, dim = 3// EOM
    macro def(u)u//
    plotMPI(Th, sol, Pk, def, real, params);
}
{
    set(N, sparams = "-ksp_view_pre -pc_type fieldsplit -fieldsplit_0_pc_type lu -fieldsplit_0_mat_mumps_cntl_3 1e-12 -fieldsplit_0_mat_mumps_icntl_24 1 -pc_fieldsplit_detect_saddle_point");
    KSPSolve(N, pF, pSol);
    if(mpirank == 0)
        pSol.resize(pSol.n - 1);
    Wh sol;
    ChangeNumbering(A, sol[], pSol, inverse = true, exchange = true);
    macro params()cmm = "Global solution #3", wait = 1, dim = 3// EOM
    macro def(u)u//
    plotMPI(Th, sol, Pk, def, real, params);
}
{
    Mat M = [[ 0 , pB'],
             [ pB, A  ]];
    set(M, sparams = "-pc_type lu");
    if(mpirank == 0) {
        real[int] tmp = F(0:F.n - 2);
        F(0) = 0;
        F(1:F.n - 1) = tmp;
    }
    real[int] x = M^-1 * F;
    if(mpirank == 0) {
        real[int] tmp = x(1:F.n - 1);
        x.resize(x.n - 1);
        x = tmp;
    }
    Wh sol;
    sol[] = x;
    macro params()cmm = "Global solution #4", wait = 1, dim = 3// EOM
    macro def(u)u//
    plotMPI(Th, sol, Pk, def, real, params);
}
