# Boundary conditions in eigenvalue problem

**URL:** <https://community.freefem.org/t/boundary-conditions-in-eigenvalue-problem/2629>\
**Category:** General Discussion\
**Created:** [July 31, 2023, 1:50pm UTC](https://community.freefem.org/t/boundary-conditions-in-eigenvalue-problem/2629 "2023-07-31T13:50:32Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![maximefio](https://avatars.discourse-cdn.com/v4/letter/m/a6a055/32.png) [@maximefio](https://community.freefem.org/u/maximefio)\
**Post date:** [July 31, 2023, 1:50pm UTC](https://community.freefem.org/t/boundary-conditions-in-eigenvalue-problem/2629/1 "2023-07-31T13:50:32Z")

</div>

Hello,

I am currently doing a 1 dimensional linear stabiity solver (for fluid mechanics). To do so I solve an eigenvalue problem as AX = lambda BX where A and B are terms from the linearized Navier-Stokes equations and of the mean flow and X is the small perturbation that is unknown (u, v, w, p).

In these matrices A and B must be added boundary conditions, for some case only Neumann or Dirichlet boundary conditions need to be imposed (case m = 0 and m\>1 in page 3). For these two cases, Freefem works perfectly well.However, for the case m=±1, the first perturbation u depends on the second perturbation v as 2u’(0)+mv’(0)=0.

 ![Capture d’écran 2023-07-31 à 15.35.15](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/a/a5b7838334d1a0fdb95af1e79cfe4010f9f35619.png)

Would there be a way to impose such kind of boundary condition if Freefem ? I have put the case that work for m = 0. In the code F stands for u, G for v, H for w and Pr for p.

real ymin = 0.;  
real ymax = 1.;  
border yline(t=ymin,ymax){x=0.;y=t;}  
int np = 1000;

meshL Th=buildmeshL(yline(np));

fespace Xh(Th, P2);  
fespace Mh(Th, P1);

Xh W=1-y^2, V=0, U=0, Re=100.0, Visco=1/Re, PrMean = 0.;  
Xh alpha=10.0, m=0.0;

fespace Vh(Th, [P2, P2, P2, P1]);  
Vh [F, G, H, Pr], [X1, X2, X3, X4];

varf opLHS([F, G, H, Pr], [X1, X2, X3, X4]) =  
int1d(Th)( -(1i_y^2/Re)dy(F)dy(X1)-((1i/Re)((m^2+1)+y^2alpha^2))FX1+y^2_alpha_W_F_X1-(1i)_(y/Re)_F_dy(X1)-(2_1i_m/Re)_G_X1+y^2_Pr_dy(X1)  
+(1i)_(-(y^2/Re)dy(G)dy(X2))-((1i)/Re)((m^2+1)+y^2alpha^2)GX2+y^2_alpha_W_G_X2-(1i)_(y/Re)_G_dy(X2)-(2_m_(1i)/Re)_F_X2+m_y_Pr_X2  
-(1i)_(y^2/Re)_dy(H)dy(X3)-(1i)(y/Re)Hdy(X3)+y^2_alpha_W_H_X3-(1i/Re)_(m^2+alpha^2_y^2)HX3+y^2_dy(W)_F_X3+1i_y^2_alpha_Pr_X3  
-X4\*(F+m_G+alpha_y_H)+y_F\*dy(X4)  
)  
+on(1, F=0, G=0)  
+on(2, F=0, G=0, H=0);

// get the RHS function  
varf opRHS([F, G, H, Pr], [X1, X2, X3, X4]) = int1d(Th)(y^2_F_X1 + y^2_G_X2 + y^2_H_X3)  
+on(1, F=0, G=0)  
+on(2, F=0, G=0, H=0);

---

<div class="post-metadata">

**Author:** ![maximefio](https://avatars.discourse-cdn.com/v4/letter/m/a6a055/32.png) [@maximefio](https://community.freefem.org/u/maximefio)\
**Post date:** [July 31, 2023, 1:55pm UTC](https://community.freefem.org/t/boundary-conditions-in-eigenvalue-problem/2629/2 "2023-07-31T13:55:47Z")

</div>

Maybe I can provide in better way the system of equations :

 ![Capture d’écran 2023-07-31 à 15.34.06](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/4/4488ffbbe7b1e727f22add12c06a8e8a33588999.png)

---

<div class="post-metadata">

**Author:** ![maximefio](https://avatars.discourse-cdn.com/v4/letter/m/a6a055/32.png) [@maximefio](https://community.freefem.org/u/maximefio)\
**Post date:** [July 31, 2023, 1:56pm UTC](https://community.freefem.org/t/boundary-conditions-in-eigenvalue-problem/2629/3 "2023-07-31T13:56:25Z")

</div>

I thank you for your help.

Maxime

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [July 31, 2023, 2:08pm UTC](https://community.freefem.org/t/boundary-conditions-in-eigenvalue-problem/2629/4 "2023-07-31T14:08:26Z")

</div>

While you can achieve what you want using Lagrange multipliers, that is not necessary.

I assume these BC came from the older works of Khorrami from the late 80’s and 90’s. Note that `u+mv=0` and `2u'+mv'=0` necessarily implies `u'=v'=0`, which is a much simpler requirement. You can impose this naturally as a Neumann BC, and, as far as I am aware, this constraint is more popular nowadays for stability analysis in cylindrical coordinates.

---

<div class="post-metadata">

**Author:** ![maximefio](https://avatars.discourse-cdn.com/v4/letter/m/a6a055/32.png) [@maximefio](https://community.freefem.org/u/maximefio)\
**Post date:** [August 1, 2023, 5:44am UTC](https://community.freefem.org/t/boundary-conditions-in-eigenvalue-problem/2629/5 "2023-08-01T05:44:44Z")

</div>

Hello Chris,

I thank you very much for your answer. This study is indeed based on the works from Khorrami as you stated. I have tried what you suggested by imposing u’=v’=0 for the case Re=2200, m=1 and alpha = 1 and I get the following most unstable mode : re(omega)=0.897, im(omega)=-0.078. The values that obtained Khorrami for these same parameters are re(omega)=0.89663 and im(omega)=−0.048114. So the real part of the most unstable mode seems to be matching but the amplification seems different. Would you know if this is what is expected ? (The modes match well both for the real and imaginary part with Khorrami for m=0 and m\>1.) Also, more generally, I would be curious to know how to impose such kind of boundary conditions that are linear conbinations of the perturbations or spatial derivatives, would you have an example of the use of the Lagrange multipliers ? I thank you very much for you help.

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [August 1, 2023, 8:25am UTC](https://community.freefem.org/t/boundary-conditions-in-eigenvalue-problem/2629/6 "2023-08-01T08:25:34Z")

</div>

I would certainly expect your results to agree with Khorrami’s. Would you mind sharing your code as an attachment?

You can see an example with Lagrange multipliers [here](https://doc.freefem.org/examples/finite-element.html#examplelagrangemultipliers) or, if using the PETSc interface, [here](https://joliv.et/FreeFem-tutorial/section_8/example5.edp.html).

---

<div class="post-metadata">

**Author:** ![maximefio](https://avatars.discourse-cdn.com/v4/letter/m/a6a055/32.png) [@maximefio](https://community.freefem.org/u/maximefio)\
**Post date:** [August 1, 2023, 9:27am UTC](https://community.freefem.org/t/boundary-conditions-in-eigenvalue-problem/2629/7 "2023-08-01T09:27:32Z")

</div>

Yes, for sure, the code is the following (I an sorry it seems that I cannot deposit attachements)

```auto
real ymin = 0.;
real ymax = 1.;
border yline(t=ymin,ymax){x=0.;y=t;}
int np = 2000;

meshL Th=buildmeshL(yline(np));

fespace Xh(Th, P2);
fespace Mh(Th, P1);
// interpolate implicitely on P2 space
Xh W=1-y^2, V=0, U=0, Re=2200.0, Visco=1/Re, PrMean = 0.;
Xh alpha=1.0, m=1.0;

fespace Vh(Th, [P2, P2, P2, P1]);
Vh<complex> [F, G, H, Pr], [X1, X2, X3, X4];

varf opLHS([F, G, H, Pr], [X1, X2, X3, X4]) =
int1d(Th)( -(1i*y^2/Re)*dy(F)*dy(X1)-((1i/Re)*((m^2+1)+y^2*alpha^2))*F*X1+y^2*alpha*W*F*X1-(1i)*(y/Re)*F*dy(X1)-(2*1i*m/Re)*G*X1+y^2*Pr*dy(X1)
           +(1i)*(-(y^2/Re)*dy(G)*dy(X2))-((1i)/Re)*((m^2+1)+y^2*alpha^2)*G*X2+y^2*alpha*W*G*X2-(1i)*(y/Re)*G*dy(X2)-(2*m*(1i)/Re)*F*X2+m*y*Pr*X2
           -(1i)*(y^2/Re)*dy(H)*dy(X3)-(1i)*(y/Re)*H*dy(X3)+y^2*alpha*W*H*X3-(1i/Re)*(m^2+alpha^2*y^2)*H*X3+y^2*dy(W)*F*X3+1i*y^2*alpha*Pr*X3
           -X4*(F+m*G+alpha*y*H)+y*F*dy(X4)
           )
           +on(1, H=0, Pr=0)
           +on(2, F=0, G=0, H=0);

varf opRHS([F, G, H, Pr], [X1, X2, X3, X4]) = int1d(Th)(y^2*F*X1 + y^2*G*X2 + y^2*H*X3)
                                              +on(1, H=0, Pr=0)
                                              +on(2, F=0, G=0, H=0);

// assemble matrices
cout << "Assemble matrices" << endl;
cout << "Matrix L"<< endl;
matrix<complex> L = opLHS(Vh, Vh, solver=GMRES, tgv=1.0e30);
cout << "Matrix B"<< endl;
matrix<complex> B = opRHS(Vh, Vh, solver=GMRES, eps=1.0e-20, tgv=1.0e30);

// save all the information
load "bfstream"
cout << "Write files "<< endl;
complex[int] C(1);
int[int] I(1),J(1);
[I,J,C]=L;
cout << "Matrix L"<< endl;
{
	ofstream fid("LIN/LNS-I.dat");
	fid.write(I);
}
{
	ofstream fid("LIN/LNS-J.dat");
	fid.write(J);
}
{
	ofstream fid("LIN/LNS-C.dat");
	fid.write(C);
}
[I,J,C]=B;
cout << "Matrix B"<< endl;
{
	ofstream fid("LIN/B-I.dat");
	fid.write(I);
}
{
	ofstream fid("LIN/B-J.dat");
	fid.write(J);
}
{
	ofstream fid("LIN/B-C.dat");
	fid.write(C);
}

// get some information about the state vector
Vh [iu1,iu2,iu3,iu4] = [0,1,2,3];
Vh [u1x,u2x,u3x,u4x] = [x,x,x,x];
Vh [u1y,u2y,u3y,u4y] = [y,y,y,y];
Vh [u1b,u2b,u3b,u4b] = [U,V,W,PrMean];

{ofstream fid("LIN/dofs.dat");
	for (int j=0; j<iu1[].n; j++)
		fid << iu1[][j] << " "<< u1x[][j] << " " << u1y[][j] << endl;
}

{ofstream fid("LIN/base.dat");
	for (int j=0; j<iu1[].n; j++)
		fid << u1b[][j] << endl;
}

```

I use Freefem to generate the matrices A and B and then I read the matrices with python and solve the eigenvalue problem with slepc5py that calls slepc. Would you want also this python script ?  
I thank you once again for your help.

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [August 1, 2023, 9:39am UTC](https://community.freefem.org/t/boundary-conditions-in-eigenvalue-problem/2629/8 "2023-08-01T09:39:16Z")

</div>

Please edit that and place the code inside

```auto
a `Preformatted text` block

```

so that I can copy/paste.

---

<div class="post-metadata">

**Author:** ![maximefio](https://avatars.discourse-cdn.com/v4/letter/m/a6a055/32.png) [@maximefio](https://community.freefem.org/u/maximefio)\
**Post date:** [August 1, 2023, 9:41am UTC](https://community.freefem.org/t/boundary-conditions-in-eigenvalue-problem/2629/9 "2023-08-01T09:41:05Z")

</div>

Yes, my apologies, I didn’t know that. This is done normally.

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [August 1, 2023, 9:45am UTC](https://community.freefem.org/t/boundary-conditions-in-eigenvalue-problem/2629/10 "2023-08-01T09:45:52Z")

</div>

It looks like you are imposing boundary conditions incorrectly. Your Dirichlet BC’s should depend on `m` and you should not enforce any BC on the pressure since you use a Taylor-Hood type FE space.

---

<div class="post-metadata">

**Author:** ![maximefio](https://avatars.discourse-cdn.com/v4/letter/m/a6a055/32.png) [@maximefio](https://community.freefem.org/u/maximefio)\
**Post date:** [August 1, 2023, 11:19am UTC](https://community.freefem.org/t/boundary-conditions-in-eigenvalue-problem/2629/11 "2023-08-01T11:19:08Z")

</div>

Ok, thank you then maybe I understood badly what you meant. I Thought that the two conditions `u+mv=0` and `2u'+mv'=0` necessarily implied `u'=v'=0`, which in the finite element approach are imposed by default. But maybe I also understood something wrong for the pressure, do I need to add something for the pressure disturbance in terms of boundary condition ? (Sorry if these are obvious questions).

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [August 1, 2023, 12:04pm UTC](https://community.freefem.org/t/boundary-conditions-in-eigenvalue-problem/2629/12 "2023-08-01T12:04:53Z")

</div>

No you do not need to impose anything special for the pressure, simply eliminate the pressure from the `on()`. I would suggest something like the following:

```auto
Xh W=1-y^2, V=0, U=0;
real Re = 2200.0, alpha=1.0; // these dont need to be defined on an FE space
int m = 1; //m is an integer (makes logical statements easier)

```

then

```auto
varf opLHS([F, G, H, Pr], [X1, X2, X3, X4])
  = int1d(Th)(
    ...
  )
  + on((abs(m) != 1)*1, F = 0, G = 0)
  + on((abs(m) > 0)*1, H = 0)
  + on(2, F=0, G=0, H=0);

```

and, of course,

```auto
varf opRHS([F, G, H, Pr], [X1, X2, X3, X4])
  = int1d(Th)(
    ...
  )
  + on((abs(m) != 1)*1, F = 0, G = 0)
  + on((abs(m) > 0)*1, H = 0)
  + on(2, F=0, G=0, H=0);

```

Also, the implemented system of equations you’ve shared is not the same as equations 7-10 above, as you must have derived the corresponding weak form. I would make sure you are certain that you didn’t mix anything up during that derivation, as it is a bit nuanced in cylindrical coordinates. I haven’t caught or checked for any errors, just thought I’d mention that. You should derive the weak form for the general Navier Stokes equations, and then apply the parallel flow assumptions to the resulting weak form.

---

<div class="post-metadata">

**Author:** ![maximefio](https://avatars.discourse-cdn.com/v4/letter/m/a6a055/32.png) [@maximefio](https://community.freefem.org/u/maximefio)\
**Post date:** [August 2, 2023, 5:35am UTC](https://community.freefem.org/t/boundary-conditions-in-eigenvalue-problem/2629/13 "2023-08-02T05:35:54Z")

</div>

hello Chris,

I thank you very much for the code, this is indeed working now, I recover the proper real and imaginary part for the most unstable mode as the one from Khorrami.

I have just one question, for the mode m=±1, normally along the axis we should impose w=0,p=0, u+mv=0 and 2u’+mv’=0, or in the code that you provided, the only condition imposed is on u=0 and v=0, could you explain me why you are doing so ?

Also, regarding the equation (7)-(10), to make them in the weak form I have transform all laplacian(perturbation) x func\_test by -grad(perturbation) x grad(func\_test) and grad(perturbation) x func\_test by -perturbation x grad(func\_test), is it what should be done ? And finally does it makes a difference to do this process of setting it in weak form the final form of the equations (i.e after linearization of the Navier-Stokes equations + perturbation decomposed as normal modes), or should this process be done right after the linearization ?

I thank you once again for your answers, it is of great help for me.

Maxime

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [August 2, 2023, 8:10am UTC](https://community.freefem.org/t/boundary-conditions-in-eigenvalue-problem/2629/14 "2023-08-02T08:10:46Z")

</div>

Hi Maxime,

The pressure BC in Khorrami’s earlier work is artificial. You can read more about that [here](https://doi.org/10.1002/fld.1650120903).

It sounds like you have done it correctly if your results match. One thing to be wary about when integrating by parts is that the complex inner product involves a conjugation, so the Fourier derivatives of the test function use the conjugate of the wavenumber. In other words, the sign of `m` and `alpha` in the test function must be flipped. This should have no effect on flows that are in the O(2) symmetry group because the spectrum is symmetric with +/- `m`, but if you break the O(2) symmetry e.g. via swirl this is no longer the case.

Best,  
Chris

---

<div class="post-metadata">

**Author:** ![Amitkaush](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/amitkaush/32/1978_2.png) [@Amitkaush](https://community.freefem.org/u/Amitkaush)\
**Post date:** [February 9, 2024, 4:16pm UTC](https://community.freefem.org/t/boundary-conditions-in-eigenvalue-problem/2629/15 "2024-02-09T16:16:01Z")

</div>

> [@cmd](#):
>
> ```auto
> 
> ```
> 
> and, of course,
> 
> ```auto
> 
> ```

hi Maxime,  
I am also working on hydrodynamic stability in cartesian coordinates. Can you share the python code to calculate eigenvector and eigenvalues.

Thanks

---

<div class="post-metadata">

**Author:** ![Amitkaush](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/amitkaush/32/1978_2.png) [@Amitkaush](https://community.freefem.org/u/Amitkaush)\
**Post date:** [February 10, 2024, 10:53am UTC](https://community.freefem.org/t/boundary-conditions-in-eigenvalue-problem/2629/16 "2024-02-10T10:53:32Z")

</div>

hi Maxime,  
I am also working on hydrodynamic stability in cartesian coordinates. Can you share the python code to calculate eigenvector and eigenvalues.

Thanks