// P.-Y. LAGREE & O. DEVAUCHELLE
// FEBRUARY 2007
// PASSAGE A GMRES

int AspectRatioH=6,n=0;

real W=1,L=W*AspectRatioH;
real slope=.05,R=120;
real beta=5,gamma=1;
real epsilon=.1,kx=2*pi/L,ky=pi/W;
real dt=.001,t=0,theta=.5;
real etam;

int npty=10,nptx=npty*AspectRatioH;

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

plot(Th,wait=0);


//fespace Vh(Th,P1,periodic=[[2,y],[4,y],[1,x],[3,x]]);
fespace Vh(Th,P2,periodic=[[2,y],[4,y]]);

Vh eta,etatest,hs,hstest;
Vh h,h2,h3;
Vh taux,tauy,normtau,charru;
Vh dxeta,dyeta;
Vh dxh,dyh;
Vh f,dxf,dyf,dxftest,dyftest,pyl,c;
real dhm;
real cpu=clock();
real epsg=1e-15;
real seed=0.123;


//h=-1+epsilon*sin(ky*y);
h=-1-epsilon*sin(kx*x)*sin(ky*y);
h=-1+epsilon/2-epsilon*exp(-(y-0.05*sin(kx*x))^2*20);
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 << 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);

//h=-1+epsilon*sin(kx*x);
//h=-1+epsilon*exp(-(2*x/L)^2-(2*y/W)^2);

problem fluid(eta,etatest,solver=GMRES,eps=epsg) = 
    - int2d(Th)(h3*dx(eta)*dx(etatest) + h3*dy(eta)*dy(etatest))
    + int2d(Th)(R*3*h2*slope*etatest*dx(eta))
    + int2d(Th)((R*h3)*dx(etatest))
    ;

problem derivation(dxf,dyf,dxftest,dyftest,solver=GMRES,eps=epsg) = 
      int2d(Th)(dxf*dxftest + dyf*dyftest)
    + int2d(Th)(f*dx(dxftest) + f*dy(dyftest))
    + int1d(Th,1,3)(-f*dxftest*N.x-f*dyftest*N.y)
    ;

problem sediment(hs,hstest,solver=GMRES,eps=epsg) =
    - int2d(Th)(h*hstest/dt) + int2d(Th)(hs*hstest/dt)
    - int2d(Th)(charru*(taux*dx(hstest)+tauy*dy(hstest)))
    + int2d(Th)(normtau*charru*(1-theta)*gamma/R*(dxh*dx(hstest)+dyh*dy(hstest)))
    + int2d(Th)(normtau*charru*theta*gamma/R*(dx(hs)*dx(hstest)+dy(hs)*dy(hstest)))
    + int1d(Th,1,3)(charru*hstest*N.y*tauy)
    - int1d(Th,1,3)(charru*hstest*gamma/R*normtau*(1-theta)*N.y*dyh)
    ;


while(t<100)
    {
    
    if(n%20==0)
        {
        plot(h,value=1,cmm="h at t = "+t+" dhm= "+dhm,fill=1,wait=0,ps="storage/t_"+t+"_R_"+R+"_slope_"+slope+".eps");
        
        ofstream file("storage/t_"+t+"_R_"+R+"_slope_"+slope+".dat");
        for(int i=0;i<Th.nt;i++)
                 {
            for(int j=0; j<3; j++)
                {
                file << Th[i][j].x << " " << Th[i][j].y << " " << h[][Vh(i,j)] << endl;
                     }
            }
        }

    n=n+1;

    h3=h*h*h;
    h2=h*h;

    fluid;
        etam=int2d(Th)(eta)/(L*W);
    eta=eta-etam;
    //plot(eta,value=1,cmm="eta before",fill=1,wait=1);
    
    f=eta;
    derivation;
    dxeta=dxf;
    dyeta=dyf;
    //plot(eta,value=1,cmm="eta after",fill=1,wait=1);

    f=h;
    derivation;
    dxh=dxf;
    dyh=dyf;

    //plot(dxeta,value=1,cmm="detax at t = "+t,fill=1,wait=1);

    taux=(slope*eta-h)*(R-dxeta)/R;
    tauy=(slope*eta-h)*(-dyeta)/R;
    //taux=-h*(R-dxeta)/R;
    //tauy=-h*(-dyeta)/R;

    normtau=sqrt(taux^2+tauy^2);
    //charru=(normtau-threshold)*(normtau > threshold);
    charru=normtau^(beta-1);

    //plot(taux,value=1,cmm="taux",fill=1,wait=1);
    //plot(tauy,value=1,cmm="tauy",fill=1,wait=1);

    //plot(normtau,value=1,cmm="normtau at t = "+t,fill=1,wait=1);
    //plot(tauy,value=1,cmm="tauy at t = "+t,fill=1,wait=1);
    //plot(charru,value=1,cmm="charru at t = "+t,fill=1,wait=1);
    //plot(dyh,value=1,cmm="dyh at t = "+t,fill=1,wait=0);

    sediment;
        dhm=(h[].max-h[].min)/2/epsilon;
      //  c=dxh*dt/(h-hs);
        //plot(c);
        cout << "t="  << t << "  >>>>>>>>>>>>>>>>>>>>> dhm =" << dhm <<    endl;
 
    h=hs;
    t=t+dt;

}


cout << " CPU time = " << clock()-cpu << endl;
/*
-- 
Olivier Devauchelle
Laboratoire de Modélisation en Mécanique
UPMC, 4 place Jussieu, case 162
75252 PARIS Cedex 05 FRANCE
TÃA
©l: +33 1 44 27 87 16

*/