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

load "PETSc"                        // PETSc plugin
macro dimension()3// EOM            // 2D or 3D
include "macro_ddm.idp"             // additional DDM functions

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

int[int] LL = [1,2, 1,2, 1,1];
mesh3 Th = cube(getARGV("-global", 10), getARGV("-global", 10), getARGV("-global", 10), [x, y, z], label = LL);    // 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) = int3d(Th)(grad(u)' * grad(v)) + int3d(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);
u[] = A^-1 * rhs;

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");
