Can you send us the code: to [email protected] or [email protected] ? Or something that reproduces the problem?
Barry > On Oct 20, 2014, at 3:31 AM, Filippo Leonardi > <[email protected]> wrote: > > 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
