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