Thank you for the answer. > Date: Fri, 17 Feb 2012 10:11:51 -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: > <CAMYG4GkKE6doSQ1FgHCSr1deZRYS0kpvLQF8daqv0au==tOgyw at mail.gmail.com> > Content-Type: text/plain; charset="iso-8859-1" > > 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. >
Indeed I would like to solve the whole linear system by a multigrid approach and not by a lu factorization. Therefore I would like to use -ksp_type richardson -pc_type mg. In this case, the preconditioned problem P^(-1) (f-A x^n) is solved exactly or it performs just a V-cycle iteration? In both cases, since I am using a one-grid multigrid (just for debugging), it should anyway provide the exact solution at the first iteration, but it is not so. > >> >>> 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. I plotted the residual and it decreases so much slowly. It seems like it is using the non-preconditioned Richardson iteration x^(n+1) = x^n + P^(-1) (f-A x^n) with P=I instead of P=A. Thank you. Armando > >> 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 >>>> >>>> >>>>
