[petsc-users] Finite difference Jacobian in TS context

[Oleksandr Koshkarov]

Dear All,

sorry for spam because of my poor PETSc knowledge (I am just starting 
with this nice framework).

I think, I figured part of it out. However, I want to point that 
src/ts/examples/tutorials/ex15.c is misleading. (or maybe it is a bug?)

in this example we have

TSGetSNES(ts,&snes);
MatCreateSNESMF(snes,&Jmf);
SNESSetJacobian(snes,Jmf,J,SNESComputeJacobianDefault,NULL);

// or this: SNESSetJacobian(snes,Jmf,J,SNESComputeJacobianDefaultColor,0);

which implies (I think) that Jacobian would be matrix free. And if one 
would use PCNONE for preconditioning the matrix would never be 
allocated. However, it seems in reality it allocates matrix.

To avoid it, I used

MatCreateSNESMF(snes,&J);

SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,0);

which seems to work fine. I am not sure I fully understand the 
difference and i have zero intuition and I also have no idea what 
happens with preconditioning in this case. If someone have some useful 
comets, please share :) (I read the relevant section in PETSc manual, 
but still not fully understanding what I should use when)

Thank you and again sorry for the spam,
Oleksandr.

On 04/28/2018 07:20 PM, Smith, Barry F. wrote:
> ~/Src/petsc/src/ts/examples/tutorials
> $ grep SNESComputeJacobianDefaultColor *.c
> ex10.c:    ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobianDefaultColor,0);CHKERRQ(ierr);
> ex15.c:      ierr = SNESSetJacobian(snes,Jmf,J,SNESComputeJacobianDefaultColor,0);CHKERRQ(ierr);
> ex17.c:    ierr = SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,0);CHKERRQ(ierr);
>
>     I don't think you need to explicitly create the MatFDColoring object.
>
>     Please take a look at ex15.c and see if you can get it working like that example. If that doesn't work let us know and we can take a closer look at it.
>
>
>     Barry
>
>
>> On Apr 28, 2018, at 8:05 PM, Oleksandr Koshkarov <olk548 at mail.usask.ca> wrote:
>>
>> Hello All,
>>
>> I hope someone can help :) I think I am doing something wrong, but cannot understand what. I have a huge time dependent system with 3d DMDA data structure and I am evolving it with explicit Runge-Kutta by using TS and basically only using "TSSetRHSFunction". Now I want to repeat it with implicit time stepper (for now Crank-Nicolson) and I am trying to provide finite difference Jacobian and I am failing miserably. I also cannot find appropriate example in PETSc tutorial (if you can point me to working example, it would be great).
>>
>> Here is my best attempt (what wrong with it?):
>>
>>    DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC,
>>             DMDA_STENCIL_STAR,
>>             NX, NY, NZ,
>>             PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE,
>>             2*3+NC*NS,
>>             1,
>>             NULL,  NULL, NULL,  &da);
>>    DMSetUp(da);
>>    DMCreateGlobalVector(da,&x);
>>    TSCreate(PETSC_COMM_WORLD,&ts);
>>    TSSetProblemType(ts,TS_NONLINEAR);
>>    TSSetRHSFunction(ts,NULL,compute_RHS,NULL);
>>    TSSetMaxTime(ts,T_FINAL);
>>    TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
>>    TSSetDM(ts,da);
>>    TSSetType(ts,TSCN); //it works with: TSSetType(ts,TSRK);
>>    set_IC(da,x);
>>    TSSetTimeStep(ts,DT);
>>    TSSetSolution(ts,x);
>>
>>    TSGetSNES(ts,&snes);
>>    SNESGetKSP(snes,&ksp);
>>    KSPGetPC(ksp,&pc);
>>    PCSetType(pc,PCNONE);
>>
>>    DMSetMatType(da,MATAIJ);
>>    DMCreateMatrix(da,&J);
>>    ISColoring iscoloring;
>>    MatFDColoring  matfdcoloring;
>>    DMCreateColoring(da,IS_COLORING_GLOBAL,&iscoloring);
>>    MatFDColoringCreate(J,iscoloring,&matfdcoloring);
>>
>>    MatFDColoringSetType(matfdcoloring,MATMFFD_DS);
>>
>> // I think I do something wrong in the following 3 lines
>>
>>    PetscErrorCode (*temp_f)(SNES,Vec,Vec,void*);
>>    SNESGetFunction(snes,NULL,&temp_f,NULL);
>>    MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))temp_f,NULL);
>>
>>    MatFDColoringSetUp(J,iscoloring,matfdcoloring);
>>
>> SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring);
>>    ISColoringDestroy(&iscoloring);
>>
>>    TSSolve(ts,x);
>>
>> Thank you,
>>
>> Oleksandr Koshkarov.
>>