//ff-mpirun -np 8 Helmholtz-2d-all.edp -ffddm_schwarz_method oras -waven 100 -ffddm_geneo_nu 60 -ffddm_geneo_threshold 0.9 -ns -wg

// 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 -ffddm_schwarz_method oras -waven 20 -ffddm_geneo_nu 20 -ffddm_geneo_threshold 0.9

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 k = getARGV("-waven",20.);
func f = 100*exp(-((x-.5)^2+(y-.5)^2)*k*50);

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 Neumann = 3;

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

int mysplit = 2;

real npplo = 8./mysplit;

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

Thc=change(Thc,refe=chlab);

mesh Th = trunc(Thc,1,split=mysplit);

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 VarfEprec(varfName, meshName, PhName)
    varf varfName(u,v) =
  int2d(meshName)(-(k^2-1i*epsilonEprec)*u*v+Grad(u)'*Grad(v))
				+ int1d(meshName,Robin)(1i*k*u*v)
				+ on(Dirichlet,u=0); // EOM

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

// k-weighted inner product
macro VarfIP(varfName, meshName, PhName)
    varf varfName(u,v) =
  int2d(meshName)(k^2*u*v+Grad(u)'*Grad(v))
        + on(Dirichlet,u=0); // EOM

// mass for DtN
macro Varfmass(varfName, meshName, PhName)
    varf varfName(u,v) =
  int1d(meshName,10)(u*v)
        + on(Dirichlet,u=0); // EOM

// tag Dirichlet dofs for DtN
macro VarfDirichlet(varfName, meshName, PhName)
    varf varfName(u,v) = on(Dirichlet,u=1); // EOM

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

ffddmbuildDmeshAug(Hmesh,Th,mpiCommWorld)

ffddmbuildDfespace(HFE,Hmesh,complex,def,init,Pk)

vsym = 0;
vtgv = -1;
vtgvelim = -1;

ffddmsetup(H,HFE,Varf,VarfOpt)

complex[int] rhs(1);

ffddmbuildrhs(H,Varfrhs,rhs)

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

HFEVhi<complex> def(u);

if (mpirank == 0) cout << endl << "ONE-LEVEL:" << endl;

u[] = HfGMRES(x0, rhs, 1.e-6, 200, "right");

Hwritesummary

if (mpirank == 0) cout << endl << "TWO-LEVEL (coarse mesh):" << endl;

ffddmcoarsemeshsetup(H,Thc,VarfEprec,VarfEprec)

u[] = HfGMRES(x0, rhs, 1.e-6, 200, "right");

Hwritesummary

if (mpirank == 0) cout << endl << "HARMONIC GENEO:" << endl;

ffddmharmonicgeneosetup(H,Varf,VarfIP,VarfIP,k^0.5)

u[] = HfGMRES(x0, rhs, 1.e-6, 200, "right");

Hwritesummary

if (mpirank == 0) cout << endl << "EXTENDED HARMONIC GENEO:" << endl;

ffddmextendedharmonicgeneosetup(H,Varf,VarfIP)

u[] = HfGMRES(x0, rhs, 1.e-6, 200, "right");

Hwritesummary

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

ffddmdtnsetup(H,Varf,Varfmass,VarfDirichlet,k)

u[] = HfGMRES(x0, rhs, 1.e-6, 200, "right");

Hwritesummary

if (mpirank == 0) cout << endl << "Hk-GENEO:" << endl;

ffddmgeneohelmholtzsetup(H,Varf,VarfIP)

u[] = HfGMRES(x0, rhs, 1.e-6, 200, "right");

Hwritesummary

ffddmplot(HFE,real(u), "Global solution");
