Resolvent operator with FF/PETSc (MatMatSolve?)

General Discussion
1

Hello FF developers,

I am interested in implementing a resolvent (input/output) analysis framework using the FreeFEM/PETSc interface.

To do this, I need to construct a matrix that is defined by a series of products and of several sparse matrices and matrix inverses. The matrix of interest is Hermitian, with a dimension of n\times n. It is defined as L=B^HR^HM_qRB.

Here, B is a m\times n matrix of 1s and 0s with m\geq n, M_q is a positive semi-definite m\times m matrix, and R is an m\times m matrix defined by an inverse as R=(i\omega M_q+J)^{-1}.

With the exception of R, I can construct the Mat objects for each of the basic components of L, and I understand how to perform the necessary MatMatMult() operations. However, I do not know how to use MUMPS to find R. It seems like this would require the MatMatSolve(), but I couldn’t find any documentation on this in FreeFEM. Has anyone encountered a similar problem and found a solution?

For context, I need to construct this matrix in order to solve the eigenvalue problem (using the SLEPc interface):

L\tilde{\mathbf{f}}=\lambda M_f\tilde{\mathbf{f}}, where M_f is a positive definite n\times n matrix.

A good overview with more detail on this subject is given here (see Section 4).

2

You cannot access R from MUMPS, nobody can, except MUMPS developers. So you’ll need to solve (i\omega M_q+J)C = B using KSPSolve() in your .edp, by first converting B to a dense Mat since KSPMatSolve() (in PETSc library) only handles dense right-hand sides. This will be extremely costly, these inversions make your L dense, and depending on which \lambda you are looking for, SLEPc may have to “invert” L, so maybe it would be best to think of an approximation of L^{-1} which would not need an explicit representation of your L.

See FreeFem-tutorial - Section 8 - example14.edp for an example of solves with multiple right-hand sides.

1 Like
3

@prj thank you for the quick reply.

You are correct that this would be very expensive. Thanks to your comment, I think I see the solution. I’ll state it briefly here in case it helps anyone in the future.

  1. Define a func which computes the action of the operator L on \tilde{\mathbf{f}} i.e. \mathbf{x}=L\tilde{\mathbf{f}}.

  2. Use this func to create the Arnoldi basis in the EVP in a matrix-free approach (Ă  la FreeFem-tutorial - Section 8 - example11.edp). This way the LU decomposition is only factored once and applied only as many times as is needed to build the Arnoldi subspace.

4

That would indeed be orders of magnitude cheaper to compute.
This is like for PCD preconditioning, cf. FreeFem-sources/oseen-2d-PETSc.edp at develop · FreeFem/FreeFem-sources · GitHub, the PCD operator is never assembled, but is computed instead by a sequence of MatMult() or KSPSolve(). In your case, you’ll need to implement:

  • MatMult() in a func as a sequence of MatMult(B, ...) + KSPSolve(R, ...) + MatMult(Mq, ...) + KSPSolveHermitianTranspose(R, ...) + MatMultHermitianTranspose(B, ...)
  • PCApply() – for L^{-1} if you are doing shift-and-invert in SLEPc for smallest \lambda – in a func (maybe an initial approximation could be simply M_q^{-1})
5

LU decomposition of what? L? You won’t have that. R or M_q? Yes, whatever you choose as a preconditioner for these two operators, PETSc is smart enough to reuse them throughout the Arnoldi iterations.

6

Thanks, yes I was referring the the LU decomposition of R. I can see that PETSc is automatically reusing that factorization.

7

Thank you for your suggestions. I believe I am close to having this working, but I am a bit confused on the required syntax.

I have defined a func complex[int] Lop(complex[int]& inPETSC) {...} to handle the sequence of MatMult()s and KSPsolve()s that define the operator L. I have tested this function independently (i.e. outside of the EPSSolve()) and can confirm that it works as intended.

My struggle right now has to do with how to get EPSSolve() to do what I want. In particular, I don’t understand how to use the precon argument used in your examples. My current code is similar to what I’ve given below. It runs without error, but the results I’m getting don’t make sense yet. This could be due to BC’s or something else, but I wondered if you see an obvious error based on my lack of experience with PETSc.

Mat<complex> L(Mf, Lop); // Mf is a Mat<complex>, Lop is a func
int nev = getARGV("-nev",1);
string EPSparams = " -eps_nev " + nev + " " +
                   " -eps_type krylovschur " +
                   " -eps_largest_real " +
                   " -st_pc_type none " +
                   " -eps_gen_hermitian ";
complex[int] val(nev);
Xh<complex>[int] deff(vec)(nev); // Xh is a fespace
int k = EPSSolve(L, Mf, vectors = vec, values = val, sparams = EPSparams);
8

I believe you don’t need -st_pc_type none with -eps_largest_real. What do you mean by “the results […] don’t make sense”? Are you getting a meaningful -eps_view? What about -eps_monitor and -eps_converged_reason?

1 Like
9

Yes, these are all meaningful. I think the issue is with my varfs/BCs. Thank you for your helpful insights as always @prj :slight_smile:

10

Hello Chris,

could you maybe share the final code?

Regards,
Tea

11

Hi there, yes. It is available via Github. What is your username?

1 Like
12

Hi Chris,

Thanks for the response!
You mean my Github username? It’s TeaVoj. You can also share the link.

13

Hi Tea,

I will let you know when the repo goes public. For now, the relevant portion of the code is below.

Cheers!

// construct matrices
M  = vMq(XMh, XMh); // Response Norm
Mf = vMf(Xh, Xh); // Forcing Norm
matrix<complex> LocPQ = vP(Xh, XMh); // Forcing/Response Correspondence
Mat<complex> PQ(M, Mf, LocPQ);

func complex[int] LHSop(complex[int]& inPETSc) {
  complex[int] temp(XMh.ndof), outPETSc(inPETSc.n);
  MatMult(PQ, inPETSc, outPETSc);
  KSPSolve(J, outPETSc, temp);
  MatMult(M, temp, outPETSc);
  KSPSolveHermitianTranspose(J, outPETSc, temp);
  MatMultHermitianTranspose(PQ, temp, outPETSc);
  return outPETSc;
}

Mat<complex> LHS(Mf, LHSop);

J  = vJ(XMh, XMh, tgv = -1); //Linear operator
set(J, sparams = KSPparams);
int k = EPSSolve(LHS, Mf, vectors = fvec, values = val,
                 sparams = "-eps_type krylovschur -eps_largest_real -eps_monitor_conv -options_left no -eps_gen_hermitian");
1 Like
14

Hi Chris,

thanks a lot. This is very helpful.

I am still having some issues with my code, it is probably due to my inexperience with FreeFem.
Could you maybe share which KSPparams you use?

15

Glad that helps. A good rule of thumb is to always start with: "-ksp_type preonly -pc_type lu"

16

Dear Chris,

I am considering using the resolvent framework, and your code is very helpful. Could you share the variational formulation of vMq, vMf and LocPQ, and the corresponding FE spaces?

17

Hi @aszaboa and @tea, I’ve opened the ff-bifbox repo to the public now. Note that the code is still under development and some things may not work perfectly yet.

You can see/test my resolvent analysis implementation by following the commands in examples/garnaud_2012.

1 Like
18

Thank you so much for sharing the code!

1 Like