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 )
Barry
> On Mar 12, 2016, at 12:04 PM, Matthew Knepley <[email protected]> wrote:
>
> On Sat, Mar 12, 2016 at 5:41 AM, Max la Cour Christensen <[email protected]> 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
