[petsc-users] SNESSetFunction and MatMFFDSetFunction
Barry Smith
bsmith at mcs.anl.gov
Fri Jan 17 15:47:51 CST 2014
On Jan 16, 2014, at 10:49 AM, Song Gao <song.gao2 at mail.mcgill.ca> wrote:
> I was looking at the example of MatMFFDSetFunction on website.
> http://www.mcs.anl.gov/petsc/petsc-current/src/snes/examples/tutorials/ex22.c.html
> I think, the line 312, the last snes should be ctx.
> 312: MatMFFDSetFunction(*A,(PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction
> ,snes);
No the code as in the example is correct. The first argument to SNESComputeFunction() is the SNES object http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/SNESComputeFunction.html#SNESComputeFunction
>
>
>
>
> On Thu, Jan 16, 2014 at 11:04 AM, Song Gao <song.gao2 at mail.mcgill.ca> wrote:
> Thank you.
>
> I have found the bug in my codes. SNESSetFunction and MatMFFDSetFunction expect functions with different interfaces, but I passed the same function to them.
>
>
> On Wed, Jan 8, 2014 at 11:52 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
> I suspect the problem is here:
>
> >
> > call MatMFFDSetBase(myJctx%mf, pet_solu_snes, PETSC_NULL_INTEGER,
> >
> > @ ierrpet)
>
> In fact I am surprised it didn’t crash at this line, since we don’t have code to handle the PETSC_NULL_INTEGER
>
> Try adding the
>
> > call SNESGetFunction(snes,f,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierrpet);
>
> before the
>
> > MatMFFDSetBase(myJctx%mf,x,f,ierrpet);
>
> that I suggested before.
>
> Does that change anything?
>
> If you still get different values here is how I would debug it next
>
> 1) 1 process
>
> 2) run each version separately in the debugger (you can use the options -start_in_debugger noxterm )
>
> 3) put a break point in MatMult(). In most debuggers you just use
>
> b MatMult
>
> 4) then type c to continue
>
> 5) when it stops in MatMult do
>
> VecView(x,0)
>
> 6) make sure both versions produce the exact same numbers (do they?)
>
> 7) then type next several times until it gets to the PetscFunctionReturn(0) line
>
> 8) do
>
> VecView(y,0)
>
> again are both answers identical? By all logic they will be different since the norms you print are different.
>
> If (6) produces the same numbers but (8) produces different ones then put a break point in MatMult_MFFD()
> and call VecView() on ctx->current_u ctx->current_f and a. For both versions they should be the same. Are they?
>
> Barry
>
>
>
> On Jan 8, 2014, at 10:47 AM, Song Gao <song.gao2 at mail.mcgill.ca> wrote:
>
> > Dear Barry,
> >
> > Thanks for reply. I basically implemented your codes. Then I have two
> > questions.
> >
> > The first is I'm working on Fortran. So I can't use MatShellSetContext to
> > set the structure. Therefore I let the variable I want to set, MyJctx, to
> > be global. Is there other way to do that?
> >
> > The second question is I did some tests. let the D vec to be zero, I
> > expect the code which I explicit set the matrix-free jacobian and the code
> > which I use runtime option -snes_mf give the same residual history. But it
> > doesn't.
> >
> > Here is the histories for
> >
> > -snes_monitor -ksp_max_it 5 -snes_converged_reason -snes_max_it 2 -ksp_converged_reason -ksp_monitor -snes_max_linear_solve_fail 300 -pc_type none -snes_view -snes_linesearch_type basic
> >
> > 0 SNES Function norm 4.272952196300e-02
> >
> > 0 KSP Residual norm 4.272952196300e-02
> > 1 KSP Residual norm 4.234712668718e-02
> > 2 KSP Residual norm 3.683301946690e-02
> >
> > 3 KSP Residual norm 3.465586805169e-02
> >
> > 4 KSP Residual norm 3.452667066800e-02
> > 5 KSP Residual norm 3.451739518719e-02
> > Linear solve did not converge due to DIVERGED_ITS iterations 5
> > 1 SNES Function norm 4.203973403992e-02
> >
> > 0 KSP Residual norm 4.203973403992e-02
> >
> > 1 KSP Residual norm 4.203070641961e-02
> > 2 KSP Residual norm 4.202387940443e-02
> > 3 KSP Residual norm 4.183739347023e-02
> > 4 KSP Residual norm 4.183629424897e-02
> >
> > 5 KSP Residual norm 4.159456024825e-02
> >
> > Linear solve did not converge due to DIVERGED_ITS iterations 5
> > 2 SNES Function norm 4.200901009970e-02
> > Nonlinear solve did not converge due to DIVERGED_MAX_IT iterations 2
> >
> >
> > Here is the histories for
> > -snes_mf -snes_monitor -ksp_max_it 5 -snes_converged_reason -snes_max_it 2 -ksp_converged_reason -ksp_monitor -snes_max_linear_solve_fail 300 -pc_type none -snes_view -snes_linesearch_type basic
> >
> >
> > 0 SNES Function norm 4.272952196300e-02
> > 0 KSP Residual norm 4.272952196300e-02
> > 1 KSP Residual norm 4.270267664569e-02
> > 2 KSP Residual norm 3.690026921954e-02
> >
> > 3 KSP Residual norm 3.681740616743e-02
> >
> > 4 KSP Residual norm 3.464377294985e-02
> > 5 KSP Residual norm 3.464376048536e-02
> > Linear solve did not converge due to DIVERGED_ITS iterations 5
> > 1 SNES Function norm 3.461633424373e-02
> >
> > 0 KSP Residual norm 3.461633424373e-02
> >
> > 1 KSP Residual norm 3.461632119472e-02
> > 2 KSP Residual norm 3.406130197963e-02
> > 3 KSP Residual norm 3.406122155751e-02
> > 4 KSP Residual norm 3.403393397001e-02
> >
> > 5 KSP Residual norm 3.403367748538e-02
> >
> > Linear solve did not converge due to DIVERGED_ITS iterations 5
> > 2 SNES Function norm 3.403367847002e-02
> > Nonlinear solve did not converge due to DIVERGED_MAX_IT iterations 2
> >
> >
> > We can see that at 0 SNES 1 KSP step, the residual norms are different. Did I do something wrong here?
> >
> > The codes are like
> >
> > type MyJContext
> > Mat mf
> > Vec D
> > Vec work
> > end type MyJContext
> >
> > c
> >
> > type(MyJContext) myJctx
> >
> > --------------------------------------------------------------------------
> > call SNESCreate(PETSC_COMM_WORLD, snes, ierpetsc)
> > call SNESSetFunction(snes, pet_rhs_snes, flowsolrhs, ctx,
> >
> > @ ierpetsc)
> >
> > c
> > call MatCreateSNESMF(snes, myJctx%mf, ierpetsc)
> > call MatMFFDSetFunction(myJctx%mf, flowsolrhs, ctx, ierpetsc)
> > call VecDuplicate(pet_solu_snes, myJctx%D, ierpetsc)
> >
> > call VecDuplicate(pet_solu_snes, myJctx%work, ierpetsc)
> >
> > call VecSet(myJctx%D, 0.0D-3, ierpetsc)
> > call MatCreateShell(PETSC_COMM_WORLD, pet_nfff, pet_nfff,
> > @ PETSC_DETERMINE, PETSC_DETERMINE, ctx, myJ, ierpetsc)
> >
> > call MatShellSetOperation(myJ, MATOP_MULT, mymultply,
> >
> > @ ierpetsc)
> > call SNESSetJacobian(snes, myJ, pet_mat_pre,
> > @ flowsoljac, ctx, ierpetsc)
> >
> > --------------------------------------------------------------------------
> >
> >
> > subroutine mymultply ( A, x, y, ierpet)
> > Mat :: A
> > Vec :: x, y
> > PetscErrorCode :: ierpet
> > c
> > call MatMult(myJctx%mf,x,y, ierpet)
> > c
> > end
> > --------------------------------------------------------------------------
> >
> >
> > subroutine flowsoljac ( snes, pet_solu_snes, pet_mat_snes,
> > @ pet_mat_pre, flag, ctxx, ierrpet )
> >
> > c explicitly assemble pet_mat_pre matrix here
> > c .........
> > c .........
> >
> > call MatMFFDSetBase(myJctx%mf, pet_solu_snes, PETSC_NULL_INTEGER,
> >
> > @ ierrpet)
> >
> > end
> >
> >
> >
> >
> > On Fri, Jan 3, 2014 at 6:46 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> >
> > Dario,
> >
> > Your discussion below (SOR, ILU(200)) seems to imply that you are providing some actual explicit representation of the Jacobian, not just doing something completely matrix free. Is this correct? But the PETSc MatMFFD() is completely matrix free, it provides only a matrix-vector product and no access to the matrix entries, hence I am slightly confused.
> >
> > If you wish to use for the Jacobian something like D + J and do it completely matrix free then rather than than monkeying with “changing the function” I would
> > use the “correct" function to compute J x using matrix free multiply and then apply the D to as an additional operation. Hence you would do something like
> >
> > typedef struct { /* data structure to store the usual matrix free matrix and the additional diagonal matrix */
> > Mat mf;
> > Vec D;
> > Vec work;
> > } MyJContext;
> >
> > MyJContext myJctx;
> >
> > MatCreateSNESMF(SNES,&myJctx.mf); /* create the usual MFFD matrix using the real nonlinear function */
> > MatMFFDSetFunction(myJctx.mf,yournonlinearfunction,nonlinearfunctionctx);
> > VecCreate(comm,&myJctx.D);
> > /* set the correct sizes for D and fill up with your diagonal matrix entries */
> > VecDuplicate(&myJctx.D,&myJCtx.work);
> > MatCreateShell(comm,…. &myJ);
> > MatShellSetOperation(myJ,MATOP_MULT, mymultiply);
> > MatShellSetContext(myJ,&myJctx);
> > SNESSetJacobian(snes,myJ,myJ, myJFunction,NULL);
> >
> > where
> >
> > PetscErrorCode mymultply(Mat A,Vec x,Vec y) /* computes y = J x + D x
> > {
> > MyJContext *myJctx;
> >
> > MatShellGetContext(A,&myJctx);
> > MatMult(myJctx->mf,x,y);
> > VecPointwiseMult(myJctx->D,x,myJctx->work);
> > VecAXPY(y,1.myJctx->work);
> > }
> >
> > and
> >
> > PetscErrorCode myJFunction(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *str,void* ctx)
> >
> > /* this is called for each new “Jacobian” to set the point at which it is computed */
> > {
> > MyJContext *myJctx;
> > Vec f;
> > MatShellGetContext(*A,&myJctx);
> > SNESGetFunction(snes,&f);
> > MatMFFDSetBase(myJctx->mf,x,f);
> >
> > /* change the D entries if they depend on the current solution etc */
> > return 0;
> > }
> >
> > Sorry now that I have typed it out it looks a bit more complicated then it really is. It does what you want but without any trickery or confusing code.
> >
> > But, of course, since it is completely matrix free you cannot use SOR with it. Of course by making D suitably large you can make it as well conditioned as you want and thus get rapid linear convergence (though that may slow down or ruin the nonlinear convergence).
> >
> > Hope this helps,
> >
> > Barry
> >
> >
> >
> >
> >
> > On Jan 3, 2014, at 12:18 PM, Dario Isola <dario.isola at newmerical.com> wrote:
> >
> > > Dear, Barry and Jed,
> > >
> > > Thanks for your replies.
> > >
> > > We understand your doubts, so let me to put our question into context. In CFD it is standard practice to solve non-linear equations of conservation for steady flows by means of a inexact Newton method. The original Jacobian matrix is modified by adding terms on the diagonal which are proportional to the Courant number and to the lumped mass matrix. This allows us to obtain two things, "relax" the solution update and increase the diagonal dominance of the matrix itself.
> > >
> > > The latter is key when simple preconditioners are adopted, in our case point Jacobi or SOR. Indeed, if the original matrix was to be used, the GMRES method would converge only on very regular meshes and only when adopting ILU preconditioners with a very high level of fill-in. As result a higher number of non-linear iterations is traded with a simpler linear system to be solved.
> > >
> > > While exploring the SNES+MF capabilities we found out that we could successfully solve the linear system only with ILU(200) or so. Of course we do not want to touch the function used to evaluate the residual, which determines the final solution. However we think that a suitable modification of the function that Petsc differences to compute the matrix vector product would allow us to obtain a behavior similar to the inexact Newton method.
> > >
> > > Best regards,
> > > Dario
> > >
> > >
> > > On 01/03/2014 12:32 PM, Song Gao wrote:
> > >>
> > >>
> > >> ---------- Forwarded message ----------
> > >> From: Jed Brown <jedbrown at mcs.anl.gov>
> > >> Date: Thu, Jan 2, 2014 at 10:20 AM
> > >> Subject: Re: [petsc-users] SNESSetFunction and MatMFFDSetFunction
> > >> To: Song Gao <song.gao2 at mail.mcgill.ca>, Barry Smith <bsmith at mcs.anl.gov>
> > >> Cc: petsc-users <petsc-users at mcs.anl.gov>
> > >>
> > >>
> > >> Song Gao <song.gao2 at mail.mcgill.ca> writes:
> > >>
> > >> > Thanks, Barry.
> > >> >
> > >> > I mean 2) providing a function that I want PETSc to difference to evaluate
> > >> > the matrix vector product.
> > >> >
> > >> > I want to make a slight modification of the matrix after PETSc evaluate the
> > >> > matrix vector product.
> > >>
> > >> Performing a matrix-vector product is not supposed to modify the matrix.
> > >> It's unlikely that you really want this.
> > >>
> > >>
> > >> On Wed, Jan 1, 2014 at 3:01 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> > >>
> > >> On Jan 1, 2014, at 11:09 AM, Song Gao <song.gao2 at mail.mcgill.ca> wrote:
> > >>
> > >> > Dear all,
> > >> >
> > >> > Happy new year!
> > >> >
> > >> > I'm using the matrix-free method to solve NS equations. I call the SNESSetFunction to set the RHS function. I think SNES also uses that RHS function to evaluate the matrix vector product.
> > >>
> > >> Yes, PETSc differences this function to evaluate the matrix vector product.
> > >>
> > >> >
> > >> > But I want to set one function to evaluate the residual, and another different function to evaluate the matrix vector product.
> > >>
> > >> Are you providing a function that
> > >>
> > >> 1) actually evaluates the matrix vector product or are you
> > >>
> > >> 2) providing a function that you want PETSc to difference to evaluate the matrix vector product?
> > >>
> > >> > How can I do that? Does MatMFFDSetFunction do this job?
> > >>
> > >> For 2) yes, but if the function you provide is different than the function provided with SNESSetFunction then the matrix-vector product obtained from differencing it will not be “correct” Jacobian for the SNESSetFunction() you are providing so I don’t see why you would do it.
> > >>
> > >> For 1) you should use MatCreateShell() and MatShellSetOperation(mat,MATOP_MULT, ….) and then pass that matrix to SNESSetJacobian(), then PETSc will use that “matrix” to do its matrix-vector products.
> > >>
> > >> Barry
> > >>
> > >>
> > >>
> > >> >
> > >> > Any suggestion is appreciated. Thank you.
> > >> >
> > >> >
> > >> > Song
> > >>
> > >
> >
> >
>
>
>
More information about the petsc-users
mailing list