Radiative feedback

[RoguePotato]

Added a radiative feedback class which modifies the "ambient" temperature of a particle to be used in the RadWS scheme. Similar implementation compared to SEREN with hdisc_single providing heating from single central object with a fixed temperature profile. sink_heating provides heating from sink accretion (not fully tested). Other regimes not implemented yet but will include: fixed temperature profiles from a binary system (hdisc_binary); combined disc heating and sink heating (hdisc_single_plus_sink_heating, hdisc_binary_plus_sink_heating).

Only continuous radiative feedback is implemented for sink heating. Episodic may come later. Goes hand in hand with RadWS scheme so branches have been merged.

[rbooth200]

Hi Anthony,

This sounds like a very useful idea, but I have a couple of suggestions / queries at the moment.

1) In general, is preferable to use enums rather than string comparison to decide between modes of operation. However, I think in this case the best solution would be to use a base class and to sub-class it for each type of behaviour.

2) The computation of temp_ambient does not look too expensive, so I'd prefer a solution that computes it on the fly rather than stores it in the particle. The CPU time in gandalf is almost always limited by memory access rather than computation, so optimal performance and scaling is normally acheived by reducing the number of loops and the amount of memory used. Would such a solution be possible? I would expect it to be.

Richard

[RoguePotato]

Hey Richard,

1. I think splitting into sub-classes will be possible, essentially DiscHeating (fixed temperature profile for the central object(s)) and SinkHeating children, where the SinkHeating will have continuous, episodic etc. children classes. Hopefully I can get something simple and understandable implemented.

2. Instead of having a part.temp_ambient, we can have the radiative feedback object take a particle for the GetAmbientTemperature function, which will be called wherever ambient temperature is used (at the moment within the EnergyRadws class). This makes the radiative feedback object a little more robust, removes a particle loop and reduces the memory footprint. If one chooses to include heating from formed sinks, the pointer to sinks within the radiative feedback object would need to be updated whenever a sink forms. This isn't much of an issue however.

Anthony

Aside: The whole issue of memory access causing slowdowns screams out for data-oriented design!

[rbooth200]

Hi Anthony,

The approach you suggest sounds reasonable. By the way, we already store a pointer to the nbody object in the locally isothermal equation of state, so I wouldn't be afraid of doing that if necessary.

Regarding data-oriented design, we are well aware of these ideas. We are now doing a lot of particle sorting and caching, which dramatically improves the memory access patterns for particle access.

The worst affected area remains the tree walk which has inherently uneven access patterns, and is often a dominant part of the cost. I expect it is possible to make significant gains in the tree walk by moving to a structure-of-arrays type approach there. However, this is by no means a quick fix, it requires restructuring the code and making trade-offs between different types of tree walks, which require different data. Lots of experimenting would be required to find the best solution.

Rich

[dhubber]

So, what exactly is the status of this branch? I read some comments on one of the related issues but wasn't sure if that meant a different branch was being worked on and this one is now 'stable'. If this pull request is all sorted and there will be no more immediate updates, then we can look to merge it soon (once we look over it one last time). Any other major changes can be submitted in a new branch/pull request.

[rbooth200]

There are still some things to do here.

First, there is an issue related to the HydroForcesPressure method, which is currently does not behave correctly. The problem is due the neighbour particles not being derived from SPH particles, and now the pressure function needs to be a virtual function, whereas before it was a template function (and virtual template functions are not allowed).

Second, I've asked Anthony to include a test, so we can easily verify if its working in the futuer.

[RoguePotato]

Here is a list of what is left to be done:

1. Remove templated HydroForcesPressure function. The particle pfactor will instead be changed to pressure and the usages of pfactor will be replaced with the full form. HydroForcesParticle will then also store invomega for GradhSph. SM2012Sph pfactor will be changed accordingly.
2. Set sink.dmdt in Sinks.cpp correctly and remove the hack currently implemented in RadiativeFB.cpp.
3. Create a cloud collapse test demonstrating the RadWS method. The form of this plot is what we will look for. The results are rather resolution dependant, but the general features should be comparable. Me and Richard agreed that a nose test would be difficult to devise.

[rbooth200]

I've just noticed that both Sinks::AccreteMassToSinks and MpiControl::UpdateSinksAfterAccretion will need to be modified to properly track the sink accretion rate, but I think that is all that is need, apart from whatever is needed in RadiativeFB.

[RoguePotato]

Okay everything is added, should be ready to merge.

[rbooth200]

I'm happy with this. David / Giovanni do you want to look over it?

[dhubber]

So, from what I understand, it might be difficult to make a nose/Travis test for these algorithms? In that case, the first question is are you satisfied that this is all working correctly? Second thing, in the absence of automatic tests, are there any other safeguards that can be implemented, such as sanity-checking asserts, to try and catch any problems if new errors are introduced in the future? It's not the end of the world if you can't think of any, but thought I will ask that now before we merge.

[RoguePotato]

The RadWS method compares well with previous implementations, cloud collapse and preliminary disc test demonstrate this. I have a lot of confidence that it is working as intended. The RadiativeFB checks out when comparing to pen and paper tests for both disc and sink heating. If I think of any tests in the future I can add them.

[RoguePotato]

Okay this is ready to merge.

[dhubber]

ok, I am happy with the changes and that the branch is ready to merge. Giovanni, since this is a rather large pull request, any further comments/suggestions yourself, or are you happy to merge now too?

[rbooth200]

I'm happy for it to go.

[giovanni-rosotti]

I'm happy for this to go. Only one last question: why is it difficult to implement an automatic test? I guess this is because there is no (semi)analytical solution. But we can still store a known working version of the end result (i.e., a snapshot) and compare it with what we have just run. Am I missing something?

[rbooth200]

It is a little difficult to come up with a reliable test. The two obvious ones are 1) cooling in a uniform box and 2) the collapse of a spherical cloud.

The difficulty with 1) is that I'm not sure how well defined the cooling of an infinite uniform box is in terms of a polytrope.

Problem 2) works quite well, but the issue is the time evolution is strongly resolution dependent at low resolution. In particular the length of time needed to reach a given central density varies considerably, as does the density-temperature relation. I think this means that if we make small changes in the scheme we could expect the answer to change, in a way that would not be obvious as to whether it's an improvement or not.

Any other suggestions? I think unit tests would be the best thing for this module, but currently they don't work.

[giovanni-rosotti]

Is the issue with 2) a real show stopper? If we are making changes to the algorithm, then we will know and we can override manually the result of the test. In all the other situations, it will still catch the case in which we have broken the algorithm by modifying something completely unrelated

[rbooth200]

Ok, I think we can do that. The test does take ~ 1minute on 1 processor (for about 4k particles IIRC). Is that going to be ok given how slow travis is?

[giovanni-rosotti]

i think we can leave with 1 minute...

[giovanni-rosotti]

sorry live

[rbooth200]

I was just trying to set up such a test, but there seems to be a problem running this from python. There is a radws_test.dat in tests/gravhydro_tests, which runs fine from the commandline, but not with python.

@giovanni-rosotti can you look into this? I have no clue what is going on.

[giovanni-rosotti]

Richard, did you clone Anthony's fork in a different folder? I also get a crash when I try do that. That's because python is actually getting the shared library from what comes first in your PYTHONPATH, i.e. most likely the other version of gandalf. If instead I add the fork as a remote to my usual folder, everything works fine even with python.

Btw, how come that the error message is in german? how did it slip into the code?!

[rbooth200]

Yeah, the python path was the issue. I've just send Anthony some tests to add that will hopefully suffice. I

Don't ask me about the german. Perhaps @RoguePotato can fix it when he adds the tests?

[RoguePotato]

Okay I've added the tests, thanks Richard. As for the german, it was legacy debug code from Paul Rohde who implemented the skeleton of the method. I've changed it to throw an exception now.

[giovanni-rosotti]

Thanks Antony, that looks fine to me. You could perhaps save a "known working snapshot" and compute the L1 error norm of rho vs r (or something similar), rather than just using the mean? A tolerance of 10% seems also a bit large (but I never played with these calculations).

[rbooth200]

I put the tolerance of 10% in, which if anything is too small given a ~ 1% increase in simulation time can produce an orders of magnitude increase in the central temperature / density.

[giovanni-rosotti]

Why can't we make the comparison at the same time?

[rbooth200]

Actually, my point is that even comparisons at the same time are not reliable. If you look at the travis tests, you'll see the runs without MPI reached the final time in ~ 250 steps, but the MPI runs stalled and doesn't get there after ~4000 steps. The test is so sensitive to small changes, even the tiny difference using MPI makes through changes in the roundoff error is enough to completely change the results at such low resolution.

I suspect that we should just forget about it until the unit tests are working, but I'm willing to try one more suggestion, which is to use an earlier time point, where the differences are hopefully smaller.

[rbooth200]

Ok, even at this time (75/99 steps without/with MPI), there are big differences between the two cases.

@RoguePotato have you confirmed that the Radws simultions are acutally working with MPI, at least at high resolution?

[RoguePotato]

The RadWS method does work with MPI. However, I'm getting large difference in the output depending on whether I'm using MPI or not. This is why the assertion is hit in the unit tests and the reason the build fails. Whether two simulations match at higher resoution with and without MPI requires some investigation.

[rbooth200]

Actually, I think there is an issue with the gpot_hydro calcultion with MPI. I'll try to explain...

When MPI is used we need to calculate the gravity from particles that are local to the processor, and those that are on distant processors. The local contribution is done normally, and included correctly, but the distant calculation is done differently.

First, we have reduced (or 'pruned') copies of the trees on other nodes. If we have enough of the tree to compute the gravitational force without exporting the particle then we use that. This is done in HydroTree::UpdateGravityExportList, so you need to add the gpot_hydro calculation to there.

In cases where we try to use parts of the tree that we did not keep, we instead export the particle to the other processor and compute the forces there. These forces then need to be collected and sent back to the original processor, which is done in MpiControl::GetExportedParticlesAccelerations. To make life easier with the export / return of partices we have introduced 'CommunicationHandler' objects that abstract this away. So all you need to do is add gpot_hydro to the communication handlers.

Those two things should be enough. Let me know if you have any questions.

[rbooth200]

This looks correct to me, but I don't think its going to fix the issue with the test. @RoguePotato: have you confirmed that the code produces the same answer for gpot on the first timestep with/without MPI? There might be tiny differences due to floating point rounding error.

[RoguePotato]

The average gpot and gpot_hydro differ at the very start of the simulation depending on the number of nodes used. This is presumably due to the different tree structure. These differences may manifest into large differences over time.

The Lombardi metric doesn't use the gravitational potential, but it does use hydrodynamical acceleration. This is currently calculated using the total acceleration minus the tree acceleration. As such, using the Lombardi metric different results are produced depending on the number of nodes used.

If we could find the hydrodynamical acceleration without using the acceleration from the tree, then the Lombardi metric should give consistent results. Obviously we require the tree for gpot calculation using the default metric.

[rbooth200]

Hhm, ok. Can you check the differences when using the brute_force tree instead of kdtree? That will avoid the problems with the tree structure.

I'll think about the a gas vs a hydro. Can you get away with the magnitude of the acceleration?

On Wed, 11 Oct 2017, 16:13 Anthony Mercer, notifications@github.com wrote:

The average gpot and gpot_hydro differ at the very start of the simulation depending on the number of nodes used. This is presumably due to the different tree structure. These differences may manifest into large differences over time.

The Lombardi metric doesn't use the gravitational potential, but it does use hydrodynamical acceleration. This is currently calculated using the total acceleration minus the tree acceleration. As such, using the Lombardi metric different results are produced depending on the number of nodes used.

If we could find the hydrodynamical acceleration without using the acceleration from the tree, then the Lombardi metric should give consistent results. Obviously we require the tree for gpot calculation using the default metric.

— You are receiving this because you commented.

Reply to this email directly, view it on GitHub https://github.com/gandalfcode/gandalf/pull/160#issuecomment-335844585, or mute the thread https://github.com/notifications/unsubscribe-auth/AOqItF9gZOfkcwLJWrZ_GEXstPCdkB9zks5srNsRgaJpZM4ObJB8 .

[RoguePotato]

For a single node, the bruteforce and kdtree results match. For two nodes, the bruteforce and kdtree do not match each other nor the single node results.

Yeah the magnitude of the hydrodynamical acceleration is what we need.

[rbooth200]

This sounds like a bug... I will dig around a bit.

[rbooth200]

I still need to work out what is going on with gravity, but there is another issue.

@RoguePotato: Are you aware that the hydro acceleration is zero at the start of the first time-step? This is because part.gamma is 1 until the energy integration object is called.

Honestly, I think the sensible thing to do here is to made the RadwsEOS object a proper EOS, that does things like compute the pressure and sound speed from rho and u using the tabulated values. Are you ok doing that? I can skype this afternoon if necessary.

[RoguePotato]

Okay I think I have something working, alleviating all problems we have been having. It's pretty much a structural rewrite so I'll commit after some general feedback. Here is the list of changes:

• Put the table of opacities/gamma/mubars into a separate class, OpacityTable.cpp which has data access helper functions.
• RadwsEos.cpp now takes a pointer to an OpacityTable, one is only created if the rad_ws eos is chosen in the params file.
• EnergyRadws.cpp now consults the RadwsEos for anything it needs (GetKappa, fcol2 etc.). This is a little ugly, lots of dereferencing e.g. eos->opacity_table->GetKappa(...).
• The gamma/mu_bar values are now found internally, meaning they can be removed from the Particle data structure. As such, initial gamma/mu_bars don't exist, instead Pressure, Temperature etc. are calculated using on-the-fly gamma and mu_bars (this removes the issue of part.gamma=1.0 at the beginning).
• The creation of the uint object for the radws method i.e. uint = new EnergyRadws(...) must specifically take a cast Radws EOS rather than the base EOS because the Radws EOS contains the OpacityTable object. Not clear to me how to get around this without giving the base EOS an OpacityTable*.
• No noticeable slowdowns (holds ones breath!).
• Consistent results for MPI and OMP! Maybe part.gamma and part.mu_bar were not being communicated across nodes?

I'm not sure if these changes are what you meant by making RadwsEOS a proper EOS Richard, but hopefully I'm on the right track. If you're happy with this, I'll clean up the code, make changes to the unit test and commit.

[rbooth200]

@RoguePotato: This is pretty much what I meant. If you push this then I can think about whether the interface needs cleaning up before we go.

I expect the MPI problem might have been part.gamma or part.mu_bar as you suggested.

[RoguePotato]

Okay the parameter stuff is done.

With regard to the functions in the RadWS equation of state class, they should be correct. The integration of specific internal energy takes into account things like molecular hydrogen dissociation etc. from the variable gamma and mu_bar which match the state changes. We also pull specific internal energy from the table when finding the equilibrium value during the du/dt integration. Up to densities 10^-2 g cm^-3 we can assume ideal gas.

I also know why the tests are failing, and it's to do with whether MPI is enabled during compilation. Without MPI, the radws_test runs for a little bit longer and produce more snapshots. During this extra time, the density and specific internal energy increase significantly. The results do match (i.e. the MPI issue is fixed), but the run-time is different. We could run the simulation up to a point where change in density and u is not as dramatic, or we could try and be more precise on the final output time.

[rbooth200]

I'll worry about the test in a bit, but we still need to fix the thermodynamics. I have a mathematica note book with the eos equations from D'Angelo & Bodenheimer (2013) I can send you if necessary.

Let's leave aside the cooling for a moment and stick to the thermodynamics, I suspect that the conversion from equilibrium temperature to internal energy is correct (from the tables or equation 37 of DWTB07). I'm aware that the EOS takes into account dissociation and ionization etc, which is why I think the expressions that we use for the adiabatic case are not valid.

The ideal gas law is fine for computing the pressure (equation 36, DWTB07), given a mean molecular weight, mu. However, if you check you will see that the internal energy you will see that P != rho u (gamma - 1), where gamma is defined by equation 26 of D'Angelo & Bodenheimer (2013). The reason for this is that some of the energy goes into dissociating H2 or ionizing H and He. I appreciate that it is possible to define a 'gamma' such that the current conversion between u and T in the eos is correct, but that gamma is manifestly not the same gamma that appears in the sound speed (which is gamma_1 in DB13).

You can see this by looking at the lower two panels of Fig. 1 from DB13: the middle panel is gamma, while the bottom panel shows the quantity that would be 1 / (mu * (gamma - 1)) if the two gammas were equal. You should be able to see that structure in middle panel that is missing from the bottom.

In general, P != rho u (gamma - 1), where gamma is the ratio of specific heats is only true if X_T and X_rho (equations 27 and 28 of DB13) are equal to 1 (i.e. mu is constant) and u = C_v T, which is only true if C_v is a constant (because C_v = du/dT, which only equals u/T if C_v is a constant.), which is not the case.

The last point, C_v = du/dT, is actually incorrect in Bodenheimer and Black (1974), which DWTB07 uses. This explains the differences between the expression for the energy associated with the excitation of rotational transitions between DWTB07 (equation 45 of DWTB07 and 25 DB13). If you do the integral of the DWTB07 with respect to temperature (using f(T) for BB74) and divide by T, you get the expression of DB13.

[distamio]

I'll worry about the test in a bit, but we still need to fix the thermodynamics. I have a mathematica note book with the eos equations from D'Angelo & Bodenheimer (2013) I can send you if necessary. Let's leave aside the cooling for a moment and stick to the thermodynamics, I suspect that the conversion from equilibrium temperature to internal energy is correct (from the tables or equation 37 of DWTB07). I'm aware that the EOS takes into account dissociation and ionization etc, which is why I think the expressions that we use for the adiabatic case are not valid. The ideal gas law is fine for computing the pressure (equation 36, DWTB07), given a mean molecular weight, mu. However, if you check you will see that the internal energy you will see that P != rho u (gamma - 1), where gamma is defined by equation 26 of D'Angelo & Bodenheimer (2013). The reason for this is that some of the energy goes into dissociating H2 or ionizing H and He. I appreciate that it is possible to define a 'gamma' such that the current conversion between u and T in the eos is correct, but that gamma is manifestly not the same gamma that appears in the sound speed (which is gamma_1 in DB13). You can see this by looking at the lower two panels of Fig. 1 from DB13: the middle panel is gamma, while the bottom panel shows the quantity that would be 1 / (mu * (gamma - 1)) if the two gammas were equal. You should be able to see that structure in middle panel that is missing from the bottom.

I agree with what you say so far. But if the ideal gas law is fine (equation 36, DWTB07), and temperature = ((gamma - 1) U mu m_hydrogen) / k (which in effect is the equation that “defines” our “gamma” from the tables, and is not the same as in BD13, i.e. it is not the adiabatic index) then pressure = (rho k T) / (mu m_hydrogen)= = (gamma - 1) rho U where again “gamma” is not the same as in DB13. What’s your opinion about this?

In general, P != rho u (gamma - 1), where gamma is the ratio of specific heats is only true if X_T and X_rho (equations 27 and 28 of DB13) are equal to 1 (i.e. mu is constant) and u = C_v T, which is only true if C_v is a constant (because C_v = du/dT, which only equals u/T if C_v is a constant.), which is not the case. The last point, C_v = du/dT, is actually incorrect in Bodenheimer and Black (1974), which DWTB07 uses. This explains the differences between the expression for the energy associated with the excitation of rotational transitions between DWTB07 (equation 45 of DWTB07 and 25 DB13). If you do the integral of the DWTB07 with respect to temperature (using f(T) for BB74) and divide by T, you get the expression of DB13.

OK, I can’t follow this point without digging into the equations, but I remember that we had found that something was incorrect in the Bodenheimer and Black (1974) approach and have corrected it. If you have already coded the BD13 method then we can do some computations and then compare the results regarding specific internal energy and pressure to see if they are compatible.

[rbooth200]

I agree with what you say so far. But if the ideal gas law is fine (equation 36, DWTB07), and temperature = ((gamma - 1) U mu m_hydrogen) / k (which in effect is the equation that “defines” our “gamma” from the tables, and is not the same as in BD13, i.e. it is not the adiabatic index) then pressure = (rho k T) / (mu m_hydrogen)= = (gamma - 1) rho U where again “gamma” is not the same as in DB13. What’s your opinion about this?

I agree that this works fine for a relationship between pressure density and temperature. But it's not the same gamma that appears in the adiabatic sound speed (equation 26 of DB13, which I'd like for the meshless). I've attached a figure showing this, the blue line is the adiabatic index, while the orange one should be similar the gamma as you defined it above, 1 + P/(rho u). I've also uploaded the mathematica script I used to generate this. There are likely to be numerous small differences to what you have...

[distamio]

Yeah of course, that's not the same gamma. Maybe we should use a different notation, e.g. gamma_eff, (eff for effective) or something like that, to avoid future misunderstandings.

[rbooth200]

Are you able to compute the gamma that appears in the sound speed easily enough and add it to the tables? If not, don't worry for now I can update the tables in the future with ones that contain the appropriate gamma.

[RoguePotato]

I'm looking into this now! I think currently the gammas are calculated using du/dT but not d(mu)/dT and d(mu)/d(rho). I think this will be a trivial addition however.

[rbooth200]

Ok, great. One of the gamma's should DB13's, the other should be 1 + P/(rho u).

I'm going to think about the tests a bit.

[RoguePotato]

Okay I've added the first adiabatic index to the table and implemented into the sound speed calculation. Here are some plots for the ratio of specific heats (calculated from the ideal gas equation) and the first adiabatic index (calculated using equation 26 of DB13):

1. Little gamma
2. Big gamma

They are made to mimic figure 1 of DM13.

[rbooth200]

Ok, I'm more or less happy for this to go. We still need to decide what we want to do with tests though.

[rbooth200]

For the tests: did you just choose some parameters that are not so sensitive?

I'm basically happy now, so I'll get Giovanni to look at it again next week