Periodic Boundary Conditions for Mixed Finite Element Spaces

General Discussion
1

The difficulty with using periodic boundary conditions on separate P2 and P1 finite element spaces in FreeFem++ is the degrees of freedom (DOF) mismatch caused by inconsistent periodic node mapping across different spaces, leading to an assertion failure.
My current code using two separate P2 and P1 spaces triggers the DOF mismatch error, and the relevant code is as follows:
mesh Th = square(cm, cm);
fespace Vh(Th, P2, periodic=[[2, y], [4, y], [1, x], [3, x]]);
fespace Wh(Th, P1, periodic=[[2, y], [4, y], [1, x], [3, x]]);

Vh u1,u2,un1,un2,du1,du2,v1,v2;
Wh rho1,rho2,p,mu1,mu2,rhon1,rhon2,drho1,drho2,dp,dmu1,dmu2,psi1, psi2,q,xi1,xi2;

real V1 = 0.3, V2 = 0.7;
rhon1 = 1.0 + 0.5sin(2pix)sin(2piy);
rhon2 = (1.0 - V1rhon1) / V2;
un1 = sin(2pix)cos(2piy);
un2 = -cos(2pix)sin(2pi*y);

The program crashes with the following error log:
– Square mesh : nb vertices =4225 , nb triangles = 8192 , nb boundary edges 256 rmdup= 0
Problem build of FEspace (2d) (may be : due to periodic Boundary condition missing ) FH
The number of DF must be 53248 and it is 54407
Assertion fail : (snbdf == NbOfDF)
line :937, in file ../femlib/FESpace.cpp

2

The following lines

int cm=10;
mesh Th = square(cm, cm);
fespace Vh(Th, P2, periodic=[[2, y], [4, y], [1, x], [3, x]]);
fespace Wh(Th, P1, periodic=[[2, y], [4, y], [1, x], [3, x]]);

Vh u1,u2,un1,un2,du1,du2,v1,v2;
Wh rho1,rho2,p,mu1,mu2,rhon1,rhon2,drho1,drho2,dp,dmu1,dmu2,psi1, psi2,q,xi1,xi2;

real V1 = 0.3, V2 = 0.7;
rhon1 = 1.0 + 0.5*sin(2*pi*x)*sin(2*pi*y);
rhon2 = (1.0 - V1*rhon1) / V2;
un1 = sin(2*pi*x)*cos(2*pi*y);
un2 = -cos(2*pi*x)*sin(2*pi*y);

terminate correctly on my computer. If your error comes at the definition of a variational formulation, can you give it explicitly so that I can test it?

3

(post deleted by author)

5

(post deleted by author)

6

with this code, after some iterations I get

  -- Solve : 
          min -3.43191e+35  max 5.29554e+27
-nan != -nan diff -nan => Sorry error in Optimization add:  (p) int2d(Th,optimize=0)(...)
 remark if you add  (..  ,   optimize=2) then  you remove this check (be careful); 
  current line = 183
Exec error : In Optimized version 
   -- number :1
Exec error : In Optimized version 
   -- number :1
 err code 8 ,  mpirank 0

It means that the iterative scheme diverges, the error does not come from FreeFem coding. You need to check your formulas and be sure that your algorithm is stable !

7

(post deleted by author)

8

Dear Sylvia,
there are the following errors in your code:
– Since you define fespace UUPh(Th, [P2,P2,P1,P1,P1,P1,P1], the unknown is a vector with 7 components. Then you declare UUPh [du1,du2,drho1,drho2,dp,dmu1,dmu2],[u1,u2,rho1,rho2,p,mu1,mu2]
which are vectors with 7 components. It follows that when you write just u1[] it gives the whole set of dof of the 7-vector. Thus you don’t have to write also u2[] (indeed u2[] is the same as u1[]: the whole set of dof of the 7-vector).
It means that instead of writing

       u1[] -= du1[]; u2[] -= du2[];
        rho1[] -= drho1[]; rho2[] -= drho2[];
        p[] -= dp[];
        mu1[] -= dmu1[]; mu2[] -= dmu2[];

you have just to write u1[] -= du1[]; The same is true for the linfty errors.
– A term with viscosity nu has a wrong sign
– 3 lines with constant terms related to \rho_i^k and \rho^ku^k have no contribution to the derivative of the functional in the Newton method. They have to be removed in that derivative.
– The invariance by adding a constant to the pressure has to be neutralized by adding a small stabilization int2d(Th)(1.e-6*p*q)

The corrected code is
simple_solver.edp (12.1 KB)

9

(post deleted by author)