Skip to content

Commit

Permalink
added intersect for magnets
Browse files Browse the repository at this point in the history
  • Loading branch information
shimwell committed Oct 8, 2024
1 parent ca09343 commit 76d6efe
Show file tree
Hide file tree
Showing 9 changed files with 91 additions and 19 deletions.
4 changes: 4 additions & 0 deletions docs/usage_spherical_tokamak.rst
Original file line number Diff line number Diff line change
Expand Up @@ -323,6 +323,7 @@ Spherical tokamak with toroidal field coils
vertical_mid_point = (600, 0),
thickness = 50,
distance = 40,
rotation_angle = 180,
with_inner_leg = True,
azimuthal_placement_angles = [0, 30, 60, 90, 120, 150, 180],
)
Expand All @@ -349,6 +350,7 @@ Spherical tokamak with toroidal field coils
r1=5,
r2=610,
azimuthal_placement_angles = [0, 30, 60, 90, 120, 150, 180],
rotation_angle = 180,
thickness = 50,
distance = 40
)
Expand Down Expand Up @@ -388,6 +390,7 @@ Spherical tokamak with toroidal field coils
thickness = 50,
distance = 40,
with_inner_leg = True,
rotation_angle = 180,
azimuthal_placement_angles = [0, 30, 60, 90, 120, 150, 180],
)
Expand Down Expand Up @@ -415,6 +418,7 @@ Spherical tokamak with toroidal field coils
r1=5,
r2=610,
azimuthal_placement_angles = [120, 150, 180],
rotation_angle = 180,
thickness = 50,
distance = 40
)
Expand Down
2 changes: 2 additions & 0 deletions docs/usage_tokamak.rst
Original file line number Diff line number Diff line change
Expand Up @@ -268,6 +268,7 @@ Tokamak with several customizations
thickness = 50,
distance = 40,
with_inner_leg = True,
rotation_angle = 180,
azimuthal_placement_angles = [0, 30, 60, 90, 120, 150, 180],
)

Expand Down Expand Up @@ -346,6 +347,7 @@ Tokamak with several customizations
thickness = 50,
distance = 40,
with_inner_leg = True,
rotation_angle = 180,
azimuthal_placement_angles = [0, 30, 60, 90, 120, 150, 180],
)
Expand Down
43 changes: 33 additions & 10 deletions examples/spherical_tokamak_from_plasma_with_tf_magnets.py
Original file line number Diff line number Diff line change
@@ -1,27 +1,49 @@
from pathlib import Path

from example_util_functions import transport_particles_on_h5m_geometry

import paramak


rotation_angle=90
tf_style_1 = paramak.toroidal_field_coil_rectangle(
horizontal_start_point = (10, 520),
vertical_mid_point = (600, 0),
thickness = 50,
distance = 40,
with_inner_leg = True,
azimuthal_placement_angles = [0, 30, 60, 90],
azimuthal_placement_angles = [0, 30, 60, 90, 120, 150, 180],
rotation_angle=rotation_angle
)

result1 = paramak.spherical_tokamak_from_plasma(
radial_build=[
(paramak.LayerType.GAP, 70),
(paramak.LayerType.SOLID, 10),
(paramak.LayerType.SOLID, 10),
(paramak.LayerType.GAP, 50),
(paramak.LayerType.PLASMA, 300),
(paramak.LayerType.GAP, 60),
(paramak.LayerType.SOLID, 10),
(paramak.LayerType.SOLID, 60),
(paramak.LayerType.SOLID, 10),
],
elongation=2.5,
rotation_angle=rotation_angle,
triangularity=0.55,
extra_cut_shapes=[tf_style_1]
)

result1.save("spherical_tokamak_from_plasma_with_rect_tf_coils.step")


tf_style_2 = paramak.toroidal_field_coil_princeton_d(
r1=5,
r2=610,
azimuthal_placement_angles = [120, 150, 180],
thickness = 50,
distance = 40
distance = 40,
azimuthal_placement_angles = [0, 30, 60, 90, 120, 150, 180],
rotation_angle=rotation_angle
)

result = paramak.spherical_tokamak_from_plasma(
result2 = paramak.spherical_tokamak_from_plasma(
radial_build=[
(paramak.LayerType.GAP, 70),
(paramak.LayerType.SOLID, 10),
Expand All @@ -34,13 +56,14 @@
(paramak.LayerType.SOLID, 10),
],
elongation=2.5,
rotation_angle=180,
rotation_angle=rotation_angle,
triangularity=0.55,
extra_cut_shapes=[tf_style_1, tf_style_2]
extra_cut_shapes=[tf_style_2]
)

result.save(f"spherical_tokamak_minimal.step")
result2.save("spherical_tokamak_from_plasma_with_prin_tf_coils.step")

# from example_util_functions import transport_particles_on_h5m_geometry
# from cad_to_dagmc import CadToDagmc
# my_model = CadToDagmc()
# material_tags = ["mat1"] * 7
Expand Down
7 changes: 5 additions & 2 deletions examples/tokamak_with_pf_tf_magnets_divertor.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,19 +2,20 @@
from example_util_functions import transport_particles_on_h5m_geometry

import paramak
from cadquery import vis, Workplane
import cadquery as cq

# makes a rectangle that overlaps the lower blanket under the plasma
# the intersection of this and the layers will form the lower divertor
points = [(300, -700), (300, 0), (400, 0), (400, -700)]
divertor_lower = Workplane('XZ', origin=(0,0,0)).polyline(points).close().revolve(180)
divertor_lower = cq.Workplane('XZ', origin=(0,0,0)).polyline(points).close().revolve(180)

# creates a toroidal
tf = paramak.toroidal_field_coil_rectangle(
horizontal_start_point = (10, 520),
vertical_mid_point = (860, 0),
thickness = 50,
distance = 40,
rotation_angle=180,
with_inner_leg = True,
azimuthal_placement_angles = [0, 30, 60, 90, 120, 150, 180],
)
Expand Down Expand Up @@ -78,6 +79,8 @@
)
my_reactor.save(f"tokamak_with_divertor.step")
print(f"Saved as tokamak_with_divertor.step")

# from cadquery import vis
# vis.show(my_reactor)

# from cad_to_dagmc import CadToDagmc
Expand Down
8 changes: 7 additions & 1 deletion src/paramak/workplanes/cutting_wedge.py
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,13 @@ def cutting_wedge(

wire = create_wire_workplane_from_points(points=points, plane=plane, origin=origin, obj=obj)

solid = wire.revolve(rotation_angle)
solid = wire.revolve(angleDegrees=rotation_angle,)
# The code can be changed to revolve it in the other direction
# solid = wire.revolve(
# angleDegrees=rotation_angle,
# axisStart=cq.Vector(0, 0),
# axisEnd=cq.Vector(0, -1)
# )
solid.name = name
solid.color = cq.Color(*color)
return solid
13 changes: 11 additions & 2 deletions src/paramak/workplanes/toroidal_field_coil_princeton_d.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,8 @@
from scipy import integrate
from scipy.optimize import minimize
from typing import List, Tuple
from ..workplanes.cutting_wedge import cutting_wedge


def _compute_inner_points(R1, R2):
"""Computes the inner curve points
Expand Down Expand Up @@ -152,6 +154,7 @@ def toroidal_field_coil_princeton_d(
r2: float = 300,
thickness: float = 30,
distance: float = 20,
rotation_angle: float = 360.0,
name: str = "toroidal_field_coil",
with_inner_leg: bool = True,
azimuthal_placement_angles: typing.Sequence[float] = [0],
Expand All @@ -169,6 +172,7 @@ def toroidal_field_coil_princeton_d(
r2 (float, optional): Outer radius of the coil. Defaults to 300.
thickness (float, optional): Thickness of the coil. Defaults to 30.
distance (float, optional): Distance to extrude the coil. Defaults to 20.
rotation_angle (float): angle of rotation in degrees, this cuts the resulting shape with a wedge. Useful for sector models.
name (str, optional): Name of the coil. Defaults to "toroidal_field_coil".
with_inner_leg (bool, optional): Whether to include the inner leg of the coil. Defaults to True.
azimuthal_placement_angles (typing.Sequence[float], optional): Angles for azimuthal placement. Defaults to [0].
Expand All @@ -192,15 +196,20 @@ def toroidal_field_coil_princeton_d(
inner_leg_connection_points=[(x,z,'straight') for x,z in inner_leg_connection_points]
# need to get square end, it appears to miss the last point in the solid, TODO fix so this append is not needed
inner_leg_connection_points.append(inner_leg_connection_points[-1])
for i in inner_leg_connection_points:
print('inner_leg_connection_points' , i)
inner_wire = create_wire_workplane_from_points(
points=inner_leg_connection_points, plane=plane, origin=origin, obj=obj
)
inner_solid = inner_wire.extrude(until=distance / 2, both=True)
inner_solid = rotate_solid(angles=azimuthal_placement_angles, solid=inner_solid)
solid = solid.union(inner_solid)

if rotation_angle < 360.:
bb=solid.val().BoundingBox()
radius = max(bb.xmax, bb.ymax)*2.1 # larger than the bounding box to ensure clean cut
height = max(bb.zmax, bb.zmin)*2.1 # larger than the bounding box to ensure clean cut
cutting_shape = cutting_wedge(height=height, radius=radius, rotation_angle=rotation_angle)
solid = solid.intersect(cutting_shape)

solid.name = name
solid.color = color
return solid
13 changes: 10 additions & 3 deletions src/paramak/workplanes/toroidal_field_coil_rectangle.py
Original file line number Diff line number Diff line change
@@ -1,13 +1,15 @@
import typing

from ..utils import create_wire_workplane_from_points, rotate_solid
from ..workplanes.cutting_wedge import cutting_wedge


def toroidal_field_coil_rectangle(
horizontal_start_point: typing.Tuple[float, float] = (20, 200),
vertical_mid_point: typing.Tuple[float, float] = (350, 0),
thickness: float = 30,
distance: float = 20,
rotation_angle: float = 360.0,
name: str = "toroidal_field_coil",
with_inner_leg: bool = True,
azimuthal_placement_angles: typing.Sequence[float] = [0],
Expand All @@ -26,9 +28,7 @@ def toroidal_field_coil_rectangle(
vertical section (cm).
thickness: the thickness of the toroidal field coil.
distance: the extrusion distance.
number_of_coils: the number of tf coils. This changes by the
azimuth_placement_angle dividing up 360 degrees by the number of
coils.
rotation_angle (float): angle of rotation in degrees, this cuts the resulting shape with a wedge. Useful for sector models.
with_inner_leg: include the inner tf leg. Defaults to True.
azimuth_start_angle: The azimuth angle to for the first TF coil which
offsets the placement of coils around the azimuthal angle
Expand Down Expand Up @@ -101,6 +101,13 @@ def toroidal_field_coil_rectangle(
inner_solid = rotate_solid(angles=azimuthal_placement_angles, solid=inner_solid)
solid = solid.union(inner_solid)

if rotation_angle < 360.:
bb=solid.val().BoundingBox()
radius = max(bb.xmax, bb.ymax)*1.1 # 10% larger than the bounding box to ensure clean cut
height = max(bb.zmax, bb.zmin)*2.1 # 10% larger than the bounding box to ensure clean cut
cutting_shape = cutting_wedge(height=height, radius=radius, rotation_angle=rotation_angle)
solid = solid.intersect(cutting_shape)

solid.name = name
solid.color = color
return solid
11 changes: 10 additions & 1 deletion tests/test_workplanes/test_toroidal_field_coil_princeton_d.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,4 +5,13 @@ def test_creation_of_inner_leg():
solid_without_inner_leg = paramak.toroidal_field_coil_princeton_d(with_inner_leg=False)
solid_with_inner_leg = paramak.toroidal_field_coil_princeton_d(with_inner_leg=True)

assert solid_without_inner_leg.val().Volume() < solid_with_inner_leg.val().Volume()
assert solid_without_inner_leg.val().Volume() < solid_with_inner_leg.val().Volume()

def test_rotation_angle():
solid_360_uncut = paramak.toroidal_field_coil_princeton_d(azimuthal_placement_angles=[0,180], rotation_angle=360)
solid_180_uncut = paramak.toroidal_field_coil_princeton_d(azimuthal_placement_angles=[0,180], rotation_angle=360)
solid_180_cut = paramak.toroidal_field_coil_princeton_d(azimuthal_placement_angles=[0,180], rotation_angle=180)

assert solid_360_uncut.val().Volume() == solid_180_uncut.val().Volume()
# checks relative volume difference
assert abs(solid_180_cut.val().Volume() - 0.5 * solid_180_uncut.val().Volume()) / (0.5 * solid_180_uncut.val().Volume()) < 0.00001
9 changes: 9 additions & 0 deletions tests/test_workplanes/test_toroidal_field_coil_rectangle.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,3 +22,12 @@ def test_volume_rotation_angle():
assert math.isclose(
half_toroidal_field_coil_rectangle.val().Volume(), toroidal_field_coil_rectangle.val().Volume() / 2
)

def test_rotation_angle():
solid_360_uncut = paramak.toroidal_field_coil_rectangle(azimuthal_placement_angles=[0,180], rotation_angle=360)
solid_180_uncut = paramak.toroidal_field_coil_rectangle(azimuthal_placement_angles=[0,180], rotation_angle=360)
solid_180_cut = paramak.toroidal_field_coil_rectangle(azimuthal_placement_angles=[0,180], rotation_angle=180)

assert solid_360_uncut.val().Volume() == solid_180_uncut.val().Volume()
# checks relative volume difference
assert abs(solid_180_cut.val().Volume() - 0.5 * solid_180_uncut.val().Volume()) / (0.5 * solid_180_uncut.val().Volume()) < 0.00001

0 comments on commit 76d6efe

Please sign in to comment.