[petsc-dev] [petsc-users] Multigrid with defect correction

[Barry Smith]

   I've looked at the code again and I'm afraid that just "adding a second DM" is not trivial. The DMKSP/SNES/TS stuff is all built around a single DM and has stuff like ->dmoriginal in it (do we have two dmoriginal with all its complications?). 

   Perhaps the entire KSPSetComputeOperators() model was misguided and I should have stuck with the "you fill up all the needed information for the levels of PCMG" before you do the multigrid setup instead of providing callbacks that can fill it them up on the fly "as needed".  We could possibly throwout all the coarsen/refine hook stuff and the DMKSP construct. Essentially we'd only be requiring SNESComputeJacobian to know about the PCMG API to provide the coarser grid operators. The coarsen/refine hook stuff seems to be there only to hide from the SNESComputeJacobian the PCMG API. At lot of complication for what real benefit besides our programming egos?

   If we threw out all this complicated stuff the user level API wouldn't change much (almost no uses directly KSPSetComputeOperators())

   Barry



> On Feb 25, 2017, at 6:26 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> 
>> 
>> On Feb 25, 2017, at 2:21 PM, Matt Landreman <matt.landreman at gmail.com> wrote:
>> 
>> Thanks everyone for the help.
>> 
>> On item 1 (using Amat different from Pmat with geometric multigrid), I tried Barry's suggestion but it did not seem to resolve the issue. For example, in ksp ex25.c, I tried adding the following lines after line 112:
>>  if (J == jac) {                                                                                                       
>>    ierr = PetscPrintf(PETSC_COMM_WORLD,"Creating a new Amat\n");                                                           
>>    ierr = DMCreateMatrix(da,&J);                                                                                          
>>    ierr = KSPSetOperators(ksp,J,jac);                                                                                      
>>  }                                                   
>>  ierr = MatShift(J,1.0); 
>> 
>> This change should set Amat (J) to be different from Pmat (jac), (specifically Amat=identity matrix), so the solution from KSP should be completely different from the original ex25.c. But viewing the solution vector, the solution is unchanged. It seems PETSc is ignoring the Amat created in this approach.
> 
>   Ahh, yes a multiple number of errors on my part.
> 
> 1) there was old code that immediately overwrote the newly set KSPSetOperators() after the call to ComputeMatrix(), hence simply removing the new user provided matrix, fixed in commit d73a0ed5842d2f684a9729e0b9c14ba0e3d554c3
> 
> 2) the default behavior of PCMG was to use the pmat on each level to compute the residual, not the mat, fixed in commit 13842ffb25245dcab28ec15282db994aa7b56c77 so even after 1) was fixed it would still be using "the wrong" matrix to compute the residuals.
> 
> 3) Even after 2) was fixed my advice still would not do what I wanted because the default setting of the matrix to compute the residual on each level (the call the PCMGSetResidual() in mg.c) occurs BEFORE any calls to ComputeMatrix() so it would still use the original matrix provided by the DM and not the new one provided by the user. 
> 
>  The best solution I see moving forward is to actually properly support providing two DMs to KSP, one responsible for the "true" operator and one responsible for the "approximate" operator used in constructing the preconditioner. Conceptually looking back this makes perfect sense and should have been done the very first day that I introduced the concept of providing a DM to the KSP. Unfortunately this requires a good amount of replumbing; wherever we do refinement/coarsening of DM we need to manage it for both, the horrible DMKSP/SNES/TS construct that Jed needed to introduce to get user callbacks working needs to handle the two DMs in some way, plus TS and SNES need to handle the two DMs. There are additional complications in managing the user callbacks for function evaluations, when they should utilize the first DM and when they should utilize the second; for Newton based TS and SNES solves presumably they should almost always be just using the higher order function, but for nonlinear multigrid and other non-Newton based nonlinear solvers there are presumably places where the lower order function should be called. The end result of the refactorization is very good however since it will then be straightforward to support many cases of the "high order discretization for operators, but low order discretization for preconditioners" paradigms, not just for DMDA but for other DMS without the usual PETSc model and code.
> 
>  Matt,
> 
>    So yes we definitely want to support what you want to do and have always intended to support it eventually, now is the time we start to add the support. I'll start a git branch shortly that begins to make it simple for your case of KSP; extending it everywhere will take more time.
> 
>  Barry
> 
> 
> 
> 
> 
> 
>> 
>> Matt K's suggestion of switching from KSP to SNES does work, allowing Amat to differ from Pmat on the finest multigrid level. (On coarser levels, it seems PETSc still forces Amat=Pmat on entry to computeMatrix). 
>> 
>> On Jed's comment, the application I have in mind is indeed a convection-dominated equation (a steady linear 3D convection-diffusion equation with smoothly varying anisotropic coefficients and recirculating convection). Gamg and hypre-boomerAMG have been working on it ok if I discretize with low-order upwind differences in Pmat and use Amat=Pmat, but I'd like higher order accuracy. Using gmres with a higher-order discretization in Amat and low-order Pmat works ok, but the number of KSP iterations required gets large as the diffusion gets small compared to convection, even with -pc_type lu. So I'm working to see if geometric mg with defect correction at each level can do better.
>> 
>> Thanks,
>> Matt Landreman
>> 
>> 
>> 
>> 
>> On Thu, Feb 23, 2017 at 5:23 PM, Jed Brown <jed at jedbrown.org> wrote:
>> Matt Landreman <matt.landreman at gmail.com> writes:
>>> 3. Is it at all sensible to do this second kind of defect correction with
>>> _algebraic_ multigrid? Perhaps Amat for each level could be formed from the
>>> high-order matrix at the fine level by the Galerkin operator R A P, after
>>> getting all the restriction matrices created by gamg for Pmat?
>> 
>> Note that defect correction is most commonly used for
>> transport-dominated processes for which the high order discretization is
>> not h-elliptic.  AMG heuristics are typically really bad for such
>> problems so stabilizing a smoother isn't really a relevant issue.  Also,
>> it is usually used for strongly nonlinear problems where AMG's setup
>> costs are likely overkill.  This is not to say that defect correction
>> AMG won't work, but there is reason to believe that it doesn't matter
>> for many important problems.