[petsc-users] Using TS

[Emil Constantinescu]

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.

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
>