Periodic Boundary Conditions

Here is an example that works with periodic boundary conditions on finite element vector spaces (thanks to Daniel LeRoux and Frederic Hecht).

Propagation of a Kelvin wave solved by Crank-Nicolson with $\theta=0.5$ and using P1 and P1nc finite elements for pressure and velocity, Coriolis force is considered also. The equations are

\begin{array}{rcl} \partial_t \vec{u}-\alpha\, \vec{u}^\perp+ \nabla p & =&0 \quad\mbox{on } \Omega,\\ \partial_t p- \nabla \cdot \vec{u} & =&0 \quad\mbox{on } \Omega. \end{array}

where $\Omega$ is a channel-like. Boundary conditions are periodic on entrance and exit and normal velocity cero on the sides of the channel.

Two different snapshots of the wave evolution

// Definition du bord du domaine
// -----------------------------
//
//                         label 3
//        (-16,4) ------------------------------ (16,4)
//                |                            |
//                |                            |
//       label 4  |                            |  label 2
//                |                            |
//                |                            |
//       (-16,-4) ------------------------------ (16,-4)
//                         label 1
//

border bottom(t=0,1){x=32.0*t-16.0; y=-4.0; label=1;};
border top(t=0,1){x=-32.0*t+16.0; y=4.0; label=3;};
border left(t=0,1){x=-16.0; y=-8.0*t+4.0; label=4;};
border right(t=0,1){x=16.0; y=8.0*t-4.0; label=2;};

// On veut imposer des conditions au bord periodiques sur les bords 2 (right) et 4 (left)
// et une condition de glissement u.n=0 (vitesse normale nulle) sur les bords 1 (bottom) et 3 (top).
// ------------------------------------------------------------------------------------------------

func perio=[[4,y],[2,y]];

//  Construction du maillage
//  ------------------------

int m=16;
int n=4;

mesh Th=buildmesh(bottom(m)+top(m)+left(n)+right(n), fixeborder=true);
plot(Th,wait=0,ps="mesh.eps");

//  Choix des parametres physiques et parametres d'adimensionalisation
//  ------------------------------------------------------------------

real omega=7.2921*1e-5;   // vitesse angulaire de la terre
real R=6.371*1e6;         // rayon de la terre (en m)
real g=9.81;              // acceleration pesanteur
real H=1.6309888;         // hauteur moyenne du fluide (equivalent depth)

real E=sqrt(2.0*R*omega/sqrt(g*H)); // en fait il s'agit de E^(1/4)
real U=sqrt(g*H);

real phi0=0.0;            // latitude
real a=E*sin(phi0);       // parametre de Coriolis = a + b * y
real b=cos(phi0);

//  Choix du schema temporel : schema de Crank Nicolson
//  ---------------------------------------------------

real dt=0.1;              // pas de temps
real theta=0.85;           // theta schema ---> Crank Nicolson (theta=0.5)
real mu=theta*dt;         // mu et nu : coefficients pour le schema temporel
real nu=(1.0-theta)*dt;

//  Choix des espaces d'approximation : P1 non conforme pour la vitesse, P1 pour la pression
//  ----------------------------------------------------------------------------------------

fespace Uh(Th,P1nc,periodic=perio);
fespace Ph(Th,P1,periodic=perio);
fespace UPh(Th,[P1nc,P1nc,P1],periodic=perio);

// u,v,p          : vitesse et pression au temps n+1
// uold,vold,pold : vitesse et pression au temps n
// uu,vv          : fonctions test pour la vitesse
// pp             : fonction  test pour la pression

//Uh u,v,uu,vv,uold,vold;
//Ph p,pp,pold;
UPh [u,v,p];
UPh [uu,vv,pp];
UPh [uold,vold,pold];
//  Ecriture du probleme
//  --------------------
// Bug in Problem with Periodic with not same FE space .
problem Kelvin([u, v, p],[ uu, vv, pp])=

int2d(Th)(u*uu - (a+b*y)*mu*v*uu + mu*dx(p)*uu
+ v*vv + (a+b*y)*mu*u*vv + mu*dy(p)*vv
+ p*pp - mu*(u*dx(pp)+v*dy(pp)))

+ int2d(Th)(-uold*uu - nu*(a+b*y)*vold*uu + nu*dx(pold)*uu
- vold*vv + nu*(a+b*y)*uold*vv + nu*dy(pold)*vv
- pold*pp - nu*(uold*dx(pp)+vold*dy(pp)))

//   Conditions no-normal flow (u.n=0) sur les bords 1 (bottom) et 3 (top)
//   Les conditions periodiques ont ete imposees sur les bords 2 (right) et 4 (left)

+ on(1,3,v=0.0); // CL : u.n=0

//  Solution initiale (propagation onde de Kelvin)
//  ----------------------------------------------

[u,v,p] =[-cos(pi*x/4.0)*exp(-y*y/2.0),0,-cos(pi*x/4.0)*exp(-y*y/2.0)];
//v=0.0;
//p=-cos(pi*x/4.0)*exp(-y*y/2.0);

//  Propagation des ondes de Kelvin : on simule 100 pas de temps
//  ------------------------------------------------------------

int imax=100;

for (int i=0; i<=imax; i++)
{
uold[]=u[];

Kelvin;

cout << "i:=";
cout << i;
cout << "\n";

cout <<i+1 << "\n ---";

plot(p, value=0, wait=0, fill=1,cmm="theta = "+theta+" dt = "+dt+" imax = "+imax);

}