[petsc-users] How to choose mat_mffd_err in JFNK
Tang, Qi
tangqi at msu.edu
Mon Aug 19 19:51:44 CDT 2019
Thanks again, Barry. I am testing more based on your suggestions.
One thing I do not understand is when I use -mat_mffd_type wp -mat_mffd_err 1e-3 -ksp_type fgmres. It computes "h" of JFNK twice sequentially in each KSP iteration. For instance, the output in info looks like
...
[0] MatMult_MFFD(): Current differencing parameter: 1.953103182078e-09
[0] MatMult_MFFD(): Current differencing parameter: 6.054449182838e-01
...
I found it called MatMult_MFFD twice in a row. Do you happen to know why petsc calls MatMult_MFFD twice when using fgmre? I also do not understand the relationship between these two h. They are very off (because apparently ||a|| is very different).
On the other hand, h is only computed once when I use gmres. The h there is clearly scaled with mat_mffd_err. That is easy to understand. However, I think I need to use fgmres because my preconditioner is not quite linear.
BTW, my snes view is:
SNES Object: () 16 MPI processes
type: newtonls
maximum iterations=20, maximum function evaluations=10000
tolerances: relative=0.0001, absolute=0., solution=0.
total number of linear solver iterations=70
total number of function evaluations=149
norm schedule ALWAYS
Eisenstat-Walker computation of KSP relative tolerance (version 3)
rtol_0=0.3, rtol_max=0.9, threshold=0.1
gamma=0.9, alpha=1.5, alpha2=1.5
SNESLineSearch Object: () 16 MPI processes
type: bt
interpolation: cubic
alpha=1.000000e-04
maxstep=1.000000e+08, minlambda=1.000000e-12
tolerances: relative=1.000000e-08, absolute=1.000000e-15, lambda=1.000000e-08
maximum iterations=40
KSP Object: () 16 MPI processes
type: fgmres
restart=50, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement
happy breakdown tolerance 1e-30
maximum iterations=10000, initial guess is zero
tolerances: relative=0.0564799, absolute=1e-50, divergence=10000.
right preconditioning
using UNPRECONDITIONED norm type for convergence test
PC Object: () 16 MPI processes
type: shell
JFNK preconditioner
No information available on the mfem::Solver
Number of preconditioners created by the factory 8
linear system matrix followed by preconditioner matrix:
Mat Object: 16 MPI processes
type: mffd
rows=787968, cols=787968
Matrix-free approximation:
err=0.001 (relative error in function evaluation)
Using wp compute h routine
Does not compute normU
Mat Object: 16 MPI processes
type: nest
rows=787968, cols=787968
Matrix object:
type=nest, rows=3, cols=3
MatNest structure:
(0,0) : type=mpiaij, rows=262656, cols=262656
(0,1) : type=mpiaij, rows=262656, cols=262656
(0,2) : NULL
(1,0) : type=mpiaij, rows=262656, cols=262656
(1,1) : type=mpiaij, rows=262656, cols=262656
(1,2) : NULL
(2,0) : type=mpiaij, rows=262656, cols=262656
(2,1) : NULL
(2,2) : type=mpiaij, rows=262656, cols=262656
________________________________
From: Smith, Barry F. <bsmith at mcs.anl.gov>
Sent: Thursday, August 15, 2019 3:30 AM
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 15, 2019, at 12:36 AM, Tang, Qi <tangqi at msu.edu> wrote:
>
> Thanks, it works. snes_mf_jorge works for me.
Great.
> It appears to compute h in every ksp.
Each matrix vector product or each KSPSolve()? From the code looks like each matrix-vector product.
>
> Without -snes_mf_jorge, it is not working. For some reason, it only computes h once, but that h is bad. My gmres residual is not decaying.
>
> Indeed, the noise in my function becomes larger when I refine the mesh. I think it makes sense as I use the same time step for different meshes (that is the goal of the preconditioning). However, even when the algorithm is working, sqrt(noise) is much less than the good mat_mffd_err I previously found (10^-6 vs 10^-3). I do not understand why.
I have no explanation. The details of the code that computes err are difficult to trace exactly; would take a while. Perhaps there is some parameter in there that is too "conservative"?
>
> Although snes_mf_jorge is working, it is very expensive, as it has to evaluate F many times when estimating h. Unfortunately, to achieve the nonlinearity, I have to assemble some operators inside my F. There seem no easy solutions.
>
> I will try to compute h multiple times without using snes_mf_jorge. But let me know if you have other suggestions. Thanks!
Yes, it seems to me computing the err less often would be a way to make the code go faster. I looked at the code more closely and noticed a couple of things.
if (ctx->jorge) {
ierr = SNESDiffParameterCompute_More(snes,ctx->data,U,a,&noise,&h);CHKERRQ(ierr);
/* Use the Brown/Saad method to compute h */
} else {
/* Compute error if desired */
ierr = SNESGetIterationNumber(snes,&iter);CHKERRQ(ierr);
if ((ctx->need_err) || ((ctx->compute_err_freq) && (ctx->compute_err_iter != iter) && (!((iter-1)%ctx->compute_err_freq)))) {
/* Use Jorge's method to compute noise */
ierr = SNESDiffParameterCompute_More(snes,ctx->data,U,a,&noise,&h);CHKERRQ(ierr);
ctx->error_rel = PetscSqrtReal(noise);
ierr = PetscInfo3(snes,"Using Jorge's noise: noise=%g, sqrt(noise)=%g, h_more=%g\n",(double)noise,(double)ctx->error_rel,(double)h);CHKERRQ(ierr);
ctx->compute_err_iter = iter;
ctx->need_err = PETSC_FALSE;
}
So if jorge is set it uses the jorge algorithm for each matrix multiple to compute a new err and h. If jorge is not set but -snes_mf_compute_err is set then it computes a new err and h depending on the parameters (you can run with -info and grep for noise to see the PetscInfo lines printed (maybe add iter to the output to see how often it is recomputed)
if ((ctx->need_err) || ((ctx->compute_err_freq) && (ctx->compute_err_iter != iter) && (!((iter-1)%ctx->compute_err_freq)))) {
so the compute_err_freq determines when the err and h are recomputed. The logic is a bit strange so I cannot decipher exactly how often it is recomputing. I'm guess at the first Newton iteration and then some number of iterations later.
You could try to rig it so that it is at every new Newton step (this would mean when computing f(x + h a) - f(x) for each new x it will recompute I think).
A more "research type" approach to try to reduce the work to a reasonable level would be to keep the same err until "things start to go bad" and then recompute it. But how to measure "when things start to go bad?" It is related to GMRES stagnating, so you could for example track the convergence of GMRES and if it is less than expected, kill the linear solve, reset the computation of err and then start the KSP at the same point as before. But this could be terribly expensive since all the steps in the most recent KSP are lost.
Another possibility is to assume that the err is valid in a certain size ball around the current x. When x + ha or the new x is outside that ball then recompute the err. But how to choose the ball size and is there a way to adjust the ball size depending on how the computations proceed. For example if everything is going great you could slowly increase the size of the ball and hope for the best; but how to detect if the ball got too big (as above but terribly expensive)?
Track the size of the err each time it is computed. If it stays about the same for a small number of times just freeze it for a while at that value? But again how to determine when it is no longer good?
Just a few wild thoughts, I really have no solid ideas on how to reduce the work requirements.
Barry
>
> Qi
>
>
> From: Smith, Barry F. <bsmith at mcs.anl.gov>
> Sent: Wednesday, August 14, 2019 10:57 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 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20190820/3ca14191/attachment-0001.html>
More information about the petsc-users
mailing list