[petsc-users] Finite difference Jacobian in TS context
Matthew Knepley
knepley at gmail.com
Mon Apr 30 05:38:57 CDT 2018
On Sun, Apr 29, 2018 at 8:29 PM, Oleksandr Koshkarov <olk548 at mail.usask.ca>
wrote:
> Thank you, it clarifies a lot :)
>
> 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.
>
>
> Yes, I was thinking to use FD matrix free Jacobian and for preconditioning
> to use Krylov method to invert the linearized Jacobian which also would be
> FD matrix free. Sounds bad?
>
It will not be scalable, unless you are using the Boundary Element Method.
PDE operators have condition number that
grows pretty rapidly with size.
> Well, from your advices I convinced to try constructing actual Jacobian
> and use it as Jacobian and matrix from which petsc would construct
> preconditioner.
>
Or the first option: Apply J matrix-free for the action in a Krylov method,
but build a perhaps simplified operator for preconditioning.
If there is no obvious simplification, then just build the whole thing.
Thanks,
Matt
> p.s. Do I reply to emails correctly? or should I replay to
> petsc-users at mcs.anl.gov only?
>
> Thank you,
> Oleksandr.
>
> On 04/29/2018 06:03 PM, Matthew Knepley wrote:
>
> 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,SNE
>>>> SComputeJacobianDefaultColor,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,mat
>>>>>> fdcoloring);
>>>>>> 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/%7Emk51/>
>
>
>
--
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/20180430/94ee157e/attachment.html>
More information about the petsc-users
mailing list