1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
|
// This a the rotating hill problem with one turn.
// First 1/2 turn is a convection equation and second 1/2 a convection diffusion
border a(t=0, 2*pi) { x = cos(t); y = sin(t); }; // the unit circle
mesh th2 = buildmesh(a(35)); // triangulates the disk
mesh th= trunc(th2,1,split=2);
fespace Vh(th,P1);
fespace Vh2(th2,P2);
Vh2 v = exp(-10*((x-0.3)^2 +(y-0.3)^2)); // initial condition
plot(v);
real dt = 0.17,t=0; // time step
Vh u1 = y, u2 = -x; // rotation velocity
int i;
Vh2 vv,vo; // work Finite element function
for ( i=0; i< 20 ; i++) {
t += dt;
vo=v;
v=convect([u1,u2],-dt,vo); // convect v by u1,u2, dt seconds, results in f
// convec(u1,u2,dt,v,v) won't work
plot(v,cmm="convection: t="+t + ", min=" + v[].min + ", max=" + v[].max,wait=0);
};
problem A(v,vv,solver=CG) = int2d(th2)(v*vv/dt + 0.01*(dx(v)*dx(vv)+dy(v)*dy(vv)) )
+ int2d(th2)(-vv*convect([u1,u2],-dt,vo)/dt) + on(a,v=0);
for ( i=0; i< 20 ; i++)
{
t += dt;
vo=v;
A; // solve le problem A
plot(v,cmm="convection& diffusive: t="+t + ", min=" + v[].min + ", max=" + v[].max);
};
plot(v,wait=1);
|
