From f1f0aaeb6b2faa3d0a12a95ee45e4cb62900b13a Mon Sep 17 00:00:00 2001 From: "Sam Ha (Gauss)" Date: Mon, 14 Oct 2024 12:23:59 +0100 Subject: [PATCH 1/2] Updating tokamak_from_plasma to base vertical_build off of inboard build rather than outboard build. When the outboard build is larger than the inboard build, the two options are: 1. Thicknesses to increase on inboard vertical runs (upper-left quadrant in +r being right, +z being up), which causes clashes with central solenoid/VV/TF coil items placed inboard of the IB shield/blanket 2. Thicknesses to increase on outboard vertical runs (upper right quadrant), which pushes out outer blanket. This pushes out the outboard build, which doesn't lead to clashes and is closer to design goals for tokamaks. The other point is that increased thickness should match where neutron doses tend to be higher, which is more pronounced on the outboard side away from the upper/lower (i.e. divertor) regions --- src/paramak/assemblies/tokamak.py | 528 +++++++++++++++--------------- 1 file changed, 266 insertions(+), 262 deletions(-) diff --git a/src/paramak/assemblies/tokamak.py b/src/paramak/assemblies/tokamak.py index 505ef437..acdbc954 100644 --- a/src/paramak/assemblies/tokamak.py +++ b/src/paramak/assemblies/tokamak.py @@ -10,291 +10,295 @@ def count_cylinder_layers(radial_build): - before_plasma = 0 - after_plasma = 0 - found_plasma = False + before_plasma = 0 + after_plasma = 0 + found_plasma = False - for item in radial_build: - if item[0] == LayerType.PLASMA: - found_plasma = True - elif item[0] == LayerType.SOLID: - if not found_plasma: - before_plasma += 1 - else: - after_plasma += 1 + for item in radial_build: + if item[0] == LayerType.PLASMA: + found_plasma = True + elif item[0] == LayerType.SOLID: + if not found_plasma: + before_plasma += 1 + else: + after_plasma += 1 - return before_plasma - after_plasma + return before_plasma - after_plasma def create_center_column_shield_cylinders(radial_build, rotation_angle, center_column_shield_height): - cylinders = [] - total_sum = 0 - layer_count = 0 + cylinders = [] + total_sum = 0 + layer_count = 0 - number_of_cylinder_layers = count_cylinder_layers(radial_build) + number_of_cylinder_layers = count_cylinder_layers(radial_build) - for index, item in enumerate(radial_build): - if item[0] == LayerType.PLASMA: - break + for index, item in enumerate(radial_build): + if item[0] == LayerType.PLASMA: + break - if item[0] == LayerType.GAP: - total_sum += item[1] - continue + if item[0] == LayerType.GAP: + total_sum += item[1] + continue - thickness = item[1] - layer_count += 1 + thickness = item[1] + layer_count += 1 - if layer_count > number_of_cylinder_layers: - break + if layer_count > number_of_cylinder_layers: + break - cylinder = center_column_shield_cylinder( - inner_radius=total_sum, - thickness=item[1], - name=f"layer_{layer_count}", - rotation_angle=rotation_angle, - height=center_column_shield_height, - ) - total_sum += thickness - cylinders.append(cylinder) - return cylinders + cylinder = center_column_shield_cylinder( + inner_radius=total_sum, + thickness=item[1], + name=f"layer_{layer_count}", + rotation_angle=rotation_angle, + height=center_column_shield_height, + ) + total_sum += thickness + cylinders.append(cylinder) + return cylinders def distance_to_plasma(radial_build, index): - distance = 0 - for item in radial_build[index + 1 :]: - if item[0] == LayerType.PLASMA: - break - distance += item[1] - return distance + distance = 0 + for item in radial_build[index + 1 :]: + if item[0] == LayerType.PLASMA: + break + distance += item[1] + return distance def create_layers_from_plasma( - radial_build, vertical_build, minor_radius, major_radius, triangularity, elongation, rotation_angle, center_column + radial_build, vertical_build, minor_radius, major_radius, triangularity, elongation, rotation_angle, center_column ): - plasma_index_rb = get_plasma_index(radial_build) - plasma_index_vb = get_plasma_index(vertical_build) - indexes_from_plamsa_to_end = len(radial_build) - plasma_index_rb - layers = [] - - cumulative_thickness_orb = 0 - cumulative_thickness_irb = 0 - cumulative_thickness_uvb = 0 - cumulative_thickness_lvb = 0 - for index_delta in range(indexes_from_plamsa_to_end): - - if radial_build[plasma_index_rb + index_delta][0] == LayerType.PLASMA: - continue - outer_layer_thickness = radial_build[plasma_index_rb + index_delta][1] - inner_layer_thickness = radial_build[plasma_index_rb - index_delta][1] - upper_layer_thickness = vertical_build[plasma_index_vb - index_delta][1] - lower_layer_thickness = vertical_build[plasma_index_vb + index_delta][1] - - if radial_build[plasma_index_rb + index_delta][0] == LayerType.GAP: - cumulative_thickness_orb += outer_layer_thickness - cumulative_thickness_irb += inner_layer_thickness - cumulative_thickness_uvb += upper_layer_thickness - cumulative_thickness_lvb += lower_layer_thickness - continue - - # build outer layer - if radial_build[plasma_index_rb + index_delta][0] == LayerType.SOLID: - outer_layer = blanket_from_plasma( - minor_radius=minor_radius, - major_radius=major_radius, - triangularity=triangularity, - elongation=elongation, - thickness=[ - upper_layer_thickness, - outer_layer_thickness, - lower_layer_thickness - ], - offset_from_plasma=[ - cumulative_thickness_uvb, - cumulative_thickness_orb, - cumulative_thickness_lvb - ], - start_angle=90, - stop_angle=-90, - rotation_angle=rotation_angle, - color=(0.5, 0.5, 0.5), - name=f"layer_{index_delta}", - allow_overlapping_shape=True, - ) - inner_layer = blanket_from_plasma( - minor_radius=minor_radius, - major_radius=major_radius, - triangularity=triangularity, - elongation=elongation, - thickness=[ - lower_layer_thickness, - inner_layer_thickness, - upper_layer_thickness, - ], - offset_from_plasma=[ - cumulative_thickness_lvb, - cumulative_thickness_irb, - cumulative_thickness_uvb, - ], - start_angle=-90, - stop_angle=-270, - rotation_angle=rotation_angle, - color=(0.5, 0.5, 0.5), - name=f"layer_{index_delta}", - allow_overlapping_shape=True, - ) - layer = outer_layer.union(inner_layer) - layers.append(layer) - # layers.append(inner_layer) - cumulative_thickness_orb += outer_layer_thickness - cumulative_thickness_irb += inner_layer_thickness - cumulative_thickness_uvb += upper_layer_thickness - cumulative_thickness_lvb += lower_layer_thickness - # build inner layer - - # union layers - - return layers + plasma_index_rb = get_plasma_index(radial_build) + plasma_index_vb = get_plasma_index(vertical_build) + indexes_from_plamsa_to_end = len(radial_build) - plasma_index_rb + layers = [] + + cumulative_thickness_orb = 0 + cumulative_thickness_irb = 0 + cumulative_thickness_uvb = 0 + cumulative_thickness_lvb = 0 + for index_delta in range(indexes_from_plamsa_to_end): + + if radial_build[plasma_index_rb + index_delta][0] == LayerType.PLASMA: + continue + outer_layer_thickness = radial_build[plasma_index_rb + index_delta][1] + inner_layer_thickness = radial_build[plasma_index_rb - index_delta][1] + upper_layer_thickness = vertical_build[plasma_index_vb - index_delta][1] + lower_layer_thickness = vertical_build[plasma_index_vb + index_delta][1] + + if radial_build[plasma_index_rb + index_delta][0] == LayerType.GAP: + cumulative_thickness_orb += outer_layer_thickness + cumulative_thickness_irb += inner_layer_thickness + cumulative_thickness_uvb += upper_layer_thickness + cumulative_thickness_lvb += lower_layer_thickness + continue + + # build outer layer + if radial_build[plasma_index_rb + index_delta][0] == LayerType.SOLID: + outer_layer = blanket_from_plasma( + minor_radius=minor_radius, + major_radius=major_radius, + triangularity=triangularity, + elongation=elongation, + thickness=[ + upper_layer_thickness, + outer_layer_thickness, + lower_layer_thickness + ], + offset_from_plasma=[ + cumulative_thickness_uvb, + cumulative_thickness_orb, + cumulative_thickness_lvb + ], + start_angle=90, + stop_angle=-90, + rotation_angle=rotation_angle, + color=(0.5, 0.5, 0.5), + name=f"layer_{index_delta}", + allow_overlapping_shape=True, + ) + inner_layer = blanket_from_plasma( + minor_radius=minor_radius, + major_radius=major_radius, + triangularity=triangularity, + elongation=elongation, + thickness=[ + lower_layer_thickness, + inner_layer_thickness, + upper_layer_thickness, + ], + offset_from_plasma=[ + cumulative_thickness_lvb, + cumulative_thickness_irb, + cumulative_thickness_uvb, + ], + start_angle=-90, + stop_angle=-270, + rotation_angle=rotation_angle, + color=(0.5, 0.5, 0.5), + name=f"layer_{index_delta}", + allow_overlapping_shape=True, + ) + layer = outer_layer.union(inner_layer) + layers.append(layer) + # layers.append(inner_layer) + cumulative_thickness_orb += outer_layer_thickness + cumulative_thickness_irb += inner_layer_thickness + cumulative_thickness_uvb += upper_layer_thickness + cumulative_thickness_lvb += lower_layer_thickness + # build inner layer + + # union layers + + return layers def tokamak_from_plasma( - radial_build: Sequence[Tuple[LayerType, float]], - elongation: float = 2.0, - triangularity: float = 0.55, - rotation_angle: float = 180.0, - extra_cut_shapes: Sequence[cq.Workplane] = [], - extra_intersect_shapes: Sequence[cq.Workplane] = [], + radial_build: Sequence[Tuple[LayerType, float]], + elongation: float = 2.0, + triangularity: float = 0.55, + rotation_angle: float = 180.0, + extra_cut_shapes: Sequence[cq.Workplane] = [], + extra_intersect_shapes: Sequence[cq.Workplane] = [], ): - inner_equatorial_point = sum_up_to_plasma(radial_build) - plasma_radial_thickness = get_plasma_value(radial_build) - outer_equatorial_point = inner_equatorial_point + plasma_radial_thickness - - # sets major radius and minor radius from equatorial_points to allow a - # radial build. This helps avoid the plasma overlapping the center - # column and other components - major_radius = (outer_equatorial_point + inner_equatorial_point) / 2 - minor_radius = major_radius - inner_equatorial_point - - # make vertical build from outer radial build - pi = get_plasma_index(radial_build) - upper_vertical_build = radial_build[pi:] - - plasma_height = 2 * minor_radius * elongation - # slice opperation reverses the list and removes the last value to avoid two plasmas - vertical_build = upper_vertical_build[::-1][:-1] + [(LayerType.PLASMA, plasma_height)] + upper_vertical_build[1:] - - return tokamak( - radial_build=radial_build, - vertical_build=vertical_build, - triangularity=triangularity, - rotation_angle=rotation_angle, - extra_cut_shapes=extra_cut_shapes, - extra_intersect_shapes=extra_intersect_shapes, - ) + inner_equatorial_point = sum_up_to_plasma(radial_build) + plasma_radial_thickness = get_plasma_value(radial_build) + outer_equatorial_point = inner_equatorial_point + plasma_radial_thickness + + # sets major radius and minor radius from equatorial_points to allow a + # radial build. This helps avoid the plasma overlapping the center + # column and other components + major_radius = (outer_equatorial_point + inner_equatorial_point) / 2 + minor_radius = major_radius - inner_equatorial_point + + # make vertical build from outer radial build + # make vertical build from outer radial build + pi = get_plasma_index(radial_build) + #upper_vertical_build = radial_build[pi:] + #changing to inner vertical build + rbi = len(radial_build)-1 - pi + upper_vertical_build = radial_build[pi-rbi:pi][::-1] + + plasma_height = 2 * minor_radius * elongation + # slice opperation reverses the list and removes the last value to avoid two plasmas + vertical_build = upper_vertical_build[::-1] + [(LayerType.PLASMA, plasma_height)] + upper_vertical_build + + return tokamak( + radial_build=radial_build, + vertical_build=vertical_build, + triangularity=triangularity, + rotation_angle=rotation_angle, + extra_cut_shapes=extra_cut_shapes, + extra_intersect_shapes=extra_intersect_shapes, + ) def tokamak( - radial_build: Union[Sequence[Sequence[Tuple[str, float]]], Sequence[Tuple[str, float]]], - vertical_build: Sequence[Tuple[str, float]], - triangularity: float = 0.55, - rotation_angle: float = 180.0, - extra_cut_shapes: Sequence[cq.Workplane] = [], - extra_intersect_shapes: Sequence[cq.Workplane] = [], + radial_build: Union[Sequence[Sequence[Tuple[str, float]]], Sequence[Tuple[str, float]]], + vertical_build: Sequence[Tuple[str, float]], + triangularity: float = 0.55, + rotation_angle: float = 180.0, + extra_cut_shapes: Sequence[cq.Workplane] = [], + extra_intersect_shapes: Sequence[cq.Workplane] = [], ): - """ - Creates a tokamak fusion reactor from a radial build and plasma parameters. - - Args: - radial_build: A list of tuples containing the radial build of the reactor. - elongation: The elongation of the plasma. Defaults to 2.0. - triangularity: The triangularity of the plasma. Defaults to 0.55. - rotation_angle: The rotation angle of the plasma. Defaults to 180.0. - extra_cut_shapes: A list of extra shapes to cut the reactor with. Defaults to []. - extra_intersect_shapes: A list of extra shapes to intersect the reactor with. Defaults to []. - - Returns: - CadQuery.Assembly: A CadQuery Assembly object representing the tokamak fusion reactor. - """ - - inner_equatorial_point = sum_up_to_plasma(radial_build) - plasma_radial_thickness = get_plasma_value(radial_build) - plasma_vertical_thickness = get_plasma_value(vertical_build) - outer_equatorial_point = inner_equatorial_point + plasma_radial_thickness - - major_radius = (outer_equatorial_point + inner_equatorial_point) / 2 - minor_radius = major_radius - inner_equatorial_point - - elongation = (plasma_vertical_thickness / 2) / minor_radius - blanket_rear_wall_end_height = sum([item[1] for item in vertical_build]) - - plasma = plasma_simplified( - major_radius=major_radius, - minor_radius=minor_radius, - elongation=elongation, - triangularity=triangularity, - rotation_angle=rotation_angle, - ) - - inner_radial_build = create_center_column_shield_cylinders( - radial_build, rotation_angle, blanket_rear_wall_end_height - ) - - blanket_layers = create_layers_from_plasma( - radial_build=radial_build, - vertical_build=vertical_build, - minor_radius=minor_radius, - major_radius=major_radius, - triangularity=triangularity, - elongation=elongation, - rotation_angle=rotation_angle, - center_column=inner_radial_build[0], # blanket_cutting_cylinder, - ) - - my_assembly = cq.Assembly() - - for i, entry in enumerate(extra_cut_shapes): - if isinstance(entry, cq.Workplane): - my_assembly.add(entry, name=f"add_extra_cut_shape_{i+1}") - else: - raise ValueError(f"extra_cut_shapes should only contain cadquery Workplanes, not {type(entry)}") - - # builds up the intersect shapes - intersect_shapes_to_cut = [] - if len(extra_intersect_shapes)>0: - all_shapes = [] - for shape in inner_radial_build + blanket_layers: - all_shapes.append(shape) - - # makes a union of the the radial build to use as a base for the intersect shapes - reactor_compound=inner_radial_build[0] - for i, entry in enumerate(inner_radial_build[1:] + blanket_layers): - reactor_compound = reactor_compound.union(entry) - - # adds the extra intersect shapes to the assembly - for i, entry in enumerate(extra_intersect_shapes): - reactor_entry_intersection = entry.intersect(reactor_compound) - intersect_shapes_to_cut.append(reactor_entry_intersection) - my_assembly.add(reactor_entry_intersection, name=f"extra_intersect_shapes_{i+1}") - - # builds just the core if there are no extra parts - if len(extra_cut_shapes) == 0 and len(intersect_shapes_to_cut) == 0: - for i, entry in enumerate(inner_radial_build): - my_assembly.add(entry, name=f"inboard_layer_{i+1})") - for i, entry in enumerate(blanket_layers): - my_assembly.add(entry, name=f"outboard_layer_{i+1})") - else: - shapes_and_components = [] - for i, entry in enumerate(inner_radial_build + blanket_layers): - for cutter in extra_cut_shapes + extra_intersect_shapes: - entry = entry.cut(cutter) - # TODO use something like this to return a list of material tags for the solids in order, as some solids get split into multiple - # for subentry in entry.objects: - # print(i, subentry) - shapes_and_components.append(entry) - - for i, entry in enumerate(shapes_and_components): - my_assembly.add(entry, name=f"layer_{i+1})") #TODO track the names of shapes, even when extra shapes are made due to splitting - - - my_assembly.add(plasma, name="plasma") - - return my_assembly + """ + Creates a tokamak fusion reactor from a radial build and plasma parameters. + + Args: + radial_build: A list of tuples containing the radial build of the reactor. + elongation: The elongation of the plasma. Defaults to 2.0. + triangularity: The triangularity of the plasma. Defaults to 0.55. + rotation_angle: The rotation angle of the plasma. Defaults to 180.0. + extra_cut_shapes: A list of extra shapes to cut the reactor with. Defaults to []. + extra_intersect_shapes: A list of extra shapes to intersect the reactor with. Defaults to []. + + Returns: + CadQuery.Assembly: A CadQuery Assembly object representing the tokamak fusion reactor. + """ + + inner_equatorial_point = sum_up_to_plasma(radial_build) + plasma_radial_thickness = get_plasma_value(radial_build) + plasma_vertical_thickness = get_plasma_value(vertical_build) + outer_equatorial_point = inner_equatorial_point + plasma_radial_thickness + + major_radius = (outer_equatorial_point + inner_equatorial_point) / 2 + minor_radius = major_radius - inner_equatorial_point + + elongation = (plasma_vertical_thickness / 2) / minor_radius + blanket_rear_wall_end_height = sum([item[1] for item in vertical_build]) + + plasma = plasma_simplified( + major_radius=major_radius, + minor_radius=minor_radius, + elongation=elongation, + triangularity=triangularity, + rotation_angle=rotation_angle, + ) + + inner_radial_build = create_center_column_shield_cylinders( + radial_build, rotation_angle, blanket_rear_wall_end_height + ) + + blanket_layers = create_layers_from_plasma( + radial_build=radial_build, + vertical_build=vertical_build, + minor_radius=minor_radius, + major_radius=major_radius, + triangularity=triangularity, + elongation=elongation, + rotation_angle=rotation_angle, + center_column=inner_radial_build[0], # blanket_cutting_cylinder, + ) + + my_assembly = cq.Assembly() + + for i, entry in enumerate(extra_cut_shapes): + if isinstance(entry, cq.Workplane): + my_assembly.add(entry, name=f"add_extra_cut_shape_{i+1}") + else: + raise ValueError(f"extra_cut_shapes should only contain cadquery Workplanes, not {type(entry)}") + + # builds up the intersect shapes + intersect_shapes_to_cut = [] + if len(extra_intersect_shapes)>0: + all_shapes = [] + for shape in inner_radial_build + blanket_layers: + all_shapes.append(shape) + + # makes a union of the the radial build to use as a base for the intersect shapes + reactor_compound=inner_radial_build[0] + for i, entry in enumerate(inner_radial_build[1:] + blanket_layers): + reactor_compound = reactor_compound.union(entry) + + # adds the extra intersect shapes to the assembly + for i, entry in enumerate(extra_intersect_shapes): + reactor_entry_intersection = entry.intersect(reactor_compound) + intersect_shapes_to_cut.append(reactor_entry_intersection) + my_assembly.add(reactor_entry_intersection, name=f"extra_intersect_shapes_{i+1}") + + # builds just the core if there are no extra parts + if len(extra_cut_shapes) == 0 and len(intersect_shapes_to_cut) == 0: + for i, entry in enumerate(inner_radial_build): + my_assembly.add(entry, name=f"inboard_layer_{i+1})") + for i, entry in enumerate(blanket_layers): + my_assembly.add(entry, name=f"outboard_layer_{i+1})") + else: + shapes_and_components = [] + for i, entry in enumerate(inner_radial_build + blanket_layers): + for cutter in extra_cut_shapes + extra_intersect_shapes: + entry = entry.cut(cutter) + # TODO use something like this to return a list of material tags for the solids in order, as some solids get split into multiple + # for subentry in entry.objects: + # print(i, subentry) + shapes_and_components.append(entry) + + for i, entry in enumerate(shapes_and_components): + my_assembly.add(entry, name=f"layer_{i+1})") #TODO track the names of shapes, even when extra shapes are made due to splitting + + + my_assembly.add(plasma, name="plasma") + + return my_assembly From a7e20bb22a341bff54ae79125a52e94cb5447c36 Mon Sep 17 00:00:00 2001 From: "Sam Ha (Gauss)" Date: Mon, 14 Oct 2024 13:41:40 +0100 Subject: [PATCH 2/2] Tidied up the indentation and entries. This implements the vertical build based on an inner radial build --- src/paramak/assemblies/tokamak.py | 529 +++++++++++++++--------------- 1 file changed, 263 insertions(+), 266 deletions(-) diff --git a/src/paramak/assemblies/tokamak.py b/src/paramak/assemblies/tokamak.py index acdbc954..d2b83970 100644 --- a/src/paramak/assemblies/tokamak.py +++ b/src/paramak/assemblies/tokamak.py @@ -10,295 +10,292 @@ def count_cylinder_layers(radial_build): - before_plasma = 0 - after_plasma = 0 - found_plasma = False + before_plasma = 0 + after_plasma = 0 + found_plasma = False - for item in radial_build: - if item[0] == LayerType.PLASMA: - found_plasma = True - elif item[0] == LayerType.SOLID: - if not found_plasma: - before_plasma += 1 - else: - after_plasma += 1 + for item in radial_build: + if item[0] == LayerType.PLASMA: + found_plasma = True + elif item[0] == LayerType.SOLID: + if not found_plasma: + before_plasma += 1 + else: + after_plasma += 1 - return before_plasma - after_plasma + return before_plasma - after_plasma def create_center_column_shield_cylinders(radial_build, rotation_angle, center_column_shield_height): - cylinders = [] - total_sum = 0 - layer_count = 0 + cylinders = [] + total_sum = 0 + layer_count = 0 - number_of_cylinder_layers = count_cylinder_layers(radial_build) + number_of_cylinder_layers = count_cylinder_layers(radial_build) - for index, item in enumerate(radial_build): - if item[0] == LayerType.PLASMA: - break + for index, item in enumerate(radial_build): + if item[0] == LayerType.PLASMA: + break - if item[0] == LayerType.GAP: - total_sum += item[1] - continue + if item[0] == LayerType.GAP: + total_sum += item[1] + continue - thickness = item[1] - layer_count += 1 + thickness = item[1] + layer_count += 1 - if layer_count > number_of_cylinder_layers: - break + if layer_count > number_of_cylinder_layers: + break - cylinder = center_column_shield_cylinder( - inner_radius=total_sum, - thickness=item[1], - name=f"layer_{layer_count}", - rotation_angle=rotation_angle, - height=center_column_shield_height, - ) - total_sum += thickness - cylinders.append(cylinder) - return cylinders + cylinder = center_column_shield_cylinder( + inner_radius=total_sum, + thickness=item[1], + name=f"layer_{layer_count}", + rotation_angle=rotation_angle, + height=center_column_shield_height, + ) + total_sum += thickness + cylinders.append(cylinder) + return cylinders def distance_to_plasma(radial_build, index): - distance = 0 - for item in radial_build[index + 1 :]: - if item[0] == LayerType.PLASMA: - break - distance += item[1] - return distance + distance = 0 + for item in radial_build[index + 1 :]: + if item[0] == LayerType.PLASMA: + break + distance += item[1] + return distance def create_layers_from_plasma( - radial_build, vertical_build, minor_radius, major_radius, triangularity, elongation, rotation_angle, center_column + radial_build, vertical_build, minor_radius, major_radius, triangularity, elongation, rotation_angle, center_column ): - plasma_index_rb = get_plasma_index(radial_build) - plasma_index_vb = get_plasma_index(vertical_build) - indexes_from_plamsa_to_end = len(radial_build) - plasma_index_rb - layers = [] - - cumulative_thickness_orb = 0 - cumulative_thickness_irb = 0 - cumulative_thickness_uvb = 0 - cumulative_thickness_lvb = 0 - for index_delta in range(indexes_from_plamsa_to_end): - - if radial_build[plasma_index_rb + index_delta][0] == LayerType.PLASMA: - continue - outer_layer_thickness = radial_build[plasma_index_rb + index_delta][1] - inner_layer_thickness = radial_build[plasma_index_rb - index_delta][1] - upper_layer_thickness = vertical_build[plasma_index_vb - index_delta][1] - lower_layer_thickness = vertical_build[plasma_index_vb + index_delta][1] - - if radial_build[plasma_index_rb + index_delta][0] == LayerType.GAP: - cumulative_thickness_orb += outer_layer_thickness - cumulative_thickness_irb += inner_layer_thickness - cumulative_thickness_uvb += upper_layer_thickness - cumulative_thickness_lvb += lower_layer_thickness - continue - - # build outer layer - if radial_build[plasma_index_rb + index_delta][0] == LayerType.SOLID: - outer_layer = blanket_from_plasma( - minor_radius=minor_radius, - major_radius=major_radius, - triangularity=triangularity, - elongation=elongation, - thickness=[ - upper_layer_thickness, - outer_layer_thickness, - lower_layer_thickness - ], - offset_from_plasma=[ - cumulative_thickness_uvb, - cumulative_thickness_orb, - cumulative_thickness_lvb - ], - start_angle=90, - stop_angle=-90, - rotation_angle=rotation_angle, - color=(0.5, 0.5, 0.5), - name=f"layer_{index_delta}", - allow_overlapping_shape=True, - ) - inner_layer = blanket_from_plasma( - minor_radius=minor_radius, - major_radius=major_radius, - triangularity=triangularity, - elongation=elongation, - thickness=[ - lower_layer_thickness, - inner_layer_thickness, - upper_layer_thickness, - ], - offset_from_plasma=[ - cumulative_thickness_lvb, - cumulative_thickness_irb, - cumulative_thickness_uvb, - ], - start_angle=-90, - stop_angle=-270, - rotation_angle=rotation_angle, - color=(0.5, 0.5, 0.5), - name=f"layer_{index_delta}", - allow_overlapping_shape=True, - ) - layer = outer_layer.union(inner_layer) - layers.append(layer) - # layers.append(inner_layer) - cumulative_thickness_orb += outer_layer_thickness - cumulative_thickness_irb += inner_layer_thickness - cumulative_thickness_uvb += upper_layer_thickness - cumulative_thickness_lvb += lower_layer_thickness - # build inner layer - - # union layers - - return layers + plasma_index_rb = get_plasma_index(radial_build) + plasma_index_vb = get_plasma_index(vertical_build) + indexes_from_plamsa_to_end = len(radial_build) - plasma_index_rb + layers = [] + + cumulative_thickness_orb = 0 + cumulative_thickness_irb = 0 + cumulative_thickness_uvb = 0 + cumulative_thickness_lvb = 0 + for index_delta in range(indexes_from_plamsa_to_end): + + if radial_build[plasma_index_rb + index_delta][0] == LayerType.PLASMA: + continue + outer_layer_thickness = radial_build[plasma_index_rb + index_delta][1] + inner_layer_thickness = radial_build[plasma_index_rb - index_delta][1] + upper_layer_thickness = vertical_build[plasma_index_vb - index_delta][1] + lower_layer_thickness = vertical_build[plasma_index_vb + index_delta][1] + + if radial_build[plasma_index_rb + index_delta][0] == LayerType.GAP: + cumulative_thickness_orb += outer_layer_thickness + cumulative_thickness_irb += inner_layer_thickness + cumulative_thickness_uvb += upper_layer_thickness + cumulative_thickness_lvb += lower_layer_thickness + continue + + # build outer layer + if radial_build[plasma_index_rb + index_delta][0] == LayerType.SOLID: + outer_layer = blanket_from_plasma( + minor_radius=minor_radius, + major_radius=major_radius, + triangularity=triangularity, + elongation=elongation, + thickness=[ + upper_layer_thickness, + outer_layer_thickness, + lower_layer_thickness + ], + offset_from_plasma=[ + cumulative_thickness_uvb, + cumulative_thickness_orb, + cumulative_thickness_lvb + ], + start_angle=90, + stop_angle=-90, + rotation_angle=rotation_angle, + color=(0.5, 0.5, 0.5), + name=f"layer_{index_delta}", + allow_overlapping_shape=True, + ) + inner_layer = blanket_from_plasma( + minor_radius=minor_radius, + major_radius=major_radius, + triangularity=triangularity, + elongation=elongation, + thickness=[ + lower_layer_thickness, + inner_layer_thickness, + upper_layer_thickness, + ], + offset_from_plasma=[ + cumulative_thickness_lvb, + cumulative_thickness_irb, + cumulative_thickness_uvb, + ], + start_angle=-90, + stop_angle=-270, + rotation_angle=rotation_angle, + color=(0.5, 0.5, 0.5), + name=f"layer_{index_delta}", + allow_overlapping_shape=True, + ) + layer = outer_layer.union(inner_layer) + layers.append(layer) + # layers.append(inner_layer) + cumulative_thickness_orb += outer_layer_thickness + cumulative_thickness_irb += inner_layer_thickness + cumulative_thickness_uvb += upper_layer_thickness + cumulative_thickness_lvb += lower_layer_thickness + # build inner layer + + # union layers + + return layers def tokamak_from_plasma( - radial_build: Sequence[Tuple[LayerType, float]], - elongation: float = 2.0, - triangularity: float = 0.55, - rotation_angle: float = 180.0, - extra_cut_shapes: Sequence[cq.Workplane] = [], - extra_intersect_shapes: Sequence[cq.Workplane] = [], + radial_build: Sequence[Tuple[LayerType, float]], + elongation: float = 2.0, + triangularity: float = 0.55, + rotation_angle: float = 180.0, + extra_cut_shapes: Sequence[cq.Workplane] = [], + extra_intersect_shapes: Sequence[cq.Workplane] = [], ): - inner_equatorial_point = sum_up_to_plasma(radial_build) - plasma_radial_thickness = get_plasma_value(radial_build) - outer_equatorial_point = inner_equatorial_point + plasma_radial_thickness - - # sets major radius and minor radius from equatorial_points to allow a - # radial build. This helps avoid the plasma overlapping the center - # column and other components - major_radius = (outer_equatorial_point + inner_equatorial_point) / 2 - minor_radius = major_radius - inner_equatorial_point - - # make vertical build from outer radial build - # make vertical build from outer radial build - pi = get_plasma_index(radial_build) - #upper_vertical_build = radial_build[pi:] - #changing to inner vertical build - rbi = len(radial_build)-1 - pi - upper_vertical_build = radial_build[pi-rbi:pi][::-1] - - plasma_height = 2 * minor_radius * elongation - # slice opperation reverses the list and removes the last value to avoid two plasmas - vertical_build = upper_vertical_build[::-1] + [(LayerType.PLASMA, plasma_height)] + upper_vertical_build - - return tokamak( - radial_build=radial_build, - vertical_build=vertical_build, - triangularity=triangularity, - rotation_angle=rotation_angle, - extra_cut_shapes=extra_cut_shapes, - extra_intersect_shapes=extra_intersect_shapes, - ) + inner_equatorial_point = sum_up_to_plasma(radial_build) + plasma_radial_thickness = get_plasma_value(radial_build) + outer_equatorial_point = inner_equatorial_point + plasma_radial_thickness + + # sets major radius and minor radius from equatorial_points to allow a + # radial build. This helps avoid the plasma overlapping the center + # column and other components + major_radius = (outer_equatorial_point + inner_equatorial_point) / 2 + minor_radius = major_radius - inner_equatorial_point + + # make vertical build from inner radial build + pi = get_plasma_index(radial_build) + rbi = len(radial_build)-1 - pi # number of unique entries in outer or inner radial build + upper_vertical_build = radial_build[pi-rbi:pi][::-1] #get the inner radial build + + plasma_height = 2 * minor_radius * elongation + # slice opperation reverses the list and removes the last value to avoid two plasmas + vertical_build = upper_vertical_build[::-1] + [(LayerType.PLASMA, plasma_height)] + upper_vertical_build + + return tokamak( + radial_build=radial_build, + vertical_build=vertical_build, + triangularity=triangularity, + rotation_angle=rotation_angle, + extra_cut_shapes=extra_cut_shapes, + extra_intersect_shapes=extra_intersect_shapes, + ) def tokamak( - radial_build: Union[Sequence[Sequence[Tuple[str, float]]], Sequence[Tuple[str, float]]], - vertical_build: Sequence[Tuple[str, float]], - triangularity: float = 0.55, - rotation_angle: float = 180.0, - extra_cut_shapes: Sequence[cq.Workplane] = [], - extra_intersect_shapes: Sequence[cq.Workplane] = [], + radial_build: Union[Sequence[Sequence[Tuple[str, float]]], Sequence[Tuple[str, float]]], + vertical_build: Sequence[Tuple[str, float]], + triangularity: float = 0.55, + rotation_angle: float = 180.0, + extra_cut_shapes: Sequence[cq.Workplane] = [], + extra_intersect_shapes: Sequence[cq.Workplane] = [], ): - """ - Creates a tokamak fusion reactor from a radial build and plasma parameters. - - Args: - radial_build: A list of tuples containing the radial build of the reactor. - elongation: The elongation of the plasma. Defaults to 2.0. - triangularity: The triangularity of the plasma. Defaults to 0.55. - rotation_angle: The rotation angle of the plasma. Defaults to 180.0. - extra_cut_shapes: A list of extra shapes to cut the reactor with. Defaults to []. - extra_intersect_shapes: A list of extra shapes to intersect the reactor with. Defaults to []. - - Returns: - CadQuery.Assembly: A CadQuery Assembly object representing the tokamak fusion reactor. - """ - - inner_equatorial_point = sum_up_to_plasma(radial_build) - plasma_radial_thickness = get_plasma_value(radial_build) - plasma_vertical_thickness = get_plasma_value(vertical_build) - outer_equatorial_point = inner_equatorial_point + plasma_radial_thickness - - major_radius = (outer_equatorial_point + inner_equatorial_point) / 2 - minor_radius = major_radius - inner_equatorial_point - - elongation = (plasma_vertical_thickness / 2) / minor_radius - blanket_rear_wall_end_height = sum([item[1] for item in vertical_build]) - - plasma = plasma_simplified( - major_radius=major_radius, - minor_radius=minor_radius, - elongation=elongation, - triangularity=triangularity, - rotation_angle=rotation_angle, - ) - - inner_radial_build = create_center_column_shield_cylinders( - radial_build, rotation_angle, blanket_rear_wall_end_height - ) - - blanket_layers = create_layers_from_plasma( - radial_build=radial_build, - vertical_build=vertical_build, - minor_radius=minor_radius, - major_radius=major_radius, - triangularity=triangularity, - elongation=elongation, - rotation_angle=rotation_angle, - center_column=inner_radial_build[0], # blanket_cutting_cylinder, - ) - - my_assembly = cq.Assembly() - - for i, entry in enumerate(extra_cut_shapes): - if isinstance(entry, cq.Workplane): - my_assembly.add(entry, name=f"add_extra_cut_shape_{i+1}") - else: - raise ValueError(f"extra_cut_shapes should only contain cadquery Workplanes, not {type(entry)}") - - # builds up the intersect shapes - intersect_shapes_to_cut = [] - if len(extra_intersect_shapes)>0: - all_shapes = [] - for shape in inner_radial_build + blanket_layers: - all_shapes.append(shape) - - # makes a union of the the radial build to use as a base for the intersect shapes - reactor_compound=inner_radial_build[0] - for i, entry in enumerate(inner_radial_build[1:] + blanket_layers): - reactor_compound = reactor_compound.union(entry) - - # adds the extra intersect shapes to the assembly - for i, entry in enumerate(extra_intersect_shapes): - reactor_entry_intersection = entry.intersect(reactor_compound) - intersect_shapes_to_cut.append(reactor_entry_intersection) - my_assembly.add(reactor_entry_intersection, name=f"extra_intersect_shapes_{i+1}") - - # builds just the core if there are no extra parts - if len(extra_cut_shapes) == 0 and len(intersect_shapes_to_cut) == 0: - for i, entry in enumerate(inner_radial_build): - my_assembly.add(entry, name=f"inboard_layer_{i+1})") - for i, entry in enumerate(blanket_layers): - my_assembly.add(entry, name=f"outboard_layer_{i+1})") - else: - shapes_and_components = [] - for i, entry in enumerate(inner_radial_build + blanket_layers): - for cutter in extra_cut_shapes + extra_intersect_shapes: - entry = entry.cut(cutter) - # TODO use something like this to return a list of material tags for the solids in order, as some solids get split into multiple - # for subentry in entry.objects: - # print(i, subentry) - shapes_and_components.append(entry) - - for i, entry in enumerate(shapes_and_components): - my_assembly.add(entry, name=f"layer_{i+1})") #TODO track the names of shapes, even when extra shapes are made due to splitting - - - my_assembly.add(plasma, name="plasma") - - return my_assembly + """ + Creates a tokamak fusion reactor from a radial build and plasma parameters. + + Args: + radial_build: A list of tuples containing the radial build of the reactor. + elongation: The elongation of the plasma. Defaults to 2.0. + triangularity: The triangularity of the plasma. Defaults to 0.55. + rotation_angle: The rotation angle of the plasma. Defaults to 180.0. + extra_cut_shapes: A list of extra shapes to cut the reactor with. Defaults to []. + extra_intersect_shapes: A list of extra shapes to intersect the reactor with. Defaults to []. + + Returns: + CadQuery.Assembly: A CadQuery Assembly object representing the tokamak fusion reactor. + """ + + inner_equatorial_point = sum_up_to_plasma(radial_build) + plasma_radial_thickness = get_plasma_value(radial_build) + plasma_vertical_thickness = get_plasma_value(vertical_build) + outer_equatorial_point = inner_equatorial_point + plasma_radial_thickness + + major_radius = (outer_equatorial_point + inner_equatorial_point) / 2 + minor_radius = major_radius - inner_equatorial_point + + elongation = (plasma_vertical_thickness / 2) / minor_radius + blanket_rear_wall_end_height = sum([item[1] for item in vertical_build]) + + plasma = plasma_simplified( + major_radius=major_radius, + minor_radius=minor_radius, + elongation=elongation, + triangularity=triangularity, + rotation_angle=rotation_angle, + ) + + inner_radial_build = create_center_column_shield_cylinders( + radial_build, rotation_angle, blanket_rear_wall_end_height + ) + + blanket_layers = create_layers_from_plasma( + radial_build=radial_build, + vertical_build=vertical_build, + minor_radius=minor_radius, + major_radius=major_radius, + triangularity=triangularity, + elongation=elongation, + rotation_angle=rotation_angle, + center_column=inner_radial_build[0], # blanket_cutting_cylinder, + ) + + my_assembly = cq.Assembly() + + for i, entry in enumerate(extra_cut_shapes): + if isinstance(entry, cq.Workplane): + my_assembly.add(entry, name=f"add_extra_cut_shape_{i+1}") + else: + raise ValueError(f"extra_cut_shapes should only contain cadquery Workplanes, not {type(entry)}") + + # builds up the intersect shapes + intersect_shapes_to_cut = [] + if len(extra_intersect_shapes)>0: + all_shapes = [] + for shape in inner_radial_build + blanket_layers: + all_shapes.append(shape) + + # makes a union of the the radial build to use as a base for the intersect shapes + reactor_compound=inner_radial_build[0] + for i, entry in enumerate(inner_radial_build[1:] + blanket_layers): + reactor_compound = reactor_compound.union(entry) + + # adds the extra intersect shapes to the assembly + for i, entry in enumerate(extra_intersect_shapes): + reactor_entry_intersection = entry.intersect(reactor_compound) + intersect_shapes_to_cut.append(reactor_entry_intersection) + my_assembly.add(reactor_entry_intersection, name=f"extra_intersect_shapes_{i+1}") + + # builds just the core if there are no extra parts + if len(extra_cut_shapes) == 0 and len(intersect_shapes_to_cut) == 0: + for i, entry in enumerate(inner_radial_build): + my_assembly.add(entry, name=f"inboard_layer_{i+1})") + for i, entry in enumerate(blanket_layers): + my_assembly.add(entry, name=f"outboard_layer_{i+1})") + else: + shapes_and_components = [] + for i, entry in enumerate(inner_radial_build + blanket_layers): + for cutter in extra_cut_shapes + extra_intersect_shapes: + entry = entry.cut(cutter) + # TODO use something like this to return a list of material tags for the solids in order, as some solids get split into multiple + # for subentry in entry.objects: + # print(i, subentry) + shapes_and_components.append(entry) + + for i, entry in enumerate(shapes_and_components): + my_assembly.add(entry, name=f"layer_{i+1})") #TODO track the names of shapes, even when extra shapes are made due to splitting + + + my_assembly.add(plasma, name="plasma") + + return my_assembly