dhubber
commented
8 years ago Okay, so first of all, I'm guessing this is implicitly about the meshless scheme since I assume 'conservative' is supposed to be in the context of a conservative FV scheme (and I'm not sure how that would relate to SPH). In that case, yes the Saitoh & Makino scheme is not conservative. I've been trying to think of ways to modify it to make it conservative but nothing satisfactory just yet.
Btw, I'm assuming this branch split off from timing since there seems to be plenty of the same code changes there? In fact, is there any different code since the above discussion talks about ideas and not about what has actually been done? Just wondering since I was thinking maybe I'm wasting time looking at this same code again :-) Actually, looking at the commit history there are a few commits more in this branch. Maybe we should (at some point soon?) accept the timing pull request and rebase this to master so we can see more clearly the code differences before delving too deep into this?
Some immediate questions/clarifications that come to mind
- You suggest a gather-based levelneib computation first and then a 'full' computation in the gradient calculation. So I assume the 'scatter' element of it is done in the gradient calculation, since obviously one is needed somewhere?
- So the difference here is there will be some do-while loop, similar to that in SphSimulation, but it will iterate over the density/gradient part, near the bottom of the main loop in the meshless?
- I'm assuming this is simply a way of fixing the Saitoh-Makino to work and can work without the need for the Springel-type timestep treewalk? If we did do a tree-walk, we would also need to include the MPI pruned trees presumably to make sure nothing funny could happen from other nodes.
Assuming I understand it all correctly, my first impression is that it could work sure. I'll have to think about it some more and look at the ordering of the main loop but seems to make sense. I have just noticed that the UpdateAllProperties is called twice in the main loop (with self-gravity) but not sure if that makes much difference here.
rbooth200
giovanni-rosotti
I'm beginning this pull request to spark discussion, rather than because it is ready to merge. It contains two methods to limit the timesteps:
The first one works by finding the largest dt such that r_ij * vsig_ij > dt for all pair in the simulation. This is done in an extra tree walk, which will need making MPI compatible later.
The second one just checks for particles that have too big timesteps compared to their neighbours and reduces their timestep at the next synchronization point. It is not exactly conservative since dQdt is used to predict dQ, but it seems to work fine for the sedov test.
In doing this work, I noticed that we don't enforce the levelneib heirachy properly at the moment. While the issue of reducing the time-step particles in the middle of their step conservatively is near impossible (although limiter 2 does is this in a non-conservative way). However, in addition we don't enforce the hierarchy at the start of the timestep either, which I think is a significant flaw. The reason that the heirachy is not enforced is that levelneib is computed before the new timesteps are evaluated, and after the new timesteps are calculated the levels shift, leading to some particles needing a smaller timestep than the one they are on because some of their neighbours are now on smaller timesteps.
I think that enforcing the hierarchy would improve the results when no extra timestep limiters are used or when the non-conservative limiters are used. This is because the fluxes are computed at the start of the timestep. Thus if we can predict that we are going to need a smaller timestep at the beginning of it we can reduce it immediately allowing a fully conservative update of these particles. With the global limiter switched on fewer particles tend to have level differences greater than requested, but some still do.
In SPH we enforce the hierarchy by iterating the whole timestep for particles on levels that are too big, but this is difficult to do correctly for the meshless because tracking the update of dQ is difficult. I have included an attemp at this by iterating on only the flux calculation but it doesn't work properly. I believe it would be much easier to iterate on the gradient calculation until the timesteps are set properly, but currently the gradients are done before the timestep calcution. Placing it after the timestep calculation has the issue that levelneib is not known in advance since we calculate it during the gradient calculation.
I would think that the most effective solution would be as follows: 1) Compute a gather based estimate of levelneib during the density calculation. 2) Compute the new timesteps. 3) Compute the full levelneib during the gradient calculation. To do this we would need to iterate on particles that had their timestep reduced to ensure it is propograted to all its neighbours. Step 3) could also be done in another tree walk to avoid recomputing the gradients.
The questions I have are as follows: How important do you think this is? Should we expend the extra cycles (probably only a few to 10 per cent) on making it more accurate? How should we implement it? I believe before the flux calculation is the only reasonable way.