[petsc-users] Multigrid as a preconditioner

[coco at dmi.unict.it]

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_hPVKOrbC452geco3Zz26SF05pGy6t7RPAq+0Zg 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?

>
>> 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?

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
>>
>>