[petsc-users] PetscFV and TS implicit
Thibault Bridel-Bertomeu
thibault.bridelbertomeu at gmail.com
Mon Aug 24 00:56:55 CDT 2020
Barry, first of all, thank you very much for your detailed answer, I keep
reading it to let it soak in - I might come back to you for more details if
you do not mind.
In the meantime, to fuel the conversation, I attach to this e-mail two pdfs
containing the pieces of the code that regard what we are discussing. In
the *timedisc.pdf, you'll find how I handle the initialization of the TS
object, and in the *petscdefs.pdf you'll find the method that calls the
TSSolve as well as the methods that are linked to the TS (the timestep
adapt, the jacobian etc ...). [Sorry for the quality, I cannot do better
than that sort of pdf ...]
Based on what is in the structured code I sent you the other day, I rewrote
the PetscJacobianFunction_JFNK. I think it should be all right, but
although it compiles, execution raises a seg fault I think when I do
ierr = TSSetIJacobian(ts, A, A, PetscIJacobian, user);
saying that A does not have the right dimensions. It is quite new, I am
still looking into where exactly the error is raised. What do you think of
this implementation though, does it look correct in your expert eyes ?
As for what we really discussed so far, it's that PetscComputePreconMatImpl
that I do not know how to implement (with the derivative of the jacobian
based on the FVM object).
I understand now that what I am showing you today might not be the right
way to go if one wants to really use the PetscFV, but I just wanted to add
those code lines to the conversation to have your feedback.
Thank you again for your help,
Thibault
Le ven. 21 août 2020 à 19:25, Barry Smith <bsmith at petsc.dev> a écrit :
>
>
> On Aug 21, 2020, at 10:58 AM, Thibault Bridel-Bertomeu <
> thibault.bridelbertomeu at gmail.com> wrote:
>
> Thank you Barry for the tip ! I’ll make sure to do that when everything is
> set.
> What I also meant is that there will not be any more direct way to set the
> preconditioner than to go through SNESSetJacobian after having assembled
> everything by hand ? Like, in my case, or in the more general case of fluid
> dynamics equations, the preconditioner is not a fun matrix to assemble,
> because for every cell the derivative of the physical flux jacobian has to
> be taken and put in the right block in the matrix - finite element style if
> you want. Is there a way to do that with Petsc methods, maybe
> short-circuiting the FEM based methods ?
>
>
> Thibault
>
> I am not sure what you mean but there are a couple of things that may
> be helpful.
>
> PCSHELL
> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCSHELL.html allows
> you to build your own preconditioner (that can and often will use one or
> more of its own Mats, and KSP or PC inside it, or even use another PETScFV
> etc to build some of the sub matrices for you if it is appropriate), this
> approach means you never need to construct a "global" PETSc matrix from
> which PETSc builds the preconditioner. But you should only do this if the
> conventional approach is not reasonable for your problem.
>
> MATNEST
> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MATNEST.html allows
> you to build a global matrix by building parts of it separately and even
> skipping parts you decide you don't need in the preconditioner.
> Conceptually it is the same as just creating a global matrix and filling up
> but the process is a bit different and something suitable for "multi
> physics" or "multi-equation" type applications.
>
> Of course what you put into PCSHELL and MATNEST will affect the
> convergence of the nonlinear solver. As Jed noted what you put in the
> "Jacobian" does not have to be directly the same mathematically as what you
> put into the TSSetI/RHSFunction with the caveat that it does have to
> appropriate spectral properties to result in a good preconditioner for the
> "true" Jacobian.
>
> Couple of other notes:
>
> The entire business of "Jacobian" matrix-free or not (with for example
> -snes_fd_operator) is tricky because as Jed noted if your finite volume
> scheme has non-differential terms such as if () tests. There is a concept
> of sub-differential for this type of thing but I know absolutely nothing
> about that and probably not worth investigating.
>
> In this situation you can avoid the "true" Jacobian completely (both for
> matrix-vector product and preconditioner) and use something else as Jed
> suggested a lower order scheme that is differentiable. This can work well
> for solving the nonlinear system or not depending on how suitable it is for
> your original "function"
>
>
> 1) In theory at least you can have the Jacobian matrix-vector product
> computed directly using DMPLEX/PETScFV infrastructure (it would apply the
> Jacobian locally matrix-free using code similar to the code that evaluates
> the FV "function". I do no know if any of this code is written, it will be
> more efficient than -snes_mf_operator that evaluates the FV "function" and
> does traditional differencing to compute the Jacobian. Again it has the
> problem of non-differentialability if the function is not differential. But
> it could be done for a different (lower order scheme) that is
> differentiable.
>
> 2) You can have PETSc compute the Jacobian explicitly coloring and from
> that build the preconditioner, this allows you to avoid the hassle of
> writing the code for the derivatives yourself. This uses finite differences
> on your function and coloring of the graph to compute many columns of the
> Jacobian simultaneously and can be pretty efficient. Again if the function
> is not differential there can be issues of what the result means and will
> it work in a nonlinear solver. SNESComputeJacobianDefaultColor
> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/SNESComputeJacobianDefaultColor.html
>
> 3) Much more outlandish is to skip Newton and Jacobians completely and use
> the full approximation scheme SNESFAS
> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNESFAS/SNESFAS.html this
> requires a grid hierarchy and appropriate way to interpolate up through the
> grid hierarchy your finite volume solutions. Probably not worth
> investigating unless you have lots of time on your hands and keen interest
> in this kind of stuff https://arxiv.org/pdf/1607.04254.pdf
>
> So to summarize, and Matt and Jed can correct my mistakes.
>
> 1) Form the full Jacobian from the original "function" using analytic
> approach use it for both the matrix-vector product and to build the
> preconditioner. Problem if full Jacobian not well defined mathematically.
> Tough to code, usually not practical.
>
> 2) Do any matrix free (any way) for the full Jacobian and
>
> a) build another "approximate" Jacobian (using any technique analytic or
> finite differences using matrix coloring on a new "lower order" "function")
> Still can have trouble if this original Jacobian is no well defined
>
> b) "write your own preconditioner" that internally can use anything in
> PETSc that approximately solves the Jacobian. Same potential problems if
> original Jacobian is not differential, plus convergence will depend on how
> good your own preconditioner approximates the inverse of the true Jacobian.
>
> 3) Use a lower Jacobian (computed anyway you want) for the matrix-vector
> product and the preconditioner. The problem of differentiability is gone
> but convergence of the nonlinear solver depends on how well lower order
> Jacobian is appropriate for the original "function"
>
> a) Form the "lower order" Jacobian analytically or with coloring and
> use for both matrix-vector product and building preconditioner. Note that
> switching between this and 2a is trivial.
>
> b) Do the "lower order" Jacobian matrix free and provide your own
> PCSHELL. Note that switching between this and 2b is trivial.
>
> Barry
>
> I would first try competing the "true" Jacobian via coloring, if that
> works and give satisfactory results (fast enough) then stop.
>
> Then I would do 2a/2b by writing my "function" using PETScFV and writing
> the "lower order function" via PETScFV and use matrix coloring to get the
> Jacobian from the second "lower order function". If this works well (either
> with 2a or 3a or both) then stop or you can compute the "lower order"
> Jacobian analytically (again using PetscFV) for a more efficient evaluation
> of the Jacobian.
>
>
>
>
> Thanks !
>
> Thibault
>
> Le ven. 21 août 2020 à 17:22, Barry Smith <bsmith at petsc.dev> a écrit :
>
>>
>>
>> On Aug 21, 2020, at 8:35 AM, Thibault Bridel-Bertomeu <
>> thibault.bridelbertomeu at gmail.com> wrote:
>>
>>
>>
>> Le ven. 21 août 2020 à 15:23, Matthew Knepley <knepley at gmail.com> a
>> écrit :
>>
>>> On Fri, Aug 21, 2020 at 9:10 AM Thibault Bridel-Bertomeu <
>>> thibault.bridelbertomeu at gmail.com> wrote:
>>>
>>>> Sorry, I sent too soon, I hit the wrong key.
>>>>
>>>> I wanted to say that context.npoints is the local number of cells.
>>>>
>>>> PetscRHSFunctionImpl allows to generate the hyperbolic part of the
>>>> right hand side.
>>>> Then we have :
>>>>
>>>> PetscErrorCode PetscIJacobian(
>>>> TS ts, /*!< Time stepping object (see PETSc TS)*/
>>>> PetscReal t, /*!< Current time */
>>>> Vec Y, /*!< Solution vector */
>>>> Vec Ydot, /*!< Time-derivative of solution vector */
>>>> PetscReal a, /*!< Shift */
>>>> Mat A, /*!< Jacobian matrix */
>>>> Mat B, /*!< Preconditioning matrix */
>>>> void *ctxt /*!< Application context */
>>>> )
>>>> {
>>>> PETScContext *context = (PETScContext*) ctxt;
>>>> HyPar *solver = context->solver;
>>>> _DECLARE_IERR_;
>>>>
>>>> PetscFunctionBegin;
>>>> solver->count_IJacobian++;
>>>> context->shift = a;
>>>> context->waqt = t;
>>>> /* Construct preconditioning matrix */
>>>> if (context->flag_use_precon) { IERR PetscComputePreconMatImpl(B,Y,
>>>> context); CHECKERR(ierr); }
>>>>
>>>> PetscFunctionReturn(0);
>>>> }
>>>>
>>>> and PetscJacobianFunction_JFNK which I bind to the matrix shell,
>>>> computes the action of the jacobian on a vector : say U0 is the state of
>>>> reference and Y the vector upon which to apply the JFNK method, then the
>>>> PetscJacobianFunction_JFNK returns shift * Y - 1/epsilon * (F(U0 +
>>>> epsilon*Y) - F(U0)) where F allows to evaluate the hyperbolic flux (shift
>>>> comes from the TS).
>>>> The preconditioning matrix I compute as an approximation to the actual
>>>> jacobian, that is shift * Identity - Derivative(dF/dU) where dF/dU is, in
>>>> each cell, a 4x4 matrix that is known exactly for the system of equations I
>>>> am solving, i.e. Euler equations. For the structured grid, I can loop on
>>>> the cells and do that 'Derivative' thing at first order by simply taking a
>>>> finite-difference like approximation with the neighboring cells,
>>>> Derivative(phi) = phi_i - phi_{i-1} and I assemble the B matrix block by
>>>> block (JFunction is the dF/dU)
>>>>
>>>> /* diagonal element */
>>>>
>>>>
>>>> for (v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pg
>>>> + v; }
>>>>
>>>>
>>>> ierr = solver->JFunction
>>>> <http://hypar.github.io/a00144.html#a25c87c7c874c7e08981688fe6c51167b>
>>>> (values,(u+nvars*p),solver->physics
>>>> <http://hypar.github.io/a00144.html#a48e81a1806b5774943fd9a26e9a190f2>
>>>> ,dir,0);
>>>>
>>>>
>>>> _ArrayScale1D_
>>>> <http://hypar.github.io/a00165.html#a1dab5f429ccd03c4a7fbb605d8ef75b6>
>>>> (values,(dxinv*iblank),(nvars*nvars));
>>>>
>>>>
>>>> ierr = MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES);
>>>> CHKERRQ(ierr);
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> /* left neighbor */
>>>>
>>>>
>>>> if (pgL >= 0) {
>>>>
>>>>
>>>> for (v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] =
>>>> nvars*pgL + v; }
>>>>
>>>>
>>>> ierr = solver->JFunction
>>>> <http://hypar.github.io/a00144.html#a25c87c7c874c7e08981688fe6c51167b>
>>>> (values,(u+nvars*pL),solver->physics
>>>> <http://hypar.github.io/a00144.html#a48e81a1806b5774943fd9a26e9a190f2>
>>>> ,dir,1);
>>>>
>>>>
>>>> _ArrayScale1D_
>>>> <http://hypar.github.io/a00165.html#a1dab5f429ccd03c4a7fbb605d8ef75b6>
>>>> (values,(-dxinv*iblank),(nvars*nvars));
>>>>
>>>>
>>>> ierr = MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES);
>>>> CHKERRQ(ierr);
>>>>
>>>>
>>>> }
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> /* right neighbor */
>>>>
>>>>
>>>> if (pgR >= 0) {
>>>>
>>>>
>>>> for (v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] =
>>>> nvars*pgR + v; }
>>>>
>>>>
>>>> ierr = solver->JFunction
>>>> <http://hypar.github.io/a00144.html#a25c87c7c874c7e08981688fe6c51167b>
>>>> (values,(u+nvars*pR),solver->physics
>>>> <http://hypar.github.io/a00144.html#a48e81a1806b5774943fd9a26e9a190f2>
>>>> ,dir,-1);
>>>>
>>>>
>>>> _ArrayScale1D_
>>>> <http://hypar.github.io/a00165.html#a1dab5f429ccd03c4a7fbb605d8ef75b6>
>>>> (values,(-dxinv*iblank),(nvars*nvars));
>>>>
>>>>
>>>> ierr = MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES);
>>>> CHKERRQ(ierr);
>>>>
>>>>
>>>> }
>>>>
>>>>
>>>>
>>>> I do not know if I am clear here ...
>>>> Anyways, I am trying to figure out how to do this shell matrix and this
>>>> preconditioner using all the FV and DMPlex artillery.
>>>>
>>>
>>> Okay, that is very clear. We should be able to get the JFNK just with
>>> -snes_mf_operator, and put the approximate J construction in
>>> DMPlexComputeJacobian_Internal().
>>> There is an FV section already, and we could just add this. I would need
>>> to understand those entries in the pointwise Riemann sense that the other
>>> stuff is now.
>>>
>>
>> Ok i had a quick look and if I understood correctly it would do the job.
>> Setting the snes-mf-operator flag would mean however that we have to go
>> through SNESSetJacobian to set the jacobian and the preconditioning matrix
>> wouldn't it ?
>>
>>
>> Thibault,
>>
>> Since the TS implicit methods end up using SNES internally the option
>> should be available to you without requiring you to be calling the SNES
>> routines directly
>>
>> Once you have finalized your approach and if for the implicit case you
>> always work in the snes mf operator mode you can hardwire
>>
>> TSGetSNES(ts,&snes);
>> SNESSetUseMatrixFree(snes,PETSC_TRUE,PETSC_FALSE);
>>
>> in your code so you don't need to always provide the option
>> -snes-mf-operator
>>
>> Barry
>>
>>
>>
>>
>> There might be calls to the Riemann solver to evaluate the dRHS / dU part
>> yes but maybe it's possible to re-use what was computed for the RHS^n ?
>> In the FV section the jacobian is set to identity which I missed before,
>> but it could explain why when I used the following :
>>
>> TSSetType(ts, TSBEULER);
>> DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM <https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/DMPlexTSComputeIFunctionFEM.html#DMPlexTSComputeIFunctionFEM>, &ctx);
>> DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM <https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/DMPlexTSComputeIJacobianFEM.html#DMPlexTSComputeIJacobianFEM>, &ctx);
>>
>> with my FV discretization nothing happened, right ?
>>
>> Thank you,
>>
>> Thibault
>>
>> Thanks,
>>>
>>> Matt
>>>
>>>
>>>> Le ven. 21 août 2020 à 14:55, Thibault Bridel-Bertomeu <
>>>> thibault.bridelbertomeu at gmail.com> a écrit :
>>>>
>>>>> Hi,
>>>>>
>>>>> Thanks Matthew and Jed for your input.
>>>>> I indeed envision an implicit solver in the sense Jed mentioned - Jiri
>>>>> Blazek's book is a nice intro to this concept.
>>>>>
>>>>> Matthew, I do not know exactly what to change right now because
>>>>> although I understand globally what the DMPlexComputeXXXX_Internal methods
>>>>> do, I cannot say for sure line by line what is happening.
>>>>> In a structured code, I have a an implicit FVM solver with PETSc but I
>>>>> do not use any of the FV structure, not even a DM - I just use C arrays
>>>>> that I transform to PETSc Vec and Mat and build my IJacobian and my
>>>>> preconditioner and gives all that to a TS and it runs. I cannot figure out
>>>>> how to do it with the FV and the DM and all the underlying "shortcuts" that
>>>>> I want to use.
>>>>>
>>>>> Here is the top method for the structured code :
>>>>>
>>>>> int total_size = context.npoints * solver->nvars
>>>>> ierr = TSSetRHSFunction(ts,PETSC_NULL,PetscRHSFunctionImpl,&context);
>>>>> CHKERRQ(ierr);
>>>>> SNES snes;
>>>>> KSP ksp;
>>>>> PC pc;
>>>>> SNESType snestype;
>>>>> ierr = TSGetSNES(ts,&snes); CHKERRQ(ierr);
>>>>> ierr = SNESGetType(snes,&snestype); CHKERRQ(ierr);
>>>>>
>>>>> flag_mat_a = 1;
>>>>> ierr = MatCreateShell(MPI_COMM_WORLD,total_size,total_size,
>>>>> PETSC_DETERMINE,
>>>>> PETSC_DETERMINE,&context,&A); CHKERRQ(ierr);
>>>>> context.jfnk_eps = 1e-7;
>>>>> ierr = PetscOptionsGetReal(NULL,NULL,"-jfnk_epsilon",&context.jfnk_eps
>>>>> ,NULL); CHKERRQ(ierr);
>>>>> ierr = MatShellSetOperation(A,MATOP_MULT,(void (*)(void))
>>>>> PetscJacobianFunction_JFNK); CHKERRQ(ierr);
>>>>> ierr = MatSetUp(A); CHKERRQ(ierr);
>>>>>
>>>>> context.flag_use_precon = 0;
>>>>> ierr = PetscOptionsGetBool(PETSC_NULL,PETSC_NULL,"-with_pc",(PetscBool
>>>>> *)(&context.flag_use_precon),PETSC_NULL); CHKERRQ(ierr);
>>>>>
>>>>> /* Set up preconditioner matrix */
>>>>> flag_mat_b = 1;
>>>>> ierr = MatCreateAIJ(MPI_COMM_WORLD,total_size,total_size,
>>>>> PETSC_DETERMINE,PETSC_DETERMINE,
>>>>> (solver->ndims*2+1)*solver->nvars,NULL,
>>>>> 2*solver->ndims*solver->nvars,NULL,&B); CHKERRQ(ierr);
>>>>> ierr = MatSetBlockSize(B,solver->nvars);
>>>>> /* Set the RHSJacobian function for TS */
>>>>> ierr = TSSetIJacobian(ts,A,B,PetscIJacobian,&context); CHKERRQ(ierr);
>>>>>
>>>>> Thibault Bridel-Bertomeu
>>>>> —
>>>>> Eng, MSc, PhD
>>>>> Research Engineer
>>>>> CEA/CESTA
>>>>> 33114 LE BARP
>>>>> Tel.: (+33)557046924
>>>>> Mob.: (+33)611025322
>>>>> Mail: thibault.bridelbertomeu at gmail.com
>>>>>
>>>>>
>>>>> Le jeu. 20 août 2020 à 18:43, Jed Brown <jed at jedbrown.org> a écrit :
>>>>>
>>>>>> Matthew Knepley <knepley at gmail.com> writes:
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> > I could never get the FVM stuff to make sense to me for implicit
>>>>>> methods.
>>>>>>
>>>>>>
>>>>>> > Here is my problem understanding. If you have an FVM method, it
>>>>>> decides
>>>>>>
>>>>>>
>>>>>> > to move "stuff" from one cell to its neighboring cells depending on
>>>>>> the
>>>>>>
>>>>>>
>>>>>> > solution to the Riemann problem on each face, which computed the
>>>>>> flux. This
>>>>>>
>>>>>>
>>>>>> > is
>>>>>>
>>>>>>
>>>>>> > fine unless the timestep is so big that material can flow through
>>>>>> into the
>>>>>>
>>>>>>
>>>>>> > cells beyond the neighbor. Then I should have considered the effect
>>>>>> of the
>>>>>>
>>>>>>
>>>>>> > Riemann problem for those interfaces. That would be in the
>>>>>> Jacobian, but I
>>>>>>
>>>>>>
>>>>>> > don't know how to compute that Jacobian. I guess you could do
>>>>>> everything
>>>>>>
>>>>>>
>>>>>> > matrix-free, but without a preconditioner it seems hard.
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> So long as we're using method of lines, the flux is just
>>>>>> instantaneous flux, not integrated over some time step. It has the same
>>>>>> meaning for implicit and explicit.
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> An explicit method would be unstable if you took such a large time
>>>>>> step (CFL) and an implicit method will not simultaneously be SSP and higher
>>>>>> than first order, but it's still a consistent discretization of the problem.
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> It's common (done in FUN3D and others) to precondition with a
>>>>>> first-order method, where gradient reconstruction/limiting is skipped.
>>>>>> That's what I'd recommend because limiting creates nasty nonlinearities and
>>>>>> the resulting discretizations lack h-ellipticity which makes them very hard
>>>>>> to solve.
>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>>
>>>>
>>>>
>>>
>>> --
>>> 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.cse.buffalo.edu/~knepley/>
>>>
>>>
>>>
>>
>>
>> --
> Thibault Bridel-Bertomeu
> —
> Eng, MSc, PhD
> Research Engineer
> CEA/CESTA
> 33114 LE BARP
> Tel.: (+33)557046924
> Mob.: (+33)611025322
> Mail: thibault.bridelbertomeu at gmail.com
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200824/915c7d67/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: petsc_fvm_part_timedisc_v2.pdf
Type: application/pdf
Size: 303982 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200824/915c7d67/attachment-0002.pdf>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: petsc_fvm_part_petscdefs_v2.pdf
Type: application/pdf
Size: 413581 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200824/915c7d67/attachment-0003.pdf>
More information about the petsc-users
mailing list