[sundials-users] CVode slowing to a crawl for no reason?

[Robert Ellis]

Hi,

Recently I’ve been working on a project in cosmology which entails solving the Boltzmann equation for 2 interacting species. The actual equations include a bunch of nasty looking integrals so I wont bother typing them out, but it boils down to a system of 4 coupled ODE’s (which correspond to the evolution of the temperatures and chemical potentials of both), parametrized by the temperature of the universe. I am using CVode which can solve my equation and gives sensible results for temperatures above ~25MeV. Below this, it slows to a crawl, and the number of steps/evaluations skyrockets for no apparent reason. The only interesting thing that occurs is one of chemical potentials reaches a maximum, but otherwise everything else is perfectly well behaved. I’ve tried adjusting the error tolerances but nothing seems to help much. There is some spiking that occurs when I zoom in (see attachment), but if I lower the tolerance sufficiently this vanishes.

I will note that I do not supply my own Jacobian due to the complexity of my system. I’m expecting this to be the culprit, but I am hoping there may be a way to work around this without providing one explicitly.

I am also not sure if this will help specify the cause of my issue, but here are the statistics I get from the computation (also in case this stands out the step size is supposed to be negative):

Before slowing down:

Time,29.81676446869958,Steps,445,Error test fails,27,NLS step fails,0,Initial step size,-6.388153926949412e-06,Last step size,-0.3377185547218121,Current step size,-0.3377185547218121,Last method order,2,Current method order,2,Stab. lim. order reductions,0,RHS fn evals,697,NLS iters,694,NLS fails,12,NLS iters per step,1.559550561797753,LS setups,76,Jac fn evals,17,LS RHS fn evals,68,Prec setup evals,0,Prec solves,0,LS iters,0,LS fails,0,Jac-times setups,0,Jac-times evals,0,LS iters per NLS iter,0,Jac evals per NLS iter,0.02449567723342939,Prec evals per NLS iter,0,Root fn evals,0

Time,28.80360880453415,Steps,448,Error test fails,27,NLS step fails,0,Initial step size,-6.388153926949412e-06,Last step size,-0.3377185547218121,Current step size,-0.3377185547218121,Last method order,2,Current method order,2,Stab. lim. order reductions,0,RHS fn evals,700,NLS iters,697,NLS fails,12,NLS iters per step,1.555803571428571,LS setups,76,Jac fn evals,17,LS RHS fn evals,68,Prec setup evals,0,Prec solves,0,LS iters,0,LS fails,0,Jac-times setups,0,Jac-times evals,0,LS iters per NLS iter,0,Jac evals per NLS iter,0.02439024390243903,Prec evals per NLS iter,0,Root fn evals,0

Time,27.95893844960629,Steps,452,Error test fails,28,NLS step fails,0,Initial step size,-6.388153926949412e-06,Last step size,-0.1689839334020173,Current step size,-0.1689839334020173,Last method order,2,Current method order,2,Stab. lim. order reductions,0,RHS fn evals,710,NLS iters,707,NLS fails,13,NLS iters per step,1.564159292035398,LS setups,78,Jac fn evals,18,LS RHS fn evals,72,Prec setup evals,0,Prec solves,0,LS iters,0,LS fails,0,Jac-times setups,0,Jac-times evals,0,LS iters per NLS iter,0,Jac evals per NLS iter,0.02545968882602546,Prec evals per NLS iter,0,Root fn evals,0

Time,26.77279710073503,Steps,457,Error test fails,28,NLS step fails,0,Initial step size,-6.388153926949412e-06,Last step size,-0.282724494022406,Current step size,-0.282724494022406,Last method order,3,Current method order,3,Stab. lim. order reductions,0,RHS fn evals,716,NLS iters,713,NLS fails,13,NLS iters per step,1.560175054704595,LS setups,79,Jac fn evals,18,LS RHS fn evals,72,Prec setup evals,0,Prec solves,0,LS iters,0,LS fails,0,Jac-times setups,0,Jac-times evals,0,LS iters per NLS iter,0,Jac evals per NLS iter,0.02524544179523142,Prec evals per NLS iter,0,Root fn evals,0

After slowing down:

Time,14.18906639590881,Steps,21795,Error test fails,1232,NLS step fails,1706,Initial step size,-6.388153926949412e-06,Last step size,-0.003183356203066493,Current step size,-0.003183356203066493,Last method order,3,Current method order,3,Stab. lim. order reductions,0,RHS fn evals,42192,NLS iters,42189,NLS fails,3477,NLS iters per step,1.935719201651755,LS setups,9090,Jac fn evals,3552,LS RHS fn evals,14208,Prec setup evals,0,Prec solves,0,LS iters,0,LS fails,0,Jac-times setups,0,Jac-times evals,0,LS iters per NLS iter,0,Jac evals per NLS iter,0.0841925620422385,Prec evals per NLS iter,0,Root fn evals,0

Time,14.08999491262614,Steps,25468,Error test fails,1555,NLS step fails,1857,Initial step size,-6.388153926949412e-06,Last step size,-1.361802565455957e-05,Current step size,-1.361802565455957e-05,Last method order,4,Current method order,4,Stab. lim. order reductions,0,RHS fn evals,48693,NLS iters,48690,NLS fails,3771,NLS iters per step,1.911810899952882,LS setups,10335,Jac fn evals,3874,LS RHS fn evals,15496,Prec setup evals,0,Prec solves,0,LS iters,0,LS fails,0,Jac-times setups,0,Jac-times evals,0,LS iters per NLS iter,0,Jac evals per NLS iter,0.07956459231875128,Prec evals per NLS iter,0,Root fn evals,0

Time,13.98997805273008,Steps,29569,Error test fails,1929,NLS step fails,1997,Initial step size,-6.388153926949412e-06,Last step size,-2.78664537136988e-05,Current step size,-2.78664537136988e-05,Last method order,4,Current method order,4,Stab. lim. order reductions,0,RHS fn evals,55764,NLS iters,55761,NLS fails,4051,NLS iters per step,1.8857925530116,LS setups,11660,Jac fn evals,4185,LS RHS fn evals,16740,Prec setup evals,0,Prec solves,0,LS iters,0,LS fails,0,Jac-times setups,0,Jac-times evals,0,LS iters per NLS iter,0,Jac evals per NLS iter,0.07505245601764675,Prec evals per NLS iter,0,Root fn evals,0

Is it likely that the Jacobian approximation is failing? Would supplying my own help avoid this slowing down? I would like to make sure before taking the time it will require for me to do so. I hope I've provided enough but would be happy to supply any other information. Any help is appreciated!

---

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

[Jason Zwolak]

Hi Robert,

I would suggest having a relatively larger tolerance for the relative tolerance. Perhaps something like 1e-4. And choosing absolute tolerance very much depends on your system of ODEs. You have to decide what the smallest meaningful value of each state variable is and then choose something like half that. The absolute tolerance is effectively the smallest meaningful value of a state variable where anything smaller is considered zero.

Having tighter tolerances will make the solver run slower. Having looser tolerances may give some jagged lines like the ones your observing in your top right figure. But jagged lines are not a problem as long as the system is still converging, which it is in your case.

Some things to look for in your equations for the different values during simulation:

1. denominator approaching zero.
2. large values for some d/dt with small values for others
3. at least one d/dt that goes quickly from a small value to a large value

If this were my code, I’d be printing out debug statements containing the state variables along with the evaluation of the right hand sides and look for something strange, like shooting off to infinity. And in particular I’d be monitoring any denominators.

Providing a Jacobian will certainly make things faster, but not necessarily fast enough to overcome the slow down you’re experiencing. Without looking closely at the equations, I personally don’t know how to tell you if that would help in your case.

Good luck!

--

Jason Zwolak

[Robert Ellis]

Jason, thank you for your reply! As you suspected, the first two of your points (and possibly the third) seem to apply. 3/4 of my equations have derivatives on the order of 1 while the other is ~1e-5, and since I obtain the derivatives by inverting two systems of linear equations each of my derivatives is proportional to the inverse of a determinant, one of which becomes very small (~1e-10). Is there anything that can be done to remedy this or is it bound to be slow regardless of what I do?

Thanks a lot.

Best,

Robert

[Balos, Cody]

Hi Robert,

 

Are you setting your absolute tolerance using the vector option? I.e., by calling CVodeSVtolerances with a vector where each entry specifies the absolute tolerance for the corresponding state variable.

 

Best,

Cody

[Robert Ellis]

Hi Cody,

I do use a vector to specify the absolute tolerances. I only use one value for relative tolerance, though (as I believe is the only supported way to do so).

Best,

Robert

[Jason Zwolak]

I concur with Cody.

Here’s the beginning of a ChatGPT chat that may help you understand the problem better and correct it.

https://chatgpt.com/share/67f6c995-b7a8-8010-a275-a327d823302f

CVODE is not a “standard numerical solver” by ChatGPT’s definition. CVODE is explicitly capable of solving stiff systems. I cannot remember all the options for CVODE. You may wish to pour through them and be sure it is configured for best results with stiff systems.

I also occasionally encounter really slow simulations and sometimes CVODE quits with an error. It’s almost always a malformed problem with constants that make the system nearly unsolvable. Like really small denominators. But sometimes the problem is that I don’t have correct absolute tolerances (or even the relative tolerance).

--

Jason Zwolak

[Alan C. Hindmarsh]

Hello Robert,

Going back to your original equations, some of your time derivatives
are determined by a linear system that becomes very nearly singular.
You are better off keeping the system that way, in its linearly implicit
form, A dy/dt = f(t,y), and solving this with IDA.  Then the division
by a small determinant goes away, and as long as the system remains
solvable even when A is singular, IDA should solve it.

In any case, the choice of tolerances is also very important, especially
the absolute tolerances when the solution values get very small.
However, a relative tolerance of 1e-4 is borderline too loose.

-Alan H

>> 1. denominator approaching zero.
>> 2. large values for some d/dt with small values for others
>> 3. at least one d/dt that goes quickly from a small value to a large

>> ------------------------------

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

>> ------------------------------

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

> ############################

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

############################

[Mehmet Erol Sanliturk]

In a matrix , when determinant approaches zero , there is a dependency of at least one

of the rows or columns on linear combinations of others .

This means that a row or column is redundant because it can be obtained from affecting rows or columns .

If you have observed or generated data about your variables , you may apply

 - Canonical correlation

 - Stepwise regression ,

 - Principal components analysis ,

 - Factor analysis , or

 - any other dimension reduction method to isolate the dependent variables .

and remove it or them from your system .

In that way , you will reach a more stable system .

After estimating the system parameters , you may generate the dependent ones during your simulations 

if they should be used .

Mehmet Erol Sanliturk

[Robert Ellis]

Hi Alan,

Thank you very much for your suggestions. It makes sense, but at the same time seems like the same problem will still be encountered. To be a bit more specific—for each of my two species I have an equation of the form:

dn_i(T_i,c_i)/dT=f(T_1,T_2,c_1,c_2,T)

dp_i(T_i,c_i)/dT=h(T_1,T_2,c_1,c_2,T)

(n_i=number and p_i=energy density, where i denotes the species. T_i is the temperature of the species, and c_i is the chemical potential of the species divided by the temperature of the species)

I truly am only interested in n_i and p_i, however I am unable to solve for them if I don’t have T_i and c_i (due to the RHS dependance on these variables). This leaves me with two approaches:

Write the above equations in terms of dT_i/dT and dc_i/dT and solve those equations, which is very straightforward

Or

Solve for n_i and p_i and invert these to find T_i and c_i. This is quite difficult—this was what I originally tried but I could not find a reliable way to invert them (as they are both integrals which I don’t believe have a formula for inversion). Even with a fairly good formula for an initial guess, Newtons method worked for larger values but failed for smaller values.

Since the second option is not viable, I am not sure how I’d be able to solve this through the approach you mentioned. Perhaps I am misunderstanding or my explanation of my problem wasn’t clear enough. If I am misunderstanding you, I would greatly appreciate if you’d be able to clarify.

Thanks again!

Best,

Robert

[Alan C. Hindmarsh]

Hi Robert,

In terms of your variables, if
 u = {all the n_i and p_i} and
 v = {all the T_i and c_i},
then your time-dependent system has the form
 du/dt = f(u,v)
 M dv/dt = g(u,v)
where matrix M becomes singular.

This is the sort of system that IDA is intended to solve.
Whether or not this system is well-posed and solvable
depends on specific properties of the function g.

For example, suppose that M is singular and that, after a suitable
change of variables, the vector v consists of subvectors v1 and v2
such that the system M dv/dt = g has the form
  d v1 / dt = g1(u,v1,v2)
             0 = g2(u,v1,v2)
If the second equation, g2 = 0, can be solved for v2 (uniquely)
as a function of u and v1, then this system is well-posed, and
IDA will solve it.  In the DAE literature it is called an
index-1 system.  You do not have to do the change of variables
or separately solve for v2.  IDA does that for you.
If M is nonsingular at time t = 0, then you can supply the
initial values of dv/dt as well as u, v, and du/dt.  But even
if M is singular at t = 0, IDA can get those initial values.

-Alan H