Thanks, it works. snes_mf_jorge works for me. It appears to compute h in every ksp.
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. 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! Qi ________________________________ From: Smith, Barry F. <[email protected]> Sent: Wednesday, August 14, 2019 10:57 PM To: Tang, Qi <[email protected]> Cc: [email protected] <[email protected]> Subject: Re: [petsc-users] How to choose mat_mffd_err in JFNK > On Aug 14, 2019, at 9:41 PM, Tang, Qi <[email protected]> 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. <[email protected]> > Sent: Tuesday, August 13, 2019 9:07 PM > To: Tang, Qi <[email protected]> > Cc: [email protected] <[email protected]> > 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 > > <[email protected]> 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
