[petsc-users] How to choose mat_mffd_err in JFNK

[Smith, Barry F.]

> 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