diff --git a/examples/example_util_functions.py b/examples/example_util_functions.py deleted file mode 100644 index d9f691c0..00000000 --- a/examples/example_util_functions.py +++ /dev/null @@ -1,102 +0,0 @@ -def transport_particles_on_h5m_geometry( - h5m_filename: str, - material_tags: list, - nuclides: list = None, - cross_sections_xml: str = None, -): - """A function for testing the geometry file with particle transport in - DAGMC OpenMC. Requires openmc and either the cross_sections_xml to be - specified. Returns the flux on each volume - - Arg: - h5m_filename: The name of the DAGMC h5m file to test - material_tags: The - nuclides: - cross_sections_xml: - - """ - - import openmc - from openmc.data import NATURAL_ABUNDANCE - - if nuclides is None: - nuclides = list(NATURAL_ABUNDANCE.keys()) - - materials = openmc.Materials() - for i, material_tag in enumerate(set(material_tags)): - # simplified material definitions have been used to keen this example minimal - mat_dag_material_tag = openmc.Material(name=material_tag) - mat_dag_material_tag.add_nuclide(nuclides[i], 1, "ao") - mat_dag_material_tag.set_density("g/cm3", 0.01) - - materials.append(mat_dag_material_tag) - - if cross_sections_xml: - openmc.config["cross_sections"] = cross_sections_xml - - else: - with open("cross_sections.xml", "w") as file: - file.write( - """ - - - - - - """ - ) - - openmc.config["cross_sections"] = "cross_sections.xml" - - dag_univ = openmc.DAGMCUniverse(filename=h5m_filename) - bound_dag_univ = dag_univ.bounded_universe() - geometry = openmc.Geometry(root=bound_dag_univ) - - # initializes a new source object - my_source = openmc.IndependentSource() - - center_of_geometry = ( - (dag_univ.bounding_box[0][0] + dag_univ.bounding_box[1][0]) / 2, - (dag_univ.bounding_box[0][1] + dag_univ.bounding_box[1][1]) / 2, - (dag_univ.bounding_box[0][2] + dag_univ.bounding_box[1][2]) / 2, - ) - # sets the location of the source which is not on a vertex - center_of_geometry_nudged = ( - center_of_geometry[0] + 0.1, - center_of_geometry[1] + 0.1, - center_of_geometry[2] + 0.1, - ) - - my_source.space = openmc.stats.Point(center_of_geometry_nudged) - # sets the direction to isotropic - my_source.angle = openmc.stats.Isotropic() - # sets the energy distribution to 100% 14MeV neutrons - my_source.energy = openmc.stats.Discrete([14e6], [1]) - - # specifies the simulation computational intensity - settings = openmc.Settings() - settings.batches = 10 - settings.particles = 10000 - settings.inactive = 0 - settings.run_mode = "fixed source" - settings.source = my_source - - # adds a tally to record the heat deposited in entire geometry - cell_tally = openmc.Tally(name="flux") - cell_tally.scores = ["flux"] - - # groups the two tallies - tallies = openmc.Tallies([cell_tally]) - - # builds the openmc model - my_model = openmc.Model(materials=materials, geometry=geometry, settings=settings, tallies=tallies) - - # starts the simulation - output_file = my_model.run() - - # loads up the output file from the simulation - statepoint = openmc.StatePoint(output_file) - - my_flux_cell_tally = statepoint.get_tally(name="flux") - - return my_flux_cell_tally.mean.flatten()[0] diff --git a/examples/spherical_tokamak_from_plasma_minimal.py b/examples/spherical_tokamak_from_plasma_minimal.py index f18917d0..48428d70 100644 --- a/examples/spherical_tokamak_from_plasma_minimal.py +++ b/examples/spherical_tokamak_from_plasma_minimal.py @@ -1,7 +1,5 @@ from pathlib import Path -from example_util_functions import transport_particles_on_h5m_geometry - import paramak my_reactor = paramak.spherical_tokamak_from_plasma( @@ -21,21 +19,3 @@ rotation_angle=180, ) my_reactor.save(f"spherical_tokamak_from_plasma_minimal.step") - - -# from cad_to_dagmc import CadToDagmc -# my_model = CadToDagmc() -# material_tags = [ -# "mat1" -# ] * 6 -# my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) -# my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) - -# h5m_filename = "dagmc.h5m" -# flux = transport_particles_on_h5m_geometry( -# h5m_filename=h5m_filename, -# material_tags=material_tags, -# nuclides=["H1"] * len(material_tags), -# cross_sections_xml="tests/cross_sections.xml", -# ) -# assert flux > 0.0 diff --git a/examples/spherical_tokamak_from_plasma_with_color.py b/examples/spherical_tokamak_from_plasma_with_color.py index bade0de2..ca186e0f 100644 --- a/examples/spherical_tokamak_from_plasma_with_color.py +++ b/examples/spherical_tokamak_from_plasma_with_color.py @@ -1,7 +1,3 @@ -from pathlib import Path - -from example_util_functions import transport_particles_on_h5m_geometry - import paramak my_reactor = paramak.spherical_tokamak_from_plasma( @@ -54,19 +50,3 @@ file_path='spherical_tokamak_from_plasma_with_colors.png' ) -# from cad_to_dagmc import CadToDagmc -# my_model = CadToDagmc() -# material_tags = [ -# "mat1" -# ] * 6 -# my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) -# my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) - -# h5m_filename = "dagmc.h5m" -# flux = transport_particles_on_h5m_geometry( -# h5m_filename=h5m_filename, -# material_tags=material_tags, -# nuclides=["H1"] * len(material_tags), -# cross_sections_xml="tests/cross_sections.xml", -# ) -# assert flux > 0.0 diff --git a/examples/spherical_tokamak_from_plasma_with_divertor.py b/examples/spherical_tokamak_from_plasma_with_divertor.py index 6d06bd2d..0d8b6a05 100644 --- a/examples/spherical_tokamak_from_plasma_with_divertor.py +++ b/examples/spherical_tokamak_from_plasma_with_divertor.py @@ -1,9 +1,5 @@ -from pathlib import Path - -from example_util_functions import transport_particles_on_h5m_geometry - import paramak -from cadquery import Workplane, vis +from cadquery import Workplane # makes a rectangle that overlaps the lower blanket under the plasma @@ -30,19 +26,3 @@ ) my_reactor.save("spherical_tokamak_from_plasma_with_divertor.step") print("written spherical_tokamak_from_plasma_with_divertor.step") - -# from cad_to_dagmc import CadToDagmc -# vis.show(my_reactor) -# my_model = CadToDagmc() -# material_tags = ["mat1"] * 21 # the two divertors split the 3 blanket layers into 9 and the magnets also splt the blanket. -# my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) -# my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) - -# h5m_filename = "dagmc.h5m" -# flux = transport_particles_on_h5m_geometry( -# h5m_filename=h5m_filename, -# material_tags=material_tags, -# nuclides=["H1"] * len(material_tags), -# cross_sections_xml="tests/cross_sections.xml", -# ) -# assert flux > 0.0 diff --git a/examples/spherical_tokamak_from_plasma_with_pf_magnets.py b/examples/spherical_tokamak_from_plasma_with_pf_magnets.py index 4fc8a21a..c4aa120b 100644 --- a/examples/spherical_tokamak_from_plasma_with_pf_magnets.py +++ b/examples/spherical_tokamak_from_plasma_with_pf_magnets.py @@ -1,5 +1,3 @@ -from example_util_functions import transport_particles_on_h5m_geometry - import paramak extra_cut_shapes = [] @@ -38,19 +36,3 @@ extra_cut_shapes=extra_cut_shapes, ) my_reactor.save(f"spherical_tokamak_from_plasma_with_pf_magnets.step") - - -# from cad_to_dagmc import CadToDagmc -# my_model = CadToDagmc() -# material_tags = ["mat1"] * 5 -# my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) -# my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) - -# h5m_filename = "dagmc.h5m" -# flux = transport_particles_on_h5m_geometry( -# h5m_filename=h5m_filename, -# material_tags=material_tags, -# nuclides=["H1"] * len(material_tags), -# cross_sections_xml="tests/cross_sections.xml", -# ) -# assert flux > 0.0 diff --git a/examples/spherical_tokamak_from_plasma_with_pf_magnets_and_divertors.py b/examples/spherical_tokamak_from_plasma_with_pf_magnets_and_divertors.py index 0923acb9..d0c72c27 100644 --- a/examples/spherical_tokamak_from_plasma_with_pf_magnets_and_divertors.py +++ b/examples/spherical_tokamak_from_plasma_with_pf_magnets_and_divertors.py @@ -1,7 +1,3 @@ -from pathlib import Path - -from example_util_functions import transport_particles_on_h5m_geometry - import paramak rotation_angle = 180 @@ -44,19 +40,3 @@ extra_cut_shapes=poloidal_field_coils, ) my_reactor.save(f"spherical_tokamak_from_plasma_with_pf_magnets_and_divertor.step") - - -# from cad_to_dagmc import CadToDagmc -# my_model = CadToDagmc() -# material_tags = ["mat1"] * 5 -# my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) -# my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) - -# h5m_filename = "dagmc.h5m" -# flux = transport_particles_on_h5m_geometry( -# h5m_filename=h5m_filename, -# material_tags=material_tags, -# nuclides=["H1"] * len(material_tags), -# cross_sections_xml="tests/cross_sections.xml", -# ) -# assert flux > 0.0 diff --git a/examples/spherical_tokamak_from_plasma_with_tf_magnets.py b/examples/spherical_tokamak_from_plasma_with_tf_magnets.py index 7d55f61d..d765bc78 100644 --- a/examples/spherical_tokamak_from_plasma_with_tf_magnets.py +++ b/examples/spherical_tokamak_from_plasma_with_tf_magnets.py @@ -61,19 +61,3 @@ ) 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 -# my_model.add_cadquery_object(cadquery_object=result, material_tags=material_tags) -# my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) - -# h5m_filename = "dagmc.h5m" -# flux = transport_particles_on_h5m_geometry( -# h5m_filename=h5m_filename, -# material_tags=material_tags, -# nuclides=["H1"] * len(material_tags), -# cross_sections_xml="tests/cross_sections.xml", -# ) -# assert flux > 0.0 diff --git a/examples/spherical_tokamak_minimal.py b/examples/spherical_tokamak_minimal.py index d1fa38c6..bd87e5cc 100644 --- a/examples/spherical_tokamak_minimal.py +++ b/examples/spherical_tokamak_minimal.py @@ -1,7 +1,3 @@ -from pathlib import Path - -from example_util_functions import transport_particles_on_h5m_geometry - import paramak my_reactor = paramak.spherical_tokamak( @@ -31,19 +27,3 @@ triangularity=-0.55, ) my_reactor.save(f"spherical_tokamak_minimal.step") - - -# from cad_to_dagmc import CadToDagmc -# my_model = CadToDagmc() -# material_tags = ["mat1"] * 6 -# my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) -# my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) - -# h5m_filename = "dagmc.h5m" -# flux = transport_particles_on_h5m_geometry( -# h5m_filename=h5m_filename, -# material_tags=material_tags, -# nuclides=["H1"] * len(material_tags), -# cross_sections_xml="tests/cross_sections.xml", -# ) -# assert flux > 0.0 diff --git a/examples/tokamak_from_plasma_animation.py b/examples/tokamak_from_plasma_animation.py index b2e397c2..f915cf3b 100644 --- a/examples/tokamak_from_plasma_animation.py +++ b/examples/tokamak_from_plasma_animation.py @@ -3,6 +3,7 @@ import cadquery_png_plugin.plugin import numpy as np import cadquery as cq +import os original_radial_build=[ (paramak.LayerType.GAP, 40), diff --git a/examples/tokamak_from_plasma_minimal.py b/examples/tokamak_from_plasma_minimal.py index 53be67b2..98001b2b 100644 --- a/examples/tokamak_from_plasma_minimal.py +++ b/examples/tokamak_from_plasma_minimal.py @@ -1,5 +1,3 @@ -from example_util_functions import transport_particles_on_h5m_geometry - import paramak my_reactor = paramak.tokamak_from_plasma( @@ -23,18 +21,3 @@ ) my_reactor.save(f"tokamak_minimal.step") print(f"Saved as tokamak_minimal.step") - -# from cad_to_dagmc import CadToDagmc -# my_model = CadToDagmc() -# material_tags = ["mat1"] * 6 # as inner and outer layers are one solid there are only 6 solids in model -# my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) -# my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) - -# h5m_filename = "dagmc.h5m" -# flux = transport_particles_on_h5m_geometry( -# h5m_filename=h5m_filename, -# material_tags=material_tags, -# nuclides=["H1"] * len(material_tags), -# cross_sections_xml="tests/cross_sections.xml", -# ) -# assert flux > 0.0 diff --git a/examples/tokamak_from_plasma_with_color.py b/examples/tokamak_from_plasma_with_color.py index 5d51f09e..c2fe6967 100644 --- a/examples/tokamak_from_plasma_with_color.py +++ b/examples/tokamak_from_plasma_with_color.py @@ -1,5 +1,3 @@ -from example_util_functions import transport_particles_on_h5m_geometry - import paramak my_reactor = paramak.tokamak_from_plasma( @@ -56,18 +54,3 @@ }, file_path='tokamak_from_plasma_with_colors.png' ) - -# from cad_to_dagmc import CadToDagmc -# my_model = CadToDagmc() -# material_tags = ["mat1"] * 6 # as inner and outer layers are one solid there are only 6 solids in model -# my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) -# my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) - -# h5m_filename = "dagmc.h5m" -# flux = transport_particles_on_h5m_geometry( -# h5m_filename=h5m_filename, -# material_tags=material_tags, -# nuclides=["H1"] * len(material_tags), -# cross_sections_xml="tests/cross_sections.xml", -# ) -# assert flux > 0.0 diff --git a/examples/tokamak_from_plasma_with_divertor.py b/examples/tokamak_from_plasma_with_divertor.py index 86f6b429..ccdd9fb6 100644 --- a/examples/tokamak_from_plasma_with_divertor.py +++ b/examples/tokamak_from_plasma_with_divertor.py @@ -1,7 +1,5 @@ -from example_util_functions import transport_particles_on_h5m_geometry -from paramak.utils import create_wire_workplane_from_points import paramak -from cadquery import vis, Workplane +from cadquery import Workplane # makes a rectangle that overlaps the lower blanket under the plasma # the intersection of this and the layers will form the lower divertor @@ -31,19 +29,4 @@ ) my_reactor.save(f"tokamak_with_divertor.step") print(f"Saved as tokamak_with_divertor.step") -# vis.show(my_reactor) -# from cad_to_dagmc import CadToDagmc -# my_model = CadToDagmc() -# material_tags = ["mat1"] * 6 # as inner and outer layers are one solid there are only 6 solids in model -# my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) -# my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) - -# h5m_filename = "dagmc.h5m" -# flux = transport_particles_on_h5m_geometry( -# h5m_filename=h5m_filename, -# material_tags=material_tags, -# nuclides=["H1"] * len(material_tags), -# cross_sections_xml="tests/cross_sections.xml", -# ) -# assert flux > 0.0 diff --git a/examples/tokamak_minimal.py b/examples/tokamak_minimal.py index 8b5ebb33..e22b26da 100644 --- a/examples/tokamak_minimal.py +++ b/examples/tokamak_minimal.py @@ -1,5 +1,3 @@ -from example_util_functions import transport_particles_on_h5m_geometry - import paramak my_reactor = paramak.tokamak( @@ -34,18 +32,3 @@ my_reactor.save(f"tokamak_minimal.step") print(f"Saved as tokamak_minimal.step") - -# from cad_to_dagmc import CadToDagmc -# my_model = CadToDagmc() -# material_tags = ["mat1"] * 6 # as inner and outer layers are one solid there are only 6 solids in model -# my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) -# my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) - -# h5m_filename = "dagmc.h5m" -# flux = transport_particles_on_h5m_geometry( -# h5m_filename=h5m_filename, -# material_tags=material_tags, -# nuclides=["H1"] * len(material_tags), -# cross_sections_xml="tests/cross_sections.xml", -# ) -# assert flux > 0.0 diff --git a/examples/tokamak_structured_mesh_simulation_with_CSG_bioshield.py b/examples/tokamak_structured_mesh_simulation_with_CSG_bioshield.py new file mode 100644 index 00000000..957c9fde --- /dev/null +++ b/examples/tokamak_structured_mesh_simulation_with_CSG_bioshield.py @@ -0,0 +1,228 @@ +""" +This script sets up and runs a neutronics simulation for a tokamak reactor +using OpenMC and Paramak. The geometry combines CSG for the simple bioshield +and CAD based (DAGMC) for the reactor. This example is minimal is terms of the +tokamak geometry, materials, neutron source. It is intended to be used as a +guide for carrying out a simulation that uses both DAGMC with CSG geometry. +This is particularly useful as it avoids using DAGMC for the entire geometry +so we can benefit from DAGMC for the more complex curved surfaces of the +tokamak and benefit from the CSG speed for the simple geometry of the bioshield. +It performs the following steps: +1. Defines the tokamak reactor geometry. +2. Removes the plasma cell as this is not needed by the neutronics simulation. +3. Converts the geometry to DAGMC format. +4. Creates simplified materials. +5. Creates a bioshield complete with floor and ceiling in CSG geometry. +6. Moves the DAGMC geometry so that it is within the CSG bioshield. +7. Adds a cylindrical mesh tally to record the heat deposited in the geometry. +8. Adds a simplified neutron source term +9. Runs the simulation. +10. Saves the tally results to a VTK file. +""" + +from pathlib import Path +import openmc +import numpy as np +import math +import paramak +import cadquery as cq +from cad_to_dagmc import CadToDagmc + +# change this path to the cross_sections.xml to the path of your cross_sections.xml +openmc.config["cross_sections"] = "/nuclear_data/cross_sections.xml" + + +# parameters for the CSG parts of the geometry. +bioshield_to_reactor_gap = 500 +bioshield_radial_thickness = 200 +floor_vertical_thickness = 200 +ceiling_vertical_thickness = 200 +reactor_to_floor_gap = 100 +reactor_to_ceiling_gap = 150 + +# creates a minimal tokamak reactor using paramak +my_reactor = paramak.tokamak_from_plasma( + radial_build=[ + (paramak.LayerType.GAP, 10), + (paramak.LayerType.SOLID, 30), # layer_1 + (paramak.LayerType.SOLID, 50), # layer_2 + (paramak.LayerType.SOLID, 10), # layer_3 + (paramak.LayerType.SOLID, 120), # layer_4 + (paramak.LayerType.SOLID, 20), # layer_5 + (paramak.LayerType.GAP, 60), + (paramak.LayerType.PLASMA, 300), # plasma + (paramak.LayerType.GAP, 60), + (paramak.LayerType.SOLID, 20), # layer_3 + (paramak.LayerType.SOLID, 120), # layer_4 + (paramak.LayerType.SOLID, 10), # layer_5 + ], + elongation=2, + triangularity=0.55, + rotation_angle=180, +) + +# removing the plasma from the CadQuery Assembly as we don't need this in the +# DAGMC mesh for the neutronics simulation. +my_reactor = my_reactor.remove(name="plasma") + + +# this prints all the names of the parts in the reactor, this is useful for +# knowing the material tags to use in the simulation +print(my_reactor.names()) + +material_tags = ["layer_1", "layer_2", "layer_3", "layer_4", "layer_5"] +my_model = CadToDagmc() + +# 6 material tags as inner and outer layers are one solid there are only 6 solids in model and plasma has been removed +my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) + +script_folder = Path(__file__).resolve().parent +h5m_filename = script_folder / "dagmc.h5m" +my_model.export_dagmc_h5m_file(filename=h5m_filename, min_mesh_size=10.0, max_mesh_size=20.0) + +# simplified material definitions have been used to keep this example minimal +mat_layer_1 = openmc.Material(name="layer_1") +mat_layer_1.add_element("Cu", 1, "ao") +mat_layer_1.set_density("g/cm3", 7) + +mat_layer_2 = openmc.Material(name="layer_2") +mat_layer_2.add_nuclide("W186", 1, "ao") +mat_layer_2.set_density("g/cm3", 0.01) + +mat_layer_3 = openmc.Material(name="layer_3") +mat_layer_3.add_nuclide("Fe56", 1, "ao") +mat_layer_3.set_density("g/cm3", 7) + +mat_layer_4 = openmc.Material(name="layer_4") +mat_layer_4.add_element("Li", 1, "ao") +mat_layer_4.set_density("g/cm3", 0.5) + +mat_layer_5 = openmc.Material(name="layer_5") +mat_layer_5.add_nuclide("Fe56", 1, "ao") +mat_layer_5.set_density("g/cm3", 7) + +mat_concrete = openmc.Material(name="concrete") +mat_concrete.add_nuclide("Fe56", 1, "ao") +mat_concrete.set_density("g/cm3", 7) + + +materials = openmc.Materials([mat_layer_1, mat_layer_2, mat_layer_3, mat_layer_4, mat_layer_5, mat_concrete]) + +# brings the DAGMC geometry into OpenMC +dag_univ = openmc.DAGMCUniverse(filename=h5m_filename) + +# gets the bounding box of the DAGMC geometry, this is used to position the CSG geometry +bbox = dag_univ.bounding_box +dagmc_geometry_offset = abs(bbox[0][2]) + floor_vertical_thickness + reactor_to_floor_gap + +dagmc_radius = max(abs(bbox[0][0]), abs(bbox[0][1]), abs(bbox[1][0]), abs(bbox[1][1])) + +# we have a 180 degree model and the YPlane acts as a reflective boundary +side_surface = openmc.YPlane(y0=0, boundary_type="reflective", surface_id=1001) +bioshield_inner_surface = openmc.ZCylinder(r=dagmc_radius + bioshield_to_reactor_gap, surface_id=1002) +ceiling_upper_surface = openmc.ZPlane( + z0=floor_vertical_thickness + + reactor_to_floor_gap + + reactor_to_ceiling_gap + + bbox.width[2] + + ceiling_vertical_thickness, + boundary_type="vacuum", + surface_id=1003, +) +ceiling_lower_surface = openmc.ZPlane( + z0=floor_vertical_thickness + reactor_to_floor_gap + reactor_to_ceiling_gap + bbox.width[2], surface_id=1004 +) +floor_upper_surface = openmc.ZPlane(z0=floor_vertical_thickness, surface_id=1005) +floor_lower_surface = openmc.ZPlane(z0=0, boundary_type="vacuum", surface_id=1006) +bioshield_outer_surface = openmc.ZCylinder( + r=dagmc_radius + bioshield_to_reactor_gap + bioshield_radial_thickness, boundary_type="vacuum", surface_id=1007 +) + +# this is a CSG region that is the inner bioshield surface +wedge_region = -bioshield_inner_surface & +floor_upper_surface & -ceiling_lower_surface & +side_surface + +# this is a CSG cell that is filled with the DAGMC geometry +dagmc_filled_cell = openmc.Cell(fill=dag_univ, cell_id=1000, region=wedge_region) +# moving the geometry so that it is within the CSG geometry bioshield. +# The CSG geometry starts at z=0, with a floor and then a gap to the reactor. +# We need to move the geometry up so that it is within the bioshield. +dagmc_filled_cell.translation = (0, 0, dagmc_geometry_offset) + + +bioshield_region = ( + -bioshield_outer_surface & +bioshield_inner_surface & +floor_upper_surface & -ceiling_lower_surface & +side_surface +) +bioshield_cell = openmc.Cell(region=bioshield_region, fill=mat_concrete, cell_id=1001) + +floor_region = -floor_upper_surface & +floor_lower_surface & -bioshield_outer_surface & +side_surface +floor_cell = openmc.Cell(region=floor_region, fill=mat_concrete, cell_id=1002) + +ceiling_region = -ceiling_upper_surface & +ceiling_lower_surface & -bioshield_outer_surface & +side_surface +ceiling_cell = openmc.Cell(region=ceiling_region, fill=mat_concrete, cell_id=1003) + +geometry = openmc.Geometry([dagmc_filled_cell, floor_cell, bioshield_cell, ceiling_cell]) + +# initializes a new source object +my_source = openmc.IndependentSource() + +# the distribution of radius is just a single value +radius = openmc.stats.Discrete([my_reactor.major_radius], [1]) + +# the distribution of source z values is just a single value +z_values = openmc.stats.Discrete([dagmc_geometry_offset], [1]) + +# the distribution of source azimuthal angles values is a uniform distribution between 0 and 2 Pi +angle = openmc.stats.Uniform(a=0.0, b=np.pi) + +# this makes the ring source using the three distributions and a radius +my_source.space = openmc.stats.CylindricalIndependent(r=radius, phi=angle, z=z_values, origin=(0.0, 0.0, 0.0)) + +# sets the direction to isotropic +my_source.angle = openmc.stats.Isotropic() + +# sets the energy distribution to a Muir distribution neutrons +my_source.energy = openmc.stats.muir(e0=14080000.0, m_rat=5.0, kt=20000.0) + + +# specifies the simulation computational intensity +settings = openmc.Settings(batches=10, particles=10000, run_mode="fixed source", source=my_source) + +# adds a tally to record the heat deposited in a mesh overlaid over entire geometry + +mesh = openmc.CylindricalMesh( + r_grid=np.linspace(0, dagmc_radius + bioshield_to_reactor_gap + bioshield_radial_thickness, 100), + z_grid=np.linspace( + 0, + floor_vertical_thickness + + reactor_to_floor_gap + + reactor_to_ceiling_gap + + bbox.width[2] + + ceiling_vertical_thickness, + 100, + ), + origin=(0, 0, 0), + phi_grid=np.linspace(0, np.pi, 100), +) +mesh_tally = openmc.Tally(name="heating") +mesh_tally.filters = [openmc.MeshFilter(mesh)] +mesh_tally.scores = ["heating"] +tallies = openmc.Tallies([mesh_tally]) + +# builds the openmc model +my_model = openmc.Model(materials=materials, geometry=geometry, settings=settings, tallies=tallies) + +# starts the simulation +output_file = my_model.run() + +# loads up the output file from the simulation +statepoint = openmc.StatePoint(output_file) + +my_heating_mesh_tally = statepoint.get_tally(name="heating") + +mesh.write_data_to_vtk( + curvilinear=True, + datasets={"mean": my_heating_mesh_tally.mean.flatten()}, + filename="structured_mesh_tally_results.vtk", +) + +print("VTK file saved to structured_mesh_tally_results.vtk") diff --git a/examples/tokamak_unstructured_mesh_simulation.py b/examples/tokamak_unstructured_mesh_simulation.py new file mode 100644 index 00000000..863de80f --- /dev/null +++ b/examples/tokamak_unstructured_mesh_simulation.py @@ -0,0 +1,166 @@ +""" +This script makes a Tokamak model and performs a neutronics simulation +using OpenMC. The model is created using the paramak package and then converted +using cad_to_dagmc package. The model is a 180 degree model so we add a wedge +shaped bounding region with reflective surfaces to make it a full torus. +This example aims to demonstrate how to use cad_to_dagmc to make a surface DAGMC +mesh and an unstructured mesh to use in an unstructured mesh tally. The script +uses minimal materials, tallies, source to keep the example concise. +""" + + +import paramak +import openmc +from cad_to_dagmc import CadToDagmc +from pathlib import Path + + +openmc.config["cross_sections"] = "/nuclear_data/cross_sections.xml" + +my_reactor = paramak.tokamak( + radial_build=[ + (paramak.LayerType.GAP, 10), # inner bore + (paramak.LayerType.SOLID, 30), + (paramak.LayerType.SOLID, 50), + (paramak.LayerType.SOLID, 10), + (paramak.LayerType.SOLID, 120), # breeder + (paramak.LayerType.SOLID, 10), # first wall + (paramak.LayerType.GAP, 60), + (paramak.LayerType.PLASMA, 300), # plasma + (paramak.LayerType.GAP, 60), + (paramak.LayerType.SOLID, 10), # first wall + (paramak.LayerType.SOLID, 120), # breeder + (paramak.LayerType.SOLID, 10), + ], + vertical_build=[ + (paramak.LayerType.SOLID, 15), + (paramak.LayerType.SOLID, 100), + (paramak.LayerType.SOLID, 10), + (paramak.LayerType.GAP, 50), + (paramak.LayerType.PLASMA, 700), # plasma + (paramak.LayerType.GAP, 60), + (paramak.LayerType.SOLID, 10), + (paramak.LayerType.SOLID, 100), + (paramak.LayerType.SOLID, 15), + ], + triangularity=0.55, + rotation_angle=180, +) + +my_reactor = my_reactor.remove(name="plasma") # removing as we don't need the plasma for this neutronics simulation + +my_model = CadToDagmc() +# as inner and outer layers are one solid there are only 6 solids in model +material_tags = ["layer_1", "layer_2", "layer_3", "layer_4", "layer_5"] +my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) +# my_model.export_dagmc_h5m_file(min_mesh_size=10.0, max_mesh_size=20.0, filename="dagmc.h5m") +# my_model.export_unstructured_mesh_file(min_mesh_size=10.0, max_mesh_size=20.0, filename="unstructured_mesh.vtk") + + +script_folder = Path(__file__).resolve().parent +h5m_filename = script_folder / "dagmc.h5m" + + +mat_layer_1 = openmc.Material(name="layer_1") +mat_layer_1.add_element("Cu", 1, "ao") +mat_layer_1.set_density("g/cm3", 7) + +mat_layer_2 = openmc.Material(name="layer_2") +mat_layer_2.add_nuclide("W186", 1, "ao") +mat_layer_2.set_density("g/cm3", 0.01) + +mat_layer_3 = openmc.Material(name="layer_3") +mat_layer_3.add_nuclide("Fe56", 1, "ao") +mat_layer_3.set_density("g/cm3", 7) + +mat_layer_4 = openmc.Material(name="layer_4") +mat_layer_4.add_element("Li", 1, "ao") +mat_layer_4.set_density("g/cm3", 0.5) + +mat_layer_5 = openmc.Material(name="layer_5") +mat_layer_5.add_nuclide("Fe56", 1, "ao") +mat_layer_5.set_density("g/cm3", 7) + + +materials = openmc.Materials([mat_layer_1, mat_layer_2, mat_layer_3, mat_layer_4, mat_layer_5]) + + +dag_univ = openmc.DAGMCUniverse(filename=h5m_filename) + +# gets the bounding box as we need this to create the correctly sized bounding cell. +bbox = dag_univ.bounding_box +dagmc_radius = max(abs(bbox[0][0]), abs(bbox[0][1]), abs(bbox[1][0]), abs(bbox[1][1])) + +cylinder_surface = openmc.ZCylinder(r=dagmc_radius, boundary_type="vacuum", surface_id=1000) +lower_z = openmc.ZPlane(bbox[0][2], boundary_type="vacuum", surface_id=1003) +upper_z = openmc.ZPlane(bbox[1][2], boundary_type="vacuum", surface_id=1004) + +# this is the surface along the side of the 180 degree model +side_surface = openmc.YPlane(y0=0, boundary_type="reflective", surface_id=1001) + +wedge_region = -cylinder_surface & +lower_z & -upper_z & +side_surface + +# bounding cell is a wedge shape filled with the DAGMC universe +bounding_cell = openmc.Cell(fill=dag_univ, cell_id=1000, region=wedge_region) + +# bound_dag_univ = dag_univ.bounded_universe() +geometry = openmc.Geometry([bounding_cell]) + +# initializes a new source object +my_source = openmc.IndependentSource() + +# the distribution of radius is just a single value +radius = openmc.stats.Discrete([my_reactor.major_radius], [1]) + +# the distribution of source z values is just a single value +z_values = openmc.stats.Discrete([0], [1]) + +# the distribution of source azimuthal angles values is a uniform distribution between 0 and 2 Pi +angle = openmc.stats.Uniform(a=0.0, b=2 * 3.14159265359) + +# this makes the ring source using the three distributions and a radius +my_source.space = openmc.stats.CylindricalIndependent(r=radius, phi=angle, z=z_values, origin=(0.0, 0.0, 0.0)) + +# sets the direction to isotropic +my_source.angle = openmc.stats.Isotropic() + +# sets the energy distribution to a Muir distribution neutrons +my_source.energy = openmc.stats.muir(e0=14080000.0, m_rat=5.0, kt=20000.0) + + +# specifies the simulation computational intensity +settings = openmc.Settings(batches=10, particles=10000, run_mode="fixed source", source=my_source) + +mesh = openmc.UnstructuredMesh(filename="unstructured_mesh.vtk", library="moab", mesh_id=1) + +# adds a tally to record the heat deposited in entire geometry +mesh_tally = openmc.Tally(name="flux") +mesh_tally.filters = [openmc.MeshFilter(mesh)] +mesh_tally.scores = ["flux"] + +# groups the two tallies +tallies = openmc.Tallies([mesh_tally]) + +# builds the openmc model +my_model = openmc.Model(materials=materials, geometry=geometry, settings=settings, tallies=tallies) + +# starts the simulation +output_file = my_model.run() + +# loads up the output file from the simulation +statepoint = openmc.StatePoint(output_file) + +mesh_tally_result = statepoint.get_tally(name="flux") + +umesh_from_sp = mesh_tally_result.find_filter(openmc.MeshFilter).mesh + +centroids = umesh_from_sp.centroids # not needed in the next release of openmc +mesh_vols = umesh_from_sp.volumes # not needed in the next release of openmc + +# writes the tally data to a VTK file for visualisation +umesh_from_sp.write_data_to_vtk( + datasets={"mean": mesh_tally_result.mean.flatten()}, + filename="unstructured_mesh_tally_results.vtk", +) + +print("VTK file saved to unstructured_mesh_tally_results.vtk") diff --git a/examples/tokamak_with_pf_magnets.py b/examples/tokamak_with_pf_magnets.py index 5e7b15f3..8aa4ed47 100644 --- a/examples/tokamak_with_pf_magnets.py +++ b/examples/tokamak_with_pf_magnets.py @@ -1,5 +1,3 @@ -from example_util_functions import transport_particles_on_h5m_geometry - import paramak extra_cut_shapes = [] @@ -53,17 +51,3 @@ my_reactor.save(f"tokamak_minimal.step") print(f"Saved as tokamak_minimal.step") -# from cad_to_dagmc import CadToDagmc -# my_model = CadToDagmc() -# material_tags = ["mat1"] * 6 # as inner and outer layers are one solid there are only 6 solids in model -# my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) -# my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) - -# h5m_filename = "dagmc.h5m" -# flux = transport_particles_on_h5m_geometry( -# h5m_filename=h5m_filename, -# material_tags=material_tags, -# nuclides=["H1"] * len(material_tags), -# cross_sections_xml="tests/cross_sections.xml", -# ) -# assert flux > 0.0 diff --git a/examples/tokamak_with_pf_tf_magnets_divertor.py b/examples/tokamak_with_pf_tf_magnets_divertor.py index 332dbbeb..96e50a9d 100644 --- a/examples/tokamak_with_pf_tf_magnets_divertor.py +++ b/examples/tokamak_with_pf_tf_magnets_divertor.py @@ -1,5 +1,3 @@ -from example_util_functions import transport_particles_on_h5m_geometry - import paramak import cadquery as cq @@ -75,21 +73,3 @@ ) 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 -# my_model = CadToDagmc() -# material_tags = ["mat1"] * 6 # as inner and outer layers are one solid there are only 6 solids in model -# my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) -# my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0) - -# h5m_filename = "dagmc.h5m" -# flux = transport_particles_on_h5m_geometry( -# h5m_filename=h5m_filename, -# material_tags=material_tags, -# nuclides=["H1"] * len(material_tags), -# cross_sections_xml="tests/cross_sections.xml", -# ) -# assert flux > 0.0