//ff-mpirun  -np 4 Richards-2d.edp -wg -ns -ffddm_disable_plots

// 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_disable_plots -nt 20

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

include "ffddm.idp"

int nt = getARGV("-nt", 100);

searchMethod = 1; // safe search algo, use brute force in case of missing point

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

real Ks=0.01,
     hg=30,
     thetas=0.3,
     eta = 6.55,
     m = 0.173,
     n = 2/(1-m);
     
     
real z0=215;
real q0=15/3600.;
real dt=60;

// $A(h) - \partial h / \partial t - div(K(h)(\nabla(h-y)) = f $ dans $ \Omega$
//  -K(h)(\nabla(h-y)). n = q0 $ sur $ \Gamma_1$
//  h = h_0$ sur $\Gamma_0$
//  + condition initial $d_d$
//   A(h) = h <0 ? C(h) : 0; 
//  


// remarque z == y 
real xmax = 300, ymax=300, x0=60, y0= 215;

border ba(t=0,ymax)   { x=0; y=ymax-t ;label=2;};  // left   
   
border bb1(t=x0,0)    { x=t; y=ymax ;label=1;};   // top     1 
border bb2(t=xmax,x0) { x=t; y=ymax ;label=2;};  // top   2   
 
border bc1(t=y0,0) { x=xmax; y=ymax-t ;label=2;};  // right     
border bc2(t=ymax,y0) { x=xmax; y=ymax-t ;label=3;};  // right  
   
border bd(t=0,xmax)   { x=t; y=0; label=4;};   // bottom

int Gamma0=3;
int Gamma1=1;

int nn=20;
int nn1=nn*x0/xmax,nn2=nn-nn1;
int ny1=nn*y0/ymax,ny2=nn-ny1;
plot(ba(nn)+bb1(nn1)+bb2(nn2)+bc1(ny1)+bc2(ny2)+bd(nn),wait=1);
mesh Th=buildmesh(ba(nn)+bb1(nn1)+bb2(nn2)+bc1(ny1)+bc2(ny2)+bd(nn));
plot(Th,wait=1);

mesh Thc = Th;

fespace Vh(Th,P1);
Vh h,v,hhh;

macro theta(h) (thetas*(1+((abs(h)/hg))^n)^(-m))//
macro dtheta(h) (m*n*thetas*(1+((abs(h)/hg))^n)^(-m-1)*(((abs(h)/hg))^(n-1))/hg)
//
macro A(h)  ( (h<=0)* dtheta(h) ) 
//
macro K(h) (Ks*((h<=0)*((theta(h)/thetas)^eta)+ (h>0)))
//

Vh hd= -y0+(ymax-y); // bof bof ????
Vh hn=hd,hh;
Vh Ahdt,Kh;

int nbiso=20;
real[int] viso(3+(75+110/2)/5);

{int k=0;
for(int i=-75;i<0;i+=5)
 viso(k++)=i;
 viso(k++)=-0.5;
 viso(k++)=0.;
 viso(k++)=0.5;
for(int i=5;i<=110;i+=5*2)
 viso(k++)=i;
}
/* 
problem Richard(h,v) =
  int2d(Th)( Ahdt * h * v+ Kh* (dx(h)*dx(v)+dy(h)*dy(v)) )
- int2d(Th)( Ahdt* hn*v - Kh* dy(v) )
- int1d(Th,Gamma1)(q0*v)
+ on(3,h=(ymax-y)-y0)
;
*/

real pena=1e10;

macro VarfPb(varfName, meshName, VhName)

  VhName mh, mhn;
  RfromVhi(hi[],VhName,mh[])
  RfromVhi(hin[],VhName,mhn[])
  VhName Ahdt = A(mh)/dt;
  VhName Kh = K(mh);

  varf varfName(h,v) =
  int2d(meshName)( Ahdt * h * v+ Kh* (dx(h)*dx(v)+dy(h)*dy(v)) )
- int2d(meshName)( Ahdt* mhn*v - Kh* dy(v) )
- int1d(meshName,Gamma1)(q0*v)
 +int1d(meshName,3)(pena*h*v)-int1d(Th,3)(pena*((ymax-y)-y0)* v)

; // EOM

plot(hn,wait=1,cmm=" hd ");

macro Rdefmplot(u)u//EOM

//macro Rwithhpddmkrylov()1//EOM

ffddmbuildDmeshAug(R,Th,mpiCommWorld)

ffddmbuildDfespace(R,R,real,def,init,Pk)

ffddmsetupinit(R,R)// 

ffddmset(R,verbosity,1)// 

RVhi hi=hd, hin=hd;

real[int] rhs;

// Richard;
// plot(hd,wait=1,cmm="hd ----");
real temps=0;
for(int ii=0;ii<nt;ii++)
{
  string scmm="h + temps "+int(temps)/3600+"h "+ ((temps)%3600/60.) + "mn ";
  for(int k=0;k<3;k++) 
  {

	ffddmsetup(R,R,VarfPb,null)
	
	ffddmbuildrhs(R,VarfPb,rhs)
	
	rhs = -rhs;
  
  //cout << " "<< Kh[].min << " " << Kh[].max << endl;
 // plot(Ahdt,fill=1,value=1,wait=1,cmm="Ahdt");
 // plot(Kh,fill=1,value=1,wait=1,cmm="Kh");
 // plot(Kh,wait=1,cmm="Kh");

	hi[] = RfGMRES(hi[], rhs, -1.e-6, 200, "right");

  // forces a refresh of the FE function h, so that it is defined on the currrent mesh Th
  // in particular, the underlying array h[] has the correct size
  h = 0.;
  
	RfromVhi(hi[],Vh,h[]);

  // plot(h,wait=0,cmm=scmm,viso=viso);
//  hhh = h <0;
  }
  if(ii%10==1) {
  Th=adaptmesh(Th,h,ratio=1.1,err=2e-2);
  broadcast(processor(0),Th);
  
  // the mesh has changed ; we create a new decomposition and obtain new local meshes
  ffddmbuildDmeshAug(R,Th,mpiCommWorld)
  ffddmbuildDfespace(R,R,real,def,init,Pk)
  // interpolate the solution at the previous time step on the new local mesh 
  hi = h;
  
//  plot(Th,h,cmm="h ",value=1,wait=1);
  }
 // plot(hhh,wait=1,cmm="h < 0");

  plot(h,wait=0,cmm=scmm,viso=viso,value=1);
//  plot(h,cmm="h ",value=1);
  hin = hi;
  temps += dt;
}

Rwritesummary
