richard.edp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

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) {
  Th=adaptmesh(Th,h,ratio=1.1);
//  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);
  hn=h;
  temps += dt;
}

cout <<  " hmax = " << hmax << endl;