Created Monte Carlo solver in C++

[ljeure]

I created and tested a Monte Carlo solver in C++. It runs mostly the same way as the python solver did, but much faster. It does not display a heat map of neutron densities. I also fixed an issue we had with sampling absorptions in the python version of the code.

I tested k-values for a geometry with homogeneous materials and reflective boundaries. I tested both a 1-group and a 2-group problem which had expected k-values of 1.5 and 1.4 respectively. Using one million neutrons, the simulation gave k-values of 1.500 and 1.40 respectively.

[geogunow]

@ljeure I have just finished my first review of this PR. Overall, it's looking pretty good. I made a few stylistic comments and we might want to also rethink the way we save fission sites.

[geogunow]

@ljeure I've finished my first review of this PR. Much more to come...but make these changes first and then I'll look at it again with fresh eyes.

[geogunow]

@ljeure what happened to all the C++ files??

[ljeure]

@geogunow I'm not sure; I added them back.

[geogunow]

I've finished my second review of the code. Once you make the requested changes, I'll try to bring in a second reviewer from either the OpenMC and OpenMOC team to give this a proper review.

[geogunow]

@friedmud has agreed to review this PR

[friedmud]

Taking a look

[friedmud]

One more general comment: typically in C++ you want your public: sections to be the first things in your class. The reason why is that that's what the outside world cares about seeing... i.e. they want to see what interface is available. The detail of what members are in the class should be at the bottom. So a class should look like this:

class Stuff
{
public:
  ...
protected:
  ...
private:
  ...
};

[ljeure]

Thank you so much for all your help!

[friedmud]

@ljeure no problem at all! Always glad to help budding object-oriented programmers!

Overall, the code looks really good. Proper use of objects and encapsulation without going overboard (always a tough thing to learn when to stop making objects and just get stuff done!).

Good job!

[friedmud]

One more thing: you might think about using rand_r() instead of urand()... and store the "state" for the random number as part of the Neutron. So to get a new random number you would do:

neutron.rand()

Also: give each neutron an ID (like 0,1,2,3,4...) and use it as the seed (or as a modification to the seed... like add it to the seed) for rand_r()

What this does is this: it gives each neutron (as identified by its ID) a unique history... regardless of what order the neutrons are handled in. This, along with using rand_r() would allow you to use OpenMP for threading this application so it could use more processors.

Let me know if that's not quite clear ;-)

[ljeure]

I've added the neutron id. What do you mean that I should use it as the seed for rand_r()?

[friedmud]

The way that it works is that you put an unsigned int _seed in your Neutron class. Set the initial value of that seed in the initialization list to be id + global_seed where global_seed is a number you can control from the input file or command-line (so that the results are always reproducible).

Then, in the constructor for Neutron call rand_r(&_seed) to get the series started. The way rand_r() works is that it stores the "state" of the random number generator in the thing you pass in as the seed. So, now add a member to Neutron called rand() that simply does return rand_r(&_seed) and use that every time you need a random number.

What this does is two fold:

1. rand_r is "re-entrant" which means it is a thread safe random number generator
2. It gives each Neutron its own unique history... so the results are reproducible... even if the Neutrons are handled "out of order" (for instance, if you do threading or use MPI to split up the number of neutrons).

This is all VERY optional... but whenever I see calls to urand() it makes me cringe a bit because of how NOT thread safe it is... and how hard it is to get reproducible results from it if you do any parallel processing!

[geogunow]

@ljeure This PR is looking really good! After a few small things are addressed it should be ready to merge.