[petsc-dev] Missing Fortran interfaces

[Abhyankar, Shrirang G.]

I am wondering if TSEvent would be of help here. TSEvent, available in
PETSc 3.6, allows one to use TS with discontinuous right hand sides, which
is what I understand is the case here with switching between 3 regions.
The switching logic is described by ³event² functions and TSEvent detects
and locates the zero crossing of such events.

A brief description of TSEvent can be found in section 6.1.7 of the manual
and there are a few examples available
$ 
<http://www.mcs.anl.gov/petsc/petsc-3.6/src/ts/examples/tutorials/ex40.c.ht
ml>PETSC_DIR/src/ts/examples/tutorials/ex40.c
$PETSC_DIR/src/ts/examples/tutorials/ex41.c
$PETSC_DIR/src/ts/examples/tests/ex22.c
$PETSC_DIR/src/ts/examples/tutorials/power_grid/ex3adj_events.c

Shri

-----Original Message-----
From: John O'Sullivan <jp.osullivan at auckland.ac.nz>
Date: Saturday, June 13, 2015 at 9:01 AM
To: barry smith <bsmith at mcs.anl.gov>
Cc: petsc-dev mailing list <petsc-dev at mcs.anl.gov>
Subject: Re: [petsc-dev] Missing Fortran interfaces

>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
>>