[petsc-users] Multiple solves with PCMG fail

[Filippo Leonardi]

Hi,

I have a very specific problem that I cannot figure out with PCMG and multiple 
solves. I got a linear system that I solve many times, with same matrix but 
different RHS. I can successfully solve the system with standard techniques, 
such as default solver or LU or PCGAMG. Even MG works if I destroy the ksp 
each time or I recompute the matrices at each time. But when I try and go to 
MG and not recomputing the matrices each time the solver fails. Any idea?

Here some detail: the setup:
ierr = KSPSetDMActive(ksp, PETSC_FALSE); CHKERRQ(ierr);
ierr = KSPSetDM(ksp,  da); CHKERRQ(ierr);
ierr = KSPSetComputeOperators(ksp, ComputeMatrix, ctx); CHKERRQ(ierr);
ierr = KSPSetComputeRHS(ksp, ComputeRHS, ctx); CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);
    
Then I solve as usual, for a large number of time steps:
ierr = KSPSolve(ksp, NULL, NULL); CHKERRQ(ierr);
ierr = KSPGetSolution(ksp, &phi);

The solver converges and does that in a reasonable number of iterations:
Linear solve converged due to CONVERGED_RTOL iterations 7
And ksp_view and ksp_monitor do not show any weird stuff.

- Weirdly enough using any solver (for instance cg+bjacobi or gamg) everything 
works (So the matrix and RHS are working fine).
- But the problem persists with Galerkin matrices (-pc_mg_galerkin) (So is not 
a ComputeMatrix problem).
- If I do:
ierr = KSPSetComputeOperators(ksp, ComputeMatrix, this); CHKERRQ(ierr);
between each solve or destroy the ksp entirely each time the solution is also 
perfect (So is not a boundary scaling or other stuff problem).
- If I run MG with only 2 levels (so just coarse) I also get a fine result.

Setup RHS is called each time as expected, setup matrix is called just once, 
also as expected.

The only thing I can think is that MG does not update some value that actually 
needs to be recomputed.

Any idea?

The solution is not that different from:
http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex29.c.html

Best,
Filippo