// P.-Y. LAGREE & O. DEVAUCHELLE
// FEBRUARY 2007
// FULLY IMPLICIT VERSION - RESCALED

int AspectRatioH=6,n=0;

real W=1,L=W*AspectRatioH;
real S=.05,R=60;
real beta=5,gamma=1;
real epsilon=.05,kx=2*pi/L,ky=pi/W;
real dt=.01,t=0;
real precision=10e-6;
real theta=0.5; 
real seed=0.123;

func real charru(real entry)	// erosion law divided by tau
	{
	return pow(entry,beta-1);
	}

func real charrubytau(real entry)	// erosion law divided by tau^2
	{
	return pow(entry,beta-2);
	}

func real dcharru(real entry)	// derivative of the above law divided by tau
	{
	return (beta-1)*pow(entry,beta-3);
	}

int npty=10,nptx=npty*AspectRatioH/2;

mesh Th=square(nptx,npty,[-L/2+L*x,-W/2+W*y]);

fespace Vh(Th,P2,periodic=[[2,y],[4,y]]);

Vh h,eta;
Vh dth,dteta;
Vh etatest,dthtest,dtetatest;

real dhm;
real cpu=clock();


h=-1-epsilon*sin(kx*x)*sin(ky*y);
epsilon=.1;
h=-1;
{
for(int i=1;i<10;i++)
{int is;
seed = (seed+pi)^3;
is=seed;
seed = (seed - is)*((-1)^is);
	cout << i<< " " << seed << " " << is << " " << clock() -cpu <<endl;
	h = h + 0.1*epsilon*seed*sin(2*pi*seed+i*pi/L*x);
	
	h = h + 0.1*epsilon*seed *sin(2*pi*seed+i*pi/W*y);
}	
}
plot(h,wait=1);


// FIRST RESOLUTION

problem flow(eta,etatest,eps=precision,solver=GMRES) =
	  int2d(Th)( h^3/(R*S)*(dx(eta)*dx(etatest) + dy(eta)*dy(etatest)) )
	+ int2d(Th)( 3*h^2*eta*dx(etatest) )
	- int2d(Th)( h^3*dx(etatest) )
	;



// THE BIG ONE : TIME EVOLUTION


// FIRST EQUATION : DEFINITION OF TAU


macro taux(eta,h)
	( eta-h + h/(R*S)*dx(eta) ) //

macro tauy(eta,h)
	( h/(R*S)*dy(eta) ) //
	
macro normtau(eta,h)
	( sqrt(taux(eta,h)^2 + tauy(eta,h)^2) ) //

macro dttaux(dteta,dth,eta,h)
	( dteta-dth + dth*dx(eta)/(R*S) + h/(R*S)*dx(dteta) ) //

macro dttauy(dteta,dth,eta,h)
	( dth*dy(eta)/(R*S) + h/(R*S)*dy(dteta) ) //

// SECOND EQUATION : FLOW

macro e2t1(dtetatest,h,dth)
	(-3*h^2*dth*dx(dtetatest)) //

macro e2t2(dtetatest,h,dth,eta)
	(6*h*dth*eta*dx(dtetatest)) //

macro e2t3(dtetatest,h,dteta)
	(3*h^2*dteta*dx(dtetatest)) //


macro e2t4(dtetatest,eta,h,dth)
	(3*h^2*dth/(R*S)*(dx(eta)*dx(dtetatest) + dy(eta)*dy(dtetatest))) //


macro e2t5(dtetatest,dteta,h)
	(h^3/(R*S)*(dx(dteta)*dx(dtetatest) + dy(dteta)*dy(dtetatest))) //


// LAST EQUATION : SEDIMENT CONSERVATION


macro e3t1(dth,dthtest)
	(-dth*dthtest) //

macro e3t2(dthtest)
	(charru(normtau(eta,h))*(dx(dthtest)*taux(eta,h) + dy(dthtest)*tauy(eta,h))) //

macro e3t3(dthtest,h)
	(-gamma/R*charru(normtau(eta,h))*normtau(eta,h)*(dx(h)*dx(dthtest) + dy(h)*dy(dthtest))) //

macro e3t4(dthtest)
	(theta*dt*dcharru(normtau(eta,h))*(dx(dthtest)*taux(eta,h) + dy(dthtest)*tauy(eta,h))*(taux(eta,h)*dttaux(dteta,dth,eta,h) + tauy(eta,h)*dttauy(dteta,dth,eta,h))) //


macro e3t5(dthtest,h)
	(-theta*dt*dcharru(normtau(eta,h))*normtau(eta,h)*gamma/R*(taux(eta,h)*dttaux(dteta,dth,eta,h) + tauy(eta,h)*dttauy(dteta,dth,eta,h))*(dx(h)*dx(dthtest) + dy(h)*dy(dthtest))) //


macro e3t6(dthtest,h)
	(-theta*dt*charrubytau(normtau(eta,h))*gamma/R*(taux(eta,h)*dttaux(dteta,dth,eta,h) + tauy(eta,h)*dttauy(dteta,dth,eta,h))*(dx(h)*dx(dthtest) + dy(h)*dy(dthtest))) //

macro e3t7(dthtest)
	(theta*dt*charru(normtau(eta,h))*(dttaux(dteta,dth,eta,h)*dx(dthtest) + dttauy(dteta,dth,eta,h)*dy(dthtest))) //

macro e3t8(dthtest,dth)
	(-theta*dt*charru(normtau(eta,h))*normtau(eta,h)*gamma/R*(dx(dth)*dx(dthtest) + dy(dth)*dy(dthtest))) //


// THE FULL PROBLEM

problem time(dth,dteta,dthtest,dtetatest,eps=precision,solver=GMRES) =

	  int2d(Th)( e2t1(dtetatest,h,dth) )
	+ int2d(Th)( e2t2(dtetatest,h,dth,eta) )
	+ int2d(Th)( e2t3(dtetatest,h,dteta) )
	+ int2d(Th)( e2t4(dtetatest,eta,h,dth) )
	+ int2d(Th)( e2t5(dtetatest,dteta,h) )

	+ int2d(Th)( e3t1(dth,dthtest) )
	+ int2d(Th)( e3t2(dthtest) )
	+ int2d(Th)( e3t3(dthtest,h) )
	+ int2d(Th)( e3t4(dthtest) )
	+ int2d(Th)( e3t5(dthtest,h) )
	+ int2d(Th)( e3t6(dthtest,h) )
	+ int2d(Th)( e3t7(dthtest) )
	+ int2d(Th)( e3t8(dthtest,dth) )

	;



// TIME RESOLUTION

// INITIAL DERIVATION

flow;

//plot(dy(h),value=1,cmm="dy(h) propre",fill=1,wait=1);
//plot(dy(h),value=1,cmm="dy(h) propre",fill=1,wait=1);

//plot(eta,value=1,cmm="initial eta",fill=1,wait=0);


// TIME LOOP

while (t<40)
	{
    	dhm=(h[].max-h[].min)/2/epsilon;	
     	//plot(h,value=1,cmm="h at t = "+t+" dhm= "+dhm,fill=1,wait=0);
       plot(h,value=1,cmm="h at t = "+t+" dhm= "+dhm,fill=1,wait=0,ps="storage/t_"+t+"_R_"+R+"_slope_"+S+".eps");

	//plot(h,value=1,cmm="h at t="+t,fill=1,wait=0);
    	//flow;

	time;

	h=h+dt*dth;
	eta=eta+dt*dteta;
    
	t=t+dt;
	cout << " CPU time = " << clock()-cpu << " t= " << t << endl;
	}


cout << " CPU time = " << clock()-cpu << endl;