fusion-energy / paramak

Create parametric 3D fusion reactor CAD and neutronics models
https://paramak.readthedocs.io/en/main/
MIT License
61 stars 18 forks source link

Full rotation shapes cause problems with cubit and openmc #144

Open cbaus opened 2 years ago

cbaus commented 2 years ago

Let me first congratulate the authors for this fantastic project and the related repos that I really came to appreciate in the past days.

About the problem When creating shapes or reactors with the rotatation_angle=360 the workflow to run in openmc via -> cad_to_h5m with cubit imprinting and openmc_dagmc_wrapper fails. I believe the issue is well known: float precision for angle resulting in 2 overlapping surfaces at the boundary of the rotation start and end of the shape. Or at least that's my guess where the problem lies. Do you experience the same

Possible solution a) I put this issue in the paramak repo even though it could fit in some other repos as well. Is there a way to create the full rotation shapes in another way, so that there are no surfaces that need to be merged? Either a flag or checking if the rotation is close to 360 degrees. b) I could also imagine that with rotation_angle slightly larger than 180 plus a reflective surface the problem might be solved. However symmetry has to be assumed and one has to be careful with the source making this not a universal solution. What do you think?

To reproduce

import paramak 
import cad_to_h5m 
shape = paramak.BlanketFP(50, -90, 270, offset_from_plasma=50.0) 

filenames = shape.export_stp("shape.stp") 

cad_to_h5m.cad_to_h5m( 
files_with_tags=[{"material_tag": "li17pb83", "cad_filename": "shape.stp"}], 
    h5m_filename='./dagmc.h5m', 
    cubit_path="/opt/Coreform-Cubit-2021.5/bin/", 
    cubit_filename='/tmp/dagmc.cub' 
) 

Errors from cubit imprint/merge

[2021-12-10 15:16:22.190] [error] In AcisModifyEngine::unite

ERROR: <stdin>, line 2. In AcisModifyEngine::unite
[2021-12-10 15:16:22.190] [error] ACIS API error number 21072 
ACIS API message = face containment cannot be determined 

ERROR: <stdin>, line 2. ACIS API error number 21072 
ACIS API message = face containment cannot be determined 
[2021-12-10 15:16:22.191] [error] UNITE failed

ERROR: <stdin>, line 2. UNITE failed
[2021-12-10 15:16:22.191] [error] UNION boolean operation failed

...
...

Auto-sizing requested with sizing factor of 5.000000
Calculating Auto Size: |                                                                                                                                                                                                                 |  0%WARNING:  Can't compute automatic size for this model 
 Try setting intervals or sizes manually if resulting sizes are not acceptable 
[2021-12-10 16:13:27.561] [error] Size must be greater than zero.

ERROR: <stdin>, line 13. Size must be greater than zero.
[2021-12-10 16:13:27.562] [error] Size must be greater than zero.

ERROR: <stdin>, line 13. Size must be greater than zero.
Journaled Command: volume all size auto factor 5

Paricles seems to get lost when running openmc. Also plotting the x-y slice is very different than when the shape is only rotated 395.5 degrees.

RuntimeError: Maximum number of lost particles has been reached.
shimwell commented 2 years ago

Super to hear the package is useful.

This is not really a solution but just a bit more info on the bug

The Paramak makes use of Cadquery and CadQuery performs a clean operation by default that might help here. I would have thought that the duplicated surface due to a 360 in the Z axis rotation gets removed

Here is a simple example with showing that we get 4 faces for a full 360 rotation

Screenshot from 2021-12-10 11-43-36

However the BlanketFP is a more complex shape and has another place where this duplicate surface can happen.

Does the error go away if the start and stop angles are changed from this

shape = paramak.BlanketFP(50, -90, 270, offset_from_plasma=50.0) 

to this

shape = paramak.BlanketFP(50, -80, 270, offset_from_plasma=50.0) 
cbaus commented 2 years ago

Yes it goes away for the proposed values. BallReactor behaves similarly. I can set the rotation up to ~359.6 degrees. If I set it larger, errors appear. As you said, it might just be a few of the volumes in the reactor giving trouble.

shimwell commented 2 years ago

I shall look into this more.

As a temporary fix there is a method of generating those reflective planes automatically with the openmc-dagmc-wrapper. So perhaps a 90 degree model wouldn't be too bad

https://github.com/fusion-energy/openmc-dagmc-wrapper/blob/3806834e22a96f5cd58d95ff5407e83f25efebae/openmc_dagmc_wrapper/Geometry.py#L18