Differences between NLTE calculation results

[yeganer]

@chvogl Brought to my attention that the NLTE result I obtained last week does not match the spectrum in figure 10 of the KerzendorfSim2014 paper very well near 3000 angstrom. Today's results don't have that problem, however the region between 4000 and 5000 angstrom does not match as well as before.

The goal of this issue is to document this behavior and provide a centralized platform for discussion, especially whether we need to look more closely at the underlying problem.

I believe this might be related to convergence criteria but I haven't confirmed anything yet.
Here are the plots (please note that for some reason the colors are reversed in both plots):

Carsus:

Last week after fixing NLTE:

[livnehra]

Hi,

I'm getting some weird results with NLTE.
I tried to follow the code, and found that in partition_function.py the data is prepared in r_ul_matrix_reshaped, but then the rates matrix is built using:
rates_matrix = r_lu_matrix + r_ul_matrix + collision_matrix
without using r_ul_matrix_reshaped.

Am I missing something?

Thanks,
Ran

[yeganer]

@livnehra Hey, I'm currently not able to test this but this might be working as intended. If it's a bug or not depends on whether r_ul_matrix_reshaped owns the data or not. AFAIK a reshape of a numpy array only create a different view of the same memory location. That means we use the reshape to more easily assign the data but when we use the data to create the rates_matrix the original view with the updated data is used (because both views point to the same memory). Please correct me when I'm wrong but I assume that's what's going on.

Nice catch! This part of the code could be written differently to better reflect what is going on. I hope this explanation helped.

[livnehra]

@yeganer Hi, of course you are correct and this is not a bug.

After some more digging, I think I may have found the issue -
I noticed that turning on NLTE for Si II changed the ion populations significantly (5x).
I made sure this was not due to temperature, n_electrons, W or other plausible effects.

Finally I noticed that this is what happens:

1. The NLTE code computes a new level_boltzmann_factor for Si II.
2. The partition_function for Si II is updated based on these (and always sums to 1.0).
3. The ion_number_density is computed based on the new partition_function.

I believe that the NLTE treatment was meant to affect only the levels without changing the ionization state. Also, using the Saha equation after messing with the partition function doesn't seem valid...

Thanks!
Ran

[livnehra]

So...

Do you think this may be a bug or am I missing something obvious?

Thanks,
Ran

[unoebauer]

@chvogl may be able to provide you with more information (and correct me if I am wrong) but for now NLTE was only implemented for the excitation balance. That means that only the level populations are calculated based on the rates equation. Ionization still has to be determined following the standard recipies, i.e. either LTE or nebular (which is an approximate NLTE treatment). Of course from a physical standpoint, you should rather use the nebular description than pure LTE.

[livnehra]

That's exactly what I understood, but I see that turning on NLTE species does change the ionization significantly.

I found that the reason for this is that the partition function is calculated by the NLTE code and then the ionization is erroneously calculated based on this new partition function.

Hoping to get some help to avoid messing with the code myself :)

[unoebauer]

Well, that's the thing. In my opinion this is actually not a bug. If you combine NLTE excitation with some approximate-NLTE (i.e. nebular) ionization scheme, I think the NLTE partition function should be used instead of an LTE one.

[livnehra]

Even if that is true, the way that the code is implemented, the Z for the particular ion selected for NLTE treatment is always identically set to 1.0, while the other Z's of that species remain at their "nebular" values.

For example, before I turn on NLTE for Si-2, the Si partition functions are:
0 10.84
1 5.83
2 1.002
3 2.00
...

After turning on NLTE for Si-2, the partition functions are:
0 10.84
1 1.00
2 1.002
3 2.00
...

I don't think that represents anything physical.

[unoebauer]

Would be good to hear @wkerzendorf's and in particular @yeganer's (who has spent a lot of time fixing NLTE) thoughts on this

[wkerzendorf]

@livnehra I think it would be great to get a better insight into your problem. Can we skype or email or something

[livnehra]

Of course!

livnehra@gmail.com

Thanks,
Ran

[livnehra]

Hi,

After correspondence with @wkerzendorf, I suggest changing partiton_function.py Line #192 to

general_level_boltzmann_factor[i].ix[species] = level_boltzmann_factor / sum(level_boltzmann_factor) * sum(general_level_boltzmann_factor[i].ix[species])

This leaves the partition function as it was.
Sorry I don't have an environment set up to branch and change the code.
Any thoughts?

Thanks,
Ran

[wkerzendorf]

@livnehra can you add the same information as you did in the email. @ssim please have a look at it.

[livnehra]

Sure,

First I tried:
general_level_boltzmann_factor[i].ix[species] = level_boltzmann_factor * g[species][0] / level_boltzmann_factor[0]

This resulted in:

1. Without NLTE:
   In [4]: model.plasma.partition_function[1][14]
   Out[4]:
   ion_number
   0 10.645095
   1 5.818883
   2 1.001149
   3 2.000014
   4 1.000000
   In [3]: model.plasma.ion_number_density[1][14]
   Out[3]:
   ion_number
   0 1.785938e+00
   1 1.798299e+07
   2 3.462116e+08
   3 4.030021e+00
   4 1.734323e-15

2. With NLTE:
   In [5]: model.plasma.partition_function[1][14]
   Out[5]:
   ion_number
   0 10.636653
   1 2.948737
   2 1.001120
   3 2.000014
   4 1.000000
   In [4]: model.plasma.ion_number_density[1][14]
   Out[4]:
   ion_number
   0 2.085735e+00
   1 1.015023e+07
   2 3.540443e+08
   3 3.537804e+00
   4 1.248089e-15

This still changes the ion_number_density drastically, when all we wanted was to change the excitation levels.

I propose instead:
general_level_boltzmann_factor[i].ix[species] = level_boltzmann_factor / sum(level_boltzmann_factor) * sum(general_level_boltzmann_factor[i].ix[species])

This gives:

2. With NLTE:
   In [2]: model.plasma.partition_function[1][14]
   Out[2]:
   ion_number
   0 10.637401
   1 5.818222
   2 1.001122
   3 2.000014
   4 1.000000
   In [3]: model.plasma.ion_number_density[1][14]
   Out[3]:
   ion_number
   0 1.971342e+00
   1 1.916908e+07
   2 3.450255e+08
   3 3.528997e+00
   4 1.278970e-15

The partition function is very similar.
There is still a small change in the ion_number_density,
but it may be due to small changes in t_rad and w:
t_rad: 8501 -> 8478 ; w: 0.4317->0.4305

What do you think?

Thanks,
Ran

[ssim]

Hello Ran - I think that this all depends on what "level_boltzmann_factor" is.

Why is is the value for "general_level_boltzmann_factor" that is being changed? Can we be a little clearer here about which physical quantity is which, and what we are trying to get right? It looks to me that "general_level_boltzmann_factor" and "level_boltzmann_factor" are both level populations here. What are they supposed to be? How are they used elsewhere in the code? There is more to getting this correct that just getting the partition function correct, right?

I think our effective definition of the partition function is

Z = g_0 * N_ion / n_0

where g_0 is the ground state statistical weight, n_0 is the ground state population and N_ion is the total ion population.

Which of those expressions does this correspond to?

If (as it seems) there is doubt about what these variables mean we should check back in the older version of the code. It's going to be extremely important that this is all consistent - so I think a clear statement of all the definitions of the variable here is needed.

[ssim]

PS. If the formula on line 58 is correct (and that is the same "general_level_boltzmann_factor" that we want here), then the equivalent thing we'd want in NLTE would (I think) be

general_level_boltzmann_factor (level_i) = g_0 * n_i / n_0 = Z / N_ion * n_i

...but I don't know if that latex formula in the code is completely correct. (If this is correct then I think it might be equivalent to the first suggestion in the post two before this one.)

[livnehra]

Perhaps,

But the way the code is implemented, the ion_number_density is computed based on the new partition_function. According to (my understanding of) the 2014 paper, the NLTE excitation mode is not supposed to (directly) affect the ion populations.

What triggered my examination of this code is that I got a ~4x factor between Si II ion population with and without NLTE. When I implemented the first suggestion there still remained a factor of almost 2x.

Maybe the solution is only to disconnect the excitation calculation from the ionization calculation when using NLTE?

[ssim]

I agree with you that it is somewhat surprising that the NLTE mode changes the partition function by as much as a factor of 2x, but I think it would be good to figure out why that is so - i.e. why is it that the partition function drops so much when the NLTE calculation is active: which levels are depopulating and why?

One could opt to keep the LTE partition functions (e.g. for testing purposes) but I don't think it's physically consistent to use a partition function in the ion calculation that is different from the one implied by the level populations. But it would be good to check for consistency with the calculations in the paper - i.e. compare to the version of the code used then: I think @wkerzendorf knows how to did out that old version for comparison?

[wkerzendorf]

@wkerzendorf is in Ghana and will be available next week

[wkerzendorf]

@livnehra check out tag v0.9 that should be used as a reference.

[livnehra]

Hi,

Glancing at v0.9.2, it seems this was once dealt with directly and not through level_boltzmann_factor (in plasma_array.py) - without affecting the ion populations.

for i in xrange(len(self.t_rads)):
relative_level_populations = np.linalg.solve(rates_matrix[:, :, i], x)
self.level_populations[i].ix[species] = relative_level_populations * self.ion_populations[i].ix[species]

[wkerzendorf]

I think it's easiest to just run it and see what the difference will be.