{
> // Direct eigenmodes plot//
> EMur=real(EigenVecu[i]); //Real part of the streamwise component direct eigenmode
> savevtk(“…Path…\EMur.vtu”, Th, EMur, dataname=DataName, order=Order);
>
> //Below plot works://
> // plotMPI(Th, EMur, P2, def1, real, cmm = “Leading StreamwiseVelocity (Real part) eigenmode”)
> }
When I open EMur.vtu in paraview I get the results from the screenshot.
---
### Author:  [@prj](https://community.freefem.org/u/prj)
#### Post date: [April 18, 2021, 8:41pm UTC](https://community.freefem.org/t/slepc-complex-eigenmodes/913/4 "2021-04-18T20:41:11Z")
Do not load `iovtk`, it’s for sequential post-processing (hence the single subdomain view). For parallel post-processing, everything is already bundled in `SLEPc-complex`.
---
### Author:  [@Robert](https://community.freefem.org/u/Robert)
#### Post date: [April 18, 2021, 8:58pm UTC](https://community.freefem.org/t/slepc-complex-eigenmodes/913/5 "2021-04-18T20:58:18Z")
Thank you, I have the full geometry now. But even in ParaVIEW I still have the partitioning lines as you can see. Is there any chance of getting rid of them?
---
### Author:  [@prj](https://community.freefem.org/u/prj)
#### Post date: [April 18, 2021, 8:59pm UTC](https://community.freefem.org/t/slepc-complex-eigenmodes/913/6 "2021-04-18T20:59:20Z")
Sure, you can add a filter of type “Threshold”, according to “Label”, and only keep positive or null values.
---
### Author:  [@Robert](https://community.freefem.org/u/Robert)
#### Post date: [April 20, 2021, 6:35pm UTC](https://community.freefem.org/t/slepc-complex-eigenmodes/913/7 "2021-04-20T18:35:21Z")
Hello again. The “Threshold” filter worked fine. Now I am facing another small problem. I want to apply a dot product between the eigenmodes plotted above and their adjoint. It should like this:
> for (i=0; i\ {
> // Sensitivity = (sqrt(abs(EigenVecu[i])^2+abs(EigenVecv[i])^2))\*(sqrt(abs(AdjEigenVecu[i])^2+abs(AdjEigenVecv[i])^2))/(int2d(Th)(EigenVecu[i]\*AdjEigenVecu[i]+EigenVecv[i]\*AdjEigenVecv[i]));
> Sensitivity = (sqrt(abs(EigenVecu[i])^2+abs(EigenVecv[i])^2))\*(sqrt(abs(AdjEigenVecu[i])^2+abs(AdjEigenVecv[i])^2));
> Sr=real(Sensitivity); //Real part of the streamwise component direct eigenmode
> Si=imag(Sensitivity);
> savevtk(“C:\Users\Results\Resultsensitivitat.vtu”, Th, [Sr, Si], dataname=DataName, order=Order);
> }
Normally, the first formula for Sensitivity is used, the one commented which has the integral on surface (Th) in its equation, but when I use it while running in parallel on 2+ cores I get results like the ones attached in the first figure, it looks like the results/the partitions are overlapped.
However, if I still use SLEPc library but with 1 core, -np 1 it works, I get the expected spatial structure but again, as I said in earlier post, the intensity seems to be scaled (I still didn’t find out if SLEPc is scaling the results in respect to something). Removing the integral when running on several cores partially solved the problem (I get the expected structure) but with very low values. Is there another way I should write the integral while running in parallel? Thank you.
---
### Author:  [@prj](https://community.freefem.org/u/prj)
#### Post date: [April 20, 2021, 7:19pm UTC](https://community.freefem.org/t/slepc-complex-eigenmodes/913/8 "2021-04-20T19:19:41Z")
The way you compute the integral in parallel is wrong. `int2d(Th)(EigenVecu[i]*AdjEigenVecu[i]+EigenVecv[i]*AdjEigenVecv[i])` is defined on each process, and they all have different values, since all `Th` are different. You can check this by simply doing a `cout` on all MPI processes. You need to do a global reduction, as well as scale the vectors so that the overlapping unknowns are counted with multiplicity.
```auto
complex wrong = int2d(Th)(EigenVecu[i]*AdjEigenVecu[i]+EigenVecv[i]*AdjEigenVecv[i]);
cout << mpirank << " " << wrong << endl;
Vh scaledVecu = EigenVecu[i];
for[i, v : the name of the Mat you use in EPSSolve.D] scaledVecu[][i] *= v;
Vh scaledVecv = EigenVecv[i][];
for[i, v : the name of the Mat you use in EPSSolve.D] scaledVecv[][i] *= v;
wrong = int2d(Th)(scaledVecu*AdjEigenVecu[i]+scaledVecv*AdjEigenVecv[i]);
complex reduce;
mpiAllReduce(wrong, reduce, mpiCommWorld, mpiSUM);
if(mpirank == 0)
cout << reduce << endl;
```
I used [this](https://github.com/FreeFem/FreeFem-sources/blob/master/examples/hpddm/laplace-2d-SLEPc-complex.edp#L156) FreeFEM example for debugging, you can copy/paste this at the end to see it working, it may then be easier for you to adapt.
```auto
diff --git a/examples/hpddm/laplace-2d-SLEPc-complex.edp b/examples/hpddm/laplace-2d-SLEPc-complex.edp
index 719d306b..77a12657 100644
--- a/examples/hpddm/laplace-2d-SLEPc-complex.edp
+++ b/examples/hpddm/laplace-2d-SLEPc-complex.edp
@@ -154,4 +154,13 @@ for(int i=0;i
### Author:  [@Robert](https://community.freefem.org/u/Robert)
#### Post date: [April 20, 2021, 7:58pm UTC](https://community.freefem.org/t/slepc-complex-eigenmodes/913/9 "2021-04-20T19:58:43Z")
Thank you, I will try this method.
---
### Author:  [@Robert](https://community.freefem.org/u/Robert)
#### Post date: [April 24, 2021, 1:09pm UTC](https://community.freefem.org/t/slepc-complex-eigenmodes/913/10 "2021-04-24T13:09:21Z")
Hello again. I tried to apply the following debug (Where Xh is my finite element space and J is the name given to Mat) and the error I get is the first photo attached:
> complex wrong = int2d(Th)(EigenVecu[i]\*AdjEigenVecu[i]+EigenVecv[i]\*AdjEigenVecv[i]);
> cout \<\< mpirank \<\< " " \<\< wrong \<\< endl;
> Xh scaledVecu = EigenVecu[i]; //EigenVecu[i] - same error
> for[i, v : J] scaledVecu[i] _= v;
> Xh scaledVecv = EigenVecv[i][];
> for[i, v : J] scaledVecv[][i] = v;
> wrong = int2d(Th)(scaledVecuAdjEigenVecu[i]+scaledVecv_AdjEigenVecv[i]);
> complex reduce;
> mpiAllReduce(wrong, reduce, mpiCommWorld, mpiSUM);
> if(mpirank == 0)
> cout \<\< reduce \<\< endl;
Another question I have is about ParaVIEW post-processing, as you can see in the second attached photo. I tried several filters (point data to cell data/point interpolate/delaunay) to interpolate/smooth the results. The results are good if I process them in MatLAB or directly from FreeFEM++, so I don’t think the mesh is a big problem. Is there another filter to apply, or a filter applied to another filter? The above filters were applied directly on the .pvd file.
Thank you very much for your help and patience
Regards,
Robert
---
### Author:  [@prj](https://community.freefem.org/u/prj)
#### Post date: [April 24, 2021, 1:32pm UTC](https://community.freefem.org/t/slepc-complex-eigenmodes/913/11 "2021-04-24T13:32:20Z")
No idea what `J` is. As you can see in my snippet above, you need the equivalent of `DistA.D`, where `DistA` is the `Mat` you pass to `EPSSolve`. What `order` are you using in `savevtk`?
---
### Author:  [@Robert](https://community.freefem.org/u/Robert)
#### Post date: [April 24, 2021, 1:44pm UTC](https://community.freefem.org/t/slepc-complex-eigenmodes/913/12 "2021-04-24T13:44:28Z")
I will check again the matrix. In savetk I use order 1, but everything is calculated with order 2. I attach a small part of the code which implies creating the finite element spaces and export data. If I change Order=[2,2…] I only get the Labels in the .pvd file.
> fespace Xh(Th,P2);
> fespace Mh(Th,P1);
> Xh EMur, EMui, EMvr, EMvi;
> Xh AdjEMur, AdjEMui, AdjEMvr, AdjEMvi;
> Mh EMpr, EMpi;
> Mh AdjEMpr, AdjEMpi;
>
> ```
> Xh Sensitivity;
> Xh Sr, Si;
>
> ```
>
> int[int] Order = [1,1,1,1,1,1,1,1,1];
> string DataName= "Velocity Pressure EigenmodesAxial EigenmodesRadial EigenmodesPressure " +
> "AdjointEigenmodesAxial AdjointEigenmodesRadial AdjointEigenmodesPressure " +
> “Sensitivity”;
> for (i=0; i\ {
>
> EMur=real(EigenVecu[i]);
> EMui=imag(EigenVecu[i]);
> EMvr=real(EigenVecv[i]);
> EMvi=imag(EigenVecv[i]);
> EMpr=real(EigenVecp[i]);
> EMpi=imag(EigenVecp[i]);
> AdjEMur=real(AdjEigenVecu[i]);
> AdjEMui=imag(AdjEigenVecu[i]);
> AdjEMvr=real(AdjEigenVecv[i]);
> AdjEMvi=imag(AdjEigenVecv[i]);
> AdjEMpr=real(AdjEigenVecp[i]);
> AdjEMpi=imag(AdjEigenVecp[i]);
> Sensitivity=…
> Sr=real(Sensitivity);
> Si=imag(Sensitivity);
> savevtk(“C:\Users\Desktop\Results\Reynolds”+Re+“\Results.vtu”, Th,
> [u, v], p,
> [EMur, EMui], [EMvr, EMvi], [EMpr, EMpi],
> [AdjEMur, AdjEMui], [AdjEMvr, AdjEMvi], [AdjEMpr, AdjEMpi],
> [Sr, Si],
> dataname=DataName, order=Order);
> }
---
### Author:  [@Robert](https://community.freefem.org/u/Robert)
#### Post date: [May 2, 2021, 5:00pm UTC](https://community.freefem.org/t/slepc-complex-eigenmodes/913/13 "2021-05-02T17:00:15Z")
Dear professor,
While I was trying to fix int2d(Th)(EigenVecu[i]\*AdjEigenVecu[i]+EigenVecv[i]\*AdjEigenVecv[i]) using the [Laplace 2d example](https://github.com/FreeFem/FreeFem-sources/blob/master/examples/hpddm/laplace-2d-SLEPc-complex.edp) I remarked the difference in writing the distributed matrices compared to [Navier-Stokes 2d SLEPc complex](https://github.com/FreeFem/FreeFem-sources/blob/master/examples/hpddm/navier-stokes-2d-SLEPc-complex.edp). I applied both ways of writing to my case which is also a flow around a cylinder. The code look like this (I didn’t write the equations as they are quite long):
> load “SLEPc-complex”
> include “macro\_ddm.idp”
> mesh Th =…
> func E = [P2, P2, P2, P1];
> fespace XXXXFES(Th,E);
>
> int[int][int] intersection; // local-to-neighbors renumbering
> real[int] D; // partition of unity
>
> //Problem definition
> varf LNSE =…
> varf B =…
>
> //Matrix creation
>
> //VARIANT I //
>
> Mat J;
> macro def(i) [i, i#A, i#B, i#C]
> macro init(i) [i, i, i, i]
> createMat(Th, J, E);
> J = LNSE(XXXXFES, XXXXFES);
> matrix B = b(XXXXFES, XXXXFES);
> Mat M(J, B);
>
> int k = EPSSolve(J, M, vectors = …, values = …, sparams = …);
>
> //VARIANT II //
>
> matrix J = LNSE(XXXXFES, XXXXFES);
> matrix B = b(XXXXFES, XXXXFES);
>
> Mat DistJ(J,intersection, D);
> Mat DistB(DistJ, B);
>
> int k = EPSSolve(DistJ, DistB, vectors = …, values = …, sparams = …);
Which is the difference between VARIANT I and VARIANT II? I am asking because the eigenvalues are the same in both cases but the eigenvectors are different. In the documentary of SLEPc it is said that the eigenvectors are normalized (I guess normalized by sqrt(a^2+b^2…) where a, b are the components of the vector). Now that I have 2 different values for the eigenvectors I am a bit confused. I looked into “macro\_ddm.idp” but I didn’t find where the eigenvectors are normalized. My questions would be:
1. What is exactly “intersection” and “D” doing to the “MatDistJ”?
2. Is it possible to eliminate the normalization of the eigenvectors (maybe from macro\_ddm.idp)?
3. Is there something I could improve in writing the files for ParaVIEW post-processing in the code posted in the last post?
Thank you,
Robert
---
### Author:  [@prj](https://community.freefem.org/u/prj)
#### Post date: [May 2, 2021, 5:28pm UTC](https://community.freefem.org/t/slepc-complex-eigenmodes/913/14 "2021-05-02T17:28:54Z")
1. Do not use variant II, it’s an old API, variant I is much easier to follow (you don’t need to guess what are those variables).
2. Do not touch `macro_ddm.idp`, it has nothing to do with the eigensolver. You can look at the SLEPc [manual](https://slepc.upv.es/documentation/slepc.pdf) to understand how the vectors are normalized and how to adjust this (I’m not sure this is possible).
3. If what you write is what you want, I guess there is no real point in trying to improve that part of the script, is there?
---
### Author:  [@Robert](https://community.freefem.org/u/Robert)
#### Post date: [May 5, 2021, 11:07am UTC](https://community.freefem.org/t/slepc-complex-eigenmodes/913/15 "2021-05-05T11:07:27Z")
Thank you for the answer. I will keep going with Variant I. About 2), I read in the manual how the vectors are normalized, I wanted to see how the norm is written in the program (especially the complex conjugate) and maybe do some simple tests. About 3) I tought there is a better way to export the results to avoid that problem in plotting, but I will just try some filter combinations in ParaVIEW.
Coming back to int2d(Th)(EigenVecu[i]\*AdjEigenVecu[i]+EigenVecv[i]\*AdjEigenVecv[i]), I attach a very simplified code with a simplified integral:
> load “SLEPc-complex”
> include “macro\_ddm.idp”
>
> mesh Th=…
>
> func E = [P2, P2, P2, P1];
>
> fespace XXXXFES(Th,E);
> fespace XFES(Th,P2);
>
> varf LNSE=…
> varf b=…
>
> Mat M, AdjM;
> macro def(j) [j, j#A, j#B, j#C]
> macro init(j) [j, j, j, j]
> createMat(Th, M, E);
>
> M = LNSE(XXXXFES, XXXXFES);
> matrix B = b(XXXXFES, XXXXFES);
>
> Mat J(M, B);
>
> int nev = 1; //nb. of eigenvalues
> int ncv = 2\*nev; //nb. of eigenvectors
> real[int] D;
>
> complex[int] EigenVal(nev); //store eigenvalues
> XXXXFES[int] EigenVecu,EigenVecv,EigenVecw,EigenVecp; //store eigenvectors
>
> int i, k;
> k = EPSSolve(M, J, values=EigenVal, vectors=EigenVecu, sparams = …);
>
> for (i=0; i\ {
>
> complex wrong = int2d(Th)(EigenVecu[i]+EigenVecv[i]);
> cout \<\< mpirank \<\< " " \<\< wrong \<\< endl;
>
> Xh scaledVecu = EigenVecu[i];
> for[i, su : M.D] scaledVecu \*= su;
>
> Xh scaledVecv = EigenVecv[i];
> for[i, su : M.D] scaledVecv \*= su;
>
> complex wrong2 = int2d(Th)(scaledVecu\*scaledVecv);
> complex reduce;
> mpiAllReduce(wrong, reduce, mpiCommWorld, mpiSUM);
> if(mpirank == 0)
> cout \<\< reduce \<\< endl;
>
> cout\<\<"wrong2 = " \<\< wrong2 \<\
> }
As you can see in the picture attached, indeed I have different values for the integral because all “Th” are different but after the global reduction the integral (redefined with wrong2) is 0 and I am unable to find the reason. I changed from your example “v” as it is defined as radial velocity in my case to “su” and scaledVecu[i] \*= v to scaledVecu \*=su due to an array error. Let me know if I should add some more information to the script.
Regards,
Robert
---
### Author:  [@prj](https://community.freefem.org/u/prj)
#### Post date: [May 5, 2021, 11:22am UTC](https://community.freefem.org/t/slepc-complex-eigenmodes/913/16 "2021-05-05T11:22:39Z")
`int2d(Th)(scaledVecu*scaledVecv)` is wrong, if you want to compute `u * v`, you need to scale either `u` or `v`, not both.
---
### Author:  [@Robert](https://community.freefem.org/u/Robert)
#### Post date: [May 5, 2021, 11:40am UTC](https://community.freefem.org/t/slepc-complex-eigenmodes/913/17 "2021-05-05T11:40:33Z")
I did the following change but I still get "wrong 2 = (0,0):
> for (i=0; i\ {
>
> complex wrong = int2d(Th)(EigenVecu[i]+EigenVecv[i]);
> cout \<\< mpirank \<\< " " \<\< wrong \<\< endl;
>
> Xh scaledVecu = EigenVecu[i];
> for[i, su : M.D] scaledVecu \*= su;
>
> //Xh scaledVecv = EigenVecv[i];
> //for[i, su : M.D] scaledVecv \*= su;
>
> complex wrong2 = int2d(Th)(scaledVecu\*EigenVecv[i]);
> complex reduce;
> mpiAllReduce(wrong, reduce, mpiCommWorld, mpiSUM);
> if(mpirank == 0)
> cout \<\< reduce \<\< endl;
>
> cout\<\<"wrong2 = " \<\< wrong2 \<\
> }
---
### Author:  [@prj](https://community.freefem.org/u/prj)
#### Post date: [May 5, 2021, 11:44am UTC](https://community.freefem.org/t/slepc-complex-eigenmodes/913/18 "2021-05-05T11:44:15Z")
Well, you didn’t used my code properly… `scaledVecu[]` should be `scaledVecu[][i]`. Just use my script, which is giving correct results, and adapt it to your needs if you don’t know how to debug your code.
---
### Author:  [@Robert](https://community.freefem.org/u/Robert)
#### Post date: [May 5, 2021, 12:05pm UTC](https://community.freefem.org/t/slepc-complex-eigenmodes/913/19 "2021-05-05T12:05:36Z")
With \* `scaledVecu[i] .\* I get the attached error. This is how I tried first time:
> for (i=0; i\ {
>
> complex wrong = int2d(Th)(EigenVecu[i]+EigenVecv[i]);
> cout \<\< mpirank \<\< " " \<\< wrong \<\< endl;
>
> Xh \< complex \> scaledVecu = EigenVecu[i]; //also tried =EigenVecu[i]
> for[i, su : M.D] scaledVecu[i] \*= su;
>
> complex wrong2 = int2d(Th)(scaledVecu\*EigenVecv[i]);
> complex reduce;
> mpiAllReduce(wrong, reduce, mpiCommWorld, mpiSUM);
> if(mpirank == 0)
> cout \<\< reduce \<\< endl;
>
> cout\<\<"wrong2 = " \<\< wrong2 \<\
> }
Line 236 being: for[i, su : M.D] scaledVecu[i] \*= su;
---
### Author:  [@prj](https://community.freefem.org/u/prj)
#### Post date: [May 5, 2021, 12:14pm UTC](https://community.freefem.org/t/slepc-complex-eigenmodes/913/20 "2021-05-05T12:14:39Z")
You are not using the proper vector then.
I’ll repeat myself.
> Just use my script, which is giving correct results, and adapt it to your needs if you don’t know how to debug your code.
[Next page](https://community.freefem.org/t/slepc-complex-eigenmodes/913.md?page=2)