Added support for random sphere packing in a conical container.

[ZoeRichter]

Description

Currently, model.pack_spheres() only supports rectangular prisms, cylinders, and spheres as container. This PR adds a _Cone subclass to _Container that allows the user to randomly pack truncated conical containers using model.pack_spheres(). This is useful for modeling pebble-bed reactors, which often use a cone in the outlet chute.

It uses the _Cylinder container as a guide because it is the closest to the cone, and I stuck as close as I could to the approach/style there as I could. The biggest deviation from the other _Container subclasses is the way it stores its limits.

The radial limit is a linear function dependent upon the z-coordinate of the pebble. It is the equation of a line that is parallel to the side of the cone, but a perpendicular distance further into the cone equal to the radius of the spheres being packed.

Two parallel lines of the form

$$\begin{aligned} ax + by + c_1 &= 0 \\\ ax + by + c_2 &= 0 \end{aligned}$$

Have a perpendicular distance, d, of

$$\begin{aligned} d &= \frac{|c_2 - c_1|}{\sqrt{a^2 + b^2}} \\\ \end{aligned}$$

With the above naming convention, the line for the edge of the cone is

$$\begin{aligned} z_{sph} &= (r_{cone} - r_{minor})\frac{2z}{r_{major} - r_{minor}} + (z0 - z) \\\ 0 &= \frac{2z}{r_{major} - r_{minor}}r_{cone} - z_{sph} - \frac{2z}{r_{major} - r_{minor}}r_{minor} + (z0 - z) \\\ \end{aligned}$$

With $d = r_{sph}$ and $m = \frac{2z}{r_{major}-r_{minor}}$, the constant for the $r_{lim}$ line can be back-calculated:

$$\begin{aligned} r_{sph} &= \frac{|c_{lim} - c_1|}{\sqrt{a^2 + b^2}} \\\ c_{lim} &= r_{sph}\sqrt{a^2 + b^2} \pm c_1 \\\ c_{lim} &= r_{sph}*\sqrt{m^2 + 1} \pm ((z0-z) - r_{minor}) \end{aligned}$$

Plugging the constant back in, solving for $r_{lim}(z_{sph})$, and simplifying, we get to the equation of the radial limit in the cone:

$$\begin{aligned} r_{lim}(z_{sph}) &= \frac{1}{m}z_{sph} - \frac{r_{sph}}{m}\sqrt{m^2 + 1} + r_{minor} - \frac{z0 - z}{m} \\\ \end{aligned}$$

Because the limit depends on the specific sphere you are checking, but I wanted to use a similar data structure to the other containers, the radial limit is broken into a 2 element list, with the first being the coefficient of $z_{sph}$, and the second being the constant in the $r_{lim}$ equation. As far as I have been able to tell, the only place this caused an issue was in updating the domain limits (line 1596 of triso.py). I got the data structure I used to work with a try/except, but there's probably a more elegant/foolproof way to do that.

I added a new section to the theory section of the docs explaining the above, and I included the discussion of the Jodrey-Tory algorithm, which is currently part of the pack_spheres() docstring, in this section as well. I have not removed the theory from the docstring, but I can if that is desired.

I also did a little test plot so I could check the sphere placement visually, in addition to the unit tests. This uses a packing fraction of 0.4, so it should be in the range that automatically uses the JT algorithm.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 15) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[ZoeRichter]

I think it is mathematically possible to use a similar approach with equations for lines and the perpendicular distance between them to work with more complex surfaces. I'm not sure of the general interest in expanding the available containers further, though - and I haven't thought as deeply about the implementation.

[paulromano]

@ZoeRichter Thanks for proposing this addition. I've also had the thought of extending pack_spheres to support a generic bounded Region. I'm going to see if I can get something working and if so, it might be better to adopt that for use cases like this unless the performance is significantly worse then a specialized version. I'll let you know how it goes!

[ZoeRichter]

Thanks Paul! A single method that can apply to a generic surface would certainly be more elegant than creating subclasses for each container.

I'd also be happy to help with developing/testing such a feature.

[MendesMichael]

Hello @ZoeRichter and @paulromano. I’m interested in packing spheres in an annular cylinder. It looks like model.pack_spheres() supports solid cylinders (and spherical shells), but not a cylindrical shell with inner/outer radii. Do you foresee supporting an annular cylinder? Would that belong here, or as a follow-on PR/issue—especially in light of the generic bounded Region idea?

[ZoeRichter]

I think there would need to be a discussion around the method going forward (individual subclasses for each shape vs a generic container that can work for an arbitrary (but realistic) region), but, I'm definitely open to collaborating to improve this feature, either within this PR, or in a new one, if that is preferred.

For what it's worth, I think the change to allow a cylindrical shell would actually be fairly straight forward - the cylinder would have a minimum r in its limits, and the checks that count the number of regions would need to allow 1 or 2 cylinders.

[ZoeRichter]

Also - I am a bit confused on why the checks are failing after changing a documentation page, but it's not failing in a documentation check. It passed all checks before I made the documentation change - is there possibly some dependency issue?

[paulromano]

Just submitted a PR for packing on generic regions in #3759! Let me know what you think