# Stokes with time stepper - PETSc
**URL:** https://community.freefem.org/t/stokes-with-time-stepper-petsc/2593
**Category:** General Discussion
**Created:** [July 10, 2023, 4:10pm UTC](https://community.freefem.org/t/stokes-with-time-stepper-petsc/2593 "2023-07-10T16:10:29Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@rirane](https://community.freefem.org/u/rirane)
#### Post date: [July 10, 2023, 4:10pm UTC](https://community.freefem.org/t/stokes-with-time-stepper-petsc/2593/1 "2023-07-10T16:10:29Z")
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](https://community.freefem.org/uploads/short-url/zZVPQ6fYW6B5zAisYBHvno2AoP3.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\*(u_v+uB_vB) + mu\*(dx(u)\*dx(v) + dy(u)_dy(v) + dx(uB)dx(vB) + dy(uB)dy(vB) )- pq1.e-6- q_div(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\*(uold_v+uoldB_vB)
> )+ 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”);
> }
---
### Author:  [@prj](https://community.freefem.org/u/prj)
#### Post date: [July 11, 2023, 4:34am UTC](https://community.freefem.org/t/stokes-with-time-stepper-petsc/2593/2 "2023-07-11T04:34:47Z")
The error is:
```auto
63 : rhs = vPbf(def(u) [u, uB, uC],Wh) error operator <6C_args>, <7E_Array>,
```
The fix is:
```auto
63 : rhs = vPbf(0,Wh);
```
Like your line 46.
---
### Author:  [@rirane](https://community.freefem.org/u/rirane)
#### Post date: [July 11, 2023, 4:54am UTC](https://community.freefem.org/t/stokes-with-time-stepper-petsc/2593/3 "2023-07-11T04:54:46Z")
Thank you for your answer!
But if I set
```auto
63 : rhs = vPbf(0,Wh);
```
my right-hand side won’t be updated with the computed velocity from the last time step.
---
### Author:  [@prj](https://community.freefem.org/u/prj)
#### Post date: [July 11, 2023, 5:05am UTC](https://community.freefem.org/t/stokes-with-time-stepper-petsc/2593/4 "2023-07-11T05:05:19Z")
Well, you need to update it yourself.
---
### Author:  [@rirane](https://community.freefem.org/u/rirane)
#### Post date: [July 11, 2023, 8:03am UTC](https://community.freefem.org/t/stokes-with-time-stepper-petsc/2593/5 "2023-07-11T08:03:36Z")
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):
```auto
// 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 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:
```auto
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?
---
### Author:  [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [July 11, 2023, 8:29am UTC](https://community.freefem.org/t/stokes-with-time-stepper-petsc/2593/6 "2023-07-11T08:29:44Z")
See here: [Finite element](https://doc.freefem.org/documentation/finite-element.html#variational-form-sparse-matrix-pde-data-vector)
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.
---
### Author:  [@rirane](https://community.freefem.org/u/rirane)
#### Post date: [July 11, 2023, 8:45am UTC](https://community.freefem.org/t/stokes-with-time-stepper-petsc/2593/7 "2023-07-11T08:45:11Z")
Thank you for the clarification. Now everything seems to work fine!
[petsc-stokes.edp](https://community.freefem.org/uploads/short-url/qT3vWNac64Kp6Q08dFeGbJ1Nsjq.edp) (1.7 KB)