Silent failure from Paramak geometry

[makeclean]

There have been a number of cases where a silent failure from the entire processing chain has lead to complete simulations being run, but with the incorrect outcome. In this case, the final geometry that is produced does not imprint correctly, despite the previous geometry differing in only one value. This failure to imprint and merge (not a failure in Cubit) is due to the two surfaces which we would consider being identical (in this case only 1 curve is different) - different in this case due to what looks to be a roundoff issue. This non-imprinted and merged surface, due to be being topologically different leads to slightly different faceting on what should be a merged surface, and thus in transport particles 'tunnel' through the other volume. When running this in an automated workflow, it would be very difficult to detect, in this case it was detected by scanning through a range of inputs, and noticing a noticble change in the 'background' result.

Its difficult to be robust to this failure, I know we've spoken about this before. I only have two possible solutions as a remedy, 1) before a real calculation is run, put a point source at the centre of the geometry and run OpenMC and record the flux in every cell, in void we should get a score in every cell 2) we can probably estimate a minimum number of merged surfaces, this can be checked after imprinting and merging.

I'll ask if the folks who produced the underlying script can share it here so you can reproduce.

[shimwell]

Interesting find and many thanks for reporting.

I'm wondering if the DAGMC overlap checker finds this error in the resulting h5m file.

conda install -c conda-forge dagmc
overlap_check dagmc.h5m -p 1000

Perhaps I should look into storing the floats differently to avoid rounding errors and use something like the Fractions class to store floats Fraction.from_float(0.3)

If anyone is able to share even just the code for the two paramak.Shape objects that break that would be useful.

Ideally we should also make this tunneling feature result in a failure during particle transport. But I understand this is difficult to implement in DAGMC

[makeclean]

With respect to the floats, you're always going to have to hand off a floating point number to OCC eventually, you're just delaying the inevitable, but it would be good to understand in this case where the problem lies. in this case it just so happens the rotation of the components is cos(pi/2) which should obviously be 0, but invariably in floating point if we try to determine it, it will be some small non-zero number e.g. 1e-16 or something.

[shimwell]

It is a bit of a hack but in the short term one could put in round(math.cos(math.pi/2), 5) instead of math.cos(math.pi/2) to get 0 returned instead of 1e-16. I guess 5 decimal places is enough

[makeclean]

Not a hack, but if you do it here, you should it for any geometric operation and consider introducing a tolerance.