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
 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
36
37
38
39
40
41
42
43
44
// Parameters
real nu = 0.1;
real T=1.;
real dt = 0.1;

// Mesh
mesh Th = square(10, 10);

// Fespace
fespace Vh(Th, P2)
Vh u, v;
Vh uu, vv;
Vh uold=0, vold=0;

fespace Qh(Th, P1);
Qh p;
Qh pp;

// Problem
problem stokes (u, v, p, uu, vv, pp)
    = int2d(Th)(
          (u*uu+v*vv)/dt
        + nu*(dx(u)*dx(uu) + dy(u)*dy(uu) + dx(v)*dx(vv) + dy(v)*dy(vv))
        - p*pp*1.e-6
        - p*(dx(uu) + dy(vv))
        - pp*(dx(u) + dy(v))
    )
    - int2d(Th)(
          (uold*uu+vold*vv)/dt
    )
    + on(1, 2, 4, u=0, v=0)
    + on(3, u=1, v=0)
    ;

// Time loop
int m, M = T/dt;
for(m = 0; m < M; m++){
    stokes;
    uold = u;
    vold = v;
}

// Plot
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
 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
fespace Xh(Th, [P2, P2, P1]);
varf aa ([u, v, p], [uu, vv, pp])
    = int2d(Th)(
          (u*uu+v*vv)/dt
        + nu*(dx(u)*dx(uu) + dy(u)*dy(uu) + dx(v)*dx(vv) + dy(v)*dy(vv))
        - p*pp*1.e-6
        - p*(dx(uu) + dy(vv))
        - pp*(dx(u) + dy(v))
    )
    + on(1, 2, 4, u=0, v=0)
    + on(3, u=1, v=0)
    ;

varf bb ([uold, vold, pold], [uu, vv, pp])
    = int2d(Th)(
          (uold*uu+vold*vv)/dt
    )
    //+ on(1, 2, 4, uold=0, vold=0)
    //+ on(3, uold=1, vold=0)
    ;

varf bcl ([uold, vold, pold], [uu, vv, pp])
    = on(1, 2, 4, uold=0, vold=0)
    + on(3, uold=1, vold=0)
    ;

matrix A = aa(Xh, Xh, solver=UMFPACK);
matrix B = bb(Xh, Xh);
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
2
3
4
5
real[int] sol(Xh.ndof), aux(Xh.ndof);
for (m = 0; m < M; m++){
    aux = B*sol; aux += b;
    sol = A^-1 * aux;
}

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