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.