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?

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 !

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)

11

In order to save a reference solution computed on a fine mesh, and later compare it to another run of your code on a rough mesh, you can follow the following sketch

// assumed known mesh Th, fespace Vh, approximate solution u in Vh

// save the mesh to a file "mesh.msh"
savemesh(Th,"mesh.msh");
// save the solution to a file "uref.d"
{
 ofstream fileref("uref.d");
 fileref << u[] << endl;
}

//-------------------------------

//recover a previously saved solution "uref.d" and its mesh "mesh.msh"
mesh Ths=readmesh("mesh.msh");
fespace Vs(Ths,....);//define the fespace exactly as it was defined when the reference solution was saved to a file
Vs us;
{
 ifstream filereadref("uref.d");
 filereadref >> us[] ;
}

// compute the error between u and the reference solution us
// the error is computed on the rough mesh Th (computation on Ths is also possible)
Vh usinterp=us;//interpolate the reference solution us defined on Ths to a function usinterp on Th
real error = sqrt(int2d(Th)((u-usinterp)^2));
12

You have to handle the fespaces in accordance to their definitions

– Case of components defined separately. When you call the array of coordinates, you have to treat separately each fespace

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,du1;
Wh rho1,drho1;
u1[]-=du1[];
u2[]-=du2[];

– Case of components defined as a vector. When you call the array of coordinates, use a single call

fespace VWh(Th, [P2,P1], periodic=[[2, y], [4, y], [1, x], [3, x]]);
VWh [u1,rho1],[du1,drho1];
u1[]-=du1[];

I think that in order to be safe with periodic spaces, you have to use the second method.
Or use the composite syntax. It means in the case of your test1.edp, write

    solve NewtonNSMS(<[du1], [du2], [drho1], [drho2], [dp], [dmu1], [dmu2]>,
                     <[v1], [v2], [psi1], [psi2], [q], [xi1], [xi2]>) =
13

The unknown variables I need to solve for are u, ρ, μ and p. Among them, u is a vector, while the other three are scalars. Why do we only update u in the Newton iteration step? I attempted to update all variables, yet it resulted in NaN values.

14

As I mentioned, when the fespace is defined as a vector, as in
fespace VWh(Th, [P2,P1]), here the vector has two components, the first in P2, the second in P1, you define (for example)
VWh [u,rho];
Then u[] represents all the coordinates of [u,rho], hence the coordinates of u AND the coordinates of rho.

16

I would recommend you to use only the final time to compute the error, and

18

The syntax error was due to a wrong definition [un1ref, un2ref,rhon1ref,rhon2ref,pnref,mun1ref,mun1ref] (mun1ref is repeated) instead of [un1ref,un2ref,rhon1ref,rhon2ref,pnref,mun1ref,mun2ref].

In order to define an array of vector finite elements, you can do as (for example)

fespace Vh(Th,[P2,P1]);
Vh [u,v],[u1,v1];
Vh[int] [ua,va](3);
[u,v]=[x+y,x-y];
[ua[0],va[0]]=[u,v];
[u1,v1]=[ua[0],va[0]];

Your code corrected with this approach:
ff.edp (17.4 KB)

19

Okay, thank you for your help. When using your previously proposed method to store solutions in files, we encountered memory problems during index access under fine mesh conditions, so I am testing a new approach.

21

The way you manage with the vector of coordinates is wrong (or at least misleading) for again the same reason:
when you write

fespace UUPhrefread(Threfread, [P2,P2,P1,P1,P1,P1,P1], periodic=[[2,y],[4,y],[1,x],[3,x]]);
UUPhrefread [u1r, u2r, rho1r, rho2r, pr, mu1r, mu2r];

the expressions u1r[] or u2r[] or rho1r[]… all refer to the same object : the vector of coordinates of the whole vector [u1r, u2r, rho1r, rho2r, pr, mu1r, mu2r]
Thus you need to write/read it only once.

However your result should be correct (writing and saving the data 7 times instead of once is useless but should not give a wrong result).
Without the full code I cannot say. Maybe the second method " store the reference solution directly into an array" is not coded correctly. I gave a model for that in the file at Periodic Boundary Conditions for Mixed Finite Element Spaces - #18 by fb77 (in the case of considering all the intermediate times).

23

I find not real mistake. Significant discrepancies come from
– You do not have the same parameters V1,V2
– You save sumMuErrSq,sumPErrSq without of with square roots
Minor corrections are in
ffall.edp (16.3 KB)
ffM.edp (16.5 KB)
(with different mesh parameters)