Convect with periodic boundaries

General Discussion
1

Hello all, I’m new here but not so new to FreeFem++.

I want to implement the convect operator on a periodic function (with a periodic velocity field) in 2d.

mesh Th=square(Nx,Ny,[Lxx,Lyy]);
fespace Den(Th,P2,periodic=[[2,y],[4,y],[1,x],[3,x]]);
Den c0=whatever0, c;
fespace Vel(Th,[P1b,P1b],periodic=[[2,y],[4,y],[1,x],[3,x]]);
Vel [u1,u2]=[whatever1,whatever2];

c = convect([u1,u2], -dt, c0);

This does not covect properly through the periodic boundaries (the density on the edges remains static). I am sure that Eulerian advection works in principle but I prefer the Semi-Langrangian route.
The workaround I am considering is:

  1. Define C0=c0, [U1,U2]=[u1,u2] on a broader mesh that extends beyond the boundaries in all directions. C0 would be periodic over Lx,Ly, same for [U1,U2].
  2. C = convect([U1,U2], -dt, C0);
  3. Restrict C to c (ie the original mesh).

Is convect supposed to work on periodic boundaries? If not, is there a more elegant implementation?

Thanks in advance

2

Yes I know this problem but it is hard to code, No bodie do the jotb Sorry,
and it is also hard to detect.

3

Thanks, Frederic. I will try the workaround I have in mind and share the code if it works.

4

Here’s a working but inefficient code that relies on interpolations. This assumes one wants to solve some PDE for psi0 with periodic BCs after applying a periodic convect.

  1. Define two meshes. Th is the periodic rectangle of interest, ETh is an extension, n is the number of elements to expand (can use n=1 if Courant<1).
mesh Th=square(Nx,Ny,[Lx*x,Ly*y]);
mesh ETh=square(Nx+2*n,Ny+2*n,[-h*n+(Lx+h*2*n)*x,-h*n+(Ly+h*2*n)*y]);
  1. FE spaces (ETh spaces are not technically periodic)
fespace Den(Th,P2,periodic=[[2,y],[4,y],[1,x],[3,x]]);
     Den psi0;
fespace Vel(Th,[P1b,P1b],periodic=[[2,y],[4,y],[1,x],[3,x]]);
     Vel [u1,u2];
fespace EDen(ETh,P2);
     EDen Epsi0, Epsi;
fespace EVel(ETh,[P1b,P1b]);
     EVel [Eu1,Eu2];
  1. Wraparound functions
func real WrapX(real X) { return ( X<0?X+Lx:(X>Lx?X-Lx:X) );}
func real WrapY(real Y) { return ( Y<0?Y+Ly:(Y>Ly?Y-Ly:Y) );}
  1. Periodic convect (psi0,[u1,u2] are explicit)
Epsi0 = psi0(WrapX(x),WrapY(y));    //extend concentration
[Eu1,Eu2] = [u1(WrapX(x),WrapY(y)), u2(WrapX(x),WrapY(y))];     //extend velocity
Epsi = convect([Eu1,Eu2], -dt, Epsi0);     //convect on ETh
psi0 = Epsi(x,y);     //restrict back to Th

This works (barring typos). However, it is highly inefficient due to the use of interpolations. So, what would be a more efficient way to implement step 4?

5

The problem is hard to code in all the case.