[petsc-users] PetscFV and TS implicit
Thibault Bridel-Bertomeu
thibault.bridelbertomeu at gmail.com
Mon Aug 24 11:38:02 CDT 2020
Thank you Barry for taking the time to go through the code !
I indeed figured out this afternoon that the function related to the
matrix-vector product is always handling global vectors. I corrected mine
so that it compiles, but I have a feeling it won't run properly without a
preconditioner.
Anyways you are right, my PetscJacobianFunction_JFNK() aims at doing some
basic finite-differencing ; user->RHS_ref is my F(U) if you see the system
as dU/dt = F(U). As for the DMGlobalToLocal() it was there mainly because I
had not realized the vectors I was manipulating were global.
I will take your advice and try with just the SNESSetUseMatrixFree.
I haven't quite fully understood what it does "under the hood" though: just
calling SNESSetUseMatrixFree(snes,PETSC_TRUE,PETSC_FALSE) before the
TSSolve call is enough to ensure that the implicit matrix is computed ?
Does it use the function we set as a RHS to build the matrix ?
To create the preconditioner I will do as you suggest too, thank you. This
matrix has to be as close as possible to the inverse of the implicit matrix
to ensure that the eigenvalues of the system are as close to 1 as possible.
Given the implicit matrix is built "automatically" thanks to the SNES
matrix free capability, can we use that matrix as a starting point to the
building of the preconditioner ? You were talking about the coloring
capabilities in PETSc, is that where it can be applied ?
Thank you so much again,
Thibault
Le lun. 24 août 2020 à 15:45, Barry Smith <bsmith at petsc.dev> a écrit :
>
> I think the attached is wrong.
>
>
> The input to the matrix vector product for the Jacobian is always global
> vectors which means on each process the dimension is not the size of the
> DMGetLocalVector() it should be the VecGetLocalSize() of the
> DMGetGlobalVector()
>
> But you may be able to skip all this and have the DM create the shell
> matrix setting it sizes appropriately and you only need to supply the MATOP
>
> DMSetMatType(dm,MATSHELL);
> DMCreateMatrix(dm,&A);
>
> In fact, I also don't understand the PetscJacobianFunction_JFKN() function
> It seems to be doing finite differencing on the
> DMPlexTSComputeRHSFunctionFVM() assuming the current function value is in
> usr->RHS_ref. How is this different than just letting PETSc/SNES used
> finite differences to do the matrix-vector product. Your code seems rather
> complicated with the DMGlobalToLocal() which I don't understand what it is
> suppose to do there.
>
> I think you can just call
>
> TSGetSNES()
> SNESSetUseMatrixFree(snes,PETSC_TRUE,PETSC_FALSE);
>
> and it will set up an internal matrix that does the finite differencing
> for you. Then you never need a shell matrix.
>
>
> Also to create the preconditioner matrix B this should work
>
> DMSetMatType(dm,MATAIJ);
> DMCreateMatrix(dm,&B);
>
> no need for you to figure out the sizes.
>
>
> Note that both A and B need to have the same dimensions on each process as
> the global vectors which I don't think your current code has.
>
>
>
> Barry
>
>
>
> On Aug 24, 2020, at 12:56 AM, Thibault Bridel-Bertomeu <
> thibault.bridelbertomeu at gmail.com> wrote:
>
> 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
>>
>>
>>
>>
>
> <petsc_fvm_part_timedisc_v2.pdf><petsc_fvm_part_petscdefs_v2.pdf>
>
>
>
--
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/934d4963/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Untitled.png
Type: image/png
Size: 120787 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200824/934d4963/attachment-0001.png>
More information about the petsc-users
mailing list