FreeFEM Documentation on GitHub

stars - forks

# Tutorial to write a transient Stokes solver in matrix form

Consider the following script to solve a time dependent Stokes problem in a cavity

 1 // Parameters
2 real nu = 0.1;
3 real T=1.;
4 real dt = 0.1;
5
6 // Mesh
7 mesh Th = square(10, 10);
8
9 // Fespace
10 fespace Vh(Th, P2)
11 Vh u, v;
12 Vh uu, vv;
13 Vh uold=0, vold=0;
14
15 fespace Qh(Th, P1);
16 Qh p;
17 Qh pp;
18
19 // Problem
20 problem stokes (u, v, p, uu, vv, pp)
21     = int2d(Th)(
22           (u*uu+v*vv)/dt
23         + nu*(dx(u)*dx(uu) + dy(u)*dy(uu) + dx(v)*dx(vv) + dy(v)*dy(vv))
24         - p*pp*1.e-6
25         - p*(dx(uu) + dy(vv))
26         - pp*(dx(u) + dy(v))
27     )
28     - int2d(Th)(
29           (uold*uu+vold*vv)/dt
30     )
31     + on(1, 2, 4, u=0, v=0)
32     + on(3, u=1, v=0)
33     ;
34
35 // Time loop
36 int m, M = T/dt;
37 for(m = 0; m < M; m++){
38     stokes;
39     uold = u;
40     vold = v;
41 }
42
43 // Plot
44 plot(p, [u, v], value=true, wait=true, cmm="t="+m*dt);


Every iteration is in fact of the form $$A[u,v,p] = B[uold,vold,pold] + b$$ where $$A,B$$ are matrices and $$b$$ is a vector containing the boundary conditions. $$A,B,b$$ are constructed by:

 1 fespace Xh(Th, [P2, P2, P1]);
2 varf aa ([u, v, p], [uu, vv, pp])
3     = int2d(Th)(
4           (u*uu+v*vv)/dt
5         + nu*(dx(u)*dx(uu) + dy(u)*dy(uu) + dx(v)*dx(vv) + dy(v)*dy(vv))
6         - p*pp*1.e-6
7         - p*(dx(uu) + dy(vv))
8         - pp*(dx(u) + dy(v))
9     )
10     + on(1, 2, 4, u=0, v=0)
11     + on(3, u=1, v=0)
12     ;
13
14 varf bb ([uold, vold, pold], [uu, vv, pp])
15     = int2d(Th)(
16           (uold*uu+vold*vv)/dt
17     )
18     //+ on(1, 2, 4, uold=0, vold=0)
19     //+ on(3, uold=1, vold=0)
20     ;
21
22 varf bcl ([uold, vold, pold], [uu, vv, pp])
23     = on(1, 2, 4, uold=0, vold=0)
24     + on(3, uold=1, vold=0)
25     ;
26
27 matrix A = aa(Xh, Xh, solver=UMFPACK);
28 matrix B = bb(Xh, Xh);
29 real[int] b = bcl(0, Xh);


Note that the boundary conditions are not specified in $$bb$$. Removing the comment // would cause the compiler to multiply the diagonal terms corresponding to a Dirichlet degree of freedom by a very large term (tgv); if so $$b$$ would not be needed, on the condition that $$uold=1$$ on boundary 3 initially. Note also that b has a tgv on the Dirichlet nodes, by construction, and so does A.

The loop will then be:

1 real[int] sol(Xh.ndof), aux(Xh.ndof);
2 for (m = 0; m < M; m++){
3     aux = B*sol; aux += b;
4     sol = A^-1 * aux;
5 }


There is yet a difficulty with the initialization of sol and with the solution from sol. For this we need a temporary vector in $$X_h$$ and here is a solution:

1 Xh [w1, w2, wp] = [uold, vold, pp];
2 sol = w1[]; //cause also the copy of w2 and wp
3 for (m = 0; m < M; m++){
4     aux = B*sol; aux += b;
5     sol = A^-1 * aux;
6 }
7 w1[]=sol; u=w1; v= w2; p=wp;
8 plot(p, [u, v], value=true, wait=true, cmm="t="+m*dt);


The freefem team agrees that the line sol=w1[]; is mysterious as it copies also w2 and wp into sol. Structured data such as vectors of $$X_h$$ here cannot be written component by component. Hence w1=u is not allowed.

Table of content