On Sat, Mar 12, 2016 at 3:19 PM, Barry Smith <[email protected]> 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 <[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
