TS: implement TS support for second-order systems

Issue #124 new
Emil Constantinescu
created an issue

Implement TS support for second-order systems to avoid using a first-order system equivalent for: 1- efficiency reasons (less memory). 2- enhanecd properties (e.g., some structural dynamics simulations need to tune the two-parameter Newmark solvers to provide extra damping for high-frequency modes).

Issue suggested by Manav Bhatia.

Comments (20)

  1. BarryFSmith

    typedef PetscErrorCode (TSIFunction2)(TS,PetscReal,Vec,Vec,Vec,Vec,void); Function pointer typedefs are not allowed in new code. Eventually, I'll remove the ones Jed put in the past.

    If we got with the "2" naming convention then it should be TS2OtherStuff, not TSOtherStuff2

  2. Lisandro Dalcin

    @BarryFSmith Please clarify: Are function typedefs still OK, right? Is that what you meant in your previous comment?.

    About the "2" naming convention, please note that in the one I used, "2" is not about a TS subtype, but a new residual/Jacobian/etc interface (in the sense rhsfunction vs ifunction). I would not say I like mine, but I don't like yours, either.

    Right now, I only have one TS implementation for 2nd. order systems (TSALPHA2, generalized alpha), however others may be added (e.g. plain Newmark algorithm). All of them should share the new API for 2nd order systems, and I think that TS2Something() may be a bit confusing, as our current name conventions would make me think that this "2" is about a particular TS subtype, while it is actually related to the kind of problem to be solved.

    And of course, remember my implementation is actually a big hack use function pointer compose/query. I think we need to provide a proper implementation by extending the TSOps structure. Do you guys agree?

  3. BarryFSmith

    I do NOT want function typedefs in new code.

    It seems you are proposing extending TS to have some types that support 2nd order systems (as the current types support first order systems)? Since the types are not interchangable, in that, for example I cannot provide TSSetIFunction() and use it with a TSALPHA2 type. It seems to me that you actually have two orthogonal abstract objects (one for 1st order systems and 1 for 2nd order systems) jammed together under the same name TS but they can share essentially no code?

    I understand your point that TS2 could be interpreted as a subclass of TS; my intention was that it is a different class than TS, not a subclass.

    One possible model is to have the abstract class TS for ALL ode solvers (which does not have TSSetIFunction or any other function setting) and then two abstract classes TS1 for first order that supports TSSetIFunction and TS2 that supports TSSetIFunction2. This would mean renaming much of TSXXX to TS1XXX. So you could have code like

    TSCreate(&ts); TSSetOrder(ts,1); TS1SetIFunction(ts,....);

    or

    TSCreate(&ts); TSSetOrder(ts,2); TS2SetIFunction(ts,....); this has different signature than TS1SetIFunction()

    note that calling TS1SetIFunction() on a ts that had TSSetOrder(2) would generate a runtime error.

    In the interest of not needing to rename everything I was thinking of just keeping the current TS name as TS and introducing a TS2 (that is not a subclass of TS) class for 2nd order systems. Since you point out that this is confusing since TS2 could be thought to be a subclass of TS we could use SOODE (terrible name) for 2nd order solvers.

    I am just throwing around ideas here hoping to trigger a brilliant suggestion from someone else on how to organize things.

  4. Emil Constantinescu reporter

    So we cannot have that IFunction has a prototype for order one and another for order two? I mean, for order two, when the system is order one we can just set u_dot_dot to zero... not that it would be smart to do so.

  5. BarryFSmith

    We could certainly. It is a question of finding a good design of interfaces/objects that is clear for users and has lowest maintenance/development cost. Would it be clearer for users if we had completely separate objects/classes for 1st and 2nd order systems or more confusing to have one object that supports both?

  6. Lisandro Dalcin

    BTW, the IFunction/IJacobian interfaces do have a meaning in my TSALPHA2 implementation, in such a case an undamped problem is assumed.

    https://bitbucket.org/dalcinl/petiga/src/default/test/Oscillator2.c?fileviewer=file-view-default#Oscillator2.c-227

    I really do not think we need to introduce at this point the "TS2" new-class concept. This is not much different to TSEULER not supporting the IFunction/IJacobian interface. Anyway, the TSSetOrder(ts,2) seems like a good idea, however the name should be better to not confuse with the method order ( in the sense of LTE order). Maybe TSSetProblemOrder(ts,2) would do?

    I think we could do I2Function,& I2Jacobian to mean "implicit, 2nd order". We still need a god name for {Set|Get}/Solution/Solve/Interpolate/EvaluateStep.

  7. Emil Constantinescu reporter

    Both orders have similar computational patterns. I think we should keep TS a single class. We can make use of

    TSSetEquationType()
    

    that we already use to distinguish among implicit/explicit ODE/DAEs to denote that the system is a second order ODE.

  8. BarryFSmith

    I am not neccessarily advocating it but we could require a VecNest for the Vec arguments in the 2nd order case and thus not need to expand the public interface much at all. So TSSolve(ts,Vec U) would be the same interface and inside for the second order case one would pull out the two nested vectors, for the values and the derivative values. The only troublesome function I see is the Jacobian that has two shift values.

  9. Lisandro Dalcin

    Mmm, I think that's overkill. I mean, the required API extension is not that large, the actual problem is that we cannot came up with good names for the new routines!!! I really think VecNest will lead to unnecesary complications (both for end-users and implementors)

  10. BarryFSmith

    Having the two TSSolve() and TSSolve2() is the most troubling to me. I could live with some of the other extra functions but this one is nasty.

    As a counter argument the unknowns for a 2nd order problem as solved are the solution and the derivative so it should just be stored in 1 vector (which contains the unknowns).

  11. Lisandro Dalcin

    @BarryFSmith @Emil Constantinescu @Shrirang Abhyankar @Mr. Hong Zhang Here you have some initial code: https://bitbucket.org/petsc/petsc/branch/dalcinl/ts-2nd-order-systems

    Some comments on the current state of things:

    • I added what I consider the minimum to make it work. I have not taken advantage of TSEquationType, that enum is not used elsewhere, need to decide about the benefits of actually using it.
    • The implementation is not yet fool-proof. There aren't any error checks for nonsensical settings. Here is where TSEquationType may prove useful.
    • I added TS2{Set|Get}Solution(), but not TS2Solve() because up to now was not strictly required, see the tutorial examples. I think we can survive without it.
    • There is no support for exact final time INTERPOLATE. This would require some extra API or extending the current one. See the following rant about it.

    I HATE the option TS_EXACTFINALTIME_INTERPOLATE, I think it is a bad API to support the feature, the current code is quite convoluted because of that. If users really care about the exact final time, I would change the implementation to STEPOVER if no adapt else MATCHSTEP, then users should do an extra call at the end of TSSolve() and then perform the interpolation. Moreover, right now TSSolve() has a lot of messy code handling the solution vectors just to support INTERPOLATE. Instead, Barry decided to give utmost importance to the exact final time option, and now users are force to set always set it or otherwise die on solve failure.

    While working on this, I started to think that our TSMonitor callback getting the solution vector does not fit well for second-order systems (same thing for TSEvent). Moreover, this is inconsistent with the way KSP and SNES do monitors, you don't get residuals or solutions on them, if you need residual/solution, you just ask KSP/SNES for them in the monitor routine. Would any of you support my proposal for refactoring TS monitors?

  12. Mr. Hong Zhang

    Mandatory usage of TSSetExactFinalTime() may look weird for users using adapt. I think it is a big deal only for fixed-step integration. Perhaps we can force users to set it only when no adapt is applied.

    @Lisandro Dalcin Regarding TSMonitor, are you suggesting to get rid of the solution vector from the argument list and let users call TSGetSolution() if the solution is needed? Based on my own experience, most of the time the solution vector is the only reason that I provide a TSMonitor. So we may end up adding TSGetSolution() inside every TSMonitor function. I am curious why the solution vector does not fit for second-order systems. Are there any special needs for second-order systems?

  13. Lisandro Dalcin

    @Mr. Hong Zhang I'm working on removal of the INTERPOLATE enum value. After that, let see if it would make sense to remove TSSetExactFinalTime() for good.

    @Mr. Hong Zhang The thing is that in second order systems you have two "solution" vectors: the usual one (let's call it displacement vector), and the one with the solution time derivative (let's velocity). Currently, TSMonitor handles the first just fine, but the second requires a call to get the other vector. Anyway, I think you gave a very strong rationale. In KSP/SNES, usually we do not care at all about intermediate residuals or solutions, just the final solution is the important stuff. In TS, most of the time we do care about the intermediate step solutions. So regarding monitoring requirements, it seems TS is a bit different from KSP/SNES. At this point, I surrender, you are right, I decline my request to refactor TSMonitor.

  14. BarryFSmith

    I don't object fundamentally to the idea of removing INTERPOLATE but I am concerned that the original user problem of users not realizing they are getting some "later time solution" and thinking that TS is broken will return. This is why the TSSetExactFinalTime() was forced on everyone. So if someone uses a method that has no time adaption and sets the final time how will you insure that they NEVER make the mistake of thinking that the solution that comes out as the end is actually at the time they requested?

  15. Emil Constantinescu reporter

    Yeah, I think so. I think we could use higher order or explicit methods as well here that are a bit more flexible like Runge-Kutta. But for now we have a working method and an interface, and that's probably a good starting point. I wish we had a good place to list for cool things to do.

  16. Log in to comment