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
 2real nu = 0.1;
 3real T=1.;
 4real dt = 0.1;
 5
 6// Mesh
 7mesh Th = square(10, 10);
 8
 9// Fespace
10fespace Vh(Th, P2)
11Vh u, v;
12Vh uu, vv;
13Vh uold=0, vold=0;
14
15fespace Qh(Th, P1);
16Qh p;
17Qh pp;
18
19// Problem
20problem 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
36int m, M = T/dt;
37for(m = 0; m < M; m++){
38    stokes;
39    uold = u;
40    vold = v;
41}
42
43// Plot
44plot(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:

 1fespace Xh(Th, [P2, P2, P1]);
 2varf 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
14varf 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
22varf 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
27matrix A = aa(Xh, Xh, solver=UMFPACK);
28matrix B = bb(Xh, Xh);
29real[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:

1real[int] sol(Xh.ndof), aux(Xh.ndof);
2for (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:

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