[petsc-users] How to choose mat_mffd_err in JFNK

[Smith, Barry F.]

> On Aug 14, 2019, at 9:41 PM, Tang, Qi <tangqi at msu.edu> wrote:
> 
> Thanks for the help, Barry. I tired both ds and wp, and again it depends on if I could find the correct parameter set. It is getting harder as I refine the mesh.
> 
> So I try to use SNESDefaultMatrixFreeCreate2, SNESMatrixFreeMult2_Private or SNESDiffParameterCompute_More in mfem. But it looks like these functions are not in linked in petscsnes.h. How could I call them?

   They may not be listed in petscsnes.h but I think they should be in the library.  You can just stick the prototypes for the functions anywhere you need them for now.


 You should be able to use

  ierr = PetscOptionsInt("-snes_mf_version","Matrix-Free routines version 1 or 2","None",snes->mf_version,&snes->mf_version,0);CHKERRQ(ierr);

then this info is passed with 

  if (snes->mf) {
    ierr = SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version);CHKERRQ(ierr);
  }

this routine has

  if (version == 1) {
    ierr = MatCreateSNESMF(snes,&J);CHKERRQ(ierr);
    ierr = MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);CHKERRQ(ierr);
    ierr = MatSetFromOptions(J);CHKERRQ(ierr);
  } else if (version == 2) {
    if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"SNESSetFunction() must be called first");
#if !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) && !defined(PETSC_USE_REAL___FP16)
    ierr = SNESDefaultMatrixFreeCreate2(snes,snes->vec_func,&J);CHKERRQ(ierr);
#else

and this routine has

  ierr = VecDuplicate(x,&mfctx->w);CHKERRQ(ierr);
  ierr = PetscObjectGetComm((PetscObject)x,&comm);CHKERRQ(ierr);
  ierr = VecGetSize(x,&n);CHKERRQ(ierr);
  ierr = VecGetLocalSize(x,&nloc);CHKERRQ(ierr);
  ierr = MatCreate(comm,J);CHKERRQ(ierr);
  ierr = MatSetSizes(*J,nloc,n,n,n);CHKERRQ(ierr);
  ierr = MatSetType(*J,MATSHELL);CHKERRQ(ierr);
  ierr = MatShellSetContext(*J,mfctx);CHKERRQ(ierr);
  ierr = MatShellSetOperation(*J,MATOP_MULT,(void (*)(void))SNESMatrixFreeMult2_Private);CHKERRQ(ierr);
  ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void (*)(void))SNESMatrixFreeDestroy2_Private);CHKERRQ(ierr);
  ierr = MatShellSetOperation(*J,MATOP_VIEW,(void (*)(void))SNESMatrixFreeView2_Private);CHKERRQ(ierr);
  ierr = MatSetUp(*J);CHKERRQ(ierr);






> 
> Also, could I call SNESDiffParameterCompute_More in snes_monitor? But it needs "a" in (F(u + ha) - F(u)) /h as the input. So I am not sure how I could properly use it. Maybe use SNESDefaultMatrixFreeCreate2 inside SNESSetJacobian would be easier to try?

If you use the flag -snes_mf_noise_file filename when it runs it will save all the noise information it computes along the way to that file (yes it is crude and doesn't match the PETSc Viewer/monitor style but it should work).

  Thus I think you can use it and get all the possible monitoring information without actually writing any code. Just

-snes_mf_version 2
-snes_mf_noise_file  filename



   Barry

> 
> Thanks again,
> Qi
> 
> From: Smith, Barry F. <bsmith at mcs.anl.gov>
> Sent: Tuesday, August 13, 2019 9:07 PM
> To: Tang, Qi <tangqi at msu.edu>
> Cc: petsc-users at mcs.anl.gov <petsc-users at mcs.anl.gov>
> Subject: Re: [petsc-users] How to choose mat_mffd_err in JFNK
>  
> 
> 
> > On Aug 13, 2019, at 7:27 PM, Tang, Qi via petsc-users <petsc-users at mcs.anl.gov> wrote:
> > 
> > Hi,
> > I am using JFNK, inexact Newton and a shell "physics-based" preconditioning to solve some multiphysics problems. I have been playing with mat_mffd_err, and it gives me some results I do not fully understand. 
> > 
> > I believe the default value of mat_mffd_err is 10^-8 for double precision, which seem too large for my problems. When I test a problem essentially still in the linear regime, I expect my converged "unpreconditioned resid norm of KSP" should be identical to "SNES Function norm" (Should I?). This is exactly what I found if I could find a good mat_mffd_err, normally between 10^-3 to 10^-5. So when it happens, the whole algorithm works as expected. When those two norms are off, the inexact Newton becomes very inefficient. For instance, it may take many ksp iterations to converge but the snes norm is only reduced slightly.
> > 
> > According to the manual, mat_mffd_err should be "square root of relative error in function evaluation". But is there a simple way to estimate it? 
> 
>   First a related note: there are two different algorithms that PETSc provides for computing the factor h using the err parameter. They can be set with -mat_mffd_type - wp or ds (see MATMFFD_WP or MATMFFD_DS) some people have better luck with one or the other for their problem. 
> 
> 
>   There is some code in PETSc to compute what is called the "noise" of the function which in theory leads to a better err value. For example if say the last four digits of your function are "noise" (that is meaningless stuff) which is possible for complicated multiphysics problems  due to round off in the function evaluations then you should use an err that is 2 digits bigger than the default (note this is just the square root of the "function epsilon" instead of the machine epsilon, because the "function epsilon" is much larger than the machine epsilon). 
> 
>    We never had a great model for how to hook the noise computation code up into SNES so it is a bit disconnected, we should revisit this. Perhaps you will find it useful and we can figure out together how to hook it in cleanly for use. 
> 
>    The code that computes and reports the noise is in the directory src/snes/interface/noise 
> 
>    You can use this version with the option -snes_mf_version 2 (version 1 is the normal behavior) 
> 
>    The code has the ability to print to a file or the screen information about the noise it is finding, maybe with command line options, you'll have to look directly at the code to see how. 
> 
>    I don't like the current setup because if you use the noise based computation of h you have to use SNESMatrixFreeMult2_Private() which is a slightly different routine for doing the actually differencing than the two MATMFFD_WP or MATMFFD_DS that are normally used. I don't know why it can't use the standard differencing routines but call the noise routine to compute h (just requires some minor code refactoring). Also I don't see an automatic way to just compute the noise of a function at a bunch of points independent of actually using the noise to compute and use h. It is sort of there in the routine SNESDiffParameterCompute_More() but that routine is not documented or user friendly. 
> 
>    I would start by rigging up calls to SNESDiffParameterCompute_More() to see what it says the noise of your function is, based on your hand determined values optimal values for err it should be roughly the square of them. Then if this makes sense you can trying using the -snes_mf_version 2  code to see if it helps the matrix free multiple behave consistently well for your problem by adaptively computing the err and using that in computing h. 
> 
>    Good luck and I'd love to hear back from you if it helps and suggestions/code for making this more integrated with SNES so it is easier to use. For example perhaps we want a function SNESComputeNoise() that uses SNESDiffParameterCompute_More() to compute and report the noise for a range of values of u. 
> 
>    Barry
> 
> 
> 
> 
> 
> > 
> > Is there anything else I could possibly tune in this context? 
> > 
> > The discretization is through mfem and I use standard H1 for my problem. 
> > 
> > Thanks,
> > Qi