[petsc-users] Multigrid as a preconditioner

[Matthew Knepley]

On Fri, Feb 17, 2012 at 6:38 AM, <coco at dmi.unict.it> wrote:

> Thank you very much for the answer, but some other doubts remain.
>
>  Date: Thu, 16 Feb 2012 13:39:11 -0600
>> From: Matthew Knepley <knepley at gmail.com>
>> Subject: Re: [petsc-users] Multigrid as a preconditioner
>> To: PETSc users list <petsc-users at mcs.anl.gov>
>> Message-ID:
>>        <CAMYG4Gmt_mn_**hPVKOrbC452geco3Zz26SF05pGy6t7**
>> RPAq+0Zg at mail.gmail.com<CAMYG4Gmt_mn_hPVKOrbC452geco3Zz26SF05pGy6t7RPAq%2B0Zg at mail.gmail.com>
>> >
>> Content-Type: text/plain; charset="iso-8859-1"
>>
>> On Thu, Feb 16, 2012 at 12:54 PM, <coco at dmi.unict.it> wrote:
>>
>>  Dear list,
>>>
>>> I would like to parallelize a multigrid code by Petsc. I do not want to
>>> use the DMMG infrastructure, since it will be replaced in the next PETSc
>>> release. Therefore I preferred to use the multigrid as a preconditioner.
>>> In
>>> practice, I use the Richardson iteration, choosing the same matrix of the
>>> linear system as a preconditioner, so that I think the Richardson
>>> iteration
>>> should converge in only one iteration, and effectively it is like solving
>>> the whole linear system by the multigrid.
>>>
>>>
>> Your understanding of the Richardson iteration is flawed. You can consult
>> Yousef Saad's book for the standard definition and anaysis.
>>
>>
> I think my explanation was not so clear. What I would like to do is to use
> a preconditioned Richardson iteration:
>
> x^(n+1) = x^n + P^(-1) (f-A x^n)
>
> Choosing P=A, I should expect to obtain the exact solution at the first
> iteration. Then, the whole linear system is solved by the preconditioner
> method that I chose. Is it what Petsc would do?
>

I am not sure what you mean by "Is it what Petsc would do?". PETSc does
what you tell it to do. If you want it
to solve in one iteration, tell it to use LU, -ksp_type richardson -pc_type
lu.


>
>>  As a first test, I tried to use a one-grid multigrid (then not a truly
>>> multigrid). I just set the coarse solver (by a standard KSPSolve), and it
>>> should be enough, because the multigrid starts already from the coarsest
>>> grid and then it does not need the smoother and the transfer operators.
>>> Unfortunately, the iteration scheme (which should be the Richardson
>>> scheme) converges with a reason=4 (KSP_CONVERGED_ITS) to a wrong
>>> solution.
>>> On the other hand, if I solve the whole problem with the standard
>>> KSPSolve
>>> (then withouth setting the multigrid as a preconditioner ...), it
>>> converges
>>> to the right solution with a reason=2.
>>>
>>>
>> Yes, give Richardson many more iterations, -ksp_max_it.
>>
>>   Matt
>>
>>
> I tried, but unfortunately nothing changed.
> Another strange phenomen is that, even with the standard KSP solver (which
> provides the right solution), if I use the flag -ksp_monitor, nothing is
> displayed. Is that an indicator of some troubles?
>

1) You must call KSPSetFromOptions() to use command line options

2) Run with -ksp_monitor -ksp_view and send the results, or there is no way
we can know what is happening.

   Matt


> Thank you in advance.
> Armando
>
>
>>  I thought that the two methods should be the exactly the same method, and
>>> I do not understand why they provide different convergence results.
>>>
>>> Here is the relevant code:
>>>
>>> // Set the matrix of the linear system
>>> Mat Mcc;
>>> ierr=MatCreate(PETSC_COMM_****WORLD,&Mcc); CHKERRQ(ierr);
>>> ierr=MatSetType(Mcc, MATMPIAIJ); CHKERRQ(ierr);
>>> ierr=MatSetSizes(Mcc,PETSC_****DECIDE,PETSC_DECIDE,1000,1000)****;
>>> CHKERRQ(ierr);
>>> ierr=setMatrix(Mcc); //It is a routine that set the values of the matrix
>>> Mcc
>>>
>>> // Set the ksp solver with the multigrid as a preconditioner
>>>  KSP ksp, KspSolver;
>>>  ierr = KSPCreate(PETSC_COMM_WORLD,&****ksp);CHKERRQ(ierr);
>>>  ierr = KSPSetType(ksp,KSPRICHARDSON);
>>>  ierr = KSPGetPC(ksp,&pc);CHKERRQ(****ierr);
>>>  ierr = PCSetType(pc,PCMG);CHKERRQ(****ierr);
>>>  ierr = PCMGSetLevels(pc,1,&PETSC_****COMM_WORLD);CHKERRQ(ierr);
>>>  ierr = PCMGSetType(pc,PC_MG_****MULTIPLICATIVE);CHKERRQ(ierr);
>>>  ierr = PCMGGetCoarseSolve(pc,&****kspCoarseSolve);CHKERRQ(ierr);
>>> ierr = KSPSetOperators(****kspCoarseSolve,Mcc,Mcc,**
>>> DIFFERENT_NONZERO_PATTERN);****CHKERRQ(ierr);
>>> ierr = KSPSetTolerances(****kspCoarseSolve,1.e-12,PETSC_**
>>> DEFAULT,PETSC_DEFAULT,PETSC_****DEFAULT);CHKERRQ(ierr);
>>> ierr = KSPSetOperators(ksp,Mcc,Mcc,****DIFFERENT_NONZERO_PATTERN);**
>>> CHKERRQ(ierr);
>>> ierr = KSPSetTolerances(ksp,1.e-12,****PETSC_DEFAULT,PETSC_DEFAULT,**
>>> PETSC_DEFAULT);CHKERRQ(ierr);
>>> ierr = KSPSolve(ksp,RHS,U);CHKERRQ(****ierr);
>>>
>>> // Solve with the standard KSPSolve
>>> KSP ksp1;
>>> ierr = KSPCreate(PETSC_COMM_WORLD,&****ksp1);CHKERRQ(ierr);
>>> ierr = KSPSetOperators(ksp1,Mcc,Mcc,****DIFFERENT_NONZERO_PATTERN);**
>>> CHKERRQ(ierr);
>>> ierr = KSPSetTolerances(ksp1,1.e-12/(****2*nn123),PETSC_DEFAULT,**
>>> PETSC_**
>>> DEFAULT,PETSC_DEFAULT);****CHKERRQ(ierr);
>>> ierr = KSPSolve(ksp1,RHS,U1);CHKERRQ(****ierr);
>>>
>>>
>>> At the end, the Vector U and U1 are different.
>>> Thank you.
>>>
>>> Best regards,
>>> Armando
>>>
>>>
>>>
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120217/d8b2bead/attachment.htm>