> On Apr 29, 2018, at 5:40 PM, Oleksandr Koshkarov <[email protected]> wrote:
>
> 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.
The above uses matrix-free to do matrix-vector products for linear system
but constructs the preconditioner by building the Jacobian via differencing and
then using that matrix to build the preconditioner.
>
> 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)
This does not build the Jacobian (only has matrix free matrix vector
products) so requires much less memory but likely for large problem sizes the
linear system solve will be slow (require many iterations) or won't converge at
all. The conditioning of the linear system depends on the exact problem you
are solving and the type of discretization you are using. There is no easy
rules that always apply but for most discretizations of PDEs the number of
iterations needed by the linear solver increases with the problem sizes. This
means for most large problems matrix-free (without building any sort of
jacobean and preconditioner) is impractical and one needs to pay the price of
using more memory to get reasonable convergence.
Barry
>
> 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 <[email protected]>
>>> 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.
>>>
>
