On Sat, Mar 12, 2016 at 8:34 PM, Emil Constantinescu <[email protected]> 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 <[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 >>> >> >> -- 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
