The paramak from fusion-energy

Silent failure from Paramak geometry

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.

Extrude rectangle shape requried

In a similar manner to the PF coils (rotated rectangles) it would be useful for a model I'm keen to build if we have an extrude rectangle.

Defined in a similar manner, center_point````, height, widthanddistanceinstead ofrotation_angle```

My use use requires the workplane and rotation_axis to be set

workplane='XY',
rotation_axis ='Z',

Unexpected blanket shape at divertors with large major radius

With a large major radius relative to the minor radius, the blanket can form a strange shape and overlap with the TF coils.
Paramak_issue_large_major_radius.pptx

immutable points

Currently the points are manipulated after input by the user. An end point the same as the first point is added and in some cases connection types are added.

The consequences are that if a variable called points are use to create a shape then the Shape.points will not be equal to the input points.

This limits the reuse of points across components.

A simple test was created that currently fails.

  def test_reuse_points(self):
      """Checks that points can be reused between shapes"""
      points = [
          (100, 0),
          (200, 0),
          (250, 50),
          (200, 100),
      ]
      test_shape = paramak.Shape(
          points=points
      )
      
      assert test_shape.points == [
          (100, 0),
          (200, 0),
          (250, 50),
          (200, 100),
      ]

The input points should not be changed and the a second processed_points can be used internally

unused import in the init.py

This is reducing the flake 8 warnings

perhaps ___all__ = [import import] could be used

as described here

https://stackoverflow.com/questions/31079047/python-pep8-class-in-init-imported-but-not-used

Conda build not installing

While trying to produce a conda package on this branch https://github.com/fusion-energy/paramak/tree/checking_dependencies

The package is not building

How to configure paramak display interface?

I open a paramak example with voila in jupyter lab as shown below:

Where to configure those check boxes and view icons on the top: Axes, Grid, Ortho, etc in paramak?

Method to rotate a shape/reactor around an axis

It could be useful to rotate a shape around an axis after construction

Currently the same result can be done by changing the azimuthal placement of the shape.

Usage

my_shape = Shape(.......)

my_rotated_shape = rotate_shape(my_shape, axis=[(0, 0, 0), (0, 0, 1)], angle=10)

Or maybe it should be a method of Shape?

my_shape = Shape(.......)

my_shape.rotate_shape(axis=[(0, 0, 0), (0, 0, 1)], angle=10)

@shimwell what do you think?

Code of conduct lacks email address

The code of conduct file requires the package maintainers to provide an email contact for enforcement. Should be a very straightforward change.

paramak/CODE_OF_CONDUCT.md

Line 59 in b5a915a

reported by contacting the project team at [INSERT EMAIL ADDRESS]. All

Paramak installation on windows 10

Are there instructions to install pymoab or Moab binary on windows 10? I got pymoab errors when I "import paramak"

Thanks a lot.

Adding precontions to type hints

When checking for types passed into functions it can be useful to ensure only positive ints or positive floats

Here is an example

https://stackoverflow.com/questions/61546201/python-typehint-int-as-positive

Perhaps a new class called PositiveInt is the way to go

Perhaps pydantic can help here https://pydantic-docs.helpmanual.io/usage/types/#constrained-types

add twine check to CI

This additional twine check command could check that the package builds

https://github.com/C-bowman/inference-tools/blob/760061e77bdc8ddeacb951196a3deb783a6ccffd/.github/workflows/tests.yml#L48

It should be in the regular CI so that it is caught before publishing action

Wrong argument type in utils.intersect_solid

paramak/paramak/utils.py

Lines 229 to 246 in bee4f04

    
           def intersect_solid(solid, intersecter): 
        
               """ 
        
               Performs a boolean intersection of a solid with another solid or iterable of 
        
               solids. 
        
               Args: 
        
                   solid Shape: The Shape that you want to intersect 
        
                   intersecter Shape: The Shape(s) that you want to be the intersecting object 
        
               Returns: 
        
                   Shape: The original shape cut with the intersecter shape(s) 
        
               """ 
        
               # Allows for multiple cuts to be applied 
        
               if isinstance(intersecter, Iterable): 
        
                   for intersecting_solid in intersecter: 
        
                       solid = solid.intersect(intersecting_solid.solid) 
        
               else: 
        
                   solid = solid.intersect(intersecter.solid) 
        
               return solid

In the code above, intersecter is a paramak.Shape but solid is paramak.Shape.solid type.

Removing old code

@RemDelaporteMathurin I was just reviewing some of the code and was wondering can we remove this code or is it useful ?

paramak/paramak/parametric_components/blanket_fp.py

Lines 262 to 322 in b91213b

    
               def create_physical_groups(self): 
        
                   """Creates the physical groups for STP files 
        
                   Returns: 
        
                       list: list of dicts containing the physical groups 
        
                   """ 
        
                   groups = [] 
        
                   number_of_volumes = 1  # only one volume 
        
                   number_of_surfaces = 2  # inner and outer 
        
                   surface_names = ["inner", "outer"] 
        
                   volumes_names = ["inside"] 
        
                   # add two cut sections if they exist 
        
                   if self.rotation_angle != 360: 
        
                       number_of_surfaces += 2 
        
                       surface_names += ["left_section", "right_section"] 
        
                       full_rotation = False 
        
                   else: 
        
                       full_rotation = True 
        
                   # add two surfaces between blanket and div if they exist 
        
                   if diff_between_angles(self.start_angle, self.stop_angle) != 0: 
        
                       number_of_surfaces += 2 
        
                       surface_names += ["inner_section", "outer_section"] 
        
                       stop_equals_start = False 
        
                   else: 
        
                       stop_equals_start = True 
        
                   # rearrange order 
        
                   new_order = [i for i in range(len(surface_names))] 
        
                   if full_rotation: 
        
                       if not stop_equals_start: 
        
                           # from ["inner", "outer", "inner_section", "outer_section"] 
        
                           # to ["inner", "inner_section", "outer", "outer_section"] 
        
                           new_order = [0, 2, 1, 3] 
        
                   else: 
        
                       if stop_equals_start: 
        
                           print( 
        
                               "Warning: If start_angle = stop_angle surfaces will not\ 
        
                                be handled correctly" 
        
                           ) 
        
                           new_order = [0, 1, 2, 3] 
        
                       else: 
        
                           # from ['inner', 'outer', 'left_section', 'right_section', 
        
                           #           'inner_section', 'outer_section'] 
        
                           # to ["inner", "inner_section", "outer", "outer_section", 
        
                           #         "left_section", "right_section"] 
        
                           new_order = [0, 4, 1, 5, 2, 3] 
        
                   surface_names = [surface_names[i] for i in new_order] 
        
                   for i in range(1, number_of_volumes + 1): 
        
                       group = {"dim": 3, "id": i, "name": volumes_names[i - 1]} 
        
                       groups.append(group) 
        
                   for i in range(1, number_of_surfaces + 1): 
        
                       group = {"dim": 2, "id": i, "name": surface_names[i - 1]} 
        
                       groups.append(group) 
        
                   self.physical_groups = groups

PoloidalFieldCoilSet ignores azimuth_placement_angle

Because PoloidalFieldCoilSet overwrites create_solid, the step rotate_solid is ignored.

Hence, the azimuth_placement_angle has no influence on this shape.

paramak/paramak/parametric_components/poloidal_field_coil_set.py

Lines 106 to 143 in bee4f04

    
               def create_solid(self): 
        
                   """Creates a 3d solid using points with straight connections 
        
                   edges, azimuth_placement_angle and rotation angle. 
        
                   individual solids in the compound can be accessed using .Solids()[i] 
        
                   where i is an int 
        
                      Returns: 
        
                         A CadQuery solid: A 3D solid volume 
        
                   """ 
        
                   iter_points = iter(self.points) 
        
                   pf_coils_set = [] 
        
                   wires = [] 
        
                   for p1, p2, p3, p4 in zip( 
        
                           iter_points, iter_points, iter_points, iter_points): 
        
                       solid = ( 
        
                           cq.Workplane(self.workplane) 
        
                           .polyline([p1[:2], p2[:2], p3[:2], p4[:2]]) 
        
                       ) 
        
                       wire = solid.close() 
        
                       wires.append(wire) 
        
                       solid = wire.revolve(self.rotation_angle) 
        
                       pf_coils_set.append(solid) 
        
                   compound = cq.Compound.makeCompound( 
        
                       [a.val() for a in pf_coils_set] 
        
                   ) 
        
                   self.wire = wires 
        
                   self.solid = compound 
        
                   return compound

Demos don't run in the CI workflows

@shimwell it seems that the workflows including jupyter-cadquery don't work anymore.

Impacted: #85 #84

Apparently your last commit 695e391 on develop didn't pass the tests either.

I shall try and rebuild everything locally

Use @types decorator to enforce arguments types

I spotted this decorator below. Thanks @EmilyBourne for sharing this repo 😄
https://github.com/EmilyBourne/Pygyro/blob/02031f96f9b05c3c523290fe349eeef9e39138ba/pygyro/advection/accelerated_advection_steps.py#L21.

We could make use of this decorator to enforce the arguments types. This could greatly reduce the amount of code in all the setters.

More doc: https://stackoverflow.com/questions/47298825/how-to-get-argument-type-hint-in-decorator

Make PF coil cases differently

the method of making cases at the moment is not ideal as it involves a connection between the innner and outer corner that results in two edges along side each other.

This causes meshing problems down the line.

Perhaps the cases can be made as a boolean cut

Workplane is not subscriptable error when rotation cadquery.Plane

It is possible to set the workplane to a string such as 'XY' or 'ZY'

It is also possible to set the workplane to an arbitrary plane using a cadquery.Plane

import cadquery as cq
workplane = cq.Plane(origin=(0, 0, 0), xDir=(-1, 1, 0), normal= (1, 1, 1))

However the Paramak code will break if the a workplane is set to a cq.Plane and a rotation is used.

This is due to the get_rotation_axis assuming that all workplanes are strings on this line

paramak/paramak/shape.py

Line 755 in 2530d7e

return rotation_axis[workplane[1]], workplane[1]

BlanketCutterParallels not working correctly

Thanks @mateczentye for spotting this issue. Incorrect BlanketCutterParallels shape being produced from input parameters.

Looks like one point is being missed when creating the shape.

InnerTfCoilsFlat inner_radius arg of 0 causes component construction to fail

Setting inner_radius=0 causes a ZeroDivisionError: float division by zero error

shape= paramak.InnerTfCoilsFlat(
    height=10,
    inner_radius=0,
    outer_radius=20,
    number_of_coils=6,
    gap_size=0,
    inner_radius_type='straight',
)

I think the component might need an if statement adding to handle this special case

preparing to move to cadquery 2.2

Currently the paramak requires cadquery 2.1

There is a really nice importBrep feature in the upcoming CadQuery 2.2 and this feature is essential for the use of Paramak geometry in the neutronics workflow. So Paramak will have to be upgraded to work with CadQuery 2.2 once it is released.

Currently the importBrep feature is available on the CadQuery master branch.

There are some API changes that appear to be in CadQuery 2.2 (currently in master branch)

This issue is aiming to note the changes so that they can be accounted for in the release of the Paramak that makes use of CadQuery 2.2

Proposal to simplify the package

So I've tracked down the cause of failing errors in #72

Basically we are trying to cram too much functionality into one package and finding conflicts between the packages

To get everything working I'm thinking we need to be more modular.

We can't have situations where cubit and moab are imported into the same python environment. It causes hdf5 conflicts and would be hard for users to resolve.

But the workflow works when the export to stp then convert to h5m and then run openmc are all separate files and are executed separately with python ... but it doesn't work as one big file.

So we could make the package more modular

make a new package called stl_to_h5m that uses moab / pymoab to convert a list of stl files and their material tags into a h5m geometry.
then remove export_h5m_with_pymoab from the paramak and remove moab as a dependency. This means the package can be installed on windows which is a nice bonus
as we already have a package cad_to_h5m which takes a list of stp files and their material tags and makes a h5m file so we can remove the export_to_h5m_with_cubit method
We could then remove the general export_to_h5m method and all the method = cubit or method = pymoab parts scattered in the code.

To make it convenient for users we would still have a neutronics_description that provides all the stp filenames and the material_tags in a dictionary form that cad_to_h5m accepts

So the files would look like this

python make_stp_files.py

import paramak

my_shape = paramak.ExtrudeStraightShape(
    points=[
        (400, 100),
        (400, 200),
        (600, 200),
        (600, 100)
        ],
    distance = 180,
)

my_shape.export_stp('steel.stp')
my_shape.export_graveyard('graveyard.stp')

python convert_stp_files_to_h5m.py

from cad_to_h5m import cad_to_h5m

cad_to_h5m(
    files_with_tags=[
        {'filename':'steel.stp', 'material_tag':'mat1'},
        {'filename':'graveyard.stp', 'material_tag':'graveyard'},
    ],
    h5m_filename='my_dagmc.h5m',
    cubit_path='/opt/Coreform-Cubit-2021.5/bin/'
)

python build_openmc_model_from_h5m.py

import openmc
import paramak_neutronics


my_source = openmc.Source()
my_source.space = openmc.stats.Point((0, 0, 0))
my_source.angle = openmc.stats.Isotropic()
my_source.energy = openmc.stats.Discrete([14e6], [1])

neutronics_model = paramak_neutronics.NeutronicsModel(
    h5m_filename='my_dagmc.h5m',
    source=my_source,
    materials={
        "mat1": "eurofer",
    },
    cell_tallies=["flux"],
    simulation_batches=5,
    simulation_particles_per_batch=1e4,
)

neutronics_model.simulate()

Any thoughts on this @RemDelaporteMathurin @billingsley-john

GitHub based CI not uploading to codecov

There is a line in one of the github actions that is broken

This line should upload to codecov

paramak/.github/workflows/ci_with_install.yml

Line 29 in 0fea938

curl -s https://codecov.io/bash

However the actions logs show that it is not working.

This is perhaps because it is running in a container and something needs changing to get it to work correctly

So currently the CircleCi tests are reporting the coverage

Automatic version number in citation file?

Currently the version number is outdated in the citation file.

changing use of trelis to cubit

Looks like we still use trelis in several places.

This should probably be updated to use cubit instead

the terminal command show use of trelis in the code

grep -rnw . -e 'trelis'

Conda build not passing tests

I'm working on a conda package and have now got the package building on the automated_conda_build branch

Unfortunately it is not currently possible to include jupyter-cadquery in the conda package meta.yml. Jupyter-cadquery depends on several packages that are themselves not on conda so it looks like it will be a significant effort to add to conda forge.

As several of the functions in the paramak make use of jupyter-cadquery some of the tests run as part of the conda package build are failing.

There is an open issue on jupyter-cadquery to add the package to conda-forge

In the meantime perhaps we should make jupyter-cadquery an optional dependency so that a conda install can be made?

Allow user to define a vector for the extrusion direction

As far as I can tell, the extruded shapes are all derived from Shape and extruded in a direction normal to the workplane and then rotated according to the definition of the rotation_axis and the azimuth_placement_angle.

However, a user frequently knows the unit vector along which they want the extrusion to occur, so it would be convenient if the various extruded shape classes could automatically calculate the rotation_axis and azumuth_placement_angle from that unit vector and the original workplane normal. Perhaps this is how path_workplane should work (?) but I can't find much documentation on that.

Even better would be a way to define this data in a number of ways:

staring point (x,y,z) and ending point (x,y,z)
unit vector of extrusion direction, reference point (x,y,z) and extents (-l1,l2)
unit vector of extrusion direction, starting point (x,y,z) and length [this is a special case of the previous where the reference point is the starting point and the extents are (0,L)]

Some componets are missing names

for example BlanketConstantThicknessArcH has no .name set

add docker instructions

We don't have very clear instructions for using the docker images.

This is particularly tricky for docker users who are not familar with the -p 8888:8888 command that is needed to share a jupyter lab connection.

For anyone building their own docker image it can be tricky to get CadQuery working in a docker image as libgl and other front ends are missing in most docker images

Perhaps we should add a section to the install instructions to help users installing inside docker and perhaps add another section to the instructions for using the existing docker images

Run example notebooks in parrallel with GA and not CircleCI to save 15+ mins at each CI job

Running the example notebooks doesn't contribute to coverage because paramak isn't imported in the test script.
CircleCI prints:

tests/test_example_components.py::TestExampleComponents::test_jupyter_notebooks_example_parametric_components PASSED [100%]
Coverage.py warning: No data was collected. (no-data-collected)

At the moment, these notebooks are run within CircleCI at the very end of the workflow (in serial).
These tests add 15 mins (without contributing to coverage).

What could be done instead is creating a new "Run Examples" workflow in Github Actions that is run in parrallel to the CircleCI build.
This would reduce the total CI workflow (Circle + GA) by 40% 🚀 and would avoid hitting the job limit of CircleCI.

adding vessel options

it would be great to add some difference vessel caps like these https://www.wermac.org/equipment/equipment_img/pvheads.gif

Updating docker images

Many of the docker images used in CI and mentioned in the docs are from the old dockerhub

I propose we migrate to using GitHub container registry instead which I have been testing out over here https://github.com/Shimwell/example_docker_to_ghcr_action

argument type hinting for filenames

filenames arguments are strings throughtout the code. Perhaps a Path object could also be accepted.

The type hinting would need changing to Union[str, PathLike[str]]

conda package is not installing

Now that the conda package builds and uploads to anaconda.org it can be downloaded with

conda install -c fusion-energy paramak

however this appears to not install (tried on new conda enviroment and latest miniconda docker image)

Add parameter to CenterColumnShieldCylinder

PR #33 really needs the docs diagram updating with the new center_point arg

Additional parameter for CenterColumnShieldCylinder

While making the CenterColumnShieldCylinder I had the need to shift the cylinder vertically on the Z axis.

This is not a supported argument at the moment and the component assumes equal distance above and below the center point of (0,0).

A vertical shift could be added or a center point option could be specified to give the extra control needed

Order of union matters?

Trying to investigate the reason of an overlap in the demo model (between the vacuum vessel and the blanket mainly) I realised that the order of the union argument list seems to matter. It shouldn't...

To reproduce simply run :

import paramak

rotation_angle = 2

blanket = paramak.RotateMixedShape(
    points=[
        (1028.5051619363035, -506.43908961374075, "spline"),
        (1192.275810836831, -281.2815003658492, "spline"),
        (1257.92101205376, -56.43827165248217, "spline"),
        (1244.8344044363996, 216.28835367966997, "spline"),
        (1153.671764932414, 400.5162013324601, "spline"),
        (988.2423861054414, 520.3252882504712, "spline"),
        (807.595124581199, 539.8030669696149, "spline"),
        (670.5280987407093, 475.18311549274034, "spline"),
        (584.8527692840913, 334.5130635036775, "spline"),
        (574.292831670967, 81.767193745908, "spline"),
        (579.3751030153657, -158.89465706473618, "spline"),
        (580.3587684368622, -363.4679185119577, "straight"),
        (686.29567682058, -355.1059282936036, "spline"),
        (657.9911865059497, -182.71061116031922, "straight"),
        (657.4125597874224, -62.3733985443065, "straight"),
        (656.6796326106211, 90.05373743597613, "straight"),
        (671.3960388185003, 294.67729656872143, "spline"),
        (725.9412508183462, 379.089387299264, "spline"),
        (804.2487333923825, 419.4532799322212, "spline"),
        (1036.1333908422225, 356.0153240651591, "spline"),
        (1095.4908483844877, 255.92292987253006, "spline"),
        (1127.4889059190507, 131.6750722070317, "spline"),
        (1128.4147086686944, -60.864467978588436, "spline"),
        (1082.0859960719372, -221.46497785651047, "spline"),
        (945.6554596218276, -418.4558632109989, "straight"),
    ],
    rotation_angle=rotation_angle,
    stp_filename="blanket.stp",
    stl_filename="blanket.stl",
    material_tag="blanket_mat",
    color=(0., 1., 0.498),
)

# SN Divertor
divertor = paramak.RotateMixedShape(
    points=[
        (580.3587684368622, -363.4679185119577, "straight"),
        (678.5035036777448, -367.1647983979668, "straight"),
        (698.2057434436014, -383.14688797319684, "straight"),
        (706.2100797165632, -415.21166249470514, "straight"),
        (698.4950568028651, -443.3154942812033, "straight"),
        (694.7436269110794, -479.42923248738805, "straight"),
        (690.8764716755882, -491.47552817037024, "straight"),
        (726.2209203989679, -495.3735987984752, "straight"),
        (745.6917094774136, -463.22080332730013, "straight"),
        (769.222529364193, -459.1341163784806, "straight"),
        (788.8861940154811, -467.093725112643, "straight"),
        (808.5691462240535, -479.0645742673395, "straight"),
        (863.866572957985, -551.0908599376131, "straight"),
        (891.3995609812453, -563.0365602495477, "straight"),
        (895.1702784303153, -530.9340624638968, "straight"),
        (883.2602118072937, -502.8931027843033, "straight"),
        (886.9923541417953, -462.76812415758457, "straight"),
        (945.6554596218276, -418.4558632109989, "straight"),
        (1028.5051619363035, -506.43908961374075, "straight"),
        (985.8893041167635, -622.9033804444106, "straight"),
        (946.7548503870296, -647.096567181423, "straight"),
        (801.5291878153039, -647.5618207725194, "straight"),
        (675.5814387491815, -575.7618746871029, "straight"),
        (612.3761135287094, -491.72701659798975, "straight"),
        (584.5345245889013, -415.60146955751554, "straight"),
    ],
    rotation_angle=rotation_angle,
    stp_filename="divertor.stp",
    stl_filename="divertor.stl",
    material_tag="divertor_mat",
    color=(1., 0.667, 0.),
)

# Vacuum vessel
inner_points = [
        (574.1267649651465, 101.10991818268735),
        (578.248235352396, 276.2517419573795),
        (576.2157294079991, 183.80759081183805),
        (574.6913499497014, 165.8825455403886),
        (592.9839034492745, 356.1350119397756),
        (618.3902277542372, 410.1718551122461),
        (660.0565996143763, 462.707847631295),
        (711.8855011965004, 503.98791778753997),
        (763.7144027786246, 528.2384308674821),
        (800.8076362638702, 535.0747131490673),
        (957.8187204685404, 529.1614255259608),
        (1009.6476220506645, 511.7304778365887),
        (1061.4765236327887, 484.3593637684569),
        (1113.3054252149127, 445.67421552338743),
        (1160.0530619360443, 396.81085941923095),
        (1196.1300424490914, 343.45306647581924),
        (1221.5363667540544, 290.06656753580216),
        (1239.3207937675284, 234.97736402844578),
        (1251.5158294339103, 177.39960361334886),
        (1259.8999164545482, 119.15248445488135),
        (1262.8978627225338, 91.48061850063345),
        (1257.7636381996099, -69.66392221833257),
        (1248.4670705173148, -117.9739684683608),
        (1236.2720348509329, -162.39122749437956),
        (1219.5038608096575, -215.31064328924913),
        (1198.1625483934886, -267.65634357273507),
        (1171.7399711163273, -320.87421415048993),
        (1140.7442554642726, -372.27816805066277),
        (1102.6347690068283, -425.0006755712611),
        (1058.9358912022924, -476.2931568156555),
        (1022.8589106892452, -523.1641842564957),
        (998.4688393564809, -582.3229466099988),
        (961.8837323573345, -633.8322193152558),
        (788.1044741113889, -641.8422854712626),
        (746.4381022512498, -626.4054742155761),
        (694.6092006691258, -595.170356965071),
        (646.3371844896964, -549.174826941633),
        (610.7683304627485, -496.43258230780384),
        (588.9188915604805, -442.35157799696844),
        (578.1466100551762, -383.4976698929381),
        (577.6991331872189, -343.8423522802541),
        (576.2157294079991, -254.96634344146696),
        (576.2157294079991, -201.97055742152952),
        (576.2157294079991, -148.97477140159208),
        (576.2157294079991, -95.97898538165487),
    ]
vac_vessel_1 = paramak.RotateStraightShape(
    points=inner_points,
    rotation_angle=rotation_angle,
    # avoid overlap between VV and blanket divertor
    union=[divertor, blanket],
    material_tag="vacuum_vessel_inner_mat",
    stp_filename='vacvesselinner.stp',
    stl_filename='vacvesselinner.stl',
)

vac_vessel_2 = paramak.RotateStraightShape(
    points=inner_points,
    rotation_angle=rotation_angle,
    # avoid overlap between VV and blanket divertor
    union=[blanket, divertor],
    material_tag="vacuum_vessel_inner_mat",
    stp_filename='vacvesselinner.stp',
    stl_filename='vacvesselinner.stl',
)

print(vac_vessel_1.volume)
# 20844066.514337115
print(vac_vessel_2.volume)
# 20844063.839064926

On the image we can see there are bits that are different but these two shapes should be exactly the same...

Read The Docs build is failing

Error message on read the docs is excess memory usage

CI trigger refractor using paths instead of paths-ignore

we currently use this

paramak/.github/workflows/ci_with_install.yml

Lines 11 to 17 in bee4f04

    
           paths-ignore: 
        
             - 'docs/**' 
        
             - '.gitignore' 
        
             - '*.md' 
        
             - 'CITATION.cff' 
        
             - 'LICENSE.txt' 
        
             - 'readthedocs.yml'

however this could be replaced with less lines

https://github.com/C-bowman/inference-tools/blob/760061e77bdc8ddeacb951196a3deb783a6ccffd/.github/workflows/tests.yml#L5-L6

Missing 3d models for several componets

We could do with 3d models for some of the missing components.

I don't think we should do every component but could do with updating the /docs/source/paramak.parametric_components.rst file for the following components:

In each case the following code will need adding

.. cadquery::
   :select: cadquery_object
   :gridsize: 0

   import paramak
   my_component = paramak.ComponentNameHere(
         ComponentArgumentsHere=Values,
         ComponentArgumentsHere=Values,
         ComponentArgumentsHere=Values,
   )

   cadquery_object = my_component.solid

Where ComponentNameHere, ComponentArgumentsHere and Values have been replaced with the component specific arguments that can be found in the examples

For an example of what is expected for each component take a look at this example PR

PRs for single components on this checklist are most welcome

Feature request: user controlled rotation center

Currently the rotation is always performed around the center of (0,0,0) it might be useful to allow the user to specify this.

The following part of code would also need changing

paramak/paramak/shape.py

Line 709 in 2530d7e

rotated_solids.append(solid.rotate(*self.get_rotation_axis()[0], angle))

Paramak docker container can't run

When running the latest paramak image, an error is raised:

bash: /home/paramak/docker-cmd.sh: No such file or directory

I suspect this is because of the entry point of the image

To reproduce
Run

docker pull ghcr.io/fusion-energy/paramak:latest
docker run ghcr.io/fusion-energy/paramak:latest

issue exporting stp files

User reported issue

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/anaconda3/envs/paramak_env/lib/python3.8/site-packages/paramak/reactor.py", line 309, in export_stp
    if len(filename) != len(self.shapes_and_components):
TypeError: object of type 'NoneType' has no len()
>>> my_reactor.export_stp()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: export_stp() missing 1 required positional argument: 'filename'

Vacuum_vessel_inner_leg cut not working

When performing a cut on the vacuum vessel inner leg component, it does not occur. this is due to the create solid method being over written in the VVIL class

don't print out warning when using blanketfp in reactor

Wondering if we can remove this warning message when using BlanketFP in a reactor model as it is confusing a few users.

update citation file

currently missing a few contributors

Seg fault when export_stl() is called on BallReactor without pf coils

import paramak
my_reactor = paramak.BallReactor(
   inner_bore_radial_thickness=10,
   inboard_tf_leg_radial_thickness=30,
   center_column_shield_radial_thickness=60,
   divertor_radial_thickness=150,
   inner_plasma_gap_radial_thickness=30,
   plasma_radial_thickness=300,
   outer_plasma_gap_radial_thickness=30,
   plasma_gap_vertical_thickness=50,
   firstwall_radial_thickness=30,
   blanket_radial_thickness=50,
   blanket_rear_wall_radial_thickness=30,
   elongation=2,
   triangularity=0.55,
   number_of_tf_coils=16,
   rotation_angle=90,
   pf_coil_case_thicknesses=[],
   pf_coil_radial_thicknesses=[],
   pf_coil_vertical_thicknesses=[],
   pf_coil_radial_position=[],
   pf_coil_vertical_position=[],
   rear_blanket_to_tf_gap=50,
   outboard_tf_coil_radial_thickness=100,
   outboard_tf_coil_poloidal_thickness=50
)
stl_filenames = my_reactor.export_stl()

results in a seg fault. Here is the trace:

Fatal Python error: Segmentation fault

Thread 0x00007f761db44700 (most recent call first):
  File "/home/.../miniconda3/lib/python3.8/threading.py", line 302 in wait
  File "/home/.../miniconda3/lib/python3.8/threading.py", line 558 in wait
  File "/home/.../miniconda3/lib/python3.8/site-packages/IPython/core/history.py", line 829 in run
  File "/home/.../miniconda3/lib/python3.8/site-packages/IPython/core/history.py", line 58 in needs_sqlite
  File "<decorator-gen-17>", line 2 in run
  File "/home/.../miniconda3/lib/python3.8/threading.py", line 932 in _bootstrap_inner
  File "/home/.../miniconda3/lib/python3.8/threading.py", line 890 in _bootstrap

Current thread 0x00007f761fddb180 (most recent call first):
  File "/home/.../.local/lib/python3.8/site-packages/cadquery/occ_impl/shapes.py", line 390 in exportStl
  File "/home/.../.local/lib/python3.8/site-packages/cadquery/occ_impl/exporters/__init__.py", line 104 in export
  File "/home/.../.local/lib/python3.8/site-packages/paramak/shape.py", line 823 in export_stl
  File "/home/.../.local/lib/python3.8/site-packages/paramak/reactor.py", line 447 in export_stl
  File "<ipython-input-2-b89d63acf246>", line 34 in <module>
  File "/home/.../miniconda3/lib/python3.8/site-packages/IPython/core/interactiveshell.py", line 3457 in run_code
  File "/home/.../miniconda3/lib/python3.8/site-packages/IPython/core/interactiveshell.py", line 3377 in run_ast_nodes
  File "/home/.../miniconda3/lib/python3.8/site-packages/IPython/core/interactiveshell.py", line 3185 in run_cell_async
  File "/home/.../miniconda3/lib/python3.8/site-packages/IPython/core/async_helpers.py", line 78 in _pseudo_sync_runner
  File "/home/.../miniconda3/lib/python3.8/site-packages/IPython/core/interactiveshell.py", line 2960 in _run_cell
  File "/home/.../miniconda3/lib/python3.8/site-packages/IPython/core/interactiveshell.py", line 2914 in run_cell
  File "/home/.../miniconda3/lib/python3.8/site-packages/IPython/terminal/interactiveshell.py", line 557 in interact
  File "/home/.../miniconda3/lib/python3.8/site-packages/IPython/terminal/interactiveshell.py", line 566 in mainloop
  File "/home/.../miniconda3/lib/python3.8/site-packages/IPython/terminal/ipapp.py", line 356 in start
  File "/home/.../.local/lib/python3.8/site-packages/traitlets/config/application.py", line 846 in launch_instance
  File "/home/.../miniconda3/lib/python3.8/site-packages/IPython/__init__.py", line 126 in start_ipython
  File "/home/.../miniconda3/bin/ipython", line 8 in <module>
Segmentation fault (core dumped)

Full rotation shapes cause problems with cubit and openmc

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.

DEMO: VV is not cut by the blanket

There is an overlap between the VV and the blanket in the DEMO model that needs fixing cause it makes the neutronics simulations fail.

	def intersect_solid(solid, intersecter):
	"""
	Performs a boolean intersection of a solid with another solid or iterable of
	solids.
	Args:
	solid Shape: The Shape that you want to intersect
	intersecter Shape: The Shape(s) that you want to be the intersecting object
	Returns:
	Shape: The original shape cut with the intersecter shape(s)
	"""

	# Allows for multiple cuts to be applied
	if isinstance(intersecter, Iterable):
	for intersecting_solid in intersecter:
	solid = solid.intersect(intersecting_solid.solid)
	else:
	solid = solid.intersect(intersecter.solid)
	return solid

	def create_physical_groups(self):
	"""Creates the physical groups for STP files

	Returns:
	list: list of dicts containing the physical groups
	"""

	groups = []
	number_of_volumes = 1 # only one volume
	number_of_surfaces = 2 # inner and outer

	surface_names = ["inner", "outer"]
	volumes_names = ["inside"]

	# add two cut sections if they exist
	if self.rotation_angle != 360:
	number_of_surfaces += 2
	surface_names += ["left_section", "right_section"]
	full_rotation = False
	else:
	full_rotation = True

	# add two surfaces between blanket and div if they exist
	if diff_between_angles(self.start_angle, self.stop_angle) != 0:
	number_of_surfaces += 2
	surface_names += ["inner_section", "outer_section"]
	stop_equals_start = False
	else:
	stop_equals_start = True

	# rearrange order
	new_order = [i for i in range(len(surface_names))]
	if full_rotation:
	if not stop_equals_start:
	# from ["inner", "outer", "inner_section", "outer_section"]

	# to ["inner", "inner_section", "outer", "outer_section"]
	new_order = [0, 2, 1, 3]
	else:
	if stop_equals_start:
	print(
	"Warning: If start_angle = stop_angle surfaces will not\
	be handled correctly"
	)
	new_order = [0, 1, 2, 3]
	else:
	# from ['inner', 'outer', 'left_section', 'right_section',
	# 'inner_section', 'outer_section']

	# to ["inner", "inner_section", "outer", "outer_section",
	# "left_section", "right_section"]
	new_order = [0, 4, 1, 5, 2, 3]
	surface_names = [surface_names[i] for i in new_order]

	for i in range(1, number_of_volumes + 1):
	group = {"dim": 3, "id": i, "name": volumes_names[i - 1]}
	groups.append(group)
	for i in range(1, number_of_surfaces + 1):
	group = {"dim": 2, "id": i, "name": surface_names[i - 1]}
	groups.append(group)
	self.physical_groups = groups

	def create_solid(self):
	"""Creates a 3d solid using points with straight connections
	edges, azimuth_placement_angle and rotation angle.
	individual solids in the compound can be accessed using .Solids()[i]
	where i is an int

	Returns:
	A CadQuery solid: A 3D solid volume
	"""

	iter_points = iter(self.points)
	pf_coils_set = []
	wires = []
	for p1, p2, p3, p4 in zip(
	iter_points, iter_points, iter_points, iter_points):

	solid = (
	cq.Workplane(self.workplane)
	.polyline([p1[:2], p2[:2], p3[:2], p4[:2]])
	)

	wire = solid.close()

	wires.append(wire)

	solid = wire.revolve(self.rotation_angle)

	pf_coils_set.append(solid)

	compound = cq.Compound.makeCompound(
	[a.val() for a in pf_coils_set]
	)

	self.wire = wires

	self.solid = compound

	return compound

	paths-ignore:
	- 'docs/**'
	- '.gitignore'
	- '*.md'
	- 'CITATION.cff'
	- 'LICENSE.txt'
	- 'readthedocs.yml'

fusion-energy / paramak Goto Github PK

paramak's Introduction

Paramak

paramak's People

Contributors

Stargazers

Watchers

Forkers

paramak's Issues

Recommend Projects

Recommend Topics

Recommend Org