[petsc-users] Using TS

[Matthew Knepley]

On Sat, Mar 12, 2016 at 8:34 PM, Emil Constantinescu <emconsta at mcs.anl.gov>
wrote:

> I also find it useful to go through one of the simple examples available
> for TS:
> http://www.mcs.anl.gov/petsc/petsc-current/src/ts/examples/tutorials/index.html
> (ex8 may be a good start).
>
> As Barry suggested, you need to implement IFunction and IJacobian. The
> argument "u" is  S_o, S_w, and p stacked together and "u_t" their
> corresponding time derivatives. The rest is calculus, but following an
> example usually helps a lot in the beginning.
>

Are you guys saying that IFunction and IJacobian are enough to do the
adjoint system as well?

  Matt


> Out of curiosity, what is the application?
>
> Emil
>
>
> On 3/12/16 3:19 PM, Barry Smith 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 )
>>
>>
>>    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/e81300c3/attachment.html>