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 <[email protected]> Date: Saturday, June 13, 2015 at 9:01 AM To: barry smith <[email protected]> Cc: petsc-dev mailing list <[email protected]> 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 [[email protected]] >Sent: Saturday, 13 June 2015 12:14 p.m. >To: John O'Sullivan >Cc: [email protected] >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 >><[email protected]> 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 >>
