enzo-project / enzo-dev

The Enzo adaptive mesh-refinement simulation code.
Other
80 stars 95 forks source link

MaximumGravityRefinement and FastSiblingLocator causes inconsistency #26

Open brittonsmith opened 10 years ago

brittonsmith commented 10 years ago

Original report by dcollins4096 (Bitbucket: dcollins4096, GitHub: dcollins4096).


The MaximumGravityRefinementLevel will cause incorrect results due to the fact that the SiblingList is not repopulated. Please stop using MaximumGravityRefinement until this is resolved.

PrepareDensityField, line 116, level = min(level, MaximumGravityRefinementLevel); and then the grid array is set from the level (which in my case was level=1, while I had MaximumRefinementLevel = 2)

It then calls

PrepareGravitatingMassField2a(Grids[grid1], grid1, SiblingList, MetaData, level, When);

where everything except SiblingList is referencing level=1. But SiblingList ws generated from level=2. Which is then passed into PrepareGravitatingMassField2a, which does some particle overlap stuff, namely calling CheckForOverlap on the GridList (on level 1) and things in the SiblingList (from level=2)

Possible solutions: -- Storing the SiblingList of MaximumGravityRefinementLevel (as a global? Passing it through the recursive call to Evolve Level)

-- Recomputing the SiblingList as it is needed

This will need testing. In the mean time, please do not use MaximumGravityRefinement, it will lead to incorrect results.

brittonsmith commented 10 years ago

Original comment by James Larrue (Bitbucket: james_larrue, ).


@dcollins4096 , in your second paragraph, did you mean to say you had MaximumGravityRefinementLevel = 1 and level=2 ?

Here is my suggestion, at a high-level, for fixing this. If the high-level concept looks fine, I can work out the details next week (back on Wednesday).

== High-Level Concept ==

1) Create an object that, if given a level number, can provide pointers to the generated grid array or created sibling list for the specified level.

2) Modify EvolveLevel and PrepareDensityField to use the object created in item 1. (There are other places that could be modified also, but see my comment below about making fewer changes, if possible. I suggest changes to other locations GenerateGridArray and CreateSiblingList be a “Phase 2”.)

After working on this for a while, I see a lot of changes that could be made if we implemented such an object, but I would rather make fewer changes, make sure things still work, then come back to do more changes, if we wish. Therefore, when I start on the details, I will not necessarily strive to implement things in the “best” way, but rather in the most “appropriate” way.

brittonsmith commented 10 years ago

Original comment by James Larrue (Bitbucket: james_larrue, ).


I implemented a solution based on my above high-level concept. It seems to be working, but I'd like to get some feedback before I make a pull request, since I expect there will be changes after the feedback.

The names I used for classes and methods are those that make sense to me, but they may not make sense to anybody else. In particular, I'd appreciate it if some people can look over the class and method names, then let me know any suggestions for changes.

My solution can be browsed here:

https://bitbucket.org/james_larrue/enzo-dev-larrue/src/9d88f9a3bd5e2cf55a63c97174c96da286046028/src/enzo/?at=week-of-code

I ended up creating 2 classes, rather than the one class referenced in the high-level concept of my previous post: EnzoDomain and GridLevel.

The EnzoDomain represents the domain used by the simulation and acts as a container for the GridLevel objects, creating and destroying them as needed. The EnzoDomain is currently implemented as a singleton, but if the current changes are all working fine, I'd like to revisit that aspect. More details about this class can be found in the code comments, or in the generated documentation I put at:

http://www.diopolis.com/enzo/apidocs/classEnzoDomain.html

The GridLevel represents a grid level. It will generate the grid hierarchy array for the grid level when needed, via the GenerateGridArray function, and save the array so it only needs to be generated once. Likewise for the sibling grid list array (created in the CreateSiblingList function). More details about this class can be found in the code comments, or in the generated documentation I put at:

http://www.diopolis.com/enzo/apidocs/classGridLevel.html

The EvolveHierarchy function was modified to initialize the EnzoDomain singleton, then destroy the EnzoDomain singleton once the main loop is complete.

The EvolveLevel function was modified to get the GridLevel object for the current level from the EnzoDomain singleton, then use the GridLevel's method to access the grid list array and sibling grid list. By using the GridLevel rather than the GenerateGridArray/CreateSiblingList functions, we will not need to call the functions a second time in PrepareDensityField. The GridLevel object will remember the values. After the main loop of EvolveLevel, in the “Clean Up” section, the current GridLevel object is destroyed, freeing up the resources it was using.

In PrepareDensityField, the variable 'level' is set to the minimum of level and MaximumGravityRefinementLevel, as before, i.e. no change here. Then the GridLevel object for 'level' is obtained from the EnzoDomain (the GridLevel object will already exist, since level <= reallevel and the GridLevel for reallevel was created by EvolveLevel, so EnzoDomain is just passing the correct GridLevel object) and both the grid hierarchy and sibling grid list are obtained from the GridLevel object. Since EvolveLevel would have previously used the grid hierarchy and the sibling grid list for GridLevels with level <= reallevel, those values will not need to be regenerated when PrepareDensityField asks for them.

My tests so far indicate that the code is running as expected, providing the results expected. The speed and memory consumption also seems to be the same as before the changes. After I incorporate any feedback, I'll make a pull request and hopefully some others can test it.

Once the solution is running correctly, I'd like to make the following further changes:

1) change hydro_rk/EvolveLevel_RK2 to have it match EvolveLevel. (I don't think Unigrid_EvolveLevel.C is used, in which case it doesn't need to be modified, but maybe it should be deleted.)

2) remove the singleton nature of EnzoDomain and instead, pass the current GridLevel as a parameter to EvolveLevel. This can reduce the number of parameters passed to EvolveLevel from 6 to 4, since GridLevel provides access to 3 of the current parameters.

3) pass the GridLevel as a parameter to PrepareDensityField, which will reduce the number of parameters passed from 5 to 2, since GridLevel provides access to 4 of the current parameters.

Please let me know if those additional changes make sense.

brittonsmith commented 10 years ago

Original comment by Nathan Goldbaum (Bitbucket: ngoldbaum, GitHub: ngoldbaum).


@james_larrue go ahead and make a pull request and mark it [WIP] (work in progress). It's much easier to discuss code changes in the context of the PR interface.