[petsc-users] Finite difference Jacobian in TS context

[Oleksandr Koshkarov]

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.