[petsc-users] PetscFV and TS implicit

[Matthew Knepley]

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.

  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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200821/d628c7be/attachment-0001.html>