[petsc-users] Questions on matrix-free GMRES implementation

[feng wang]

Dear All,

I am new to petsc and trying to implement a matrix-free GMRES. I have assembled an approximate Jacobian matrix just for preconditioning. After reading some previous questions on this topic, my approach is:

the matrix-free matrix is created as:

      ierr = MatCreateMFFD(*A_COMM_WORLD, iqe*blocksize, iqe*blocksize, PETSC_DETERMINE, PETSC_DETERMINE, &petsc_A_mf); CHKERRQ(ierr);
      ierr = MatMFFDSetFunction(petsc_A_mf, FormFunction_mf, this); CHKERRQ(ierr);

KSP linear operator is set up as:

      ierr = KSPSetOperators(petsc_ksp, petsc_A_mf, petsc_A_pre); CHKERRQ(ierr); //petsc_A_pre is my assembled pre-conditioning matrix

Before calling KSPSolve, I do:

         ierr = MatMFFDSetBase(petsc_A_mf, petsc_csv, petsc_baserhs); CHKERRQ(ierr); //petsc_csv is the flow states, petsc_baserhs is the pre-computed right hand side

The call back function is defined as:

   PetscErrorCode cFdDomain::FormFunction_mf(void *ctx, Vec in_vec, Vec out_vec)
  {
      PetscErrorCode ierr;
      cFdDomain *user_ctx;

      cout << "FormFunction_mf called\n";

      //in_vec: flow states
      //out_vec: right hand side + diagonal contributions from CFL number

      user_ctx = (cFdDomain*)ctx;

      //get perturbed conservative variables from petsc
      user_ctx->petsc_getcsv(in_vec);

      //get new right side
      user_ctx->petsc_fd_rhs();

      //set new right hand side to the output vector
      user_ctx->petsc_setrhs(out_vec);

      ierr = 0;
      return ierr;
  }

The linear system I am solving is (J+D)x=RHS. J is the Jacobian matrix.  D is a diagonal matrix and it is used to stabilise the solution at the start but reduced gradually when the solution moves on to recover Newton's method. I add D*x to the true right side when non-linear function is computed to work out finite difference Jacobian, so when finite difference is used, it actually computes (J+D)*dx.

The code runs but diverges in the end. If I don't do matrix-free and use my approximate Jacobian matrix, GMRES  works. So something is wrong with my matrix-free implementation. Have I missed something in my implementation? Besides,  is there a way to check if the finite difference Jacobian matrix is computed correctly in a matrix-free implementation?

Thanks for your help in advance.
Feng
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210311/57896ac3/attachment-0001.html>