Implementation of Triply Periodic Minimal Surfaces (TPMS) in OpenMC #3292
Conversation
…ve - nothing is stable so far
…t everything in tpms file.
…ithm by adding a max_range parameter. Modified the way cells compute min distance to surface. Added a example of a high cst TPMS that needed that modification.
…acxing algorithm var computations. improved API user friendly routines.
There was a problem hiding this comment.
Here are some preliminary comments on your code. Regardless of the utility of the functionality being introduced, I think there’s room for improvement in how the data is structured. In particular, some existing OpenMC functionality could be leveraged—such as using the XML-to-xtensor reader instead of building custom string-to-matrix readers.
I also recommend reviewing the developer style guide, especially regarding function naming. The use of camel case is inconsistent with our convention, and keeping things uniform improves readability and maintainability.
Lastly, have you considered using mesh-based geometries for these exotic surfaces? Is there a measurable speed advantage to the current approach? Given enough complexity, mesh-based ray tracing might actually be faster than handling intricate ray-surface intersections.
| # boost library | ||
| #=============================================================================== | ||
|
|
||
| find_package(Boost 1.85.0 REQUIRED) |
There was a problem hiding this comment.
It appears you have included an implementation of bisection root-finding here, which I assume is why you had Boost here in the first place. The code doesn't need the Boost library any more, right?
| @@ -52,6 +56,7 @@ class Surface { | |||
| //! \return true if the point is on the "positive" side of the surface and | |||
| //! false otherwise. | |||
| bool sense(Position r, Direction u) const; | |||
| virtual bool is_tpms() const { return false; }; | |||
There was a problem hiding this comment.
Ideally, the base class shouldn't have to have any methods like this. I believe this is not in line with the "dependency inversion principle" which, if followed, means you'll have cleaner object oriented code. do you know of any other way to accomplish this without letting the base class know we're dealing with a TPMS surface?
| @@ -76,7 +81,15 @@ class Surface { | |||
| //! \param u The direction of the ray. | |||
| //! \param coincident A hint to the code that the given point should lie | |||
| //! exactly on the surface. | |||
| virtual double distance(Position r, Direction u, bool coincident) const = 0; | |||
There was a problem hiding this comment.
We should probably leave this as =0. That prevents programmer mistakes and ensures that base classes implement this method, which will surely be unique for each surface type.
| @@ -327,6 +340,55 @@ class SurfaceQuadric : public CSGSurface { | |||
| double A_, B_, C_, D_, E_, F_, G_, H_, J_, K_; | |||
| }; | |||
|
|
|||
| //============================================================================== | |||
| // TPMS surfaces | |||
| //============================================================================== | |||
There was a problem hiding this comment.
"TPMS", to me, means "tire pressure monitoring system". There should probably be some documentation explaining what this is.
| std::string surface_type; | ||
|
|
||
| private: | ||
| std::list<std::string> tpms_types = {"Schwarz_P", "Gyroid", "Diamond"}; |
There was a problem hiding this comment.
std::array<std::string, 3> is the better choice here. std::list is a linked list under the hood; you probably don't need to use a linked list here.
| std::vector<double> x_grid; | ||
| std::vector<double> y_grid; | ||
| std::vector<double> z_grid; | ||
| std::vector<std::vector<std::vector<double>>> matrix; |
There was a problem hiding this comment.
you should probably use an xtensor object here instead. vectors of vectors make sense only if the array is jagged, which, considering this is named matrix, I assume it is not.
| double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)}; | ||
| // TPMS : compute the distance in the end to lower computation time | ||
| const auto& surface = model::surfaces[abs(token) - 1]; | ||
| if (surface->is_tpms()) // Add token to tpms_token vector if SurfaceTPMS |
There was a problem hiding this comment.
can you explain why TPMS are treated as a special case in distance calculations? IMO, if there is some generalized logic here like ray tracing using a fast lower bound distance, that could be included as a general feature in the surface class. For example, I have been thinking about implementing a helix surface, where a lower bound on the ray trace distance could also be helpful for determining the intersection.
| } | ||
|
|
||
| // Function to convert a string to a 3D matrix of doubles | ||
| std::vector<std::vector<std::vector<double>>> convertStringTo3DMatrix( |
There was a problem hiding this comment.
this should already be covered in the utilities in the code, to read to an xtensor object. however I may be wrong.
| } | ||
|
|
||
| // Function to convert a string to a matrix of doubles | ||
| std::vector<std::vector<double>> convertStringToMatrix(const std::string& input) |
There was a problem hiding this comment.
check the naming convention for functions.
https://docs.openmc.org/en/latest/devguide/styleguide.html#naming
| double SurfaceTPMS::evaluate(Position r) const | ||
| { | ||
| double value; | ||
| if (surface_type == "Schwarz_P") { |
There was a problem hiding this comment.
I think you should just implement each as its own class, and not switch based on string comparison in a run-time class like this. building the object on each evaluation, too, might be computationally wasteful depending on what's happening in the constructor.
Description
Hello,
This is my first contribution to OpenMC. I did my best to fit the guidelines but it might not be perfect. I am very open to critics so I can be more efficient in the future. Sorry in advance for potential headaches as it is a big feature.
I am publishing a paper soon that will dive in the technical details of this pull request. Here is a reference to provide some context on TPMS:
https://www.sciencedirect.com/science/article/pii/S014919702300330X
This development introduces two new Surfaces objects:
TPMSwhich represent a TPMS surface which can be "Schwarz-P", "Gyroid" or "Diamond" . It has got apitchand anisovaluecoefficient alongside 12 other coefficients who define the translation and rotation of the surface.FunctionTPMSwhich represent a TPMS whosepitchandisovalueare not constant, but interpolated from a 3D matrix.To support the distance calculation, an algorithm of ray-tracing has been developed in
tpms.cpp. Comments explain the philosophy behind the implementation. This has been tested and verificated against a Serpent + CAD files workflow. The very workflow used in the article above.The implementation have been made so it will be easy in the future to implement other kind of TPMS or Function. Helper functions to perform all necessary operations have been implemented. a bisection function has been implemented, directly taken from the Boost library.
Fixes # (issue)
Initial implementation. No fixes.
Checklist