Re: [petsc-users] Finite difference Jacobian in TS context

[Oleksandr Koshkarov]

Hello all,

Sorry for stupid question, but do "SNESComputeJacobianDefault" or 
"SNESComputeJacobianDefaultColor" make Jacobian matrix free? (I was 
thinking so from the start... but now it seems I am totally wrong)

Best regards,

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 <olk...@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.