> 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
