For the record, I did as Barry suggested, and it worked well for our problem.
Thank you very much. On Fri, Jan 17, 2014 at 4:47 PM, Barry Smith <[email protected]> wrote: > > On Jan 16, 2014, at 10:49 AM, Song Gao <[email protected]> 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 <[email protected]> > 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 <[email protected]> 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 <[email protected]> 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 <[email protected]> > 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 <[email protected]> > 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 <[email protected]> > > > >> Date: Thu, Jan 2, 2014 at 10:20 AM > > > >> Subject: Re: [petsc-users] SNESSetFunction and MatMFFDSetFunction > > > >> To: Song Gao <[email protected]>, Barry Smith < > [email protected]> > > > >> Cc: petsc-users <[email protected]> > > > >> > > > >> > > > >> Song Gao <[email protected]> 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 <[email protected]> > wrote: > > > >> > > > >> On Jan 1, 2014, at 11:09 AM, Song Gao <[email protected]> > 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 > > > >> > > > > > > > > > > > > > > > > > >
