concerns about how we report and save cooling/luminosity

[jhmatthews]

Current Status: Unsolved concern.
Version discovered in: Python 76 but seems to be as far back as 58.
Brief statement of problem: The cooling we report on screen, and save, is not the same one we use to minimise the temperature
How to reproduce: Any model.
Discussion:

This was discovered when testing macro atoms, as a difference was observed in the line and bound-free luminosities between versions 76 and 58.

Here are printouts in three different runs (spanning both macro, non macro and 58, 76), of the sum over lum_rad in all cells before, and also the line luminosity and total luminosity in a specific cell of interest.

yso_76_macro.diag:  
wind call pre: Wind luminosity  1.84e+36 cell 460 lum lines 1.52e+33 cell 460 lum rad 1.62e+33
wind call post: Wind luminosity  1.99e+36 cell 460 lum lines 1.70e+33 cell 460 lum rad 1.81e+33

yso_58_macro.diag:
wind call pre: Wind luminosity  1.84e+36 cell 460 lum lines 1.53e+33 cell 460 lum rad 1.62e+33
wind call post: Wind luminosity  4.00e+36 cell 460 lum lines 4.73e+33 cell 460 lum rad 4.98e+33

yso_76_nonmacro.diag:
wind call pre: Wind luminosity  1.33e+34 cell 460 lum lines 2.31e+30 cell 460 lum rad 1.14e+31
wind call post: Wind luminosity  1.69e+34 cell 460 lum lines 7.90e+30 cell 460 lum rad 1.65e+31

So a few fairly obvious observations from this:

• directly before the call to wind_luminosity in wind_updates2d.c, the luminosities in cells are different. I think this is wrong, or am I misunderstanding the routine?
• regardless of whether that is a problem, in 58 and 76 the figures agree before the call and disagree after...so something is definitely up here
• this also happens in non macro mode.

So why does a call to wind_luminosity in wind_update lead to a different calculation of luminosity?

Here's why.

We start with an initial guess of temperature and level populations(matom)/ion abundances (non-matom)

T_0, P_0

after photon transport we having a heating and cooling rate

H_1, C_1

These are then used to obtain a temperature solution, which means our four values of interest are

H_1, C_1, T_1, P_0

After we have our new temperature, we work out a new set of level populations / ion abundances (calling either lucy, nabular_concentrations or macro_pops depending on mode)

H_1, C_1, T_1, P_1

Here, our heating and cooling are consistent with the temperature but not really with the level populations / ion abundances. However, when we call wind_luminosity the cooling rate is calculated using the altered level populations.

H_1, C_2, T_1, P_1

This means that the cooling we report on screen, and save in the windsave file is not the same one we use to minimise the temperature

However, I think there's a slightly bigger discussion to be had here, which is the question of whether we can improve the scheme. JM and CK had a discussion about this, and wonder if we can iterate using the new level populations so that the temperature, heating, cooling and level populations / ion abundances are mutually consistent. Is this possible?

With our current scheme, it is not possible to report values which are all consisten with each other. This isn't really a problem as long as we deal in converged models.

I propose a few things

• a fix which reports C_1 rather than C_2
• a discussion about the above point about an iteration scheme inside each cycle with provides four values which are consistent with each other, and whether this is possible
• people should bear in mind that in models where some cells are not converged there may be values reported in the cells which are not consistent, even if we report C_1 rather than C_2.

Note that there is a bug report to follow about what looks to be a bug in macro_pops.

Old red herring issue- convergence, particularly in macro atom mode

Edit: I had some concerns about our solution for T_e, but this turned out to be a red herring. There is a loop in calc_te in ionization.c which does the following before every T_e solution.

for (n = 0; n < 190; n++)
    {
      temp = temptmin + (n * dt);
      zemtemp = zero_emit (temp);
    }

Don't really know why it's there...but I think that clears up issue 2! I suggest we get rid of that as I can't see what purpose it serves unless someone is still using it for diagnostics.

[kslong]

Issue 1:

I am doubtful that something is really wrong here, but it is worth investigating. Wind update changes the temperature and abundances in the wind to match the cooling int he previous cycle so the real question is are you going in the right direction to get close to equilibrium. The only time I would think these would all end up the same is after the wind had converged.

Issue 2:
I agree that the bottom plot looks quite odd. There should be limits on how much the temperature can move in a single iteration and it looks like that is what is happening. The rise to 200,000 degrees should also not be happening. Note that there could be some kind of issue having to do with the fact that H starts to recombine about 10,000 K, and this is where a YSO wind is going to end up.

(JM edit: red herring bug)

—
Reply to this email directly or view it on GitHubhttps://github.com/agnwinds/python/issues/36.

[jhmatthews]

• Are we happy with how we currently report our temperatures and heating and cooling? i.e., are we ok with the fact that the temperature we print out is not the same as that used to find the cooling of the wind?
• Do we feel comfortable with our current method? Should we be doing an iteration wi

If we are happy, I can close. This links with the macor pops issue in that there are a number of areas where it feels like we could try and converge on certain values better within a given ionization cycle (e.g. macro populations, trying to make heating cooling temp and ion fractions mutual consisten within a non converged wind model - not sure if this is even possible...)

[jhmatthews]

I'm closing this, it is kind of the least of our worries at the moment but we should bear it in mind whenever we compare non converged temperatures and heating and cooling

[Higginbottom]

This is a minor concern, but one I'd like to quickly close.
If we call zero_emit once more after we have selected a new value of t_e, then the cooling terms will be correct for the actual t_e in the cell rather than just the last value of t_e used in the minimising routine. I propose simply calling zero_emit once more after the lines in one_shot (around line440 in 77a_dev) where t_e is finally set.
I'm quite sure it won't change anything about how the code runs, but in trying to understand if anything is going wrong, it makes things a little more consistent.

[jhmatthews]

@Higginbottom did you make the above change here?

[Higginbottom]

No - and I also haven't put in a minimum temperature. This should be done since the fix to recomb_cooling allows the plasma to cool down forever in some cases. Before the recomb_cooling would equal zero below a certain temperature - but this was due to a bug. So the fix is good, but the results are not good. I think I logged a bug for this...

[jhmatthews]

This specific issue was discussed at the meeting in July 2014. Closed with caveat that one just needs to be careful whenever comparing non-converged heating/cooling/temperatures, as they may not be consistent.