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

[Smith, Barry F.]

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