[petsc-users] Using TS

[Matthew Knepley]

On Sat, Mar 12, 2016 at 3:19 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:

>
>    This is only a starting point, Jed and Emil can fix my mistakes and
> provide additional details.
>
>     In your case you will not provide a TSSetRHSFunction and
> TSSetRHSJacobian since everything should be treated implicitly as a DAE.
>
>     First move everything in the three equations to the left side and then
> differentiate through the \partial/\partial t so that it only applies to
> the S_o, S_w, and p. For example for the first equation using the product
> rule twice you get something like
>
>     \phi( p ) \rho_o( p ) \partial S_o/ \partial t + phi( p ) S_o
> \partial \rho_o( p )  \partial t  + \rho_o( p ) S_o \partial \phi( p )
> \partial t - F_o = 0
>
>     \phi( p ) \rho_o( p ) \partial S_o/ \partial t  + phi( p ) S_o
> \rho_o'(p)  \partial p \partial t + \rho_o( p ) S_o \phi'( p ) \partial  p
>  \partial t  - F_o = 0
>
> The two vector arguments to your IFunction are exactly the S_o, S_w, and p
> and \partial S_o/ \partial t ,  \partial S_w/ \partial t, and  \partial p/
> \partial t so it is immediate to code up your IFunction once you have the
> analytic form above
>
> For the IJacobian and the "shift business" just remember that dF/dU means
> take the derivative of the IFunction with respect to S_o, S_w, and p
> treating the \partial S_o/ \partial t ,  \partial S_w/ \partial t, and
> \partial p/ \partial t as if they were independent of S_o, S_w, and p. For
> the dF/dU_t that means taking the derivate with respect to the \partial
> S_o/ \partial t ,  \partial S_w/ \partial t, and  \partial p/ \partial t
> treating the S_o, S_w, and p as independent of \partial S_o/ \partial t ,
> \partial S_w/ \partial t, and  \partial p/ \partial t.   Then you just need
> to form the sum of the two parts with the a "shift" scaling dF/dU +
> a*dF/dU_t
>
> For the third equation everything is easy. dF/dS_o = 1 dF/dS_w = 1 dF/dp =
> 0  dF/d (\partial S_o)/\partial t = 0  (\partial S_w)/\partial t = 0
> (\partial p)/\partial t = 0
>
> Computations for the first two equations are messy but straightforward.
> For example for the first equation dF/dS_o = phi( p ) \rho_o'(p)  \partial
> p \partial t  + \rho_o( p ) \phi'( p ) \partial  p + dF_o/dS_o  and dF/d
> (\partial S_o)/\partial t) = \phi( p ) \rho_o( p )


Max and Stefan,

You can see me trying to do this in a more generic sense here:


https://bitbucket.org/petsc/petsc/src/f0b116324093eeda71fbbb2872e1bdb3483ad365/src/snes/utils/dmplexsnes.c?at=master&fileviewer=file-view-default#dmplexsnes.c-1678

This is intended to work with both FEM and FVM, which makes it look a
little messier
than I would like right now.

  Thanks,

    Matt


>
>   Barry
>
> > On Mar 12, 2016, at 12:04 PM, Matthew Knepley <knepley at gmail.com> wrote:
> >
> > On Sat, Mar 12, 2016 at 5:41 AM, Max la Cour Christensen <mlcch at dtu.dk>
> wrote:
> >
> > Hi guys,
> >
> > We are making preparations to implement adjoint based optimisation in
> our in-house oil and gas reservoir simulator. Currently our code uses
> PETSc's DMPlex, Vec, Mat, KSP and PC. We are still not using SNES and TS,
> but instead we have our own backward Euler and Newton-Raphson
> implementation. Due to the upcoming implementation of adjoints, we are
> considering changing the code and begin using TS and SNES.
> >
> > After examining the PETSc manual and examples, we are still not
> completely clear on how to apply TS to our system of PDEs. In a simplified
> formulation, it can be written as:
> >
> > \partial( \phi( p ) \rho_o( p ) S_o )/ \partial t = F_o(p,S)
> > \partial( \phi( p ) \rho_w( p ) S_w )/ \partial t = F_w(p,S)
> > S_o + S_w = 1,
> >
> > where p is the pressure,
> > \phi( p ) is a porosity function depending on pressure,
> > \rho_x( p ) is a density function depending on pressure,
> > S_o is the saturation of oil,
> > S_g is the saturation of gas,
> > t is time,
> > F_x(p,S) is a function containing fluxes and source terms. The primary
> variables are p, S_o and S_w.
> >
> > We are using a lowest order Finite Volume discretisation.
> >
> > Now for implementing this in TS (with the prospect of later using
> TSAdjoint), we are not sure if we need all of the functions:
> TSSetIFunction, TSSetRHSFunction, TSSetIJacobian and TSSetRHSJacobian and
> what parts of the equations go where. Especially we are unsure of how to
> use the concept of a shifted jacobian (TSSetIJacobian).
> >
> > Any advice you could provide will be highly appreciated.
> >
> > Barry and Emil,
> >
> > I am also interested in this, since I don't know how to do it.
> >
> >   Thanks,
> >
> >      Matt
> >
> > Many thanks,
> > Max la Cour Christensen
> > PhD student, Technical University of Denmark
> >
> >
> >
> > --
> > 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
>
>


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