[petsc-users] Multiple solves with PCMG fail

[Filippo Leonardi]

Actually sorry, I have figured out the issues.

The first iterate in the for cycle (int i = 0) sets the vector to be zero, 
which, for what I understand from the discussion above, is not handled 
correctly by PCMG (and only him) (because I guess the answer should still be 
computable).

What I did in my real code was to solve an equation and, strangely enough, for 
some initial data the first rhs would be zero, exactly as in the code I 
provided (and I did not thing of that trough). Skipping the first solve of the 
system I get rid of the problem and things **seems** to b working correctly.

What is even weirder is that in debug mode I get the mentioned EigenValue 
problem, but in Optimized mode the code runs through everything, but at the 
end the solutions of the whole process is completely wrong.

Now things seem to be solved,
thanks.

On Tuesday 21 October 2014 21:27:06 you wrote:
>   There are a few issues with your code
> 
> 1) If you use KSPSetComputeRHS() then the function you provide is called
> immediately before the solve automatically for every solve. Thus your code
> with the vector operations here
> 
>    ierr = VecSet(ctx.rhs, (float) i * mx); CHKERRQ(ierr);
>     ierr = VecScale(ctx.rhs,  rank); CHKERRQ(ierr);
>     ierr = KSPSolve(ksp, NULL, NULL); CHKERRQ(ierr);
> does nothing. And each linear solve ends up being with a zero right hand
> side.
> 
> 2) -pc_type mg doesn’t give enough information for the multigrid to know how
> many levels to use. In fact it is always using one hence it is not
> multigrid. You  need to use -p_mg_levels 2 or 3 or whatever to get it to
> use more than one level.
> 
>   Barry
> 
> > On Oct 21, 2014, at 1:50 AM, Filippo Leonardi
> > <filippo.leonardi at sam.math.ethz.ch> wrote:
> > 
> > Here,
> > 
> > run this e.g. with
> > mpirun -np 4 ./test1 -draw_pause -1 -pc_type mg -pc_mg_galerkin
> > 
> > In Debug mode it crashes, in standard mode it fails. Using Petsc current.
> > 
> > Run without mg and works!
> > 
> > Hope it's not me being stupid and doing something terribly wrong!
> > 
> > Best,
> > Filippo
> > 
> > On Monday 20 October 2014 07:29:04 Barry Smith wrote:
> >>  Can you send us the code: to petsc-maint at mcs.anl.gov or
> >> 
> >> petsc-users at mcs.anl.gov ? Or something that reproduces the problem?
> >> 
> >>  Barry
> >>  
> >>> On Oct 20, 2014, at 3:31 AM, Filippo Leonardi
> >>> <filippo.leonardi at sam.math.ethz.ch> 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/tutorial
> >>> s/
> >>> ex29.c.html
> >>> 
> >>> Best,
> >>> Filippo
> > 
> > <test1.cpp>