Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Formal Integral: including electron scattering #759

Merged
merged 18 commits into from
Aug 1, 2017

Conversation

unoebauer
Copy link
Contributor

@unoebauer unoebauer commented Jul 4, 2017

The goal of this PR is the inclusion of the electron scattering into the Formal Integral scheme, introduced by @yeganer in PR #740.

For this purpose, the mean intensity in the blue and red wings of each line transition have to be calculated (see Eqs. 23, 24 and 28 in Lucy 1999) and passed to the Formal Integral integrator. The most obvious place for this calculation is alongside the determination of the attenuated line source functions att_S_ul in the routine make_source_function in formal_integral.py.

The intensities have to be passed to the formal integrator (alongside att_S_ul). In the integrator, the recursive stepping through the line list has to be adjusted to include the electron scattering contribution (see Eqs. 26, 27 and 28 in Lucy 1999). For this purpose, the electron scattering optical depth accumulated along the constant p paths has to be calculated.

Finally, unit tests for these additions have to be devised, the sanity checks revised and the integration with the configuration system checked.

Note: Even after this change, the Formal Integral will only work with the downbranching scheme for now.

Milestones:

  • calculate Jblues_lu and Jreds_lu
  • calculate electron scattering optical depths on integration segments
  • adjust recursive integration along z for const-p rays
  • properly account for the edges of the line list
  • properly treat e-scattering effect across shell interfaces (not only record dtau but dtau * Jkkp)
  • add unit tests
  • verification with single and multiline setups
  • update documentation

* in preparation for the inclusion of electron scattering into the
  formal integral formalism, the Jblues are calculated
* the calculation occurs in the same routine in which att_S_ul is
  prepared
* WARNING: somehow I broke the integration
* Jred and Jblue are now passed to the integrator
* some preparatory steps for the e-scattering optical depths have been
  made
@unoebauer
Copy link
Contributor Author

Somehow, I broke the integration with the last commit. Any ideas why, @yeganer?

@@ -258,11 +270,23 @@ _formal_integral(
++patt_S_ul)
{
if (*pline < nu_end)
// Calculate e-scattering optical depth to grid cell boundary
zend = storage->time_explosion / C_INV * (1. - nu_end / nu);
dtau += (zend - zstart) * escat_op;
break;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This break get's called every loop because the if statement is not followed by a scope {}

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah - my bad! Coding too much in python lately...

I_nu_r[p_idx] = I_nu_b[p_idx] * (*pexp_tau) + *patt_S_ul;
// In the absence of electron scattering, simple recursion as in Lucy 1991
// TODO: e-scattering: Replace by Lucy 1999, Eq. 27
I_nu_b[p_idx] = I_nu_r[p_idx];

// reset e-scattering opacity
dtau = 0;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure why you reset dtau every iteration.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Because dtau measures the accumulated e-scattering optical depth between consecutive resonances. I still have to modify the calculation of I_nu_b, where this quantity is then used. Afterwards it has to be reset.

@@ -258,11 +270,23 @@ _formal_integral(
++patt_S_ul)
{
if (*pline < nu_end)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a second loop end condition because it's not easy to do that nicely

// Initialize pointers for inner loop
pline = storage->line_list_nu + idx_nu_start;
pexp_tau = exp_tau + offset + idx_nu_start;
patt_S_ul = att_S_ul + offset + idx_nu_start;
pJred_lu = Jred_lu + offset + idx_nu_start;
pJblue_lu = Jblue_lu + offset + idx_nu_start;

for (;pline < storage->line_list_nu + size_line;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe I can change pline < storage->line_list_nu + size_line to
(pline < storage->line_list_nu + size_line) & (*pline < nu_end)

Alternatively one could add idx_nu_end and then change the loop condition to pline < storage->line_list_nu + idx_nu_end.

* these modifications do not change the accurate performance of the
  formal integral method when e-scattering is deactivated
* HOWEVER: some safety checks are still missing! Use with caution!
@unoebauer
Copy link
Contributor Author

Seems to work pretty well already:

formal_integral

@wkerzendorf
Copy link
Member

@unoebauer that looks very good!


double R_ph = storage->r_inner[0];
double R_max = storage->r_outer[size_shell - 1];
double pp[N];
double *exp_tau = calloc(size_tau, sizeof(double));
#pragma omp parallel firstprivate(L, exp_tau)
#ifdef WITHOPENMP
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For the formal integral openmp parallelization I chose a different approach. By providing a dummy header, we can limit #ifdef WITHOPENMP to the header part of the document and I would strongly advise against including it anywhere in the body.

I know this implementation is basically a copy/paste of my progress from cmontecarlo.c but I think we are now able to provide a more readable solution.


double R_ph = storage->r_inner[0];
double R_max = storage->r_outer[size_shell - 1];
double pp[N];
double *exp_tau = calloc(size_tau, sizeof(double));
#pragma omp parallel firstprivate(L, exp_tau)
#ifdef WITHOPENMP
#pragma omp parallel firstprivate(L, exp_tau, finished_nus)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we make finished_nus shared we can make the bottom printing much nicer

@unoebauer
Copy link
Contributor Author

Unfortunately, and despite the excellent results obtained with tardis_example, the current formal integral scheme goes completely mad in a simple pure hydrogen atmosphere, in which only the Halpha is active. This is independent of whether electron scattering is active or not.

I'll have a closer look at it tomorrow.

@yeganer
Copy link
Contributor

yeganer commented Jul 10, 2017

@unoebauer Does this also apply to the version in the current master or only to this PR?

@yeganer
Copy link
Contributor

yeganer commented Jul 10, 2017

Okay, I think I found two bugs.

  1. We use line_search to determine idx_nu_start. The output of line_search is the index of the next line to the red, or no_of_lines if the value is redder than the reddest line. That means, *pline becomes a pointer to uninitialized memory and ALL 2D arrays point to the first line in the next shell.
    We should really think about ending the line list with a 0.

  2. If there is no line interaction AT ALL for one p ray, then I_nu_r never get's updated, that means the intensity is 0 and not related to the value we initialized I_nu_b with.

}
}
I_nu[p_idx] *= p;
I_nu_r[p_idx] *= p;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If nu is redder than the reddest line I_nu_r[p_idx] will always be 0. So we have to account for the fact, that there may not be a line interaction. I propose to initialize I_nu_r as it physically does make sense in my opinion. I_nu_r is the intensity after a interaction. We assume a photosphere and emit radiation there. I think it makes sense to say "After they leave the photosphere there are at the red wing of some transition that happens at the boundary of the photosphere".
For the case that the ray does not intersect the photosphere, we initialize to 0. Here again I think we can apply the argument that we are at the red wing of a hypothetical transition after which our intensity is 0. Then we propagate through matter and apply the electron scattering.

Does this interpretation make sense? It would at least solve the problem and make the code a little bit easier (only Jkk has to be calculated in the first block)

What's left is to decide what happens in terms of electron scattering if there is no resonance.

Swapping the initialization to I_nu_r and changing the code in first does produce the correct spectrum for disabled electron scattering.

@unoebauer
Copy link
Contributor Author

I've now reintroduced I_nu and thus avoided the problem of an uninitialised integration in case of no resonances on the ray. With this change the simple H-atmosphere problem works again when electron scattering is switched off. Once activated, I still get a terrible match between formal integral and virtual packet spectrum.

The problem with the line_list_search is still present, but as I will detail in the next comment, the problems in the simple H-atmosphere with electron scattering are more fundamental and cannot be resolved with the current information we are extracting from the Monte Carlo simulation.

Pragmatically, however, the current state of the Formal Integral simulation should work for all realistic supernova calculations (the formal integral matches the virtual spectrum very well in the tardis example setup, both with and without e-scattering). Since we will typically have a very large line list once a realistic supernova ejecta is set up, we will over the entire spectrum not reach the boundaries of this list and thus never run into the numerical problems associated with the line_list_search and also not into the conceptual problems outline below.

@unoebauer
Copy link
Contributor Author

unoebauer commented Jul 12, 2017

In order to estimate the effect of electron scattering on the various rays, we need to know the strength of the diffuse radiation field at the LF frequency of the frequency bin we are interested in. Lucy 1999 estimates the diffuse radiation field by considering the mean intensity in the blue and red wings of the next and previous lines. This is justified as long as the current LF frequency is not too far from the frequencies of the lines.

It is however a problem for applications like the H-atmosphere with only H-alpha. Here, most parts of the spectrum are far beyond the frequency of this line and the diffuse radiation field cannot be accurately estimated by considering the intensity in the blue or red wing of the H-alpha line.

To correctly calculate the e-scattering contribution in such cases a J_nu estimator would have to be recorded in the Monte Carlo simulation. In principle this is possible (@chvogl has done this already in his Type IIP branch) but associated with additional effort.

@unoebauer unoebauer requested a review from chvogl July 12, 2017 15:17
@unoebauer
Copy link
Contributor Author

@wkerzendorf, @yeganer, @chvogl: please have a look at this PR. I think this is as far as we can go without introducing additional estimators.

// conditions) is not in Lucy 1999; should be
// re-examined carefully
escat_contrib += (zend - zstart) * escat_op * (*pJblue_lu - I_nu[p_idx]) ;
I_nu[p_idx] = I_nu[p_idx] + escat_contrib;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you want you can move this out of the if statement.

// Account for e-scattering, c.f. Eqs 27, 28 in Lucy 1999
Jkkp = 0.5 * (*pJred_lu + *pJblue_lu);
escat_contrib += (zend - zstart) * escat_op * (Jkkp - I_nu[p_idx]) ;
I_nu[p_idx] = I_nu[p_idx] + escat_contrib;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you want you can move this out of the if statement.

@@ -124,4 +128,9 @@ def make_source_function(self):
att_S_ul = ( wave * (q_ul * e_dot_u) * t / (4*np.pi) )

result = pd.DataFrame(att_S_ul.as_matrix(), index=transitions.transition_line_id.values)
return result.ix[atomic_data.lines.index.values].as_matrix(), e_dot_u
att_S_ul = result.ix[atomic_data.lines.index.values].as_matrix()
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

if you just do .values that is maybe clearer that you just want the array and not as a real matrix.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

actually, this is on you :) You've implemented this ages ago (SOCIS 2016) and I've just shuffled things aroun

@wkerzendorf wkerzendorf merged commit 96a94b5 into tardis-sn:master Aug 1, 2017
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

4 participants