<html><head><meta http-equiv="Content-Type" content="text/html charset=iso-8859-1"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space;" class=""><br class=""><div><blockquote type="cite" class=""><div class="">On Mar 15, 2016, at 6:33 AM, Max la Cour Christensen <<a href="mailto:mlcch@dtu.dk" class="">mlcch@dtu.dk</a>> wrote:</div><br class="Apple-interchange-newline"><div class="">
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1" class="">
<div style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; " class="">
<div class=""><br class="">
</div>
<div class="">Barry, Matt, Emil, thanks for the explanations! They are very useful.</div>
<div class=""><br class="">
</div>
<div class=""><b class="">Emil:</b> The application is the flow of oil and gas in the subsurface. The code uses a model of the rock in the subsurface and given a setup of wells, it computes a prediction for the production of oil and gas. It can also be used to model CO2 storage
in the subsurface.</div>
<div class=""><br class="">
</div>
<div class="">I think we will attempt to switch to TS. We are working with different formulations of the equations. In the following, I have summarised how I think IFunction and IJacobian should be implemented for the formulation of the most interest to us. If you can
bare with me and confirm this that would be great. The equations are:</div>
<div class=""><br class="">
</div>
<div class="">dm_o/dt = f_o(m,p)</div>
<div class="">dm_g/dt = f_g(m,p)</div>
<div class="">S_o(m,p) + S_g(m,p) = 1,</div>
<div class=""><br class="">
</div>
<div class="">where</div>
<div class=""><br class="">
</div>
<div class="">m_o is the mass of oil</div>
<div class="">m_g is the mass of gas</div>
<div class="">t is time,</div>
<div class="">f(m_x,p) contains fluxes for mass transfer between adjacent cells and source terms for wells</div>
<div class="">S_o is the oil saturation</div>
<div class="">S_g is the gas saturation</div>
<div class="">p is the pressure.</div>
<div class=""><br class="">
</div>
<div class="">In this formulation, the primary variables are m_o, m_g and p. </div>
<div class=""><br class="">
</div>
<div class="">The above looks simpler than it is. The code to evaluate these things is around 90k lines of Fortran.</div>
<div class=""><br class="">
</div>
<div class="">If I understood correctly, IFunction and IJacobian would implement the following:</div>
<div class=""><br class="">
</div>
<div class=""><b class="">IFunction: </b></div>
<div class=""><br class="">
</div>
<div class="">f_o(m,p) - dm_o/dt</div>
<div class="">f_g(m,p) - dm_g/dt</div>
<div class="">S_o(m,p) + S_g(m,p) - 1,</div>
<div class=""><br class="">
</div>
<div class="">where dm_o/dt and dm_g/dt is given by the time derivative (u_t in PETSc) from IFunction()</div>
<div class=""><br class="">
</div>
<div class=""><b class="">IJacobian:</b></div>
<div class=""><br class="">
</div>
<div class="">First 2 equations: \partial f(m,p) / \partial m - a and \partial f(m,p) / \partial p</div>
<div class="">Last equation: \partial (S_o(m,p)+S_g(m,p)-1) / \partial m and \partial (S_o(m,p)+S_g(m,p)-1) / \partial p,</div>
<div class=""><br class="">
</div>
<div class="">where a comes from PETSc IJacobian().</div>
<div class=""><br class="">
</div>
<div class="">Is this correct usage?</div>
<div class=""><br class="">
</div>
<div class=""><b class="">Emil: </b>Now for the adjoint implementation, we will need to provide a function for TSAdjointSetRHSJacobian, which computes the derivative of our system wrt. our control variables (control variables being the p in the PETSc manual) and some additional
cost stuff: TSSetCostGradients() and possibly TSSetCostIntegrand(). Is this correct? Our cost function would typically be to maximize profit of an oil reservoir. In a simple way, this could be to maximize the number of oil barrels produced over the lifetime
of the reservoir (e.g. 50 years). In the equations above, the number of oil barrels is contained in the source terms for the wells.</div></div></div></blockquote><div><br class=""></div><div>If your cost function depends on the intermediate states of the simulation, e.g. in an integral form </div><div><br class=""></div><div>\int F(m_o,m_g,p;c) dt, </div><div><br class=""></div>where c is your control variable, then the integrand F() needs to be provided with TSSetCostIntegrand() so that PETSc can evaluate the cost function and its sensitivities to the control variables. If the cost function depends only on the final states, e.g. the values of m_o,m_g,p at the end of simulation, this function is not needed. By the way, is the number of oil barrels represented in integer or real values?</div><div><br class=""></div><div>Hong </div><div><br class=""><blockquote type="cite" class=""><div class=""><div style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; " class="">
<div class="">Many thanks,</div>
<div class="">Max</div>
<div class=""><br class="">
</div>
<br class="">
<div class="">
<div class="">On Mar 13, 2016, at 4:02 AM, Emil Constantinescu <<a href="mailto:emconsta@mcs.anl.gov" class="">emconsta@mcs.anl.gov</a>></div>
<div class=""> wrote:</div>
<br class="Apple-interchange-newline">
<blockquote type="cite" class="">On 3/12/16 8:37 PM, Matthew Knepley wrote:<br class="">
<blockquote type="cite" class="">On Sat, Mar 12, 2016 at 8:34 PM, Emil Constantinescu<br class="">
<<a href="mailto:emconsta@mcs.anl.gov" class="">emconsta@mcs.anl.gov</a> <<a href="mailto:emconsta@mcs.anl.gov" class="">mailto:emconsta@mcs.anl.gov</a>>> wrote:<br class="">
<br class="">
I also find it useful to go through one of the simple examples<br class="">
available for TS:<br class="">
<a href="http://www.mcs.anl.gov/petsc/petsc-current/src/ts/examples/tutorials/index.html" class="">http://www.mcs.anl.gov/petsc/petsc-current/src/ts/examples/tutorials/index.html</a><br class="">
(ex8 may be a good start).<br class="">
<br class="">
As Barry suggested, you need to implement IFunction and IJacobian.<br class="">
The argument "u" is S_o, S_w, and p stacked together and "u_t"<br class="">
their corresponding time derivatives. The rest is calculus, but<br class="">
following an example usually helps a lot in the beginning.<br class="">
<br class="">
<br class="">
Are you guys saying that IFunction and IJacobian are enough to do the<br class="">
adjoint system as well?<br class="">
<br class="">
</blockquote>
<br class="">
Pretty much yes, but it depends on the cost function. This is the beauty of discrete adjoints - if you have the Jacobian (transpose, done internally through KSP) you're done. You need IJacobian for sure to do the backward propagation. If you have that, the
rest is usually trivial. Mr. Hong Zhang (my postdoc) set up quite a few simple examples.<br class="">
<br class="">
Emil<br class="">
<br class="">
<blockquote type="cite" class=""> Matt<br class="">
<br class="">
Out of curiosity, what is the application?<br class="">
<br class="">
Emil<br class="">
<br class="">
<br class="">
On 3/12/16 3:19 PM, Barry Smith wrote:<br class="">
<br class="">
<br class="">
This is only a starting point, Jed and Emil can fix my<br class="">
mistakes and provide additional details.<br class="">
<br class="">
<br class="">
In your case you will not provide a TSSetRHSFunction and<br class="">
TSSetRHSJacobian since everything should be treated implicitly<br class="">
as a DAE.<br class="">
<br class="">
First move everything in the three equations to the left<br class="">
side and then differentiate through the \partial/\partial t so<br class="">
that it only applies to the S_o, S_w, and p. For example for the<br class="">
first equation using the product rule twice you get something like<br class="">
<br class="">
\phi( p ) \rho_o( p ) \partial S_o/ \partial t + phi( p )<br class="">
S_o \partial \rho_o( p ) \partial t + \rho_o( p ) S_o<br class="">
\partial \phi( p ) \partial t - F_o = 0<br class="">
<br class="">
\phi( p ) \rho_o( p ) \partial S_o/ \partial t + phi( p )<br class="">
S_o \rho_o'(p) \partial p \partial t + \rho_o( p ) S_o \phi'( p<br class="">
) \partial p \partial t - F_o = 0<br class="">
<br class="">
The two vector arguments to your IFunction are exactly the S_o,<br class="">
S_w, and p and \partial S_o/ \partial t , \partial S_w/<br class="">
\partial t, and \partial p/ \partial t so it is immediate to<br class="">
code up your IFunction once you have the analytic form above<br class="">
<br class="">
For the IJacobian and the "shift business" just remember that<br class="">
dF/dU means take the derivative of the IFunction with respect to<br class="">
S_o, S_w, and p treating the \partial S_o/ \partial t ,<br class="">
\partial S_w/ \partial t, and \partial p/ \partial t as if they<br class="">
were independent of S_o, S_w, and p. For the dF/dU_t that means<br class="">
taking the derivate with respect to the \partial S_o/ \partial t<br class="">
, \partial S_w/ \partial t, and \partial p/ \partial t<br class="">
treating the S_o, S_w, and p as independent of \partial S_o/<br class="">
\partial t , \partial S_w/ \partial t, and \partial p/<br class="">
\partial t. Then you just need to form the sum of the two<br class="">
parts with the a "shift" scaling dF/dU + a*dF/dU_t<br class="">
<br class="">
For the third equation everything is easy. dF/dS_o = 1 dF/dS_w =<br class="">
1 dF/dp = 0 dF/d (\partial S_o)/\partial t = 0 (\partial<br class="">
S_w)/\partial t = 0 (\partial p)/\partial t = 0<br class="">
<br class="">
Computations for the first two equations are messy but<br class="">
straightforward. For example for the first equation dF/dS_o =<br class="">
phi( p ) \rho_o'(p) \partial p \partial t + \rho_o( p ) \phi'(<br class="">
p ) \partial p + dF_o/dS_o and dF/d (\partial S_o)/\partial t)<br class="">
= \phi( p ) \rho_o( p )<br class="">
<br class="">
<br class="">
Barry<br class="">
<br class="">
On Mar 12, 2016, at 12:04 PM, Matthew Knepley<br class="">
<<a href="mailto:knepley@gmail.com" class="">knepley@gmail.com</a> <<a href="mailto:knepley@gmail.com" class="">mailto:knepley@gmail.com</a>>> wrote:<br class="">
<br class="">
On Sat, Mar 12, 2016 at 5:41 AM, Max la Cour Christensen<br class="">
<<a href="mailto:mlcch@dtu.dk" class="">mlcch@dtu.dk</a> <<a href="mailto:mlcch@dtu.dk" class="">mailto:mlcch@dtu.dk</a>>> wrote:<br class="">
<br class="">
Hi guys,<br class="">
<br class="">
We are making preparations to implement adjoint based<br class="">
optimisation in our in-house oil and gas reservoir<br class="">
simulator. Currently our code uses PETSc's DMPlex, Vec, Mat,<br class="">
KSP and PC. We are still not using SNES and TS, but instead<br class="">
we have our own backward Euler and Newton-Raphson<br class="">
implementation. Due to the upcoming implementation of<br class="">
adjoints, we are considering changing the code and begin<br class="">
using TS and SNES.<br class="">
<br class="">
After examining the PETSc manual and examples, we are still<br class="">
not completely clear on how to apply TS to our system of<br class="">
PDEs. In a simplified formulation, it can be written as:<br class="">
<br class="">
\partial( \phi( p ) \rho_o( p ) S_o )/ \partial t = F_o(p,S)<br class="">
\partial( \phi( p ) \rho_w( p ) S_w )/ \partial t = F_w(p,S)<br class="">
S_o + S_w = 1,<br class="">
<br class="">
where p is the pressure,<br class="">
\phi( p ) is a porosity function depending on pressure,<br class="">
\rho_x( p ) is a density function depending on pressure,<br class="">
S_o is the saturation of oil,<br class="">
S_g is the saturation of gas,<br class="">
t is time,<br class="">
F_x(p,S) is a function containing fluxes and source terms.<br class="">
The primary variables are p, S_o and S_w.<br class="">
<br class="">
We are using a lowest order Finite Volume discretisation.<br class="">
<br class="">
Now for implementing this in TS (with the prospect of later<br class="">
using TSAdjoint), we are not sure if we need all of the<br class="">
functions: TSSetIFunction, TSSetRHSFunction, TSSetIJacobian<br class="">
and TSSetRHSJacobian and what parts of the equations go<br class="">
where. Especially we are unsure of how to use the concept of<br class="">
a shifted jacobian (TSSetIJacobian).<br class="">
<br class="">
Any advice you could provide will be highly appreciated.<br class="">
<br class="">
Barry and Emil,<br class="">
<br class="">
I am also interested in this, since I don't know how to do it.<br class="">
<br class="">
Thanks,<br class="">
<br class="">
Matt<br class="">
<br class="">
Many thanks,<br class="">
Max la Cour Christensen<br class="">
PhD student, Technical University of Denmark<br class="">
<br class="">
<br class="">
<br class="">
--<br class="">
What most experimenters take for granted before they begin<br class="">
their experiments is infinitely more interesting than any<br class="">
results to which their experiments lead.<br class="">
-- Norbert Wiener<br class="">
<br class="">
<br class="">
<br class="">
<br class="">
<br class="">
--<br class="">
What most experimenters take for granted before they begin their<br class="">
experiments is infinitely more interesting than any results to which<br class="">
their experiments lead.<br class="">
-- Norbert Wiener<br class="">
</blockquote>
</blockquote>
</div>
<br class="">
</div>
</div></blockquote></div><br class=""></body></html>