//ff-mpirun -np 4 Helmholtz-2d-HPDDM-BGMRES.edp -wg -raspart -ffddm_schwarz_method oras -waven 30

// If you have openmpi you may need to add the option --oversubscribe to allow more processes than the number of cores available on your computer

// for the make check:
// NBPROC 4
// PARAM -raspart -ffddm_schwarz_method oras -waven 20

macro dimension 2// EOM            // 2D or 3D

include "ffddm.idp"

macro def(i)i// EOM                         // scalar field definition
macro init(i)i// EOM                        // scalar field initialization
func Pk = P2;

real xf = 0.5;
real yf = 0.5;

real k = getARGV("-waven",40.);
func f = 100*exp(-((x-xf)^2+(y-yf)^2)*k*10);

real lambda = 2*pi/k;

real epsilonprob = 0;
real be = getARGV("-betaEprec",1.);
real epsilonEprec = 0;//k^be;

int Dirichlet = 1;
int Robin = 2;

//int[int] chlab=[1, Robin, 2, Robin, 3, Robin, 4, Robin];
int[int] chlab=[1, Dirichlet, 2, Dirichlet, 3, Dirichlet, 4, Dirichlet];

real npplo = 15;

mesh Th = square(npplo*1./lambda,npplo*1./lambda);    // global mesh

Th=change(Th,refe=chlab);

macro Grad(u) [dx(u),dy(u)] // EOM

macro Varf(varfName, meshName, PhName)
    varf varfName(u,v) =
  int2d(meshName)(-(k^2-1i*epsilonprob)*u*v+Grad(u)'*Grad(v))
				+ int1d(meshName,Robin)(1i*k*u*v)
				+ on(Dirichlet,u=0); // EOM

// for the preconditioner
macro VarfOpt(varfName, meshName, PhName)
    varf varfName(u,v) =
  int2d(meshName)(-(k^2-1i*epsilonEprec)*u*v+Grad(u)'*Grad(v))
				+ int1d(meshName,Robin,10)(1i*k*u*v)
				+ on(Dirichlet,u=0); // EOM

macro Varfrhs(varfName, meshName, PhName)
    varf varfName(u,v) =
                -int2d(meshName)(f*v)
                + on(Dirichlet,u=0); // EOM

ffddmbuildDmesh(Hmesh,Th,mpiCommWorld)
ffddmbuildDfespace(HFE,Hmesh,complex,def,init,Pk)

vsym = 2; // symmetric, not hermitian
vtgv = 1.e+30;
vtgvelim = 1.e+30;

macro Hwithhpddm()1//

ffddmsetup(H,HFE,Varf,VarfOpt)

int nbrhs = 10;

complex[int] brhs(nbrhs*HFEVhi.ndof);
complex[int] bu(nbrhs*HFEVhi.ndof);

complex[int] rhs(0);
for (int i=0; i<nbrhs; i++) {
      xf = 0.4+i*0.4/nbrhs;
      ffddmbuildrhs(H,Varfrhs,rhs)
      brhs(i*HFEVhi.ndof:(i+1)*HFEVhi.ndof-1) = rhs;
}

complex[int] x0(brhs.n);
x0 = 0;

if (mpirank == 0) cout << "Pseudo-Block GMRES :" << endl;

bu = HfGMRES(x0, brhs, 1.e-6, 200, "right");

Hwritesummary

if (mpirank == 0) cout << endl << "Block GMRES :" << endl;

set(HhpddmOP,sparams="-hpddm_H_krylov_method bgmres");

bu = HfGMRES(x0, brhs, 1.e-6, 200, "right");

Hwritesummary

HFEVhi<complex> def(u);
u[] = bu(0:HFEVhi.ndof-1); 

ffddmplot(HFE,real(u), "u[0]");
