[petsc-users] Finite difference Jacobian in TS context
Matthew Knepley
knepley at gmail.com
Sun Apr 29 19:03:46 CDT 2018
On Sun, Apr 29, 2018 at 7:44 PM, Oleksandr Koshkarov <olk548 at mail.usask.ca>
wrote:
> Thank you. A little clarification:
>
> 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.
>>
> So SNESComputeJacobianDefaultColor will use memory (construct Jacobian)
> even if preconditioning is set to PCNONE and J=NULL? (that is what I saw in
> my example)
>
1) SNESComputeJacobianDefaultColor() creates the Jacobian, fullstop.
2) When you pass J=NULL, the Jacobian is created automatically by DMDA
3) The PC has no influence on the assembly process.
> 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.
>>
> So the following command ignores precoditioning?
Ignore is the wrong word. How do you form a preconditioner from an object
with no values in it?
> Can I still use something like PCSHELL?
>
What exactly would you do in your shell? If you are only using the action
of an operator, it
usually equivalent to some Krylov method.
> SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,0);
>
> p.s. for my problems if probably be unrealistic to construct Jacobian
> (state size will start from N > 1000^4).
>
The number of dofs does not tell us anything. You would need to know the
sparsity. People regularly solve
problems with billions of unknowns.
Thanks,
Matt
> Thank you,
>
> Oleksandr.
>
>
> On 04/29/2018 05:26 PM, Smith, Barry F. wrote:
>
>>
>> On Apr 29, 2018, at 5:40 PM, Oleksandr Koshkarov <olk548 at mail.usask.ca>
>>> 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,SNESC
>>>> omputeJacobianDefaultColor,0);CHKERRQ(ierr);
>>>> ex15.c: ierr = SNESSetJacobian(snes,Jmf,J,SNE
>>>> SComputeJacobianDefaultColor,0);CHKERRQ(ierr);
>>>> ex17.c: ierr = SNESSetJacobian(snes,J,J,SNESC
>>>> omputeJacobianDefaultColor,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.
>>>>>
>>>>>
>
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/~mk51/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180429/6393acc2/attachment.html>
More information about the petsc-users
mailing list