How to simulate accurately with states approaching zero?

Hello,

I'm struggling with a simulation of a network of chemical reactions (CSTR reactor model). As some chemical species concentrations approach zero, I had issues with negative values crashing my simulation. Naturally, I opted for an integrator limited and set the lower limit to zero. This worked fine at first, but now I have issues with the simulation running slowly as the model got more complex. The issues stems from solver resets caused by block changes. This makes the solver reset at every timestep, which is the cause for the slow simulation as in understand. Any tricks of using Saturation block or a Compare-to-Zero block just shift the problem around, causing Jacobian updates at every step which takes just as long.

Now the root cause (the negative concentration values) seem to be caused not by a faulty model but rather by integration errors. I suspect this as I checked the value of the integrator inputs and outputs. At times, always positive inputs resulted in negativ integration result. I found someone on Stackoverflow with a similar issues but no resolve.

I suppose this could be solved by a smaller absolute tolerance, but with low states commonly in the e-25 range and some states reaching as low as e-230, that's not feasible either. Now this is something that should be occuring commonly with states approaching zero, isn't it?

I would be grateful for any hints on solving this issue.