[sundials-users] Inverting a SUNDIALS SUNMartix/ Mass Matrix setup error in ARKODE

[Anirudh Jonnalagadda]

Hi everyone,

I'm trying to solve a coupled system of steady, non-linear ODEs of the form R(S) = 0. While going through the KINSOL documentation, I came across the note that considering the time-dependant ODE dS/dt = - M^{-1} R(S), where M is the approximation to the Jacobian Matrix corresponding to R(S), may be suitable for nonlinear, stiff systems. I've currently been able to implement dS/dt = -R(S) but I am not sure how I can invert the SUNMatrix_Dense Jacobian representation.

I also tried avoiding the matrix inversion by using ARKode, but I am getting the following errors:

1. [ARKode::ARKStep ERROR]  arkStep_Init
     The mass matrix solver's setup routine failed.
   
   
2. [ARKode ERROR]  arkInitialSetup
     Error in initialization of time stepper module

Has anyone come across such a situation before? Any pointers will be very helpful!

Best,

Anirudh

---

To unsubscribe from the SUNDIALS-USERS list: write to: mailto:SUNDIALS-USERS-...@LISTSERV.LLNL.GOV

[Hindmarsh, Alan Carleton]

Hello Anirudh,

It sounds like you have not completed your

setup with the ability to solve linear systems

M x = b. This is necessary in order to get the

initial value of dS/dt  (from M Sdot(0) = -R(S(0)).

If ARKode is failing for some other reason, you

might try IDA on the system M dS/dt = -R(S).

But that still requires you to supply the initial

dS/dt, or (using IDACalcIC) the ability to solve

linear systems M x = b.

-Alan H

---

From: sundials-users <sundial...@llnl.gov> on behalf of Anirudh Jonnalagadda <000028a06e27442...@LISTSERV.LLNL.GOV>
Sent: Tuesday, April 25, 2023 1:36 AM
To: sundials-users <sundial...@llnl.gov>
Subject: [sundials-users] Inverting a SUNDIALS SUNMartix/ Mass Matrix setup error in ARKODE

 

[Anirudh Jonnalagadda]

Hi Dr. Hindmarsh, and everyone else,

Thanks for getting back to me and my apologies for my long delay in replying. It turns out I was doing things COMPLETELY wrong. I have since corrected my implementations and got things running with ARKode (I essentially have been able to tackle the problem using KLU and then iterate!)

However I have run into a more problematic situation where I have been trying to solve my problem using the wrong approach. My DAE system involves a constraint which is supposed to intrinsically adjust itself to a match boundary condition (I'm essentially trying to build up to the dynamic optimization work done here but using AMReX, PelePhysics and sundials). Is it possible to solve such a problem via M ds/dt = -R(s) using arkode?

I was trying to start off with a much simpler case (essentially a hydrodynamic situation with only mass, momentum and energy conservation) and was able to set up a very simple python based solution with the levenberg-marquadt method using numpy (I wanted to confirm that my RHS function and Jacobian implementations were correct). I then tried to set up the same case using Kinsol. With a direct sparse linear solver using KLU, I get the KIN_MXNEWT_5X_EXCEEDED error. I also tried using GMRES without preconditioning - that failed with KIN_LINSOLV_NO_RECOVERY. Currently, I am trying to turn preconditioning on - I already have a CSR representation of the Jacobian and am trying to figure out how the pSolve method needs to be defined.

Any suggestions on how I could proceed will be really helpful :)

Thanks!

[Reynolds, Daniel]

Hello Anirudh,

I'm not sure that I completely understand what you're attempting here, but let me answer some at-least-hopefully-related questions.

While ARKODE is designed to solve ODE IVPs that include a potentially time-dependent mass matrix, e.g., M(t) ds/dt = f1(s,t) + f2(s,t), where f1 is treated implicitly and f2 explicitly, it currently requires that the matrix M(t) be nonsingular.  Thus, it is not currently designed to support DAEs (that would typically have singular M), although that kind of extension is under consideration.  If you truly have a DAE, then IDA is the best option in SUNDIALS.

Turning to your KINSOL issues, it is unclear to me what you implemented in Python, so I cannot comment on how that implementation differs from your use of KINSOL.  

However, if your CSR implementation of the Jacobian is correct (and a direct solve is not too costly with respect to memory or computational time), then the KLU solve should be more robust than GMRES (even with a preconditioner).  However, given that KLU fails for you, then I question (a) the accuracy of the CSR Jacobian, and/or (b) the initial guess provided to KINSOL.  If your Python implementation uses the same initial guess as you supply to KINSOL, then that would leave only the CSR Jacobian.  Is it possible to coarsen the problem somewhat and try to use the dense linear solver with difference-quotient Jacobian approximation (i.e., you just don't provide a Jacobian construction routine)?  This would help isolate whether the issue is with your CSR Jacobian, or something else.

On the other hand, I'm perplexed by the error that you get when trying non-preconditioned GMRES.  The KIN_LINSOLV_NO_RECOVERY error should indicate a "recoverable" failure in the linear solver where the residual is not reduced, and that KINSOL cannot remedy within the nonlinear solve.  However, without preconditioning there are very few "recoverable" but non-residual-reduced failure modes from our GMRES solver.  If you want to pursue the GMRES option further, then immediately after it fails can you call the function KINGetLastLinFlag, and share the last value returned by the linear solver?

Daniel R. Reynolds (he/him)

Professor and Department Chair

Mathematics @ SMU

214-768-4339

https://people.smu.edu/dreynolds/

---

From: sundial...@llnl.gov <sundial...@llnl.gov> on behalf of Anirudh Jonnalagadda <jonnalagad...@GMAIL.COM>
Sent: Sunday, July 9, 2023 12:13
To: SUNDIAL...@LISTSERV.LLNL.GOV <SUNDIAL...@LISTSERV.LLNL.GOV>
Subject: Re: [sundials-users] Inverting a SUNDIALS SUNMartix/ Mass Matrix setup error in ARKODE

 

[EXTERNAL SENDER]