Icsilcc

[dhubber]

Although this branch still needs some revision and testing before merging, I've made this pull request as a place to comment and discuss anything on the branch. Contains :

• A total refactoring of the initial conditions class, with each individual IC in its own class. This is principally to make it easier to use regularisation/relaxation algorithms to generate more complex initial conditions (i.e. which are not simply uniform or spherically symmetric)
• Ewald gravity added to all N-body classes. Not fully tested (suspect I've forgotten to add code in one or two places) but most of the work is done.
• Added some new initial conditions, such as ICs for SILCC like simulations and filaments
• Added random supernova driver, for testing, and code for SILCC supernova driver for future

[giovanni-rosotti]

I just wrote a series of comments; I haven't tested the nbody yet.

I gave a look at travis, and the failure is coming from the sedov test. Could it be that you have left the supernova switched on?

[dhubber]

I think the reason is that the IC class is trying to do far too much on its own. I feel like the base IC class should be doing just the bare minimum: It should provide some utlity functions for the derived classes such as AddRandomSphere, but no more. I'm not even convinced that functions like AddBinaryStar and GenerateTurbulentVelocityField belong here either.

There are one or two functions that I've moved to other classes (e.g. the Lane-Emden functions are moved to the PolytropeIc class) but some, like the two you mention (AddBinaryStar and GenerateTurbulentVelocityField) are needed by more than one IC class so should (at least for now) live in the IC class.

I am not convinced that CalculateMassTable should be part of this object and I'm pretty convinced that ContactDiscontinutiy, GaussianRing and SedovBlastWave should not.

If it were to live anywhere, then this is the place it should be (where else would it go without making some other class family?). However, it's now been deleted anyway since it's not needed now the Monte-Carlo + regularisation ICs are working okay.

As I mentioned in a previous comment to Giovanni, I've now moved the GaussianRing code to its own class and file. But I'm not sure what you mean about the other two since they're already part of their own class (unless I forgot some bit of code somewhere?).

Also, is there code to generate a ndim-dimensional glass in here somewhere?

In a way, but not in an ideal way at all. Hopefully we can chat about this when we meet soon since it's the last major issue of the IC component of this branch.

[dhubber]

So, I've run the two proposed tests for the N-body Ewald gravity. The first one (a uniform line of stars) works fine and has zero acceleration so that's great. I also got the Jeans test with stars to work. However, that reveals a problem since the accelerations seem to be inverted with some weird discontinuity or something. Basically, I think there's a wrong sign somewhere so will get onto tracking it down. But once it's found and I can get it to work with this test, then I think that will be that for the N-body Ewald.

[dhubber]

Results for the Jeans test for the N-body case using the Leapfrog-KDK. Seems to work pretty well I think so I'm quite happy. The case I was mentioning earlier was for the Hermite scheme so there is perhaps a problem there that I want to track down first. Once that's done, and I've cleaned the new code up, I'll push everything later!

[giovanni-rosotti]

Tha's great. Can you make this into a test? Presumably it doesn't take long to run

On Wed, Dec 14, 2016 at 4:11 PM, dhubber notifications@github.com wrote:

[image: ewaldnbody] https://cloud.githubusercontent.com/assets/4923605/21189778/1233e71c-c220-11e6-886e-81322830d143.png

Results for the Jeans test for the N-body case using the Leapfrog-KDK. Seems to work pretty well I think so I'm quite happy. The case I was mentioning earlier was for the Hermite scheme so there is perhaps a problem there that I want to track down first. Once that's done, and I've cleaned the new code up, I'll push everything later!

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/gandalfcode/gandalf/pull/114#issuecomment-267076006, or mute the thread https://github.com/notifications/unsubscribe-auth/ABdSgoDuayFDSN4jjwk0UkcrGm_DmlbXks5rIBUsgaJpZM4LFchZ .

[dhubber]

Yes, I had the same idea to make it into a test (since I'd already done 80% of the work for that). I've added it to the tests/nbody_tests sub-directory (called nbodyjeans.dat with corresponding python script). I've still not technically tested this with both gas and stars but the Ewald now works with both separately so that element of this branch is now complete I believe (finally).

[giovanni-rosotti]

Cool. So now it's Richard's turn to refactor the IC class. David, I'll turn the python script you've added into a nose test so that it can be run automatically by travis. Is there anything else left in this branch to do?

[rbooth200]

Can you point me to where the E-wald correction is done in th N-body force calculation? I remember having some corcerns about it that I'd like to check

[rbooth200]

I'm not convinced that we will get the Ewald correction working with the Nbody Hermite integration without a fair chunk of work. The reason is I don't see a way to compute the \dot{a} term from the Ewald correction properly. This needs to be done carefully since some of the terms in the Ewald sum will change sign. We'd need to implement a new function for that...

[dhubber]

Is there anything else left in this branch to do?

Like I mentioned, it would be nice to test hydro + nbody with Ewald to be sure it all works, but that's not essential right now. So there's only really Richard's refactoring plus whatever I have to do later to 'fix' it and then it'll be ready.

Can you point me to where the E-wald correction is done in th N-body force calculation? I remember having some corcerns about it that I'd like to check

Unfortunately, it's done in several places (and I have to admit there's some code repetition that should be eliminated in truth); every N-body child class in principle can have its own 'CalculateDirectGravForces', 'CalculateDirectSmoothedGravForces', 'CalculateDirectHydroForces' and 'CalculateAllStartupQuantities' which all take simbox and ewald as parameters to compute the Ewald corrections.

I remember you mentioning about some of them not having 1st time derivative terms computed (with the Ewald correction). This is partly because I don't know what this term is and if it's possible to determine simply (it might be in the literature but I'm not aware of it). Luckily it is of less importance due to the 1/r^3 nature of the time derivative but would still be good to do it if possible.

[rbooth200]

I remember you mentioning about some of them not having 1st time derivative terms computed (with the Ewald correction). This is partly because I don't know what this term is and if it's possible to determine simply (it might be in the literature but I'm not aware of it). Luckily it is of less importance due to the 1/r^3 nature of the time derivative but would still be good to do it if possible.

I'm not really willing to accept this solution to be honest. I'd much prefer it if we just disallow Hermite integration with periodic systems.

[rbooth200]

So I've made some changes to the particle regularization that should make the interface a bit more flexible. At this stage I've only implemented a bare bones version, with a single example of how it should be used, based on the Silcc ICs. The new implementation means that it should now work with the meshless as well as SPH.

David, could you test? Having never used it, I'm not sure what to expect...

[dhubber]

Great, thanks. I'll have a look at it asap!

[dhubber]

Richard, a few questions about the regularisation

1) Atm, these changes essentially just encapsulate all of the regularisation functionality in a new class and then this is called immediately after the IC-Generator step (but before the PostInitialConditionsSetup)? So we still need to add in the step of setting the particle properties (e.g. velocity, internal energy, etc..) after the regularisation step. But this probably can be done quite easily now immediately after the regularisation itself (inside the if-statement block).

2) For every IC that we wish to use the 'regularizer' (btw, I will un-americanise that later ;-) ), I'm guessing we need to create its own ICRegularizer class including the density function? I'm guessing this also means that currently the GetValue function is not being used (although it or something similar will be needed to set particle properties as discussed in Q1).

3) One of your aims was to be able to generate a glass easily to be used for making ICs. I'm guessing it should be easy to do here but we should probably make some 'AddGlass' function similar to the other random-particle/lattice functions for ease of use?

4) "ParticleRegularizer(simparams, simbox)(hydro, neib, nbody, *reg_func) ;" What the hell kind of Satanic C++ witchcraft is this? :-) It seems like it's maybe calling a constructor but not actually making an object while also calling an operator??

Anyway, in the meantime, assuming I understand it, I'll have a play around checking it's working okay.

[giovanni-rosotti]

While you two guys have fun with c++, just an update on why this is crashing on travis. The issue comes from MPI in combination with pure Nbody simulations after I've added the nbody jeans test. I've fixed it in the fix_stars branch, which however is based on neibmanager. I can't replicate the fix as it's in process parameters which was heavily modified, so we would need to modify it again. Shall I skip the offending test in this branch (only when using MPI) for now? Once this has been merged I can then rebase fix_stars onto master and reactivate it there...

[dhubber]

Sure, if you're confident you know what the problem is and can fix it later (presumably once we've merged/rebased the various branches into master) then we can skip it for now to get this branch done asap!

[giovanni-rosotti]

Cool, I'll do that asap then. Just to warn you, merging this will be a nightmare as there are zillion of conflicts... This branch basically modifies every single line of code in GANDALF!

[dhubber]

Merry Christmas to all of us then! >.<

[giovanni-rosotti]

Just to let you know, I've now fixed the tests...

[rbooth200]

I've rebased this branch. Wasn't too painful...

[dhubber]

Richard, I said a couple of days ago I was happy with the regularisation classes, but on some reflection I have a few (quite big) concerns that I'd like to address now before continuing :

(a) We're splitting a single IC over two classes (the IC class itself and the regularisation function class) which is perhaps not the best idea. This is not the end of the world for these simple cases, but if we're doing something more complex like loading in ICs from a file (e.g. from grid/AMR ICs) then it's a problem. The functions that create (e.g. Generate) and define the ICs (i.e. GetValue, Density, etc..) should really all be in one class, the IC class. (b) The IC class is one of the few classes, maybe the only class, that we should expect users (the ones who want to get their hands dirty with C++ that is) to create their own versions. Therefore it should be made as simple as possible and I feel this current way is overly complex. I don't mind the complexity being hidden away so the (casual) user does not see it or need to bother with it, but in this form they must make their own inherited versions of these various classes, functions, etc.. so they are exposed to it directly.

My basic question/suggestion is can this be simply modified to use the IC object instead of making a brand new RegulariserFunction class? The GetValue/Density functions will be in the IC class as before and perhaps the IC class can be passed in instead.

Btw, on the subject of how the current regularisation works/performs, I only had to fix one small bug (essentially it goes wrong if h is not initialised to a non-zero value) but otherwise it works fine so there are no other issues on that front.

[rbooth200]

Hi David,

Firstly, we should do away with the global GetValue functions. I genuinely believe that a solution that relies on these cannot be generic, and will cause problems later. Any class that does something useful must not be forced to work with them.

With regards to making life easy for people to write their own IC class, I personally don't mind making them write a second. If they are already writing one derived class, then surely they have the skills to write a second one?

But I do accept that there will be cases (IcSilcc may be one of them) where the RegularizerFunction might be better off calling a density function on the IC class. In this case, subclassing DefaultRegularizerFunctionBase to an "ICDefaultRegularizerFunction", which just forwards the density call to an IC->density(part) would work fine.

One option to make this easy for user to implement is use a pointer to an IC base class, and then call IC->density virtually. However, this isn't good practice since most IC classes don't need a density function, which spreads like a cancer if it's declared in the base the class. So I think a template ICDefaultRegularizerFunction would work nicely, which people could use.

I will push an example for you to peruse...

Rich

[rbooth200]

Ok, Have a look at the changes that I made to DefaultRegularizerFunctionBase ans SilccIc. Hopefully that is simple enough?

[dhubber]

Firstly, we should do away with the global GetValue functions. I genuinely believe that a solution that relies on these cannot be generic, and will cause problems later. Any class that does something useful must not be forced to work with them.

Okay, I still do not understand this argument at all. What problems exactly?

With regards to making life easy for people to write their own IC class, I personally don't mind making them write a second. If they are already writing one derived class, then surely they have the skills to write a second one?

Because of the reasons I outlined, this is one case where easy of use must trump over design philosophy. I do not want us to force users to become C++ gurus in order to use the code, because it is far more likely they will (a) give up, or (b) send us all e-mails about how to do it. And I don't want either scenario tbh! This has to be made as simple as possible!

[dhubber]

ok, just saw you made another commit as I was posting. Will have a peek at that now!

[rbooth200]

Re. the 'Ease of Use' argument.

I think is a good reason to create an AddUnifromGlass or InsertUniformGlass function. It will have to be slightly different to the AddCubicLattice functions though, since as we discussed it might be near impossible to add multiple glasses this way.

[dhubber]

ok, that's definitely an improvement from the previous version from the general users point of view.

Before getting into that too much, I'm now just realising, is your objection to the GetValue function that it can return everything (i.e. rho, v, u, etc..)? You keep saying 'global' (but not actually explaining your objection fully enough), but is it simply that you'd prefer them to be separated into say GetDensity, GetVelocity, GetInternalEnergy, etc. functions in the Ic (or whatever) class?

[rbooth200]

No, my objection is related to trying to use it to set important particle quantities. As with the discussion for CreateParticles, this function doesn't easily allow users to ensure that the particles are valid - they have no way of knowing what the need to set. While one could come up with a solution to the latter problem of knowing what needs to be set, it's not the best approach. There are just too many edge cases that can crop up this way.

It would be much better to have a ResetParticleProperties function on the IC class that handles everything. After all, this even supports users trying to do things we haven't even thought of, since they could do just about anything they want in that function. This way it's up to the IC writer to ensure the particle is valid, but then again they are the person who knows for sure what is valid for their IC.

On Sun, 18 Dec 2016 18:04 dhubber, notifications@github.com wrote:

ok, that's definitely an improvement from the previous version from the general users point of view.

Before getting into that too much, I'm now just realising, is your objection to the GetValue function that it can return everything (i.e. rho, v, u, etc..)? You keep saying 'global' (but not actually explaining your objection fully enough), but is it simply that you'd prefer them to be separated into say GetDensity, GetVelocity, GetInternalEnergy, etc. functions in the Ic (or whatever) class?

— You are receiving this because you commented.

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

[dhubber]

No, my objection is related to trying to use it to set important particle quantities. As with the discussion for CreateParticles, this function doesn't easily allow users to ensure that the particles are valid - they have no way of knowing what the need to set. While one could come up with a solution to the latter problem of knowing what needs to be set, it's not the best approach. There are just too many edge cases that can crop up this way.

It would be much better to have a ResetParticleProperties function on the IC class that handles everything. After all, this even supports users trying to do things we haven't even thought of, since they could do just about anything they want in that function. This way it's up to the IC writer to ensure the particle is valid, but then again they are the person who knows for sure what is valid for their IC.

Okay, I partly see what you're saying now (at long last, lol) in terms of objections to a GetValue function. But (as you say in the second comment) it's up to the user to decide how to properly construct the ICs and there's no failsafe method that prevents them doing something stupid. And more importantly (as least for my concerns about ease-of-use), this extra layer of regularisation classes seems unnecessarily more complex than having just a (Get)Density function (that is used by the MC and regularisation scheme) and a (Re)SetParticleProperties function in each Ic derived class. Limiting the Ic class to having just those two functions would seem a much better compromise solution to me than this extra class layer.

Anyway, for now I will have a look at this (Re)SetParticleProperties function idea and see how it works in terms of with and without regularisation.

[rbooth200]

So yes, i agree with you that just these two methods on the IC class is the best way forward.

The reason for the extra layer, as i see it, is that we might want to do particle regularisation in another part of the code, e.g. when injecting particles. That's why I wanted to separate the regularisation from the IC class. The cleanest way to do this is with an extra layer.

On Sun, 18 Dec 2016 18:33 dhubber, notifications@github.com wrote:

No, my objection is related to trying to use it to set important particle quantities. As with the discussion for CreateParticles, this function doesn't easily allow users to ensure that the particles are valid - they have no way of knowing what the need to set. While one could come up with a solution to the latter problem of knowing what needs to be set, it's not the best approach. There are just too many edge cases that can crop up this way.

It would be much better to have a ResetParticleProperties function on the IC class that handles everything. After all, this even supports users trying to do things we haven't even thought of, since they could do just about anything they want in that function. This way it's up to the IC writer to ensure the particle is valid, but then again they are the person who knows for sure what is valid for their IC.

Okay, I partly see what you're saying now (at long last, lol) in terms of objections to a GetValue function. But (as you say in the second comment) it's up to the user to decide how to properly construct the ICs and there's no failsafe method that prevents them doing something stupid. And more importantly (as least for my concerns about ease-of-use), this extra layer of regularisation classes seems unnecessarily more complex than having just a (Get)Density function (that is used by the MC and regularisation scheme) and a (Re)SetParticleProperties function in each Ic derived class. Limiting the Ic class to having just those two functions would seem a much better compromise solution to me than this extra class layer.

Anyway, for now I will have a look at this (Re)SetParticleProperties function idea and see how it works in terms of with and without regularisation.

— You are receiving this because you commented.

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

[dhubber]

ok, still not 100% convinced about the extra class layer (at least how it's done atm; need to think more about it) but it's not super important right this minute. I just tried adding what we talked about (the 2 functions alone in the Ic class). Seemed pretty straight-forward so far. Still need to try it on a better example though (since these ones had zero-velocity ICs with no variations).

The other thing I need to try next is to add the Smoothed Density Ics, since I'll need those pretty soon. Just need to think how this all clips together in the new classes.

[rbooth200]

Looks good.

The only complaint I have at the moment is that the regularization only works with one particle type. We really want this to be able to simulataneously regularize different density distributions. Can you add that in? I think it means changing the interface to GetDensity, and having the Regularizer function return zero for neighbours that are not included in the density mask.

[dhubber]

ok, so I just pushed both a GetDensity with particle types and GetSmoothedDensity options. I just realised my GetDensity with types is not really done correctly (since it doesn't have the particle type itself) so going to re-do it now! Definitely worth trying some set of initial conditions (maybe a dust test?) to be sure this is all working fine!

[rbooth200]

As a way to test this, can we do a random cube type test? Instead we'd just use the regularization routines to do the regularization and not take any timesteps (or only 1 perhaps). Afterwards we'd check that the standard deviation of the density is below some threshold.

Dust could be included as a second fluid to regularize in this test and we should check that it works with the meshless, perhaps both the smoothed and unsmoothed varients too?

[dhubber]

As a way to test this, can we do a random cube type test? Instead we'd just use the regularization routines to do the regularization and not take any timesteps (or only 1 perhaps). Afterwards we'd check that the standard deviation of the density is below some threshold.

Okay, I assume you're talking from the point of view of having an automated test here? That perhaps would be sensible sure. I'm actually not sure what happens if Nstepsmax = 0. Perhaps it will only run the Setup and not bother with the MainLoop at all?

[rbooth200]

I'm actually not sure what happens if Nstepsmax = 0.

I don't know either, but Nstepmax = 1 would be fine too.

[dhubber]

Why are we using an int here to decide the type? Surely it should be an enum ptype?

Because I did try to use the enum originally but it was complaining for some reason so switched to int for an easier life. Actually we were mixing enums and ints before anyway in various places so should stick with one of the two (preferably enums sure, but ints are not end of the world).

[rbooth200]

I agree its not the end of the world. Having had a look, the problem is because ptype has gone and we only have parttype left.

I'd suggest that we rename parttype -> ptype and make a consciencious effort to change int to enum ptype as we go about our business in the future?

[dhubber]

I'd suggest that we rename parttype -> ptype

No!! Part of the issue (at least that I was worried about) was I was using ptype often as a variable (both before and after these changes) and one of the enums was also called ptype. So, now I've used 'ptype' as the variable name (both in the Particle class and in various subroutines that need it) and the enum is now called parttype. There's still some messy variable naming I saw in other parts of the code but I tried to get this bit consistent at least. It's possible some confusion with the naming is what caused my compiling issues that made me use ints. If we get this all consistent, I'm sure changing the ints back to enums won't be a problem.

[rbooth200]

Ok, so we should just conscienciously change int ptype to enum parttype ptype...

[rbooth200]

When you're happy that this is working as intended, it can go. I'm more or less happy with it as it is.

[dhubber]

Closing this pull request (and will soon delete this branch) since all commits are in final_paper_tests (which has additional fixes).