//  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;
int[int] n2o;
macro ThN2O()n2o//
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)) + int2d(Th)(v) + on(1, u = 0.0);
real[int] rhs = vPb(0, Wh);

set(A, sparams = "-ksp_view");
Wh<real> u;                         // local solution

A = vPb(Wh, Wh);
real memory = PetscMemoryGetCurrentUsage();
u[] = A^-1 * rhs;
memory = PetscMemoryGetCurrentUsage() - memory;
if(mpirank == 0)
    cout << memory << " bytes of memory in usage" << endl;

real[int] err = A * u[];            // global matrix-vector product
real[int] transpose = A' * u[];
exchange(A, rhs, scaled = true);
err -= rhs;

macro def(u)u//
plotMPI(Th, u, Pk, def, real, cmm = "Global solution");
u[] = err;
plotMPI(Th, u, Pk, def, real, cmm = "Global residual");

Wh<real> Rb[1];
Rb[0] = 1;
set(A, sparams = "-pc_type gamg -ksp_type gmres -ksp_max_it 200", nearnullspace = Rb);
u[] = 0.0;
u[] = A^-1 * rhs;
plotMPI(Th, u, Pk, def, real, cmm = "Global solution");
{
    mesh ThGlobal = square(getARGV("-global", 40), getARGV("-global", 40)); // global mesh
    mesh ThSplit = trunc(ThGlobal, 1, split = ThRefinementFactor);
    fespace WhGlobal(ThSplit, Pk);
    WhGlobal sol, solReduced;
    int[int] rest = restrict(Wh, WhGlobal, n2o);
    {
        real[int] tmp;
        ChangeNumbering(A, u[], tmp);
        ChangeNumbering(A, u[], tmp, inverse = true, exchange = false);
    }
    sol[] = 0;
    for[i, v : rest] sol[][v] = u[][i];
    mpiReduce(sol[], solReduced[], processor(0), mpiSUM);
    if(mpirank == 0)
        plot(solReduced, cmm = "Reconstructed global solution");
}
