### Author:  [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [August 18, 2023, 1:13pm UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/5 "2023-08-18T13:13:26Z")
One thing I can definitely tell you off the bat is that your mixed finite element space is not inf-sup stable. You have to use an inf-sup stable FEspace – something like Taylor–Hood (`[P2,P2,P2,P1]`) or MINI (`[P1b,P1b,P1b,P1]`) elements. You may also have to algebraically simplify your weak form (i.e. enforce that `x*(1/x)=1`) so that you dont divide everything by zero along the coordinate axis.
---
### Author:  [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [August 18, 2023, 1:17pm UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/6 "2023-08-18T13:17:45Z")
Also, it does not look like you integrated the viscous term by parts correctly. Note that you are integrating by parts using vectorial test and trial functions – you cannot treat the x, r, and theta momentum equations separately when integrating by parts.
---
### Author:  [@Robert](https://community.freefem.org/u/Robert)
#### Post date: [August 18, 2023, 2:40pm UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/7 "2023-08-18T14:40:58Z")
Thank you, I changed to P2 elements type and enforced the term x\*(1/x)=1, huge differences indeed. As with the mathematical part, I surely miss something. Do you refer to the viscous terms of second order partial derivative, highlighted in red?
Regards,
Robert
---
### Author:  [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [August 18, 2023, 2:59pm UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/8 "2023-08-18T14:59:28Z")
Not just the parts highlighted in red – the entire viscous term. For example, see slide 9 [here](https://www.ljll.math.upmc.fr/~frey/cours/UPMC/NM491/P.%20Frey,%20chap3.pdf). Note that in your case, the gradient and divergence operators must be the correct ones for a cylindrical coordinate system.
---
### Author:  [@Robert](https://community.freefem.org/u/Robert)
#### Post date: [August 18, 2023, 4:22pm UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/9 "2023-08-18T16:22:25Z")
Thank you, Chris for your patience. So I took separately only the viscous term and this is what I’ve got.
I guess there is still something misunderstood in the equations or a typo in the code
Regards,
Robert
---
### Author:  [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [August 18, 2023, 4:41pm UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/10 "2023-08-18T16:41:02Z")
I think the main issue is that you are still integrating by parts incorrectly. You have to integrate the VECTOR Laplacian by parts. What you are doing doesn’t make sense since you are treating the viscous term as three separate instances of a scalar Laplacian.
---
### Author:  [@Robert](https://community.freefem.org/u/Robert)
#### Post date: [August 18, 2023, 5:51pm UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/11 "2023-08-18T17:51:51Z")
Wouldn’t that be simply
where the term involving the partial derivative in respect to azimuthal direction is null?
---
### Author:  [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [August 18, 2023, 6:03pm UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/12 "2023-08-18T18:03:35Z")
No, you are integrating a vector-valued function by parts with respect to a vector-valued test function.
---
### Author:  [@Robert](https://community.freefem.org/u/Robert)
#### Post date: [August 18, 2023, 6:17pm UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/13 "2023-08-18T18:17:15Z")
Thank you for your guidance, but I still miss something to get the correct weak formulation. Is it possible to show me the solution?
Regards,
Robert
---
### Author:  [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [August 19, 2023, 9:38am UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/14 "2023-08-19T09:38:21Z")
Here is one way to write it for the axisymnmetric case. Here, x and y are along the axial and radial directions, respectively. [vax,vra,vth] is the velocity test function and [uax, ura, uth] is the velocity trial function.
```auto
1/Re*(y*(dx(vax)*dx(uax) + dy(vax)*dy(uax) + dx(vra)*dx(ura) + dy(vra)*dy(ura) + dx(vth)*dx(uth) + dy(vth)*dy(uth)) + (vth*uth + vra*ura)/y)
```
---
### Author:  [@Robert](https://community.freefem.org/u/Robert)
#### Post date: [August 19, 2023, 2:07pm UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/15 "2023-08-19T14:07:47Z")
Thank you for your time!
So, basically the problem is the following:
And I wrote the following script:
```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[ue,ve,we,pe]; //Eigenmodes: uaxial, uradial, utheta
XXXXXFES[hy,hx,hz,q]; //Test functions, axial, radial, tangential directions
XFES EMur, EMui;
func E = [P2, P2, P2, P1];
// u = Axial Velocity;
// v = Radial velocity;
// w = Tangential Velocity;
// p = Static Pressure;
{
ifstream file("Axial.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(u, wait=true);
// plot(v, wait=true);
// plot(w, wait=true);
// plot(p, wait=true);
real Re = 500; // Reynolds number
int alpha = 1; // Wave number
real shiftr = 0; // Real part of the shift value
real shifti = 0.2; // Imaginary part of the shift value
complex shift = shiftr+shifti*1i; // Shift value (eigenvalue)
real eps=1e-6; // Convergence criteria
int nev=3; // Number of wanted eigenvalues
int ncv=10*nev; // Number of eigenvectors
int[int] Order = [1,1];
string DataName = "DEMu DEMv ";
varf LNSE ([ue,ve,we,pe],[hy,hx,hz,q])
=int2d(Th)
(
-((v*dx(ve)+ve*dx(v)+u*dy(ve)+ue*dy(v))*hx*x+(w*ve*1i*alpha-2*w*we)*hx)
-((v*dx(we)+ve*dx(w)+u*dy(we)+ue*dy(w))*hz*x+(v*we+ve*w+w*we*1i*alpha)*hz)
-((v*dx(ue)+ve*dx(u)+u*dy(ue)+ue*dy(u))*hy*x+(w*ue*1i*alpha)*hy)
+pe*dx(hx)*x+pe*dx(hy)*x-pe/x*1i*alpha*hz*x
-(1/Re)*(x*(dx(ve)*dx(hx)+dy(ve)*dy(hx)+dx(we)*dx(hz)+dy(we)*dy(hz)+dx(ue)*dx(hy)+dy(ue)*dy(hy)))
-(1/Re)*(we*hz+ve*hx)/x
-shift*(ue*hy+ve*hx+we*hz)*x
)
+on(7,12,11, ue=0., ve=0., we=0.) //7, 12, 11 -- inlet and walls
+on(10, ue=0., pe=0.); //symmetry axis, pressure needed?? pe=0.
varf b([ue,ve,we,pe],[hy,hx,hz,q])
=int2d(Th)((ue*hy+ve*hx+we*hz)*x);
matrix OP = LNSE(XXXXXFES,XXXXXFES, solver=sparsesolver);
matrix B = b(XXXXXFES,XXXXXFES, solver=CG);
complex[int] EigenVal(nev);
XXXXXFES[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
### Author:  [@marchywka](https://community.freefem.org/u/marchywka)
#### Post date: [August 19, 2023, 2:55pm UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/16 "2023-08-19T14:55:32Z")
Did you put the edp in the rar file? The copy from screen seems to drop asterisks
and make other problems on my Chrome.
---
### Author:  [@Robert](https://community.freefem.org/u/Robert)
#### Post date: [August 19, 2023, 3:02pm UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/17 "2023-08-19T15:02:49Z")
Hello,
No, I attach a new .rar with the .edp file too
Regards,
Robert
[Test.rar](https://community.freefem.org/uploads/short-url/3NLSQhk4iCnPFn4Ia8yuXmXVplO.rar) (488.5 KB)
---
### Author:  [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [August 19, 2023, 9:44pm UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/18 "2023-08-19T21:44:39Z")
I wonder if there is an issue with your varf and the singularity at r=0. Try replacing the `1/y` with `1/(y+(y==0))`.
Note that one could probably find a better way to avoid dividing by zero.
---
### Author:  [@Robert](https://community.freefem.org/u/Robert)
#### Post date: [August 22, 2023, 8:36am UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/19 "2023-08-22T08:36:23Z")
Thank you very much, Chris. It seems that did the trick, it is working now. The results are however different than what I expected, I will keep looking to the equations/boundary conditions/perform some mesh sensitivity analysis to find the problem.
Regards,
Robert
---
### Author:  [@marchywka](https://community.freefem.org/u/marchywka)
#### Post date: [August 22, 2023, 10:49am UTC](https://community.freefem.org/t/get-the-weak-form-help/2660/20 "2023-08-22T10:49:27Z")
What results did you expect "? Do you have a paper to cite?
AFAICT that change was just intended to determine if the divide by
zero was causing the immediate problem. It may require some reformulation
to get a useful result. On a quick look of the literature this seems
to be a bit of an open topic. Apparently there are some frequency
domain approaches. These may help with some singularities
but are still difficult with non-linear terms.
[Next page](https://community.freefem.org/t/get-the-weak-form-help/2660.md?page=2)