Incorporate PETSc's TSEvent functionality to handle discrete events

[abhyshr]

GridPACK's dynamics simulation consists of a lot of discrete events, such as fault on/off, line tripping, generator/exciter/turbine-governor limiters, relays, etc.. To handle this correctly with the DAE implicit integration approach implemented in GridPACK, PETSc's TSEvent functionality needs to incorporated within GridPACK. This will involve adding the TSEvent API to the DAE solver, adding methods to DAE models to set event and postevent functions, and assembling/mapping the events from bus, generator, etc., models to the DAE solver interface.

[wperkins]

@abhyshr, I have some questions:

• Do "events" only (or primarily) come from individual network component states?
• If events are generated outside of network components, what does that look like?
  • Can an "event" generated by a network component be solely defined and handled by that network component? In other words, can a bus or branch generate an "event" value based solely on it's own state and that of it's immediate neighbors?
  • If yes, computation of "event" values by a bus/branch only needs to access it's own part of the solution vector, which can be managed.
  • If no, things get really complicated.

[wperkins]

If events can be checked and handled solely by a network component, I suggest an implementation something like this:

• The DAE solver would (optionally) have an event manager.
• Event descriptions would be identified when the problem is set up (matrix/vector sizing).
• A bus/branch needs to identify if it has some event(s) that need to be monitored and how many values are involved.
• If so, a bus/branch produces a single "event checker/handler" thing. This class would have some methods:
  • to identify how many event values it produces/checks
  • to produce and fill an array with those event values
    • This would need to read solution vector values in some form.
    • this would require local index information from the mapper
  • to take action after an event happens
    • this would need to be able to modify solution vector values
    • this would require local information from the mapper
  • access the component's part of the solution vector
• The DAE event manager might do the following:
  • collect and store event checkers from buses and branches
  • maybe correlate local solution vector indexes with the event checker
  • construct appropriate event descriptions for PETSc TS
  • provide a single event checker routine for PETSc TS
    • Just loops through the local bus/branch checkers and fills an event value vector at the appropriate indexes
  • provide a post event handler to PETSc TS
    • Using event index, the appropriate bus/branch post event handler is called.

Thoughts?

[abhyshr]

• Do "events" only (or primarily) come from individual network component states?

Primarily, they would be come from network components, e.g. some exciter state hitting its limit, or some bus voltage falls below certain threshold. These are events for network components only. There can also be "system" or "sub-system" level events, e.g. the total line flow over a set of lines exceeds some threshold or the average frequency deviation for a bunch of generators goes beyond a certain threshold.

• If events are generated outside of network components, what does that look like?

The system or sub-system level events are one example of events generated outside of network components. There may be more types as well. I can only think of these as of now.

• Can an "event" generated by a network component be solely defined and handled by that network component? In other words, can a bus or branch generate an "event" value based solely on it's own state and that of it's immediate neighbors?

Yes, the network component would generate the event (and handle post an event) based on its own state and the states from neighbors it is connected with.

• If yes, computation of "event" values by a bus/branch only needs to access it's own part of the solution vector, which can be managed.
• If no, things get really complicated.

Yes for component-based events. No, for system-based events.

[bjpalmer]

So what does this thing look like? If we have a state vector X and its time derivative dX/dt, is the event handle also a vector of some type or is it a function f(X) of the current state? Is it a vector function of the current state? Also, does the event handle get invoked at each time step?

[abhyshr]

So what does this thing look like? If we have a state vector X and its time derivative dX/dt, is the event handle also a vector of some type or is it a function f(X) of the current state?

See TSSetEventHandler. It defines the PETSc API for setting events and the API for event functions

Is it a vector function of the current state? Also, does the event handle get invoked at each time step?

Yes. The event function is invoked after every time-step. If there is a post-event processing function also get, then the event function is called afterr the post-event function to reset the event values (since they may be changed because the state has changed).

The way the event handling is as follows:

• User sets the number of events, event handling function, and other information via TSSetEventHandler
• PETSc calls the event handling function after every time-step. If it finds a sign-change in the event value at a certain time, it goes into a root-finding mode where it goes back one time-step, shrinks the time-step, does the DAE solve, and calls the event function again. It keeps on doing this till the event-value is zero (or within threshold).
• Once it finds the event value zero-crossing, it calls the user post-event callback function, resyncs the time-step (for fixed time-step) and then keeps chugging along in time.

[bjpalmer]

I've looked at the TSSetEventHandle documentation and it mostly tells you how to notify PETSc about where your event handler is located (as far as I can tell). It doesn't tell you what the event handler itself is suppose to look like. From what I can gather, it looks like the event handler is some sort of vector function where the vector are all the variables that could potentially trigger and event, I'n our case, it looks more like what we have is a vector of functions, each of which could potentially trigger an event. Am I missing something here?

[abhyshr]

@wperkins I like what you have sketched out. One additional note, an important one, for post-event handling : Typically, for faults/line-switching and other events, an extra nonlinear solve on the entire system is done after the event occurs (post-event). For DAE equations are of the form

dx/dt = f(x,y)
 0     = g(x,y)

The post-event nonlinear solve equations are

g(y) = 0, x = xpre

i.e., the post-event only solves for the algebraic variables y and holds the differential variables x at their pre-event values xpre. Alternately, one can also do the following to avoid resizing and use the same DAE vectors/matrices.

0  =  x - xpre
0  =  g(x,y)

[abhyshr]

I've looked at the TSSetEventHandle documentation and it mostly tells you how to notify PETSc about where your event handler is located (as far as I can tell). It doesn't tell you what the event handler itself is suppose to look like.

Are you asking about the prototype for the event handling and postevent handling routines?

PetscErrorCode PetscEventHandler(TS ts,PetscReal t,Vec U,PetscScalar fvalue[],void* ctx)
PetscErrorCode PostEvent(TS ts,PetscInt nevents_zero,PetscInt events_zero[],PetscReal t,Vec U,PetscBool forwardsolve,void* ctx)

This is given in the TSSetEventHandler documentation.

From what I can gather, it looks like the event handler is some sort of vector function where the vector are all the variables that could potentially trigger and event, I'n our case, it looks more like what we have is a vector of functions, each of which could potentially trigger an event. Am I missing something here?

Yes, that's correct. We have many events, so the event function is a vector.

[bjpalmer]

So if I understand the PetscEventHandler interface, the vector U is the current state of the system and the array fvalue are the corresponding values of the event functions? Does the nevents parameter in the the TSSetEventHandler function refer to the dimension of fvalue?

[wperkins]

So if I understand the PetscEventHandler interface, the vector U is the current state of the system and the array fvalue are the corresponding values of the event functions? Does the nevents parameter in the the TSSetEventHandler function refer to the dimension of fvalue?

I believe the answer is yes to both questions.

[bjpalmer]

Does the fvalue array represent a replicated data structure across the entire solver domain or does it just represent the function values associated with a particular processor?

Either way, one way to handle this (if I understand this correctly) would be to create a framework level function that pushes the current state of the system onto the network, updates the network to reflect current state and then queries each of the network elements to find the current value of the event functions. The framework function then gathers all the current values and stores them in the fvalue array. The user is only required to implement a function Event on each bus (assuming we are only paying attention to buses for the moment) that returns a vector of values from the bus representing whatever event functions the bus is responsible for. The framework then organizes the values into a larger vector that is passed out of the event handler function.

I think you could do something similar on the post event side.

[abhyshr]

@wperkins : See this example from PETSc. It solves the power grid DAE equations and has both time-based and state-based events. The time-based events are triggered on reaching the fault-on and fault-off times. The state-based events are triggered when reaching the thresholds.

This is a slightly involved, but typical case, where each state-based event is a set of two complementary events. The first one checks whether the state has reached its limit (hold condition) if it is within limits, and the second one checks whether the derivative has reached 0 (release condition) if the state is held at its limit. The application uses a flag to signal whether the state is within limits or otherwise and this is updated in the postevent function.

[wperkins]

@abhyshr and @bjpalmer , The feature/ds-dae branch of my repository is the GridPACK feature/ds-dae brought up-to-date with my master branch, which has DAE Events and is up to date with the GridPACK master branch. I have recoded the fault in dsim so that DAESolver::Events are used. The application runs to completion (10s). I have not checked the output, but the residuals are very similar both before and after the fault.
There are some things that need to be worked out. I've assumed that the DAESolver::Event is created local to a processor and can only compute and react to events generated by locally stored parts of the network. I've assumed that a specialized DAESolver::EventManager performs any collective operations needed after all local events are handled. It may be that an explicitly collective event class is necessary that would compute event values and handle events collectively. I suspect we'll need to work that out in order to make this application parallel.
I've attached output from the application in my repository (dsim.txt) and the original in the GridPACK repositiory (dsim.orig.txt).
Please check over the code changes (SHA: 69aadc9f2a5) and see if it makes any sense.
dsim.orig.txt dsim.txt

[abhyshr]

@wperkins : Thanks for working on this. I looked at the changes you've made and they look good to me. I agree with you that in parallel, there needs to be some reduction operation when doing the resolve.

[abhyshr]

@wperkins I am getting this error when running the DAE application.

bash-3.2$ ./dsim.x

GridPACK math module configured on 1 processors
DSim: Finished Reading data files case9.raw and case9_GENROU_ESST1A_WSIEG1.dyr
0: I have 9 buses and 9 branches
DSim:Finished partitioning network
DSim:Finished setting up factory
DSim:Finished setting up mappers
DSim:Finished setting up DAE solver
DSim:Set up completed
start solving:
0 TS dt 0.01 time 0.
    0 SNES Function norm 7.074816342825e-04 
    1 SNES Function norm 2.405869430116e-12 
[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Object is in wrong state
[0]PETSC ERROR:  Vec is locked read only, argument # 1
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.8.4, Mar, 24, 2018 
[0]PETSC ERROR: ./dsim.x on a arch-gridpack-2 named WE37252 by abhy245 Tue Jul  7 11:24:47 2020
[0]PETSC ERROR: Configure options --download-metis --download-mpich --download-parmetis --download-suitesparse --download-superlu --download-superlu_dist --with-cc=gcc --with-clanguage=cxx --with-cxx=g++ --with-fc=/usr/local/gfortran/bin/gfortran --with-scalar-type=complex --with-shared-libraries=0 --with-x=0 PETSC_ARCH=arch-gridpack-2
[0]PETSC ERROR: #1 VecGetArray() line 1576 in /Users/abhy245/software/petsc-3.8.4/src/vec/vec/interface/rvector.c
[0]PETSC ERROR: #2 VecSetValues_Seq() line 697 in /Users/abhy245/software/petsc-3.8.4/src/vec/vec/impls/seq/bvec2.c
[0]PETSC ERROR: #3 VecSetValues() line 854 in /Users/abhy245/software/petsc-3.8.4/src/vec/vec/interface/rvector.c
[0]PETSC ERROR: #4 p_setOrAddElements() line 190 in /Users/abhy245/Documents/software/GridPACK/src/math/petsc/petsc_vector_implementation.hpp
libc++abi.dylib: terminating with uncaught exception of type gridpack::math::PETScException: PETSc error (0): Error detected in C PETSc
Abort trap: 6

Are you trying to insert/add values to the solution vector after the time-step is completed?