> On Feb 22, 2017, at 4:53 AM, Matt Landreman <[email protected]> wrote:
>
> Hi, PETSc developers,
>
> 1. Consider solving a linear problem with -pc_type mg, assembling the matrix
> using KSPSetComputeOperators(ksp,ComputeMatrix,&user) and MatSetValuesStencil
> as in ksp ex25.c. I’d like Pmat to be different than Amat (a stable
> low-order discretization and a high-order discretization, respectively.
Presumably the high-order discretization will be "matrix-free" and you'd
like to use a MATSHELL for it?
> ) However the two Mat arguments that PETSc passes to ComputeMatrix appear to
> be the same object, so Pmat and Amat are forced to be the same. (For example
> in ex25, J==jac). How do I allow Amat and Pmat to differ?
It looks like you are using a DMDA and having the DMDA create the matrices
for all levels instead of setting them yourself. The problem comes about
because of this code:
PetscErrorCode KSPSetUp(KSP ksp)
{
PetscErrorCode ierr;
.....
ierr = DMCreateMatrix(ksp->dm,&A);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
ierr = PetscObjectDereference((PetscObject)A);CHKERRQ(ierr);
The easiest way for you to work around this is when you enter you
ComputeMatrix() if the two matrices are the same then create your shell matrix
amat and call KSPSetOperators(ksp,amat,pmat thus "replacing the first use of
the pmat with you new shell matrix. Then just "fill up" the pmat as usual and
do what you need for you shell matrix amat.
The next time ComputeMatrix() is called for this level the two matrices will
already be different so you just set the entries in the pmat and do what you
need for the amat shell matrix.
This approach is "tacky" but allows you do exactly what you want.
Barry
>
> 2. I suppose if the above issue were resolved, the resulting gmres/mg
> iteration would amount to doing defect correction on only the finest grid.
No, I think
> Another defect correction method is to use a different Mat for smoothing and
> residual computation at every level. (Stable low-order for smoothing,
> high-order for residual.) This latter approach is recommended on page 103 of
> Brandt’s multigrid guide
> (http://www.wisdom.weizmann.ac.il/~achi/classics.pdf). What would be the best
> way to do this second kind of defect correction in PETSc mg? Perhaps loop
> over every multigrid level and call
> PCMGSetResidual(pc,l,PCMGResidualDefault,Amat) on each level, after calling
> ComputeMatrix to build the appropriately-sized Amat?
You could use PCMGSetResidual but I don't think it is necessary.
>
>
> 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?
>
> Thanks,
> Matt Landreman
