# Get the weak form help

**URL:** <https://community.freefem.org/t/get-the-weak-form-help/2660>\
**Category:** General Discussion\
**Created:** [August 17, 2023, 7:21pm UTC](https://community.freefem.org/t/get-the-weak-form-help/2660 "2023-08-17T19:21:19Z")\
**Posts on this page:** 6\
**Page:** 2

<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 22, 2023, 11:17am UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/21 "2023-08-22T11:17:53Z")

</div>

Yes, this was only a test. A better solution would be to just overwrite the singular terms using the `tgv=-1` trick. You will also need to copy the BC from `LNSE` into the varf for `b` to get this to work properly.

```auto
matrix<complex> OP = LNSE(XXXXXFES,XXXXXFES, tgv = -1, solver=sparsesolver);
matrix<complex> B = b(XXXXXFES,XXXXXFES, tgv = -10, solver=CG);

```

I also wonder what the `alpha` in your equations represents. If it is the azimuthal wavenumber it should be set to zero since your formulation is written only for the axisymmetric case.

---

<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 22, 2023, 12:17pm UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/22 "2023-08-22T12:17:08Z")

</div>

If you need to do non-axisymmetric calculations and the coordinate singularity turns out to be an issue, you might consider using the approach taken by [Meliga et al. (Phys. Fluids, 2011),](https://hal.science/hal-01067912v1/document) which uses a premultiplication by the radial coordinate to eliminate the singularity. (Note that another factor of r is introduced through the volume integral in cylindrical coordinates). You can find a variant of the resulting weak form [here](https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/nonlinear-dynamics-of-fully-developed-swirling-jets/FE28A6F2254BAE7A97058457F4C9E3E5#eqn4). For simplicity, you can just neglect the boundary integrals in your case.

---

<div class="post-metadata">

**Author:** ![Robert](https://avatars.discourse-cdn.com/v4/letter/r/e95f7d/32.png) [@Robert](https://community.freefem.org/u/Robert)\
**Post date:** [August 23, 2023, 10:48am UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/23 "2023-08-23T10:48:01Z")

</div>

Hello, sorry for the late response,

I do have some papers, here [Paper\_1](https://polytechnique.hal.science/hal-01002591/file/13259357.pdf) and [Paper\_2](https://web.stanford.edu/group/ctr/Summer/SP16/Reports/III_Turbulence_and_Transition_Physics/64_Rigas.pdf).  
So firstly I did a 3D simulation in another program, I performed a time and azimuthally average on the flow components and interpolated the results on this 2D axisymmetric geometry. Now I have to do the linear global stability analysis in 2D.  
Indeed alpha is the azimuthal wave number, I set it to 1 as it seems to be the same in these papers.

Regards,  
Robert

---

<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 23, 2023, 1:38pm UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/24 "2023-08-23T13:38:22Z")

</div>

In that case, you will have to include all of the non-axisymmetric terms in the governing equations.

---

<div class="post-metadata">

**Author:** ![Robert](https://avatars.discourse-cdn.com/v4/letter/r/e95f7d/32.png) [@Robert](https://community.freefem.org/u/Robert)\
**Post date:** [September 2, 2023, 12:04pm UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/25 "2023-09-02T12:04:23Z")

</div>

Thank you, Chris

After studying more the problem I reached to this form of the code:

```auto
load "UMFPACK64"
include "macro_ddm.idp"		
load "gmsh"					
load "iovtk"

mesh Th=gmshload("Test.msh");
plot(Th, wait=true);

fespace XFES(Th, P1);
fespace XMFES(Th, P2);
fespace XXXXXFES(Th,[P2,P2,P2,P1]); 	

		XFES Axial, Radial, Tangential, Pressure;
		XXXXXFES[u,v,w,p];
		XXXXXFES<complex>[uy,ux,uz,up]; //Eigenmodes: uaxial, uradial, utheta
		XXXXXFES<complex>[hy,hx,hz,hp]; //Test functions, axial, radial, tangential directions

		XFES EMur, EMui;
func E = [P2, P2, P2, P1];				

// uy = Axial Velocity;
// ux = Radial velocity;
// uz = Tangential Velocity;
// up = Static Pressure;

{
ifstream file("DAxial.txt");
file >> Axial[];
}
{
ifstream file("Radial.txt");
file >> Radial[];
}
{
ifstream file("Tangential.txt");
file >> Tangential[];
}
{
ifstream file("Pressure.txt");
file >> Pressure[];
}

// Combine the files above in one finite element space - XXXXXFES instead of using XFES for each of them
[u,v,w,p] = [Axial, Radial, Tangential, Pressure]; //Time and azimuthal averaged flow components
plot(u, wait=true);
plot(v, wait=true);
plot(w, wait=true);

// plot(uy, wait=true);
// plot(ux, wait=true);
// plot(uz, wait=true);
// plot(up, wait=true);

real Re = 500; // Reynolds number                       
int m = 1; // Wave number
real shiftr = 0; // Real part of the shift value
real shifti = 0.2; // Imaginary part of the shift value
real nu=1/Re;
complex shift = shiftr+shifti*1i; // Shift value (eigenvalue)     
real eps=1e-12; // Convergence criteria 
int nev=5; // Number of wanted eigenvalues
int ncv=10*nev; // Number of eigenvectors
// real epss = 1e-12;					   
int[int] Order = [1,1];
string DataName = "DEMu DEMv ";
macro JJJ x // multiply by "r" which is "x"

macro div(im,u) (dy(u#y)+dx(u#x)+u#x/x+im/x*u#z)// macro for divergence 
macro Grad(im,u) [
					[dy(u#y), dx(u#y), im/x*u#y], 
					[dy(u#x), dx(u#x), im/x*u#x-u#z/x],
					[dy(u#z), dx(u#z), im/x*u#z+u#x/x]
				 ] // macro for vecocity gradient tensor
macro D(im,u) [	
				[dy(u#y), .5*(dy(u#x)+dx(u#y)), .5*(im/x*u#y+dy(u#z))], 
				[.5*(dy(u#x)+dx(u#y)), dx(u#x), .5*(im/x*u#x-u#z/x+dy(u#z))],
				[.5*(im/x*u#y+dy(u#z)), .5*(im/x*u#x-u#z/x+dx(u#z)), im/x*u#z+u#x/x]
				] // macro for rate-of-deformation tensor

varf LNSE ([uy,ux,uz,up],[hy,hx,hz,hp])								
=int2d(Th)
		  (
		  (
		  -(v*dx(ux)+ux*dx(v)+u*dy(ux)+uy*dy(v)-2*w*uz/(x+(x==0))+w*ux*1i*m/(x+(x==0)))*hx
		  -(v*dx(uz)+ux*dx(w)+u*dy(uz)+uy*dy(w)+v*uz/(x+(x==0))+ux*w/(x+(x==0))+w*uz*1i*m/(x+(x==0)))*hz	  
		  -(v*dx(uy)+ux*dx(u)+u*dy(uy)+uy*dy(u)+w*uy*1i*m/(x+(x==0)))*hy	
		  +up*div(-1i*m,h)
		  +div(1i*m,u)*hp
		  -2*nu*(D(1i*m,u):D(-1i*m,h))
		  -shift*(uy*hy+ux*hx+uz*hz)
		  +eps*up*hp
		  )*JJJ
		  )	  
		+on(7,12,11, uy=0., ux=0., uz=0.) //7, 12, 11 -- inlet and walls
		+on(10, uy=0., pe=0.);			           
			
varf b([uy,ux,uz,up],[hy,hx,hz,hp])
=int2d(Th)((uy*hy+ux*hx+uz*hz)*JJJ);

matrix<complex> OP = LNSE(XXXXXFES,XXXXXFES, solver=sparsesolver);
matrix<complex> B = b(XXXXXFES,XXXXXFES, solver=CG);
complex[int] EigenVal(nev); 		
XXXXXFES<complex>[int] [EigenVecu,EigenVecv,EigenVecw,EigenVecp](nev);

int k;
int i;
k = EigenValue(OP,B,sigma=shift,value=EigenVal,nev=nev,vector=EigenVecu,ncv=ncv,tol=eps);
// PLOTS//
for (i=0; i<nev; i++)
{
 EMur=real(EigenVecu[i]); //Real part of the streamwise component direct eigenmode
 EMui=imag(EigenVecu[i]); //Imaginary part of the streamwise component direct eigenmode
savevtk("EigenModes"+i+".vtu", Th, [EMur,EMui], dataname = DataName, order = Order);
}

```

However, I see in this example [here](https://gitlab.com/stabfem/StabFem/-/blob/master/SOURCES_FREEFEM/StabAxi.edp?ref_type=heads) a boundary condition written as an integral on line on the symmetry axis, which is labeled with “6” in the code (see lines 186-187).

\*I changed the wave number name from _alpha_ to _m_.  
Normally, in such cases I should impose on:  
Inlet and walls: velocity components of the eigenmodes equal to 0;  
Symmetry axis: axial velocity component and pressure equal to 0;

Is there some boundary condition that I am missing?

Best regards,  
Rober

---

<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:** [September 2, 2023, 12:59pm UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/26 "2023-09-02T12:59:20Z")

</div>

I’m not sure that I understand your question. The integral in that code is basically equivalent to `on()` with `tgv=1e30`. Also, you should definitely not be imposing any boundary conditions on the pressure.

[Previous page](https://community.freefem.org/t/get-the-weak-form-help/2660.md?page=1)