<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">Message: 4<br>
Date: Sun, 14 Jun 2015 16:47:02 +0000<br>
From: "Abhyankar, Shrirang G." <<a href="mailto:abhyshr@anl.gov">abhyshr@anl.gov</a>><br>
To: John O'Sullivan <<a href="mailto:jp.osullivan@auckland.ac.nz">jp.osullivan@auckland.ac.nz</a>>, "Smith, Barry F."<br>
<<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>><br>
Cc: "<a href="mailto:petsc-dev@mcs.anl.gov">petsc-dev@mcs.anl.gov</a>" <<a href="mailto:petsc-dev@mcs.anl.gov">petsc-dev@mcs.anl.gov</a>><br>
Subject: Re: [petsc-dev] Missing Fortran interfaces<br>
Message-ID: <<a href="mailto:D1A2F17C.230C9%25abhyshr@mcs.anl.gov">D1A2F17C.230C9%abhyshr@mcs.anl.gov</a>><br>
Content-Type: text/plain; charset="Windows-1252"<br>
<br>
I am wondering if TSEvent would be of help here. TSEvent, available in<br>
PETSc 3.6, allows one to use TS with discontinuous right hand sides, which<br>
is what I understand is the case here with switching between 3 regions.<br>
The switching logic is described by ?event? functions and TSEvent detects<br>
and locates the zero crossing of such events.<br></blockquote><div><br></div><div>I had a similar thought about event functions; glad to see they're implemented in PETSc 3.6!</div><div><br></div><div>Although I can't say I'm familiar with the details of PFLOTRAN or Dr. O'Sullivan's applications, I would think that the switching function here would be crossing the spinodal manifold, due to the discussion of phase changes. In a single-phase system, it's usually when the Gibbs energy is at an inflection point, so I'd guess that more generally, you would have to determine when the Hessian of the Gibbs free energy is singular.</div><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<br>
A brief description of TSEvent can be found in section 6.1.7 of the manual<br>
and there are a few examples available<br>
$<br>
<<a href="http://www.mcs.anl.gov/petsc/petsc-3.6/src/ts/examples/tutorials/ex40.c.ht
ml" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/petsc-3.6/src/ts/examples/tutorials/ex40.c.ht<br>
ml</a>>PETSC_DIR/src/ts/examples/tutorials/ex40.c<br>
$PETSC_DIR/src/ts/examples/tutorials/ex41.c<br>
$PETSC_DIR/src/ts/examples/tests/ex22.c<br>
$PETSC_DIR/src/ts/examples/tutorials/power_grid/ex3adj_events.c<br>
<br>
Shri<br>
<br>
-----Original Message-----<br>
From: John O'Sullivan <<a href="mailto:jp.osullivan@auckland.ac.nz">jp.osullivan@auckland.ac.nz</a>><br>
Date: Saturday, June 13, 2015 at 9:01 AM<br>
To: barry smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>><br>
Cc: petsc-dev mailing list <<a href="mailto:petsc-dev@mcs.anl.gov">petsc-dev@mcs.anl.gov</a>><br>
Subject: Re: [petsc-dev] Missing Fortran interfaces<br>
<br>
>Hi Barry,<br>
><br>
>At this stage we are planning to implement the same variable switching<br>
>algorithm that TOUGH2 uses. In the future we'll be working on different<br>
>algorithms but hopefully they will work with much the same Petsc<br>
>architecture.<br>
><br>
>Having spent 2-3 years working on the TOUGH2 phase change/variable<br>
>switching algorithm its probably easiest if I explain what (AU)TOUGH2<br>
>does.<br>
><br>
>I should say that we've had a lot of success improving the stability and<br>
>convergence of the TOUGH2 algorithm in the last couple of years and our<br>
>version of TOUGH2 (AUTOUGH2) is very reliable now. In fact the code we<br>
>are developing with Petsc is initially going to be based on the AUTOUGH2<br>
>algorithm but will be leveraging the power of Petsc in terms of<br>
>parallelisation and linear solvers.<br>
><br>
>Let's consider a simple pure water example. Sorry if this is too<br>
>simplified or too confusing...<br>
><br>
>There are 3 regions: liquid, gas and 2-phase.<br>
><br>
>As Peter said there are lots of options for primary variables in each<br>
>region but lets stick to the TOUGH2 EOS1 choices. In the liquid and gas<br>
>regions the primary variables are pressure and temperature which are<br>
>independent. In the 2-phase region primary variables are pressure and gas<br>
>saturation.<br>
><br>
><br>
>Here's the basic algorithm for a time step. In TOUGH2 it is working to<br>
>determine the increment to the primary variables as a result of the time<br>
>step.:<br>
><br>
>a. Check the primary variables to see if a phase change has occurred<br>
>b. Change variables if required and update their increment<br>
>c. Calculate the residuals for each block<br>
>e. Calculate derivatives for each block to construct the Jacobian<br>
>f. If the residuals are all lower than the convergence criteria then<br>
>update the primary variables using their increments and move on to the<br>
>next time step.<br>
>g. If not then solve the system for an update to the increment and return<br>
>to (a) if the number of Newton steps is less than the maximum otherwise<br>
>reduce the timestep size and start again at (a).<br>
><br>
>In TOUGH2 you we don't iterate on the variable change because the<br>
>adaptive timestepping takes care of this by reducing the timestep size<br>
>until an acceptable variable change takes place.<br>
><br>
>What you're really after is the detail of what happens in (b). I have the<br>
>arcane Fortran77 TOUGH2 code in front of me so I'll try to summarise it:<br>
><br>
>Let Xn be the current primary variables and DXn+1 be their increment for<br>
>the timestep.<br>
><br>
>b.1. Liquid to 2-phase<br>
>- detected by comparing the saturation pressure (Ps(Tc)) at the current<br>
>temperature (Tc) with the current pressure (Pc)<br>
>- if Pc < Ps then change to 2-phase (gas phase evolves)<br>
>- X: (Pn, Tn) => X: (Pc, Sc)<br>
>- Sc = small (1e-6)<br>
>- Pc = Ps(Tc)<br>
>- DXn+1 = (Pc-Pn, Sc-Tn)<br>
><br>
><br>
>b.2. Gas to 2-phase<br>
>- detected by comparing the saturation pressure (Ps(Tc)) at the current<br>
>temperature (Tc) with the current pressure (Pc)<br>
>- if Pc > Ps then change to 2-phase (liquid phase evolves)<br>
>- X: (Pn, Tn) => X: (Pc, Sc)<br>
>- Sc = big (1.0-1e-6)<br>
>- Pc = Ps(Tc)<br>
>- DXn+1 = (Pc-Pn, Sc-Tn)<br>
><br>
><br>
>b.3. 2-phase to Liquid<br>
>- detected by comparing the current gas saturation (Sc) to 0.0<br>
>- if Sc < 0.0 then change to liquid (gas phase disappears)<br>
>- X: (Pn, Sn) => X: (Pc, Tc)<br>
>- Tc = Ts(Pc-1) (the saturation temperature<br>
>from the previous nonlinear solve iteration)<br>
>- Pc = Ps(Tc)*(1.0 + small) (just above the saturation pressure)<br>
>- DXn+1 = (Pc-Pn, Tc-Sn)<br>
><br>
><br>
>b.4. 2-phase to Gas<br>
>- detected by comparing the current gas saturation (Sc) to 1.0<br>
>- if Sc > 1.0 then change to gas (liquid phase disappears)<br>
>- X: (Pn, Sn) => X: (Pc, Tc)<br>
>- Tc = Ts(Pc-1) (the saturation temperature<br>
>from the previous nonlinear solve iteration)<br>
>- Pc = Ps(Tc)*(1.0 - small) (just below the saturation pressure)<br>
>- DXn+1 = (Pc-Pn, Tc-Sn)<br>
><br>
><br>
>In TOUGH2 the primary variables are used to calculate secondary variables<br>
>and then both are fed into the balance equations used in the nonlinear<br>
>solve. The components of the Jacobian are constructed explicitly using a<br>
>finite difference approach. This uses the primary variables with a tiny<br>
>increment and calculates the corresponding secondary variables.<br>
><br>
>The key to successful variable swapping in TOUGH2 is ensuring that the<br>
>primary variables and the secondary variables are continuous across a<br>
>variable swap.<br>
><br>
>Also in TOUGH2 the balance equations themselves are the same in all<br>
>regions (2-phase equations) so no logic is required to change them when<br>
>variables are switched. All that happens is that the X vector corresponds<br>
>to different terms in the same equations.<br>
><br>
>ie F(P,T,S,rho,nu,h?.) = 0<br>
><br>
>One thing to note is that during each iteration of the nonlinear solve<br>
>TOUGH2 is solving for the update to the increment DXn+1 so that<br>
>F(Xn+1,Xn) = 0. As b.1-b.4 above show DXn+1 gets recalculated to<br>
>accommodate each variable switch and the nonlinear solve subsequently<br>
>starts adjusting DXn+1 for the new variable set.<br>
><br>
>How I was planning to implement it in Petsc was basically the same:<br>
><br>
>In the SNESSetUpdate function:<br>
>1. Check which region each block is in (using the logic in b.1-b.4)<br>
>2. If the region is different than the previous iteration then change<br>
>variables by updating X<br>
><br>
>In the RHSFuntion:<br>
>1. Calculate all the fluid properties (equivalent to primary and<br>
>secondary variables in TOUGH2) using X. This will be done differently for<br>
>each region<br>
>2. Calculate and return the residuals of the balance equations using the<br>
>fluid properties.<br>
><br>
>There are lots of unanswered questions though. As I said in TOUGH2 we<br>
>manually adjust DXn+1 where in Petsc this will be internal to SNESSolve<br>
>so I don?t know:<br>
><br>
>a. if we need to adjust it<br>
>b. if we do, then how we would.<br>
><br>
>Also for calculating the Jacobian we were in the first instance going to<br>
>leave this up to Petsc so I?m sure of what the implications of the<br>
>variable switching will be. There are some subtle things about the<br>
>balance equations that may mean that the Jacobian needs to be done<br>
>explicitly. As I said before the key to good phase changing in TOUGH2 is<br>
>continuous fluid properties across the variable switch. However the<br>
>gradients of these properties which feed into the coefficients of the<br>
>Jacobian are almost always not continuous. This obviously will affect the<br>
>topography of the objective function so I don?t know what it would do to<br>
>the Petsc SNESSolve? In TOUGH2 essentially the timestep is reduced to the<br>
>point that the gradients look sufficiently continuous for the nonlinear<br>
>solve to converge.<br>
><br>
>Finally there are lots of situations where calculations of the fluid<br>
>properties etc fail because the newton step has taken you outside the<br>
>range of certain equations' validity. In these circumstances we need to<br>
>terminate the SNESSolve prematurely and reduce the timestep. I don't know<br>
>how this should be done in Petsc?<br>
><br>
>My approach was going to be very much an engineers approach. I was going<br>
>to try things with the tools we have and see if it worked. If not, I was<br>
>going to try to work out why not? We have a lot of good benchmarking<br>
>tests now for phase changing and variable switching from TOUGH2 so any<br>
>solution we come up with will be put thoroughly through its trials.<br>
><br>
>OK, sorry for the essay. I hope this helps explains where we?re at and<br>
>what we?re trying to achieve.<br>
><br>
>Thanks very much for all your help.<br>
><br>
>Cheers<br>
>John<br>
><br>
>--<br>
>Dr John O'Sullivan<br>
>Lecturer<br>
>Department of Engineering Science<br>
>University of Auckland, New Zealand<br>
>email: jp.osullivan at <a href="http://auckland.ac.nz" rel="noreferrer" target="_blank">auckland.ac.nz</a><br>
>tel: +64 (0)9 923 85353<br>
><br>
>________________________________________<br>
>From: Barry Smith [<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>]<br>
>Sent: Saturday, 13 June 2015 12:14 p.m.<br>
>To: John O'Sullivan<br>
>Cc: <a href="mailto:petsc-dev@mcs.anl.gov">petsc-dev@mcs.anl.gov</a><br>
>Subject: Re: [petsc-dev] Missing Fortran interfaces<br>
><br>
> John,<br>
><br>
> We can provide any Fortran stubs you need but before we do could you<br>
>please outline to us (in algebraic notation; not in reference to PETSc<br>
>code) how you plan to handle your "variable switching". That is, what is<br>
>the problem being solved, what are the constraints, when do you need to<br>
>switch variables etc?<br>
><br>
>I would say everyone who has tried to do variable switching in PETSc<br>
>before has failed because we (the PETSc developers) did not help them<br>
>enough and understand what is needed in the Newton solver. I suspect we<br>
>need to provide a bit more functionality in PETSc to make variable<br>
>switching a success and I'd like to understand what that is and add it to<br>
>the PETSc algorithms. Also if you know of any references that actually<br>
>explain variable switching in detail, that would be great, all I've found<br>
>are vague hand-wavy discussions that leave too many important details<br>
>unexplained.<br></blockquote><div><br></div><div>I've had similar thoughts about how to incorporate equation switching in a reduced order modeling context. There is much discussion of how to switch among multiple reduced order models on-the-fly while solving an ODE (or an optimization problem, or a nonlinear system of algebraic equations), but I think the solver infrastructure (possibly also the theory) really lags in this area because I don't think anyone has demonstrated compelling results with good performance. (I should note that Zahr and Farhat have some interesting results for the optimization problem case in <a href="http://arxiv.org/pdf/1407.7618v1.pdf">http://arxiv.org/pdf/1407.7618v1.pdf</a>, but because the reduced order model is applied to the PDE-solve in a nested-analysis-and-design approach, by treating the PDE state variables as an implicit function of the design variables, the "variable-switching" issue can be sidestepped somewhat.) If variable switching gets sorted out for phase changes, there could be potentially interesting research questions to sort out in other areas as well. </div><div><br></div><div>-- <br><div class="gmail_signature"><div dir="ltr">Geoffrey Oxberry, Ph.D., E.I.T.<br></div></div></div><div><a href="mailto:goxberry@gmail.com" target="_blank">goxberry@gmail.com</a> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
><br>
> Thanks<br>
><br>
> Barry<br>
><br>
>> On Jun 12, 2015, at 5:54 PM, John O'Sullivan<br>
>><<a href="mailto:jp.osullivan@auckland.ac.nz">jp.osullivan@auckland.ac.nz</a>> wrote:<br>
>><br>
>> Hi,<br>
>><br>
>> We're working on changing variables during an SNES solve so that we can<br>
>>handle phase changes.<br>
>><br>
>> I'm looking into using SNESSetUpdate but there doesn't seem to be a<br>
>>fortran interface. I get an error:<br>
>><br>
>> Undefined symbols for architecture x86_64:<br>
>> "_snessetupdate_", referenced from:<br>
>> _MAIN__ in phaseChange.o<br>
>><br>
>> I've looked in the interface files and can't find it either. I tried<br>
>>writing one myself but things get a bit too complicated to guess...<br>
>><br>
>> Can an interface be added please?<br>
>><br>
>> Is this the best to do variable switching? If not where would be the<br>
>>right place and are there Fortran interfaces? We will be changing from<br>
>>Pressure, Temperature to Pressure, Saturation etc. But for now I'm just<br>
>>working on a simple example code.<br>
>><br>
>> Cheers<br>
>> John<br>
>> --<br>
>> Dr John O'Sullivan<br>
>> Lecturer<br>
>> Department of Engineering Science<br>
>> University of Auckland, New Zealand<br>
>> email:<br>
>> jp.osullivan at <a href="http://auckland.ac.nz" rel="noreferrer" target="_blank">auckland.ac.nz</a><br>
>><br>
>> tel: +64 (0)9 923 85353<br>
>></blockquote></div><div class="gmail_signature"><div dir="ltr"><span></span></div></div>
</div></div>