sorry for my late reply.

@yongxing thank you for pointing this out. I had already noticed that the code doesnâ€™t work well for larger time steps, but since I use small time steps, I didnâ€™t investigate it. maybe, itâ€™s a better idea to switch to newton method and a SNES solver for NS equations, and I should update the code to work on that basis.

I used such technique (putting `convect`

inside `varf`

) because I has seen a wide range of FF examples for solving NS writing it in that way (like this, this, and this to name a few), but apparently it has the problem @prj mentioned (at least in parallel).

for clarifying this to people seeing this later, I prepared these minimal examples, one with `convect`

in `varf`

and one using `convectParallel`

. they produce the same results in small time steps, but the results will be different in larger time steps.

first example:

```
// using convect in the varf statement
load "PETSc"
macro dimension()2// EOM
include "macro_ddm.idp"
real dt = 0.02;
int n = 200;
mesh Mesh = square(n, n);
fespace SpaceP1(Mesh, P1);
SpaceP1 c, cold;
c = exp(-5*((x-0.3)^2 + (y-0.3)^2));
Mat A;
createMat(Mesh, A, P1)
c = c;
// SpaceP1 u1 = 0.3, u2 = 0.3;
SpaceP1 u1 = y, u2 = x;
varf convection(c, v) = intN(Mesh)(c/dt*v)
+ intN(Mesh)(convect([u1, u2], -dt, cold)*v/dt);
matrix Loc;
real[int] b(SpaceP1.ndof);
Loc = convection(SpaceP1, SpaceP1);
A = Loc;
set(A, sparams="");
for (int i=0; i < 50; i++)
{
if (mpirank == 0)
cout << "step " + i << endl;
cold = c;
b = convection(0, SpaceP1);
c[] = A^-1 * b;
savevtk("vel-1.vtu", Mesh, c, dataname="c", append = i?true:false);
}
```

second example:

```
// using convectParallel
load "PETSc"
macro dimension()2// EOM
include "macro_ddm.idp"
real dt = 0.02;
int n = 200;
mesh Mesh = square(n, n);
fespace SpaceP1(Mesh, P1);
SpaceP1 c;
c = exp(-5*((x-0.3)^2 + (y-0.3)^2));
buildDmesh(Mesh);
c = c;
// SpaceP1 u1 = 0.3, u2 = 0.3;
SpaceP1 u1 = y, u2 = x;
for (int i=0; i < 50; i++)
{
if (mpirank == 0)
cout << "step " + i << endl;
convectParallel(Mesh, [u1, u2], -dt, c, 20)
savevtk("vel-2.vtu", Mesh, c, dataname="c", append = i?true:false);
}
```