# Resolvent operator with FF/PETSc (MatMatSolve?)

**URL:** <https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603>\
**Category:** General Discussion\
**Created:** [March 16, 2022, 4:05pm UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603 "2022-03-16T16:05:45Z")\
**Posts on this page:** 5\
**Page:** 2

<div class="post-metadata">

**Author:** ![clthu](https://avatars.discourse-cdn.com/v4/letter/c/f04885/32.png) [@clthu](https://community.freefem.org/u/clthu)\
**Post date:** [August 9, 2024, 3:14am UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/21 "2024-08-09T03:14:21Z")

</div>

OK, I understand. I am not so familiar with tgv before.

Thanks very much for your quick reply!

---

<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:** [August 22, 2025, 10:57am UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/22 "2025-08-22T10:57:59Z")

</div>

Hi Chris,

I went through your code. Can you pinpoint the file where you calculate the inverse of R=(iωMq+J)−1. I am not that much experienced with PETSc and I know this is too much to ask, but if you can help it will be helpful.

thanks

---

<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, 2025, 11:12am UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/23 "2025-08-22T11:12:33Z")

</div>

See [here](https://github.com/cmd8/ff-bifbox/blob/56500ca27c77b0d4b3f388603978c05dfb608d2b/rslvcompute.edp#L114C1-L121C2). Note that we never compute the resolvent operator explicitly. We only compute its action on a vector.

---

<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:** [September 15, 2025, 12:25pm UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/24 "2025-09-15T12:25:32Z")

</div>

[resol.edp](https://community.freefem.org/uploads/short-url/w6MLmEZ2KZ51jY6qjdZXlDEKjJy.edp) (6.6 KB)  
Hi Chris,  
I went through your code and adopted the portion where you indirectly form the LHS operator in my code.

I have written a code to do resolvent analyis of Blasius boundary layer (code is attached).

I am able to match the gains value with those reported in the literature. I am interested in plotting the forcing modes against the domain in MATLAB (1D in this case i.e., along the y axis). Can you suggest a way to do this, since the values of y in the Y.txt are not in ascending manner?

Here is the the snippet from my code where I am saving my forcing modes in “.txt” file.

```auto
          int k = EPSSolve(LHS, Mf, values = val, vectors =fv,
                 sparams = "-eps_type krylovschur -eps_largest_real -eps_monitor_conv -options_left no -eps_gen_hermitian");
                
  
                 
                 real[int] gains = sqrt(val.re);
                 
  
  VhC<complex>[int] defu(uvec)(val.n);
  
  ChangeNumbering(Mf, fv[0][], gm);
      MatMult(PQ, gm, qm);
      KSPSolve(J, qm, qm);
      ChangeNumbering(Mf, fv[0][], gm, inverse = true);
      ChangeNumbering(J, uvec[0][], qm, inverse = true);
  
  fespace V1(Thg,P2);
V1 Y = y;   

 
   Xhg<complex> deff(tempf), deff(fmg);
    XMhg<complex> defu(tempu), defu(umg);
    GH<complex> deff(fm);
    VhC<complex> defu(um);
      fm[] = fv[0][];
      um[] = uvec[0][];
      tempu[](restu) = um[];
      mpiAllReduce(tempu[], umg[], mpiCommWorld, mpiSUM);
      tempf[](restf) = fm[];
      mpiAllReduce(tempf[], fmg[], mpiCommWorld, mpiSUM);
      
      if(mpirank==0){
      cout<<gains[0]<<endl;
      
     
       ofstream ofile("forcing"+kz+".txt");
       ofile << fmg[] << endl;
       ofstream ofile1("response"+kz+".txt");
       ofile1 << umg[]<<endl;
       ofstream ofile2("Y.txt");
       ofile2<< Y[]<<endl;}

```

Your suggestion will be helpful. Also thank you for your code, it was helpful for writing my own code.

---

<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 15, 2025, 1:17pm UTC](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603/25 "2025-09-15T13:17:12Z")

</div>

I think there are a lot of ways to get what you want. It would be easy to get this with ff-bifbox via ParaView output. Is there a reason why you’re starting from scratch? You can also look at the file I/O examples in the FreeFEM documentation. Note that there is no need to have the values of y in ascending order, as this can be done in a post-processing step.

[Previous page](https://community.freefem.org/t/resolvent-operator-with-ff-petsc-matmatsolve/1603.md?page=1)