Convect( ) in parallel

General Discussion
6

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);
}
1 Like