[petsc-dev] Missing Fortran interfaces

[John O'Sullivan]

Hi all,

Yes the standard versions of TOUGH2 are well know for getting stuck iterating on phase changes. However by being extremely carefully about the continuity of the fluid properties we have largely overcome this problem. The code still iterates but once the time step is reduced sufficiently, the block changes phase correctly and the time step immediately returns up to its previous level.

And yes we also see the phase changes occurring over volumes including large number of blocks in our 3D simulations. But one of the advantages of using the adaptive timestepping to do the phase change iterating is that all the blocks are essentially being iterated at the same time. By the time the time step is small enough for the most problematic block to change phase, all the other blocks have successfully changed phase.

To be honest though, while we run large simulations from a geothermal reservoir model point of view, they're quite small compared to other applications (max ~200k blocks). This is because the parallel version of TOUGH2 was never really finished properly except for the nuclear storage application. So while our simulations contain a wide range of block sizes, including some quite small ones, they're not true multi-scale simulations.

That's why we're working on a Petsc code capable of the same sorts of simulations.

I'm optimistic that our adapted TOUGH2 phase change algorithm will work for much larger simulations also long as we're careful about the implementation.

Cheers
John


--
Dr John O'Sullivan
Lecturer
Department of Engineering Science
University of Auckland, New Zealand
email: jp.osullivan at auckland.ac.nz
tel: +64 (0)9 923 85353

________________________________________
From: Peter Lichtner [peter.lichtner at gmail.com]
Sent: Sunday, 14 June 2015 11:45 a.m.
To: John O'Sullivan
Cc: Barry Smith; petsc-dev at mcs.anl.gov
Subject: Re: [petsc-dev] Missing Fortran interfaces

One more comment: with pflotran we have really not been able to make use of HPC with multiphase simulations. The problem seems to be that as the grid size is reduced so also is the time step. This issue is related to the dimensionality of the problem: for 1D phase changes occur at a single cell during a time step roughly; for 2 and 3D problems phases changes occur over lines and areas, and surfaces and volumes, respectively. As the grid spacing is reduced a greater number of cells undergo a phase change within a time step. There is also the problem of oscillations where on one iteration a phase change occurs and on the next iteration it is changed back to the original state etc.

I would be curious if John observes similar problems.
…Peter

> On Jun 13, 2015, at 08:01, John O'Sullivan <jp.osullivan at auckland.ac.nz> wrote:
>
> Hi Barry,
>
> At this stage we are planning to implement the same variable switching algorithm that TOUGH2 uses. In the future we'll be working on different algorithms but hopefully they will work with much the same Petsc architecture.
>
> Having spent 2-3 years working on the TOUGH2 phase change/variable switching algorithm its probably easiest if I explain what (AU)TOUGH2 does.
>
> I should say that we've had a lot of success improving the stability and convergence of the TOUGH2 algorithm in the last couple of years and our version of TOUGH2 (AUTOUGH2) is very reliable now. In fact the code we are developing with Petsc is initially going to be based on the AUTOUGH2 algorithm but will be leveraging the power of Petsc in terms of parallelisation and linear solvers.
>
> Let's consider a simple pure water example. Sorry if this is too simplified or too confusing...
>
> There are 3 regions: liquid, gas and 2-phase.
>
> As Peter said there are lots of options for primary variables in each region but lets stick to the TOUGH2 EOS1 choices. In the liquid and gas regions the primary variables are pressure and temperature which are independent. In the 2-phase region primary variables are pressure and gas saturation.
>
>
> Here's the basic algorithm for a time step. In TOUGH2 it is working to determine the increment to the primary variables as a result of the time step.:
>
> a. Check the primary variables to see if a phase change has occurred
> b. Change variables if required and update their increment
> c. Calculate the residuals for each block
> e. Calculate derivatives for each block to construct the Jacobian
> f.  If the residuals are all lower than the convergence criteria then update the primary variables using their increments and move on to the next time step.
> g. If not then solve the system for an update to the increment and return to (a) if the number of Newton steps is less than the maximum otherwise reduce the timestep size and start again at (a).
>
> In TOUGH2 you we don't iterate on the variable change because the adaptive timestepping takes care of this by reducing the timestep size until an acceptable variable change takes place.
>
> What you're really after is the detail of what happens in (b). I have the arcane Fortran77 TOUGH2 code in front of me so I'll try to summarise it:
>
> Let Xn be the current primary variables and DXn+1 be their increment for the timestep.
>
> b.1. Liquid to 2-phase
> - detected by comparing the saturation pressure (Ps(Tc)) at the current temperature (Tc) with the current pressure (Pc)
> - if Pc < Ps then change to 2-phase (gas phase evolves)
> - X: (Pn, Tn) => X: (Pc, Sc)
> - Sc = small (1e-6)
> - Pc = Ps(Tc)
> - DXn+1 = (Pc-Pn, Sc-Tn)
>
>
> b.2. Gas to 2-phase
> - detected by comparing the saturation pressure (Ps(Tc)) at the current temperature (Tc) with the current pressure (Pc)
> - if Pc > Ps then change to 2-phase (liquid phase evolves)
> - X: (Pn, Tn) => X: (Pc, Sc)
> - Sc = big (1.0-1e-6)
> - Pc = Ps(Tc)
> - DXn+1 = (Pc-Pn, Sc-Tn)
>
>
> b.3. 2-phase to Liquid
> - detected by comparing the current gas saturation (Sc) to 0.0
> - if Sc < 0.0 then change to liquid (gas phase disappears)
> - X: (Pn, Sn) => X: (Pc, Tc)
> - Tc = Ts(Pc-1)                            (the saturation temperature from the previous nonlinear solve iteration)
> - Pc = Ps(Tc)*(1.0 + small)      (just above the saturation pressure)
> - DXn+1 = (Pc-Pn, Tc-Sn)
>
>
> b.4. 2-phase to Gas
> - detected by comparing the current gas saturation (Sc) to 1.0
> - if Sc > 1.0 then change to gas (liquid phase disappears)
> - X: (Pn, Sn) => X: (Pc, Tc)
> - Tc = Ts(Pc-1)                            (the saturation temperature from the previous nonlinear solve iteration)
> - Pc = Ps(Tc)*(1.0 - small)      (just below the saturation pressure)
> - DXn+1 = (Pc-Pn, Tc-Sn)
>
>
> In TOUGH2 the primary variables are used to calculate secondary variables and then both are fed into the balance equations used in the nonlinear solve. The components of the Jacobian are constructed explicitly using a finite difference approach. This uses the primary variables with a tiny increment and calculates the corresponding secondary variables.
>
> The key to successful variable swapping in TOUGH2 is ensuring that the primary variables and the secondary variables are continuous across a variable swap.
>
> Also in TOUGH2 the balance equations themselves are the same in all regions (2-phase equations) so no logic is required to change them when variables are switched. All that happens is that the X vector corresponds to different terms in the same equations.
>
> ie    F(P,T,S,rho,nu,h….) = 0
>
> One thing to note is that during each iteration of the nonlinear solve TOUGH2 is solving for the update to the increment DXn+1 so that F(Xn+1,Xn) = 0. As b.1-b.4 above show DXn+1 gets recalculated to accommodate each variable switch and the nonlinear solve subsequently starts adjusting DXn+1 for the new variable set.
>
> How I was planning to implement it in Petsc was basically the same:
>
> In the SNESSetUpdate function:
> 1.    Check which region each block is in (using the logic in b.1-b.4)
> 2.    If the region is different than the previous iteration then change variables by updating X
>
> In the RHSFuntion:
> 1.    Calculate all the fluid properties (equivalent to primary and secondary variables in TOUGH2) using X. This will be done differently for each region
> 2.    Calculate and return the residuals of the balance equations using the fluid properties.
>
> There are lots of unanswered questions though. As I said in TOUGH2 we manually adjust DXn+1 where in Petsc this will be internal to SNESSolve so I don’t know:
>
> a.    if we need to adjust it
> b.    if we do, then how we would.
>
> Also for calculating the Jacobian we were in the first instance going to leave this up to Petsc so I’m sure of what the implications of the variable switching will be. There are some subtle things about the balance equations that may mean that the Jacobian needs to be done explicitly. As I said before the key to good phase changing in TOUGH2 is continuous fluid properties across the variable switch. However the gradients of these properties which feed into the coefficients of the Jacobian are almost always not continuous. This obviously will affect the topography of the objective function so I don’t know what it would do to the Petsc SNESSolve? In TOUGH2 essentially the timestep is reduced to the point that the gradients look sufficiently continuous for the nonlinear solve to converge.
>
> Finally there are lots of situations where calculations of the fluid properties etc fail because the newton step has taken you outside the range of certain equations' validity. In these circumstances we need to terminate the SNESSolve prematurely and reduce the timestep. I don't know how this should be done in Petsc?
>
> My approach was going to be very much an engineers approach. I was going to try things with the tools we have and see if it worked. If not, I was going to try to work out why not… We have a lot of good benchmarking tests now for phase changing and variable switching from TOUGH2 so any solution we come up with will be put thoroughly through its trials.
>
> OK, sorry for the essay. I hope this helps explains where we’re at and what we’re trying to achieve.
>
> Thanks very much for all your help.
>
> Cheers
> John
>
> --
> Dr John O'Sullivan
> Lecturer
> Department of Engineering Science
> University of Auckland, New Zealand
> email: jp.osullivan at auckland.ac.nz
> tel: +64 (0)9 923 85353
>
> ________________________________________
> From: Barry Smith [bsmith at mcs.anl.gov]
> Sent: Saturday, 13 June 2015 12:14 p.m.
> To: John O'Sullivan
> Cc: petsc-dev at mcs.anl.gov
> Subject: Re: [petsc-dev] Missing Fortran interfaces
>
>  John,
>
>   We can provide any Fortran stubs you need but before we do could you please outline to us (in algebraic notation; not in reference to PETSc code) how you plan to handle your "variable switching". That is, what is the problem being solved, what are the constraints, when do you need to switch variables etc?
>
> I would say everyone who has tried to do variable switching in PETSc before has failed because we (the PETSc developers) did not help them enough and  understand what is needed in the Newton solver. I suspect we need to provide a bit more functionality in PETSc to make variable switching a success and I'd like to understand what that is and add it to the PETSc algorithms. Also if you know of any references that actually explain variable switching in detail, that would be great, all I've found are vague hand-wavy discussions that leave too many important details unexplained.
>
>  Thanks
>
>  Barry
>
>> On Jun 12, 2015, at 5:54 PM, John O'Sullivan <jp.osullivan at auckland.ac.nz> wrote:
>>
>> Hi,
>>
>> We're working on changing variables during an SNES solve so that we can handle phase changes.
>>
>> I'm looking into using SNESSetUpdate but there doesn't seem to be a fortran interface. I get an error:
>>
>> Undefined symbols for architecture x86_64:
>>  "_snessetupdate_", referenced from:
>>      _MAIN__ in phaseChange.o
>>
>> I've looked in the interface files and can't find it either. I tried writing one myself but things get a bit too complicated to guess...
>>
>> Can an interface be added please?
>>
>> Is this the best to do variable switching? If not where would be the right place and are there Fortran interfaces? We will be changing from Pressure, Temperature to Pressure, Saturation etc. But for now I'm just working on a simple example code.
>>
>> Cheers
>> John
>> --
>> Dr John O'Sullivan
>> Lecturer
>> Department of Engineering Science
>> University of Auckland, New Zealand
>> email:
>> jp.osullivan at auckland.ac.nz
>>
>> tel: +64 (0)9 923 85353
>>
>