[petsc-users] Multiple solves with PCMG fail
Filippo Leonardi
filippo.leonardi at sam.math.ethz.ch
Wed Oct 22 07:36:26 CDT 2014
Granted that I solved my problem by **not** solving the system if the Rhs is
zero. Here's some result:
1) Using Multigrid with debug crashes (see below)
2) Using MG without debug doesn't crash, but gives wrong solution
3) defaults or gamg or wathever*works both with and without debug
Now weirder:
4) changing line 145 from
for(int i = 1; i < N; ++i) {
to
for(int i = 0; i < N; ++i) {
i.e. avoiding broken case 1) above (rhs never 0) makes 1), 2) and 3) work
perfectly (and any combination thereof)
(You can look at the differences in the pictures given by Draw to see that the
solution is more or less wrong.)
Details:
Release build with no debug:
mpirun -np 4 ./test1 -pc_type mg -pc_mg_levels 5 -draw_pause -1 -ksp_monitor -
ksp_monitor_true_residual
0 KSP Residual norm 0.000000000000e+00
0 KSP preconditioned resid norm 0.000000000000e+00 true resid norm
0.000000000000e+00 ||r(i)||/||b|| -nan
0 KSP Residual norm 3.964678029172e+02
0 KSP preconditioned resid norm 3.964678029172e+02 true resid norm
1.118033988750e+00 ||r(i)||/||b|| 1.000000000000e+00
1 KSP Residual norm 6.290716828453e-13
1 KSP preconditioned resid norm 6.290716828453e-13 true resid norm
3.788861141557e+00 ||r(i)||/||b|| 3.388860427931e+00
0 KSP Residual norm 7.929356058344e+02
0 KSP preconditioned resid norm 7.929356058344e+02 true resid norm
2.236067977500e+00 ||r(i)||/||b|| 1.000000000000e+00
1 KSP Residual norm 1.258143365691e-12
1 KSP preconditioned resid norm 1.258143365691e-12 true resid norm
7.577722283114e+00 ||r(i)||/||b|| 3.388860427931e+00
i.e. I don't see the residual to be zero each time.
Debug build:
mpirun -np 4 ./test1 -pc_type mg -pc_mg_levels 5 -draw_pause -1 -ksp_monitor -
ksp_monitor_true_residual
[1]PETSC ERROR: [0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Null argument, when expecting valid pointer
[0]PETSC ERROR: Null Pointer: Parameter # 3
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for
trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.5.0, Jun, 30, 2014
[0]PETSC ERROR: ./test1 on a arch-linux2-cxx-debug named Besikovitch-II by
filippo Wed Oct 22 14:20:46 2014
[0]PETSC ERROR: Configure options --with-clanguage=cxx --download-scalapack --
download-mumps --download-parmetis --download-metis --download-fftw
[0]PETSC ERROR: #1 KSPComputeEigenvalues() line 119 in
/home/filippo/Workspace/petsc-3.5.0/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #2 KSPChebyshevComputeExtremeEigenvalues_Private() line 328 in
/home/filippo/Workspace/petsc-3.5.0/src/ksp/ksp/impls/cheby/cheby.c
[0]PETSC ERROR: #3 KSPSolve_Chebyshev() line 373 in
/home/filippo/Workspace/petsc-3.5.0/src/ksp/ksp/impls/cheby/cheby.c
[0]PETSC ERROR: #4 KSPSolve() line 459 in
/home/filippo/Workspace/petsc-3.5.0/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #5 PCMGMCycle_Private() line 19 in
/home/filippo/Workspace/petsc-3.5.0/src/ksp/pc/impls/mg/mg.c
Debug build again with defaults:
mpirun -np 4 ./test1 -draw_pause -1 -ksp_monitor -ksp_monitor_true_residual
0 KSP Residual norm 0.000000000000e+00
0 KSP preconditioned resid norm 0.000000000000e+00 true resid norm
0.000000000000e+00 ||r(i)||/||b|| -nan
0 KSP Residual norm 1.839181443725e+00
0 KSP preconditioned resid norm 1.839181443725e+00 true resid norm
1.118033988750e+00 ||r(i)||/||b|| 1.000000000000e+00
1 KSP Residual norm 1.767908128952e+00
1 KSP preconditioned resid norm 1.767908128952e+00 true resid norm
1.182712814433e+00 ||r(i)||/||b|| 1.057850500373e+00
[...]
82 KSP Residual norm 1.642117972681e-05
82 KSP preconditioned resid norm 1.642117972681e-05 true resid norm
1.556660502826e-05 ||r(i)||/||b|| 1.392319480883e-05
[...]
0 KSP Residual norm 3.678362887450e+00
0 KSP preconditioned resid norm 3.678362887450e+00 true resid norm
2.236067977500e+00 ||r(i)||/||b|| 1.000000000000e+00
1 KSP Residual norm 3.535816257904e+00
1 KSP preconditioned resid norm 3.535816257904e+00 true resid norm
2.365425628867e+00 ||r(i)||/||b|| 1.057850500373e+00
On Wednesday 22 October 2014 07:12:50 you wrote:
> > On Oct 22, 2014, at 1:09 AM, Leonardi Filippo
> > <filippo.leonardi at sam.math.ethz.ch> wrote:
> >
> > For 1): I must be missing something then: in main i set and scale ctx.rhs,
> > then I call KSPSolve that calls ComputeRhs that further scales ctx.rhs
> > and copy it to b, which is then used by the solve. To my understanding I
> > never override ctx.rhs.
> Oh yes. I missed that.
>
> What I observed was when I ran your code each solve had an initial
> residual norm of exactly zero. I’m not sure why that happens.
> > For 2): that was a small typo, problem is present no matter what number of
> > levels I use.
> >
> > Thanks,
> > Filippo
> >
> >> On 22-ott-2014, at 04:27, "Barry Smith" <bsmith at mcs.anl.gov> 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/tutori
> >>>>> als/
> >>>>> ex29.c.html
> >>>>>
> >>>>> Best,
> >>>>> Filippo
> >>>
> >>> <test1.cpp>
--
Filippo Leonardi
SAM ETH Zürich
(HG J 45) Rämistrasse 101, Zürich (CH)
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ETHZ.vcf
Type: text/vcard
Size: 594 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20141022/aaf09087/attachment.bin>
More information about the petsc-users
mailing list