[petsc-users] Computing residual norm in KSPFGMRESCycle()

[Dave Lee]

Thanks Barry, good to know that MatMult() is necessary to get a fresh
residual estimation.

The only thing I don't want to "see" all of the equations is the residual
computation. Replacing VecNorm() with a custom version seems to be fine
everywhere else, but the UpdateHessenberg() call is problematic. I will
have a think about Matt's idea of customising the convergence test.

Thanks again, Dave.

On Fri, Jul 5, 2019 at 3:24 PM Smith, Barry F. <bsmith at mcs.anl.gov> wrote:

>
>
> > On Jul 4, 2019, at 11:22 PM, Dave Lee <davelee2804 at gmail.com> wrote:
> >
> > Hi Matt and Barry,
> >
> > thanks for the good ideas.
> >
> > I wasn't aware of the function KSPSetConvergenceTest(), thanks for
> alerting me to this Matt. Unfortunately I really do need to exclude a
> couple of specific degrees of freedom from the norm, and as such I want to
> make sure I can replicate the norm as computed within UpdateHessenberg()
> before I replace this.
> >
> > Apologies Barry, there was a stupid bug in that code, it should have
> read:
> >     KSPFGMRESUpdateHessenberg(ksp,loc_it,hapend,&res_norm);
> >     KSPFGMRESBuildSoln(fgmres->nrs,ksp->vec_sol,tmp_sol,ksp,loc_it+1);
> >     VecAXPY(tmp_sol,-1.0,ksp->vec_rhs);
> >     VecNorm(tmp_sol,NORM_2,&res_tmp);
>
>   This won't work.
> >
> > Which does not give the same result as the norm for UpdateHessenberg().
> Unfortunately I want to avoid using
> > KSPBuildResidual(ksp,NULL,NULL,&resid);
> > as this calls KSP_MatMult(), which will in turn call my snes residual
> assembly routine. This routine in turn calls a Navier Stokes solver to run
> for a fixed (long) time period, and so is expensive.
>
>    If you want a "freshly" computed residual you have no choice but to
> call MatMult()
> >
> > Replacing VecNorm() with my own version is indeed what I'm doing
> elsewhere, I'm just not sure of how to replace the UpdateHessenberg() norm
> right now.
>
>    It sounds like you don't want FGMRES to "see" some of your equations at
> all? You could do this by taking the result of your true matrix-vector
> product and zeroing the values for the equation you don't want FGMRES to
> "see". This will introduce a null space for your linear operator but that
> may not matter. Of course since FGMRES doesn't "see" these components it
> will make no effort to make them small so the residual for those equations
> will be whatever which means that the solution to your linear system is not
> the solution to the original system.
>
>   Barry
>
>
> >
> > Cheers, Dave.
> >
> > On Fri, Jul 5, 2019 at 12:59 PM Smith, Barry F. <bsmith at mcs.anl.gov>
> wrote:
> >
> >
> > > On Jul 4, 2019, at 7:29 PM, Dave Lee via petsc-users <
> petsc-users at mcs.anl.gov> wrote:
> > >
> > > Hi PETSc,
> > >
> > > I have a problem for which I need to exclude certain degrees of
> freedom from the evaluation of the norm as used to monitor the residual in
> FGMRES (don't ask!). Therefore I would like to replace the existing norm
> evaluation with my own version.
> > >
> > > Currently this is done within:
> > > KSPFGMRESUpdateHessenberg(ksp,loc_it,hapend,&res_norm)
> > > as:
> > > *res = PetscAbsScalar(*RS(it+1));
> > >
> > > Firstly, I would like to make sure I can replicate the existing
> residual norm. I tried to do this as:
> > > KSPFGMRESUpdateHessenberg(ksp,loc_it,hapend,&res_norm);
> > > KSPFGMRESBuildSoln(fgmres->nrs,ksp >vec_sol,tmp_sol,ksp,loc_it);
> >
> >    I don't understand the lines of code below. I don't see how they
> could compute the norm of the residual.
> >
> > From KSPMonitorTrueResidualNorm() you can compute the unpreconditioned
> residual with
> >
> >   ierr = KSPBuildResidual(ksp,NULL,NULL,&resid);CHKERRQ(ierr);
> >   ierr = VecNorm(resid,NORM_2,&truenorm);CHKERRQ(ierr);
> >   ierr = VecDestroy(&resid);CHKERRQ(ierr);
> >
> > since FGMRES always uses right preconditioning (and hence the
> unpreconditioned residual) this should match the value computed internally
> by FGMRES (subject to roundoff differences)
> >
> > To exclude certain values in the residual from the norm calculation you
> could write a custom "VecNorm" norm that skipped those entries.
> >
> >   Barry
> >
> >
> >
> >
> >
> > > VecNorm(tmp_sol,NORM_2,&res_tmp);
> > > VecNorm(ksp->vec_rhs,NORM_2,&res_rhs);
> > > res_norm = fabs(res_tmp-res_rhs);
> > >
> > > But this doesn't match the norm that comes out of UpdateHessenberg.
> > >
> > > Any ideas on how I can re-create the norm derived from
> UpdateHessenberg as an expansion over the Kyrlov vectors in FGMRES?
> > >
> > > Cheers, Dave.
> >
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20190705/b0c9d6c8/attachment.html>