//ff-mpirun -np 8 natural_convection.edp -raspart -ffddm_schwarz_method oras -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 -raspart -ffddm_schwarz_method oras -Ra 1e6 -global 20

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

include "ffddm.idp"

searchMethod = 1;

macro def(i)[i, i#B, i#C, i#D]// EOM     // vector field definition
macro init(i)[i, i, i, i]// EOM        // vector field initialization
macro grad(u)[dx(u), dy(u)]// EOM   // two-dimensional gradient
macro Grad(u1,u2)[dx(u1),dy(u1),dx(u2),dy(u2)]//
macro div(u)(dx(u) + dy(u#B))// EOM
macro UgradU(uw,u)[[uw,uw#B]'*grad(u),[uw,uw#B]'*grad(u#B)]//
func Pk = [P2, P2, P1, P1];             // finite element space 

real Ra=getARGV("-Ra", 1.e7);
real Pr=0.71;
real Rey=sqrt(Ra/Pr);
real IRe=1./Rey, IRa=Ra/(Rey*Rey*Pr), IPr=1./(Rey*Pr) ;

real Thot=0.5, Tcold=-0.5;

real eps = 1e-7; // penalization in Navier-Stokes

int newtonMax = 10; 

real coef;

int np = getARGV("-global", 40);

int mysplit = 4;

mesh Thc = square(np/mysplit, np/mysplit, [x, y]);
mesh ThGlobal = trunc(Thc,1,split=mysplit);

//macro NCFEdefplot(u)u#D//

ffddmbuildDmeshAug(NCmesh,ThGlobal,mpiCommWorld)
ffddmbuildDfespaceAug(NCFE,NCmesh,real,def,init,Pk)

//ffddmbuildAug(NC,ThGlobal,real,def,init,Pk,mpiCommWorld)

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

macro NCwithhpddmkrylov()1//

ffddmsetupinit(NC,NCFE)

if (NCmeshpCS > 1) {
	int[int] nupart(Thc.nt);
	nupart=0;
	if (mpirank == 0)
		scotch(nupart, Thc, NCmeshpCS);
	if (mpirank < NCmeshpCS) {
		broadcast(processor(0,NCmeshcommCS),nupart);
		Thc=change(Thc,fregion=nupart[nuTriangle]);
	}
}

ffddmcoarsemeshsetupinit(NC,Thc)

NCVhCoarse [uG,uGB,pG,TG];

ffddmset(NC,verbosity,1)

NCFEVhi [ui,uiB,pi,Ti], [upi,upiB,ppi,Tpi];

real[int] rhs;

macro Varf(varfName, meshName, VhName)
	VhName [uG,uGB,pG,TG];
	NCfromVhi(ui[],VhName,uG[])
	
    varf varfName([uw, uwB, pw, Tw], [v, vB, q, W]) =
     int2d(meshName)(- div(uw)*q -div(v)*pw - eps*pw*q
         + UgradU(uG,uw)'*[v,vB]
         + UgradU(uw,uG)'*[v,vB]
         + ( Grad(uw,uwB)'*Grad(v,vB))*IRe
         - coef*IRa*Tw*vB
         -[uG,uGB]'*grad(W)*Tw
         -[uw,uwB]'*grad(W)*TG
         + grad(Tw)'*grad(W)*IPr)
    + on(1,2,3,4, uw=0,uwB=0)
    + on(2,Tw=Tcold) + on(4,Tw=Thot);
// EOM

macro VarfCS(varfName, meshName, VhName)
    varf varfName([uw, uwB, pw, Tw], [v, vB, q, W]) =
     int2d(meshName,mpirank)(- div(uw)*q -div(v)*pw - eps*pw*q
         + UgradU(uG,uw)'*[v,vB]
         + UgradU(uw,uG)'*[v,vB]
         + ( Grad(uw,uwB)'*Grad(v,vB))*IRe
         - coef*IRa*Tw*vB
         -[uG,uGB]'*grad(W)*Tw
         -[uw,uwB]'*grad(W)*TG
         + grad(Tw)'*grad(W)*IPr)
    + on(1,2,3,4, uw=0,uwB=0)
    + on(2,Tw=Tcold) + on(4,Tw=Thot);
// EOM

macro VarfOpt(varfName, meshName, VhName)
    varf varfName([uw, uwB, pw, Tw], [v, vB, q, W]) =
     int2d(meshName)(- div(uw)*q -div(v)*pw - eps*pw*q
         + UgradU(ui,uw)'*[v,vB]
         + UgradU(uw,ui)'*[v,vB]
         + ( Grad(uw,uwB)'*Grad(v,vB))*IRe
         - coef*IRa*Tw*vB
         -[ui,uiB]'*grad(W)*Tw
         -[uw,uwB]'*grad(W)*Ti
         + grad(Tw)'*grad(W)*IPr)
    + int1d(meshName,10)(20./hTriangle*IRe*[uw,uwB]'*[v,vB]+20*1./hTriangle*[pw, Tw]'*[q, W])
    + on(1,2,3,4, uw=0,uwB=0)
    + on(2,Tw=Tcold) + on(4,Tw=Thot);
// EOM
   
macro Varfrhs(varfName, meshName, VhName)
	varf varfName([uw, uwB, pw, Tw], [v, vB, q, W]) =
     int2d(meshName)(
         UgradU(ui,ui)'*[v,vB]       
         - [ui,uiB]'*grad(W)*Ti       
     )
     + on(1,2,3,4, uw=0,uwB=0)
     + on(2,Tw=Tcold) + on(4,Tw=Thot);
// EOM

int nstep= max((int(log10(Ra)-4)^2), 3)*max(log10(Ra)-7.,1.);
real tol = 1e-6;

[ui,uiB,pi,Ti]=[0,0,0,Thot+x*(Tcold-Thot)];
nstep *=2;

int iter=0;

for(int step=1; step <= nstep; ++step)
{
   iter++;

    real err=1e100;
    real errp = err;
    coef = (step/real(nstep))^3; // continuation sur le coef en temperature ...
    if (mpirank == 0)
      cout<< "---------------  continuation coef*Ra=" << coef*Ra << endl;
    real tolk= coef<1 ? 1e-2 : tol;

   // Iteration de Newton ---------------------------
   for(int niter=0 ; niter < newtonMax ; niter++){
   	
	 NCfromVhi(ui[],NCVhCoarse,uG[])

	 ffddmsetup(NC,NCFE,Varf,VarfOpt)
	 
	 vtgv = 1.e+30/NCmeshpCS;	 
	 ffddmcoarsemeshsetup(NC,Thc,VarfCS,null) 
	 vtgv = -1;
	 
     ffddmbuildrhs(NC,Varfrhs,rhs)

     ui[] = NCfGMRES(ui[], rhs, 1.e-6, 200, "left");
     
     upi[] -= ui[];
     real reduce = upi[].linfty;
     upi[]  = ui[];
     mpiAllReduce(reduce, err, mpiCommWorld, mpiMAX);

     if(mpirank == 0)
	   cout << "niter = " << niter << " err =  " << err << endl;

     if(err > errp*100) break;
     if(err < tolk) break;
   }
   
   ffddmplot(NCFE,Ti, "Global solution");

}

NCwritesummary
