// NBPROC 50
// PARAM -ns -3  -k 2 -gmres 3
/*
  a first true parallele example first freefem++ 
  Ok up to 200 proc for a Poisson equation.. 
  See the Doc for full explaination

  F Hecht Dec. 2010. 
  -------------------
usage :
ff-mpirun [mpi parameter]  DDM-Schwarz-Lap-2dd.edp  [-glut ffglut]  [-n N] [-k K]  [-d D] [-ns] [-gmres [0|1]
 argument: 
   -glut ffglut : to see graphicaly the process
   -n N:  set the mesh cube split NxNxN
   -d D:  set debug flag D must be one for mpiplot 
   -k K:  to refined by K all  elemnt
   -ns: reomove script dump
   -gmres 0   : use iterative schwarz algo.  
          1   :  Algo GMRES on residu of schwarz algo.
          2   :  DDM GMRES 
          3   :  DDM GMRES with coarse grid preconditionner (Good one)  
*/
//
//  I don't seem to have the MPICG module available.
//
load "MPICG"  load "medit"  load "metis"
include "getARGV.idp"
include "MPIplot.idp"
include "DDM-Schwarz-macro.idp"
include "AddLayer2d.idp"
include "DDM-funcs-v2.idp"


searchMethod=1; // more safe seach algo (warning can be very expensive in case lot of ouside point) 
// 0 by default

assert(version >=3.11);
real[int] ttt(10);int ittt=0;
macro settt {ttt[ittt++]=mpiWtime();}//


verbosity=getARGV("-vv",0);
int vdebug=getARGV("-d",1);
int ksplit=getARGV("-k",1);
int nloc = getARGV("-n",10);
string sff=getARGV("-p,","");
int gmres=getARGV("-gmres",0); 
int nC = getARGV("-N" ,max(nloc/10,5)); 
int sizeoverlaps=1; // size of overlap
bool RAS=1;  // Global Variable ..



if(mpirank==0 && verbosity)
  cout << " vdebug: " << vdebug << " kspilt "<< ksplit << " nloc "<< nloc << " sff "<< sff <<"."<< endl;

string sPk="P2-2gd";     
func Pk=P2;

func bool  plotMPIall(mesh &Th,real[int] & u,string  cm)
{ PLOTMPIALL(mesh,Pk, Th, u,{ cmm=cm,nbiso=20,fill=1,dim=3,value=1}); return 1;}

mpiComm comm(mpiCommWorld,0,0);// trick : make a no split mpiWorld 


int iiipart= mpiRank(comm); // current partition number 


if(iiipart==0)  cout << " Final N=" << ksplit*nloc << " nloc =" << nloc << " split =" << ksplit <<  endl;
int[int] l111=[1,2,2,2]; 
settt 

mesh Thg=square(nloc,nloc*5,[x,5*y],label=l111);
mesh ThC=square(nC,nC*5,[x,5*y],label=l111);//   Caarse Mesh
fespace VhC(ThC,P1); // of the coarse problem.. 


BuildPartitioning(sizeoverlaps,mesh,Thg,Thi,aThij,RAS,pii,jpart,comm,vdebug)

if(ksplit>1)
{
 for(int jp=0;jp<jpart.n;++jp)
  aThij[jp]  = trunc(aThij[jp],1,split=ksplit);
 Thi =   trunc(Thi,1,split=ksplit);
}

BuildTransferMat(iiipart,mesh,Pk,1,[0],
                 Thi,Whi,Whij,Thij,aThij,Usend,Vrecv,jpart,vdebug)

Whi ui,vi; 



/* the definition of the Problem .... */
func G=1.; /* ok  */
func F=1.; /* ok  */
macro grad(u) [dx(u),dy(u)] //
// warning for Dir. BC. the last win 
varf vPb(U,V)= int2d(Thi)(grad(U)'*grad(V)) + int2d(Thi)(F*V) + on(10,U=0)+on(1,U=G) ; //');// for emacs
varf vPbC(U,V)= int2d(ThC)(grad(U)'*grad(V))  +on(1,U=0) ; //');// for emacs
varf vPbon10(U,V)=on(10,U=1)+on(1,U=0);

varf vPBC(U,V)=on(1,U=G);


real[int] onG10 = vPbon10(0,Whi); // on 1 
real[int] Bi=vPb(0,Whi);


matrix Ai = vPb(Whi,Whi,solver=sparsesolver); 

DMMDeffuncAndGlobals(Lap2,comm,jpart,Whi,Vhc,1,Ai,vPbC,onG10,Pii,Usend,Vrecv,[0])

Lap2CheckUpdate(); 
  
Whi u=0,v;
 

u[]=vPBC(0,Whi,tgv=1); 
real eps=1e-10;
Lap2DDMSolver(Bi,u,v,gmres,eps,vdebug)


real errg =1,umaxg;

{ 
  real umax = u[].max,umaxg;
  real[int] aa=[umax], bb(1);
  mpiAllReduce(aa,bb,comm,mpiMAX);
  errg=bb[0];
  if(iiipart==0)
    cout << " umax global  = " << bb[0] << " Wtime = " << (ttt[ittt-1]-ttt[ittt-2])  << " s " <<  " " << Lap2kiter <<  endl;
}

Lap2Saveff(sff,eps,ksplit,nloc,sizeoverlaps);