Stokes with time stepper - PETSc

General Discussion
1

Hi,
I am new to FreeFem and PETSc and need some help. I want to solve a time dependent stokes equation, where the righthandside is updated with the old velocity uold in each iteration step. Until the time loop everything works fine but then I get the error "Error line number 63, in file petsc-stokes.edp, before token ". So something is wrong with the rhs in the loop. Do I need to set the value of the rhs in the loop in a different way? Does anybody has an idea how to fix this problem/script?

petsc-stokes.edp (1.7 KB)

// run with MPI: ff-mpirun -np 4 script.edp

load “PETSc” // PETSc plugin
macro dimension()2// EOM // 2D or 3D
include “macro_ddm.idp” // additional DDM functions
macro def(i)[i, i#B, i#C]// EOM // vector field definition
macro init(i)[i, i, i]// EOM // vector field initialization
macro grad(u)[dx(u), dy(u)]// EOM // two-dimensional gradient

macro div(u) (dx(u) + dy(u#B)) //EOM
func Pk = [P2, P2, P1]; // finite element space

//Parameters
real dt = 0.1;
real tau = 1/dt;
real T = 1.;
real mu = 0.1;
int n = 10;

mesh Th = square(getARGV(“-global”,n),getARGV(“-global”,n));

Mat A;
macro ThRefinementFactor()getARGV(“-split”, 0)//
MatCreate(Th, A, Pk);

//Fespace
fespace Wh(Th, Pk);

//variational formultation
varf vPbA([u, uB, p], [v, vB, q]) = int2d(Th)(
tau*(uv+uBvB) + mu*(dx(u)*dx(v) + dy(u)dy(v) + dx(uB)dx(vB) + dy(uB)dy(vB) )- pq1.e-6- qdiv(u)
)+ on(1, 2, 4, u=0, uB=0)+ on(3, u=1, uB=0)
;

//variational formualtion for the righthandside
varf vPbf([uold, uoldB, p], [v, vB, q]) = int2d(Th)(
tau*(uoldv+uoldBvB)
)+ on(1, 2, 4, uold=0, uoldB=0)+ on(3, uold=1, uoldB=0)
;
real[int] rhs = vPbf(0,Wh);

set(A, sparams = “-pc_type lu”);
Wh def(u);

A = vPbA(Wh, Wh);
u = A^-1 * rhs;

macro def2(u)[u, u#B] //EOM
macro def1(u)u //EOM
plotMPI(Th, def2(u), [P2,P2], def2, real, cmm = “Global velocity”);

// Time loop
int M = T/dt;
for(int m = 0; m < M; m++){

rhs = vPbf(def(u),Wh);

set(A, sparams = “-pc_type lu”);
Wh def(u);

A = vPbA(Wh, Wh);
u = A^-1 * rhs;

plotMPI(Th, def2(u), [P2,P2], def2, real, cmm = “Global velocity”);
}

2

The error is:

   63 : 	rhs = vPbf(def(u)  [u, uB, uC],Wh) error operator   <6C_args>, <7E_Array>, <PP5v_fes>

The fix is:

   63 : 	rhs = vPbf(0,Wh);

Like your line 46.

3

Thank you for your answer!
But if I set

   63 : 	rhs = vPbf(0,Wh);

my right-hand side won’t be updated with the computed velocity from the last time step.

4

Well, you need to update it yourself.

5

Okay thanks. Is there any documentation why rhs has to be set to vPbf(0,Wh) and can’t be something like vPbf([uold,uoldB,p],Wh)?

I am trying to update it myself for the first iteration step (so without the loop):

 //  run with MPI:  ff-mpirun -np 4 script.edp

load "PETSc"                        // PETSc plugin
macro dimension()2// EOM            // 2D or 3D
include "macro_ddm.idp"             // additional DDM functions

macro def(i)[i, i#B, i#C]// EOM     // vector field definition
macro init(i)[i, i, i]// EOM        // vector field initialization
macro grad(u)[dx(u), dy(u)]// EOM   // two-dimensional gradient
macro div(u) (dx(u) + dy(u#B)) //EOM
func Pk = [P2, P2, P1];             // finite element space

//Parameters
real dt = 0.1;
real tau = 1/dt;
real T = 1.;
real mu = 0.1;
int n = 10;
real uold = 0.;
real uoldB = 0.;




mesh Th = square(getARGV("-global",n),getARGV("-global",n));

Mat A;
macro ThRefinementFactor()getARGV("-split", 0)//
MatCreate(Th, A, Pk);

//Fespace
fespace Wh(Th, Pk);

//variational formultation for the left-hand side
varf vPbA([u, uB, p], [v, vB, q]) = int2d(Th)(
tau*(u*v+uB*vB) + mu*(dx(u)*dx(v) + dy(u)*dy(v) + dx(uB)*dx(vB) + dy(uB)*dy(vB) )
- p*q*1.e-6
- q*div(u)
)
+ on(1, 2, 4, u=0, uB=0)
+ on(3, u=1, uB=0)
;




//variational formualtion for the right-hand side
varf vPbf([u, uB, p], [v, vB, q]) = int2d(Th)(
tau*( (u+uold)*v + (uB+uoldB)*vB )
)
+ on(1, 2, 4, u=0, uB=0)
+ on(3, u=1, uB=0)
;
real[int] rhs = vPbf(0,Wh);


set(A, sparams = "-pc_type lu");
Wh<real> def(u);

A = vPbA(Wh, Wh);
u[] = A^-1 * rhs;

macro def2(u)[u, u#B] //EOM
macro def1(u)u //EOM
plotMPI(Th, def2(u), [P2,P2], def2, real, cmm = "Global velocity");

I get a similar error:

 Error line number 49, in file petsc-stokes.edp, before  token

If I define vPbf without uold and uoldB everything works fine. Now uold and uoldB are real just parameters like tau. Do I need to define uold and uoldB in a different way?

6

See here: Finite element

Note that vPbf(0,Wh) does not set [uold,uoldB,p] to zero. It just tells FreeFEM to evaluate the expression as a linear form rather than a bilinear one.

7

Thank you for the clarification. Now everything seems to work fine!
petsc-stokes.edp (1.7 KB)