CVODES Error -4

[AS]

Hi All,

My Cantera code is crashing-- details below:

***********************************************************************
CanteraError thrown by CVodesIntegrator:
CVodes error encountered. Error code: -4
At t = 0.642546 and h = 2.37699e-09, the corrector convergence test failed repeatedly or with |h| = hmin.

Components with largest weighted error estimates:
3: -85203.2
7: 2261.79
4: 2251.67
1: 535.26
2: 518.651
5: -420.552
6: -1.83677e-13
0: -6.09082e-14

I found some previous threads about this:

https://groups.google.com/forum/#!searchin/cantera-users/cvodes$20-4|sort:relevance/cantera-users/Lzl89B-dVM0/c3JLJx1m-OYJ
and
https://groups.google.com/forum/#!searchin/cantera-users/cvodes$20-4|sort:relevance/cantera-users/5lhfTcpRQbk/FJEqMPbe9s0J

but these fixes did not help. I tried to change the tolerances and the amount of tries that the solver takes. Attached is my code. I wrote the mechanism based on information from a text. Any help would be greatly appreciated!

Thanks!!

[Kyle B]

AS:

Both this and your previous post seem eerily similar to the trouble I am currently having with my PSR code (just posted a moment ago regarding solution error in 2.3.0a2).  

Kyle

[AS]

Hi Kyle,

I just read over your thread.

They do seem pretty similiar. I have been able to run mine with the GRi 3.0 and it seems to work okay. When I switch to other mechanisms (smaller mechanisms) it doesn't work.

I've tried changing tolerances and attempts before it crashes, but it always seems to find another point to crash. My first question involves an extremely long code, and I can't get it uncouple, which is infuriating. The second one is a simple PSR where the other mechanisms work, but not this one. May be related.

Have you got a good solution to yours? Have you tried running with the code I posted above? Does it crash for you, too? I don't have anybody here who can help, but I really need to solve this problem. Hopefully somebody seems my flag and can find the issue!

[Kyle B]

AS:

I've been running with GRI 3.0 as well as a jet fuel mechanism, and ended up with integrator problems in both cases; no solution as yet. Tolerances don't seem to make a difference, nor does changing the max_err_test_fails value.  As a check, what version of Cantera are you using and from what platform?

I'm unfortunately not set up to run from C++, and it would not be straightforward to give that computer this capability, so I won't be able to test your code.  However I would be interested to see what happens if you switch it to use .advance instead of .step.  My code seemed to fail less when using step, and has had a lot of trouble using .advance, which is very weird since .advance is basically repeated calls to .step.  

Kyle

[AS]

Hi Kyle,

I just switched to advance from step-- it crashes even sooner.

Not totally sure where to go from here...

[Ray Speth]

Hi,

I am not able to replicate the error you mention using either the current stable or development versions of Cantera. Instead, I find that the integration fails at t=0.187997 s, at a point where the interaction between the fixed inlet / outlet mass flow rates and the exceedingly large wall expansion rate coefficient causes the volume to approach zero. If you want to represent a constant pressure WSR, then you should do so by using a MassFlowController for the inlet and a PressureController linked to that. I think you will have better luck with that approach.

Regards,

Ray

[AS]

As always, thanks Ray.

I am confused why you can't reproduce this error, though... This is what I get:

CanteraError thrown by CVodesIntegrator:
CVodes error encountered. Error code: -4

At t = 0.36875 and h = 5.13086e-09, the corrector convergence test failed repeatedly or with |h| = hmin.

Components with largest weighted error estimates:

3: -59955.4
7: 889.532
4: 887.465
1: 348.224
2: 336.033
5: -286.102
0: -4.18195e-14
6: 3.54014e-14
***********************************************************************

I am running on Linux Fedora 20 using Cantera 2.2.0. Is there something I am doing wrong? Why can't we see the same thing?

Secondly, I am confused why a pressure controller would be appropriate instead of a wall, shouldn't it do the same thing? I don't know much about the pressure controllers, but I'm not sure that a pressure controller is appropriate here. The example posted code above may be misleading, since in my real case I have:

reservoir ----> WSR ----> WSR ----> WSR ----> ...... ----> reservoir

and I want the pressure between all of the WSRs constant, and I want to dictate all of the mass flows between them individually. From what I've read, pressure controllers are designed to be used on reactor outlets, but what if the outlet is an inlet to another reactor and I want that reactor to be at the same pressure AND have a mass flow associated with it that I can set?

I'm not sure how this would look, but I can do some research if you think this is still appropriate...

Thanks a ton for the help!!!!

[Kyle B]

AS:

Other may have more of an idea, but couldn't you use IdealGasConstPressureReactor?  That includes a variable volume as part of the formulation and might avoid some of the issues you're having with the wall definition.

Kyle

[AS]

Hi Kyle,

This is true-- I believe I can. I'm not sure what the difference is though. Shouldn't using:

IdealGasConsPressureReactor, IdealGasReactor + wall, IdealGasReactor + PressureController, IdealGasReactor + Valve, etc.

all give the same answer? I'm going for computational speed and have tens or hundreds of these reactors all joined together in various architectures based on what my code dictates. I had a conversation with somebody who told me that speed will be gained by using IdealGasReactor + wall, but I'm afraid now that I was led awry.

Also, I still want to know why my code is not crashing for Ray, but does crash for me! If anybody has any ideas, it would be very helpful!

Additionally, I have  noticed something very peculiar: My code works well with larger mechanisms than smaller (eg: GRI vs. single step). I try to stick with small mechanisms to get the code to run fast, but only larger mechanisms run without failure (so far). For example, the code I posted above was run with my single step mechanism and crashes. When I run it with GRI 3.0, I can't get it to crash, even if I try!

-AS

[Ray Speth]

AS,

These different options can result in a constant pressure reactor, but the equations solved for each are different. The differences can be both in terms of the "physical" system being solved, or in terms of the properties of the equations.

`IdealGasConstPressureReactor` and IdealGasReactor + Wall with infinite expansion rate coefficient (`K`) are equivalent in that they both keep the reactor pressure constant by changing the volume. However, using the IdealGasReactor+Wall method results in an equation with a time constant that scales inversely with `K`. If you set `K` too large, then it will make the system unnecessarily stiff and result in solver failures. In addition, both of these options result in reactor where the volume changes. If you're trying to control the residence time in each part of a chain of reactors, this is presumably not what you want since in this case your steady-state solution will be path-dependent.

IdealGasReactor + either PressureController or Valve will keep the volume of the reactor constant. As with using a Wall, increasing the valve coefficient can generate unnecessary stiffness. One of the benefits of PressureController is that in the steady state, the error in the pressure will be zero, regardless of the valve coefficient, whereas a Valve will always have some pressure differential, so in the latter case you are more tempted to use a high value. In the case of a chain of WSRs, you should use a MFC for the first, and then the rest should be PressureControllers linked to the first MFC.

Regards,

Ray

[AS]

Ray, thanks for answering! I want to clarify something if you don't mind!

I am going for speed  + robustness (as fast as possible while NOT crashing!!!!!) I think I misled you before. I don't have a simple chain of reactors in a row, I have an entire matrix (attached) with various reservoirs + reactors which are interconnected (interconnected right now with MFCs and walls)...

It sounds like you think I should use a mass flow controller from each reservoir and then pressure controllers between reactors which are linked to one of the many reservoirs? I just want to clarify that this is the best way to do it. Remember, each reactors sits at the same and constant pressure while I want to control the mass flow between all of the arrows. I'm not sure how this would look with pressure controllers instead of MFCs, but if you say this is the way to do it then I can research it and figure it out. What would confuse me is which pressure controller links to which mass flow controller. I just want your input that this should be a MFC from each res to each reactor and then pressure controllers between instead of something like a constant pressure reactor.

Please see the figure and let me know if you don't mind. This will be extremely helpful to me!

Thanks a ton, Ray!

-AS

[Bryan W. Weber]

Hi AS,

Maybe this is a misunderstanding, but your system can never reach steady state, because the lower left and middle reactors don't have outlets. Is that intentional, or just because this is a sample case?

Regards,

Bryan

[AS]

Hi Bryan,

This is just a simple case. My fault for making it "busy". I was not careful. each reactor will have at least one outlet so that yes, it can reach steady state. I was not careful when I made this figure! My apologies.

In general, there are nxn reactors in a matrix and each reactor does indeed have at least one outlet. Again, my fault. Each reactor can have 1-3 inlets and 1-3 outlets, but number of inlets and number of outlets add to four ALWAYS, for each reactor.

With this being said (and I'm sorry again), what is the best thing to do in this situation, via my above question.

Thanks Bryan and Ray for your help! I really appreciate it, because I'm really struggling!

-AS

[Ray Speth]

Hi AS,

Are all of the mass flow rates constant (not necessarily equal)? In that case, what you might want to do is link each neighboring reactor with a MFC and a pair of valves, one going in each direction (each valve only permits flow in one direction). When setting the valve coefficients, recall that this value is mass flow in kg/s per Pa of pressure drop -- even a valve coefficient of 1.0 is quite large. In the steady state, the pressure difference will go to zero even with a finite valve coefficient, so long as the net mass flow rate through the MFCs is zero for each reactor.

Your situation is a bit more complicated than what PressureController can handle on its own, but this arrangement should provide the same effect.

Regards,

Ray

[AS]

Hi Ray, thanks for getting back to me!!

To complicate things more, the mass flow rates are unfortunately NOT constant. They change as the simulation progresses.

Does this change your answer? Does it look like the best way to do this is still one MFC and two valves (one in each direction) between reactors? Are walls still out of the question (right now they are making things crash?) Doesn't the valve still make the equations stiff?

What's the best way to handle this? Please help!

Thanks a ton, Ray. This is a really useful conversation, I appreciate it!

[AS]

Also, is just a an IdealGasConstPressureReactor out of the question, also?

[Ray Speth]

AS,

If the mass flow rates are not constant, it will just require more care to make sure that the sum of mass flow rates in and out of of each reactor adds up to zero at all times. Other than that, I think my suggestion still stands.

Valves can add stiffness, just like the walls. You need to exercise some judgement when picking the valve coefficients so that the resulting time constants are short enough that you get the constant pressure behavior you want but not so short as to introduce any additional stiffness.

I think using Walls or IdealGasConstPressureReactor will not do what you want, because the volume of each reactor will change with time, which (I presume) does not represent the system that you're trying to model.

Regards,

Ray

[AS]

Hi Ray,

Thanks. And yes-- I do ensure that mass flow out minus mass flow in for all reactors = 0 at all times (to some tolerance, usually less than ~1e-4).

I've haven't had luck with using walls or const pressure reactors. I'll try the valves.

Thanks for the insight. I'm going to fool around a bit with this and see if it helps.

Just to be clear, I made a diagram and I want to know if this is what you mean. I think it is, but I'm struggling, and just want to make sure your suggestion is what I am going to try. This shows just one reactor in top left of the network...

Thanks a ton!!!!!!

-AS

[AS]

Here, I just set the MFC to have the correct mass flow rates (again ensuring that all out-in = 0), and don't put anything through the valves, just let them control the pressure, and i control the mass flow rate through the MFCs?

[AS]

My apologies for the overwhelming amount of posts, I just wanted to give some more information again:

I running the case where walls sit between, and changed the expansion coefficient from 1.0e6 to 1.0e4. It ran later before it crashed. won't this just happen if I use valves, too. I mean, does it matter which flow controllers I use? As long as the coefficneients are stiff, won't it crash either way? what is special about switching to the valve vs. the wall?

[Ray Speth]

AS,

Yes, that picture is essentially what I was suggesting. One modification that I might suggest would be for the places where you have Reactors connected to Reservoirs, only keep the Valve that flows in the same direction as the imposed mass flow. Otherwise, you might get flow coming in from your outlet reservoirs.

Yes, a high wall expansion coefficient and a high valve coefficient both introduce stiffness. There are two reasons to use valves instead of walls: first, it preserves the geometry of your reactor network. With walls, your reactor volumes will change in some way, which will not give you the expected residence time within each reactor for the specified mass flow rates. Second is that a moving wall can drive the reactor volume to zero, in which case the integrator is essentially guaranteed to fail (e.g. even something as simple as calculating the density will fail when the mass and volume of the reactor both go to zero).

Regards,

Ray

[AS]

Got it, Ray. Thanks a lot for the help.

I started putting in valves to the main code. Seems to run. Looks like it takes a while to get to the steady state solution, though. I have to explore this more.

Thanks a ton!!

-AS

[AS]

Hi again, Ray,

So I did some studies. It turns out that you only really need valves between the reactors (and one reservoir) in order to keep the pressure constant. This makes sense to me. (Attached figure)

I did some studies, and if you put valves on the front and rear reservoirs, it runs slowly in the beginning, but then speeds up. Either way, runtime averages out to about the same.

It still fails sometimes (err code -3), but I can try and catch this. It definitely doesn't fail as easily.

I do have another question though. Not sure if I need to open up another topic or if it can go here. I am going to post it below, so sorry if it screws things up-- but it is still related to what I have:

-----------------------------------------

It seems like it takes forever to get to steady state. I've seen tons of topics here and other places that compare chemkin and cantera reactors. Chemkin has the ability to solve algebraic equations to get to SS. I think my problem is that I just have so many things interacting, and it takes a while for all reactors to exchange mass/radicals/etc. to react to get to a SS solution. This may be a longshot, but is there any trick I can plan to get me here faster for a large system?? I ran out to 3s. (which took hours for just a small system of about 10 reactors), and it still wasn't at SS. It looked like it was because temperature wasn't changing drastically (but still changing a bit nonetheless). I am really trying to capture the SS heat release, but the enthalpies and net production rates keep changing, even when it looks like the T is relatively constant. Any ideas?

Thanks a lot again for all the help, and all of the continued help!

[AS]

I realize that what I wrote may not have been explicit:

The temperature stays "relatively" constant, but the "heat release" looks like it is increasing continuously, and I'm not sure why.

The code I use to calculate heat release is what I believe everybody uses:

heatRel = sum(enthalpies * netProdRate) * vol * T * R * netTime,

which should give the actual heat release in J in a given time.

I start with a high T at the beginning to ensure things react quickly. They do react, but it takes a really long time for the heat release rate to get "high"-- ie. what I believe the actual heat release to be. Does this make sense? Is the above formula correct? Is it a steady state issue? Anybody have any ideas why the heat release continuously changes?

[Ray Speth]

AS,

I’m not sure about your formula. What should go to a constant value at steady state is the heat release rate (in J/s), which is just:

sum(gas.partial_molar_enthalpies * gas.net_production_rates)

I’m not sure why it takes so long for your system to reach steady state. How large is your reaction mechanism? Are you still working with the single-step mechanism that you posted before? It definitely shouldn’t take hours to reach a steady-state solution with that mechanism and a system of 10 reactors. With a larger mechanism and system of reactors, things can get pretty slowly, as Cantera uses a dense solver for the stiff equations which scales as O(N^3) where N is the total number of state variables.

Your idea of only having one pair of valves connected to the reservoirs makes perfect sense to me.

Regards,
Ray

[AS]

Hey Ray,

Thanks for the info.

I did actually switch reaction schemes, and I am now using a 2 step mechanism (2 reactions and six species, but only 5 of them react since N2 only dilutes). It's pretty close to the single step I was using before, but I couldn't get those heat releases to match. This mechanism I got online and there was a paper on it, so I'm assuming it will be a bit better than the one I created. The mechanism is the 2 step kerosene luche mechanism (2S_KERO_LUCHE.cti).

The network I ran had 16 reactors. It probably took a few hours before the temperatures and heat releases stopped changing enough that I was comfortable to consider this as steady state. Perhaps this is just normal. Does this time frame sound reasonable?

Also, it does not seem to be very stable. I typically get error code -3 unless I net.setMaxErrTestFails(), which I usually use 10. Is this because of something I'm doing? Or Is it the nature of my network?

Thanks a ton!!

[Ray Speth]

Hi AS,

With 6 species * 16 reactors, I'd expect the convergence time to be no more than a couple of minutes. Without a concrete example, it's hard to determine what the cause of the slowdown is, though. I would look at the values you are using for the valve coefficients as a starting point, and also make sure you aren't running a copy of Cantera where compiler optimizations have been disabled (they should be on by default).

Regards,

Ray

[AS]

Thanks, Ray, I will check the valve coefficients! I've just left them all constant and this point, perhaps they need to be all different.... or just all lower.

I'm looking in my /cantera-2.2.0/ directory, but can't find where to look to determine if the optimization flag was on during compilation/install. Can you tell me where to look?

Thanks!
-AS

[Ray Speth]

AS,

Unfortunately, I don't think there's a way to see the compilation options that were used for an installed copy of Cantera. I'll have to think about whether there's a good way to include this information with the installation in future versions.

Regards,

Ray

[AS]

Crud! okay....

I took the cantera tarball that I have and made a new directory to reinstall. Note the directions located here:

http://www.cantera.org/docs/sphinx/html/install.html#ubuntu

for whatever reason *looks* like it has optimizations turned off in the example:

scons build optimize=n blas_lapack_libs=blas,lapack prefix=/opt/cantera

I put optimize=y and reran. It didn't seem to really make a difference, though. I may very well still not have optimizations enabled, and just not know what I am doing!

Also, all of my valves have coeffients of 1, which should be alright. I fooled around with them a bit, but this seems to make it run fastest. Not sure what else could be happening. I want a high pressure coefficient to keep the pressure constant, but I also don't want ti introduct stuffness...

Is there a good way to calcualte a "good" coefficient so that the pressure is constant but not very stiff??

-AS 

[AS]

Something else to note, and not sure if it helps.

The mass flow rates to the reactors change quite frequently, in addition to the volumes being quite low (~0.001 m^3). Does this generally create more stiffness??

-AS