[petsc-users] PetscFV and TS implicit
Thibault Bridel-Bertomeu
thibault.bridelbertomeu at gmail.com
Mon Aug 24 14:20:49 CDT 2020
Good evening everyone,
Thanks Barry for your answer.
Le lun. 24 août 2020 à 18:51, Barry Smith <bsmith at petsc.dev> a écrit :
>
>
> On Aug 24, 2020, at 11:38 AM, Thibault Bridel-Bertomeu <
> thibault.bridelbertomeu at gmail.com> wrote:
>
> 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 ?
>
>
> All it does is "replace" the A matrix with one automatically created
> for the job using MatCreateMFFD(). It does not touch the B matrix, it does
> not build the matrix but yes if does use the function to provide to do the
> differencing.
>
OK, thank you. This MFFD Matrix is then called by the TS to construct the
linear system that will be solved to advance the system of equations, right
?
>
> 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 ?
>
>
> No the MatrixFree doesn't build a matrix, it can only do matrix-vector
> products with differencing.
>
My bad, wrong word. Yes of course it's all matrix-free hence it's just a
functional, however maybe the inner mechanisms can be accessed and used for
the preconditioner ?
> You were talking about the coloring capabilities in PETSc, is that where
> it can be applied ?
>
>
> Yes you can use that. See MatFDColoringCreate() but since you are using
> a DM in theory you can use -snes_fd_color and PETSc will manage everything
> for you so you don't have to write any code for Jacobians at all. Again it
> uses your function to do differences using coloring to be efficient to
> build the Jacobian for you.
>
I read a bit about the coloring you are mentioning. As I understand it, it
is another option to have a matrix-free Jacobian behavior during the
Newton-Krylov iterations, right ? Either we use the SNESSetUseMatrixFree()
alone, then it works using "basic" finite-differencing, or we use the
SNESSetUseMatrixFree + MatFDColoringCreate
& SNESComputeJacobianDefaultColor as an option to SNESSetJacobian to access
the finite-differencing based on coloring. Is that right ?
Then if i come back to my preconditioner problem ... once you have set-up
the implicit matrix with one or the other aforementioned matrix-free ways,
how would you go around setting up the preconditioner ? In a matrix-free
way too, or rather as a real matrix that we assemble ourselves this time,
as you seemed to mean with the previous MatAij DMCreateMatrix ?
Sorry if it seems like I am nagging, but I would really like to understand
how to manipulate the matrix-free methods and structures in PETSc to run a
time-implicit finite volume computation, it's so promising !
Thanks again,
Thibault
> Barry
>
> Internally it uses SNESComputeJacobianDefaultColor() if you are interested
> in what it does.
>
>
>
>
>
> 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.
>>
>> <Untitled.png>
>>
>> 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/51771727/attachment-0001.html>
More information about the petsc-users
mailing list