Richard Problem

Snapshots for the time evolution of the solution

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;
real hmax=0;

// $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;};  // gauche

border bb1(t=x0,0)    { x=t; y=ymax ;label=1;};   // haut     1
border bb2(t=xmax,x0) { x=t; y=ymax ;label=2;};  // haut   2

border bc1(t=y0,0) { x=xmax; y=ymax-t ;label=2;};  // droite
border bc2(t=ymax,y0) { x=xmax; y=ymax-t ;label=3;};  // droite

border bd(t=0,xmax)   { x=t; y=0; label=4;};   // bas

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);

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;
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)
+int1d(Th,3)(pena*h*v)-int1d(Th,3)(pena*((ymax-y)-y0)* v)

;

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

Ahdt=0; Kh=1;
// Richard;
// plot(hd,wait=1,cmm="hd ----");
real temps=0;
for(int i=0;i<100;i++)
{
string scmm="h + temps "+int(temps)/3600+"h "+ ((temps)%3600/60.) + "mn ";
for(int k=0;k<3;k++)
{
Kh=K(h);
Ahdt=A(h)/dt;

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");
Richard;
cout << " h: min " << h[].min << " max  " << h[].max <<endl;
hmax=h[].max ;
// plot(h,wait=0,cmm=scmm,viso=viso);
//  hhh = h <0;
}
if(i%10==1) {