Dust

[rbooth200]

So it's time to think about merging the dust back into the master branch. I have been keeping the dust branch up to date with the master, so this should (currently be straightforward).

The changes are pretty extensive, affecting the tree and grad-h sph, as well as some changes to boundary conditions, IO, and even the meshless (which was very minor). I have modified the ICs to add a couple of new test problems and added dust to some of the others as well (but in a not in systematic way). I have also written some documentation. I've pulled changes from some different places, such as the userguide branch, but I don't envisage any direct problems,

The current issues that I'm aware about are:

• interfacing with python
  • sinks

I haven't done anything regarding sinks, but I'm sure that there will be issues that need fixing. Some help/guidance on this would be much appreciated.

[giovanni-rosotti]

Thank you Richard for your passionate work about this! I can confirm that the code compiles and the dust tests seem to be running fine. Given the size of the request it's a bit difficult for me to do more than this - although I will try to go through the code as well.

• A problem I noticed is that problems without periodic boundaries (I think, at least, I tried freefall and bossbodenheimer in test folder + bossbodenheimer and evrard in dust_tests) crash with the message "Error : looking for incorrect timing block : SPH_ALL_PERIODIC_FORCES". Guess some timing code call is not correct, I'll give a look to it later.

About the userguide, it is largely fine.

• There are a few typos here and there but I can correct those.
• About the papers you reference, I guess we should then include a bibliography?
• What I would like expanded is the section about the tests. If I just go in the tests/dust folder one is a bit confused by finding only parameter files with unfriendly names. Maybe we can:
  • add a couple of images to the userguide (a shocktube or something)
  • add a couple of python scripts in the folder that do some plots - helps in reducing the shock. Currently there are two scripts, but one doesn't do anything and the other one wants to load ADSOD1 - which is not produced by any of the files in the folder.
  • should we setup a radial drift problem in a disc?
• a few scripts (unrelated to dust) have gandalf changed to gandalf_dust - which will prevent them from running for many people...
• until we (I) have fixed the python, how do we distinguish gas from dust particles? it's just the position in the output, is that right? this probably needs to be said in the userguide. And I need to know it anyway for fixing python
• about sinks: David it's best if you comment about that. At least, definitely for the smooth accretion. For standard accretion, I'm not sure why anything should go wrong... I'm sure funny things will happen near the sink where your solution is wrong anyway, but I don't think it should crash. Am I too naive?

I'm sure other things will come to my mind, but this is a long message already so I'll stop now!

PS with hindsight, we should have rebased onto master rather than merging - in this way the list of changes would include only your changes and not also the other ones... but now it's too late. Perhaps we should give it as a policy for future pull requests

[rbooth200]

I've pushed some changes that should fix the names (I'd forgotten that I was playing around with them). I've also fixed the timing names, which arose because I got rid of a load of functions in GradhSphTree.cpp.

The particles are currently identified within the code via ptype, and are stored in the snapshot in order of type, so the type can be worked out from where they are in the file (which is what splash does). I haven't tested the read to see if it's working in gandalf yet...

I'll work on some of the rest.

[dhubber]

Just found a few EndTiming calls that didn't match in GradhSphTree so fixed them. Otherwise, seemed to be working fine for me for the gas-only tests

[dhubber]

A few more comments :

• Are you talking about the issues of how the dust will couple with sinks? I guess that is something that needs to be addressed sure. It shouldn't be too difficult, especially if the dust particle mass is much smaller than the surrounding gas.
• If we do start adding implementation details and therefore referencing papers, then sure I think we will need some kind of bibliography.
• Regarding the parameters for the dust physics; currently they are in the userguide in its own section. Should this not be together with the other parameters? And if not, and if we decide it's cleaner to keep the parameters in the section that is describing the physics, then should we do the same for other algorithms (e.g. tree, sinks, periodic boundaries, etc..). This would of course lead us to expanding greatly on the different parts of GANDALF, not in itself a bad thing (and would improve the userguide further no doubt) but would obviously be a lot more work.
• Regarding the dust particles in python, we would probably need to modify the snapshot class to contain arbitrary types of particles. I would recommend this would either lead to having a single particle array with particles identified by type with an itype variable, OR each snapshot object would need to allocate memory for each individual type that is used. In either case, it would make sense to remove the 'special' star arrays and have a unified approach to this.

[rbooth200]

Thanks for fixing the end timing calls.

In direct reply to your points:

• Sinks: I'm mean both the accretion of dust particles and the creation of sinks. It would we worth doing some testing to make sure that sink creation isn't affected by the dust. More urgently for disc people would be whether accretion happens correctly. I would think that it might not for the smoothed accretion, but should be ok for the vacuum cleaner sinks... Your wisdom on this issue would be useful.
• Documentation: I completely agree about integrating the dust parameters etc with the rest of the documentation, but perhaps we should do this one we've collected the material together? Content is the most important part, so we should get that sorted. Personally I think that working with this in one section is easiest until we're happy with the content.

Similarly, I didn't add a bibliography since there wasn't one, but it would make sense to use one if we go down that road. I added the description of the algorithms since I thought it would be useful, but we don't have to keep. I do think that for non-standard things like dust, radiation transport etc. we should include references (where possible) so that people can follow it up.

• Particle types in python: Giovanni is going to look into doing this properly...

[rbooth200]

I just fixed a bug that crept into the non-periodic tree forces. All tests seem to be working fine for me now

[rbooth200]

The dust particle read/write in C++ should now be working for both seren formats (needs testing). Column format is not supported at the moment. What should we do about this?

[giovanni-rosotti]

I have modified the python library to handle also dust particle types (which actually required changing a lot of C++ code and only two lines of python code...). I tried using it and it seems fine to me. While doing so I also noticed that the render functions have the "sph" type hard-coded. I guess we actually want to pass that as a parameter as happens in the other plotting functions?

[rbooth200]

That would be good, with a default to sph?

[rbooth200]

Also, it might make sense to support multiple types in the render, although if this is more work it could probably be delayed.

[giovanni-rosotti]

What do you mean with multiple types? Rendering the total density (gas+dust)?

[rbooth200]

Exactly

[rbooth200]

I've just added a note about sinks & dust in the user guide. It would be good if someone could check this.

[rbooth200]

The merge appears successful: Travis checks are passed and test problems look fine.

[giovanni-rosotti]

Great - I think we are ready to merge - I've also added the possibility of specifying the type to render. I will also add multiple types as Richard suggests but it takes a bit longer. I think it's actually best to create a separate issue for that

[rbooth200]

I'm happy with suggestion. There are simple enough work arounds in the mean time, (e.g. the density can be summed).

[dhubber]

Okay, I've run a few more (non-dust) tests to make sure the basic hydro + gravity is still working and it seems fine. I also run a couple with mirror boundaries and seems to work so I'm happy to pull this if everyone else is. At some point, we'll have to do some tests with other particle types (e.g. cdm) but for now, this branch is supposed to be for dust so the others can be done later.

[dhubber]

Although before we pull, I just noticed the warning that this branch has conflicts (presumably problems merging with master now we've merged other pull requests). I guess we should resolve these first.

[rbooth200]

I'll try to rebase the conflicts

[rbooth200]

This build failure is due to asserts in the openmp region. Should this be an issue?

[rbooth200]

It also appears that MPI is now broken in the dust branch. I don't know when this happened though...

[dhubber]

Remove the asserts in the openmp sections. I did the same mistake (which I`m guessing you copied) but then I deleted mine. Not sure about the mpi though. Did you test the dust with mpi before? And also, is this pure mpi or hybrid?

[rbooth200]

Ok, thanks. I assumed asserts in OpenMP should be fine. I'll fix it though.

I've found the bug in the MPI, it's related to the 2nd dust implementation only and is due to the MPI ghosts. It should be ok to fix tonight.

[rbooth200]

I've pushed a fix for the MPI and the asserts in the openmp + mpi build. I've run a few dust and pure gas tests of differing nature and all seems to be working fine in serial, openmp and mpi. Let me know if you have any issues...

[rbooth200]

If you find conflicts when pulling the changes, it is due to the rebase that I did. Seems like this a bad idea on branches that have already been pushed. The best thing to do is to delete your dust branch and pull from the repository again.