Fixes for chemeq integrator #214
Conversation
|
I've tested this by running the AGORA disk, along with GM1 for 64 steps using the metal cooling + superbubble before and after this fix was applied. The rho-T distributions look pretty much identical, although the total thermal energy comes out slightly different (by ~1e-3). I'm not sure if this is a problem, as I also see a small but measurable difference in the thermal energy when comparing the chemeq and NR integrator on the base branch. I should note, all comparisons are made after f01bb6d was applied to both branches, since the temperature distribution obviously comes out different otherwise. |
| q[i] = 0.0; | ||
| d[i] = 0.0; | ||
| y0[i] = y[i]; | ||
| y[i] = max(y[i], ymin[i]); |
There was a problem hiding this comment.
It is a bad idea to delete this line: if y[i] is zero, you can infinities later on.
| temp = -fabs(ascr-d[i])*scr2; | ||
| /* If the species is already at the minimum, disregard | ||
| destruction when calculating step size */ | ||
| if (y[i] == ymin[i]) temp = 0.0; |
There was a problem hiding this comment.
I think we might want to keep this line.
|
I'll try putting those two lines back in. I'm running into NaN temperatures occasionally with Michael's simulation. This might be the cause. |
stiff.c
Outdated
| */ | ||
| y[i] = max(ys[i] + dt*scrarray[i], ymin[i]); | ||
| y[i] = ys[i] + dt*scrarray[i]; | ||
| if (i == 0 && y[0] < ymin[0]) return; |
There was a problem hiding this comment.
Do you mean "if (iter == 0 && ...". Testing for "i" looks redundant here.
There was a problem hiding this comment.
This is meant to skip the y[0] < ymin[0] comparison when the for loop is iterating over the abundances i=1...4
It would probably make more sense to just place the check after the for loop. I imagine the compiler might be doing this during optimization anyways.
* Recalculate EMin after calling clAbunds
2 participants
The metal and H2 cooling modules behave strangely as particles approach dCoolingTmin. These problems become especially apparent when the superbubble feedback is enabled.
This PR fixes two issues with the cooling code at low temperature:
Before the StiffStep runs, the minimum energy is calculated using fixed values for the abundances. This creates a temperature offset between the edge of the cooling table and where the integrator actually stops. To fix this, the current abundances for a gas particle are now used to calculate EMin.
The chemeq version of StiffStep tries to prevent particles from falling off the cooling table bounds during integration. This seems to force particles near the bounds to take extremely tiny integration steps, which kills performance. To fix this, the integrator is now allowed to drop below EMin, but will immediately exit upon doing so. This probably means solutions for particles near dCoolingTmin are slightly wrong, but it also gives a 5-10x boost for performance when superbubble feedback is enabled.