Resolvent operator with FF/PETSc (MatMatSolve?)

General Discussion
21

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

Thanks very much for your quick reply!

22

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

23

See here. Note that we never compute the resolvent operator explicitly. We only compute its action on a vector.

24

resol.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.

          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.

25

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.