changewait;precise;
/*
border(1,0,4,30) begin x:=t; y:=0; end;
border(2,0,1,11) begin x:=4; y:=t; end;
border(3,0,4,30) begin x:=4-t; y:=1; end;
border(4,0,1,11) begin x:=0; y:=1-t; end;
border(0,0,4,30) begin x:=4-t; y:=0.1; end;

buildmesh(1500);
*/
ro =1;
u=2.9;
v=0;
p=1/1.4;
rou=ro*u;
rov=ro*v;
h = p/0.4 + rou*u /2 + rov*v/2;
pu=p*u; pv=p*v;

dt:= 4/(2.9*100);

h3:=1.52819/0.4+ 0.5*1.69997
*(2.61934*2.61934+0.50633*0.50633);
h4:=1/(1.4*0.4) + 2.9*2.9/2;
i:=1;

iter(200) begin
eps = 0.5*dt*dt*p/ro;
D=1-dx(u)*dt-dy(v)*dt;
solve(ro,i) begin
onbdy(3) ro=1.69997;
onbdy(4) ro =1;
pde(ro) id(ro)-laplace(ro)*eps= D*convect(ro,u,v,dt);
end;
plot(ro); i:=-1;
solve(rou,i) begin
onbdy(3) rou=1.69997*2.61934;
onbdy(4) rou =2.9;
pde(ro) id(rou)-laplace(rou)*eps= D*convect(rou,u,v,dt)
-dt*dx(p);
end;

solve(rov,i) begin
onbdy(3) rov=-1.69997*0.50633;
onbdy(4) rov =0;
pde(ro) id(rov) -laplace(rov)*eps= D*convect(rov,u,v,dt)
-dt*dy(p);
end;
 
solve(h,i) begin
onbdy(3) h=h3;
onbdy(4) h =h4;
pde(ro) id(h)-laplace(h)*eps= D*convect(h,u,v,dt)
-dt*(dx(pu)+dy(pv));
end;

u = rou/ro; 
v = rov/ro;
p = 0.4*(h- rou*u/2 - rov*v/2);
pu = p*u;
pv = p*v;
v=v -2*v* (y<0.05);
/*
D=1-dx(u)*dt-dy(v)*dt;
v=v -100*v*one(y-0.01);
ro= D*convect(ro,u,v,dt);
rou= D*convect(rou,u,v,dt)-dt*dx(p);
rov= D*convect(rov,u,v,dt)-dt*dy(p);
h= D*convect(h,u,v,dt)
        -dt*(dx(pu)+dy(pv));
plot(ro);
*/
end;