[petsc-users] PetscFV and TS implicit

Fri Aug 21 08:10:00 CDT 2020

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.

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