diff --git a/README.md b/README.md index 67de6bc..c0bf169 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,11 @@ # paramak-neutronics -Adds support for neutronics simulations to the paramak package + +Adds support for neutronics simulations to the paramak package. + +Install with Pip + +```bash +pip install paramak-neutronics +``` + +For examples see the [examples folder](https://github.com/fusion-energy/paramak-neutronics/tree/main/examples) diff --git a/examples/ball_reactor.ipynb b/examples/ball_reactor.ipynb new file mode 100644 index 0000000..be32cff --- /dev/null +++ b/examples/ball_reactor.ipynb @@ -0,0 +1,108 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This example makes a reactor geometry and a neutronics model. A homogenised material made of enriched lithium lead and eurofer is being used as the blanket material for this simulation in order to demonstrate the use of more complex materials." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import neutronics_material_maker as nmm\n", + "import openmc\n", + "import paramak\n", + "\n", + "# makes the 3d geometry\n", + "my_reactor = paramak.BallReactor(\n", + " inner_bore_radial_thickness=1,\n", + " inboard_tf_leg_radial_thickness=30,\n", + " center_column_shield_radial_thickness=60,\n", + " divertor_radial_thickness=50,\n", + " inner_plasma_gap_radial_thickness=30,\n", + " plasma_radial_thickness=300,\n", + " outer_plasma_gap_radial_thickness=30,\n", + " firstwall_radial_thickness=3,\n", + " blanket_radial_thickness=100,\n", + " blanket_rear_wall_radial_thickness=3,\n", + " elongation=2.75,\n", + " triangularity=0.5,\n", + " number_of_tf_coils=16,\n", + " rotation_angle=359.9, # when using trelis method this can be set to 360\n", + ")\n", + "\n", + "# method is set to Trelis or Cubit by default to avoid overlaps in the geometry\n", + "# pymoab is used as it is open source and can be tested in the CI\n", + "# if you have Trelis or Cubit then this line can be deleted\n", + "my_reactor.method='pymoab'\n", + "\n", + "# makes a homogenised material for the blanket from lithium lead and\n", + "# eurofer\n", + "blanket_material = nmm.Material.from_mixture(\n", + " fracs=[0.8, 0.2],\n", + " materials=[\n", + " nmm.Material.from_library(\n", + " name='Pb842Li158',\n", + " enrichment=90,\n", + " temperature=500),\n", + " nmm.Material.from_library(name='eurofer')\n", + " ])\n", + "\n", + "source = openmc.Source()\n", + "# sets the location of the source to x=0 y=0 z=0\n", + "source.space = openmc.stats.Point((my_reactor.major_radius, 0, 0))\n", + "# sets the direction to isotropic\n", + "source.angle = openmc.stats.Isotropic()\n", + "# sets the energy distribution to 100% 14MeV neutrons\n", + "source.energy = openmc.stats.Discrete([14e6], [1])\n", + "\n", + "# makes the neutronics material\n", + "neutronics_model = paramak.NeutronicsModel(\n", + " geometry=my_reactor,\n", + " source=source,\n", + " materials={\n", + " 'inboard_tf_coils_mat': 'copper',\n", + " 'center_column_shield_mat': 'WC',\n", + " 'divertor_mat': 'eurofer',\n", + " 'firstwall_mat': 'eurofer',\n", + " 'blanket_mat': blanket_material, # use of homogenised material\n", + " 'blanket_rear_wall_mat': 'eurofer'},\n", + " cell_tallies=['TBR'],\n", + " simulation_batches=2,\n", + " simulation_particles_per_batch=10, # this will need increasing to obtain accurate results\n", + ")\n", + "\n", + "# starts the neutronics simulation\n", + "neutronics_model.simulate()\n", + "\n", + "# prints the results to screen\n", + "print('TBR', neutronics_model.results['TBR'])" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/ball_reactor_minimal.ipynb b/examples/ball_reactor_minimal.ipynb new file mode 100644 index 0000000..e3da518 --- /dev/null +++ b/examples/ball_reactor_minimal.ipynb @@ -0,0 +1,101 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This is a minimal example that obtains the TBR (Tritium Breeding Ratio) for a parametric ball reactor" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import openmc\n", + "import paramak\n", + "\n", + "\n", + "# makes the 3d geometry from input parameters\n", + "my_reactor = paramak.BallReactor(\n", + " inner_bore_radial_thickness=50,\n", + " inboard_tf_leg_radial_thickness=200,\n", + " center_column_shield_radial_thickness=50,\n", + " divertor_radial_thickness=50,\n", + " inner_plasma_gap_radial_thickness=50,\n", + " plasma_radial_thickness=100,\n", + " outer_plasma_gap_radial_thickness=50,\n", + " firstwall_radial_thickness=1,\n", + " blanket_radial_thickness=100,\n", + " blanket_rear_wall_radial_thickness=10,\n", + " elongation=2,\n", + " triangularity=0.55,\n", + " number_of_tf_coils=16,\n", + " rotation_angle=359.9, # when using trelis method this can be set to 360\n", + ")\n", + "\n", + "# method is set to Trelis or Cubit by default to avoid overlaps in the geometry\n", + "# pymoab is used as it is open source and can be tested in the CI\n", + "# if you have Trelis or Cubit then this line can be deleted\n", + "my_reactor.method='pymoab'\n", + "\n", + "source = openmc.Source()\n", + "# sets the location of the source to x=0 y=0 z=0\n", + "source.space = openmc.stats.Point((my_reactor.major_radius, 0, 0))\n", + "# sets the direction to isotropic\n", + "source.angle = openmc.stats.Isotropic()\n", + "# sets the energy distribution to 100% 14MeV neutrons\n", + "source.energy = openmc.stats.Discrete([14e6], [1])\n", + "\n", + "# makes the neutronics model from the geometry and material allocations\n", + "neutronics_model = paramak.NeutronicsModel(\n", + " geometry=my_reactor,\n", + " source=source,\n", + " materials={\n", + " 'inboard_tf_coils_mat': 'eurofer',\n", + " 'center_column_shield_mat': 'eurofer',\n", + " 'divertor_mat': 'eurofer',\n", + " 'firstwall_mat': 'eurofer',\n", + " 'blanket_rear_wall_mat': 'eurofer',\n", + " 'blanket_mat': 'Li4SiO4'},\n", + " cell_tallies=['TBR', 'heating'],\n", + " simulation_batches=2,\n", + " simulation_particles_per_batch=10, # this will need increasing to obtain accurate results\n", + ")\n", + "\n", + "# simulate the neutronics model\n", + "neutronics_model.simulate()\n", + "print(neutronics_model.results)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/ball_reactor_source_plot.ipynb b/examples/ball_reactor_source_plot.ipynb new file mode 100644 index 0000000..0226baf --- /dev/null +++ b/examples/ball_reactor_source_plot.ipynb @@ -0,0 +1,153 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This example makes a reactor geometry and a neutronics model. A homogenised material made of enriched lithium lead and eurofer is being used as the blanket material for this simulation in order to demonstrate the use of more complex materials." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import neutronics_material_maker as nmm\n", + "import openmc\n", + "import paramak\n", + "\n", + "\n", + "# makes the 3d geometry\n", + "my_reactor = paramak.BallReactor(\n", + " inner_bore_radial_thickness=10,\n", + " inboard_tf_leg_radial_thickness=30,\n", + " center_column_shield_radial_thickness=60,\n", + " divertor_radial_thickness=150,\n", + " inner_plasma_gap_radial_thickness=30,\n", + " plasma_radial_thickness=300,\n", + " outer_plasma_gap_radial_thickness=30,\n", + " firstwall_radial_thickness=30,\n", + " blanket_radial_thickness=50,\n", + " blanket_rear_wall_radial_thickness=30,\n", + " elongation=2,\n", + " triangularity=0.55,\n", + " outboard_tf_coil_radial_thickness=100,\n", + " outboard_tf_coil_poloidal_thickness=50,\n", + " rotation_angle=90,\n", + ")\n", + "\n", + "# method is set to Trelis or Cubit by default to avoid overlaps in the geometry\n", + "# pymoab is used as it is open source and can be tested in the CI\n", + "# if you have Trelis or Cubit then this line can be deleted\n", + "my_reactor.method='pymoab'\n", + "\n", + "my_reactor.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "source = openmc.Source()\n", + "radius = openmc.stats.Discrete([300], [1]) # at the major radius\n", + "z_values = openmc.stats.Discrete([0], [1]) # middle of the plasma (vertically)\n", + "angle = openmc.stats.Uniform(a=0., b=2 * 3.14159265359) # 360 degrees\n", + "source.space = openmc.stats.CylindricalIndependent(\n", + " r=radius,\n", + " phi=angle,\n", + " z=z_values,\n", + " origin=(0.0, 0.0, 0.0)\n", + ")\n", + "source.angle = openmc.stats.Isotropic()\n", + "source.energy = openmc.stats.Discrete([14_060_000], [1])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# makes the neutronics material\n", + "neutronics_model = paramak.NeutronicsModel(\n", + " geometry=my_reactor,\n", + " source=source,\n", + " materials={\n", + " 'inboard_tf_coils_mat': 'copper',\n", + " 'center_column_shield_mat': 'WC',\n", + " 'divertor_mat': 'eurofer',\n", + " 'firstwall_mat': 'eurofer',\n", + " 'blanket_mat': 'eurofer',\n", + " 'blanket_rear_wall_mat': 'eurofer'},\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# exports the geometry and source in 2d (RZ) viewplane where R stands for\n", + "# radius\n", + "neutronics_model.export_html(\n", + " filename='2d_source.html',\n", + " view_plane='RZ',\n", + " number_of_source_particles=100\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# exports the geometry and source in 3d (XYZ) viewplane\n", + "neutronics_model.export_html(\n", + " filename='3d_source.html',\n", + " view_plane='XYZ',\n", + " number_of_source_particles=100\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/center_column_study_reactor.ipynb b/examples/center_column_study_reactor.ipynb new file mode 100644 index 0000000..d5b3af0 --- /dev/null +++ b/examples/center_column_study_reactor.ipynb @@ -0,0 +1,113 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\"\"\"This example makes a reactor geometry and a neutronics model, the addition\n", + "of a Python for loop allow a parameter sweep of the neutronics results for\n", + "different geometries. The distance between the plasma and center column is\n", + "varried while simulating the impact on the heat depositied in the center\n", + "column.\"\"\"\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import openmc\n", + "import paramak\n", + "\n", + "\n", + "total_heats_in_MW = []\n", + "plasma_to_center_column_gaps = []\n", + "\n", + "# this will take a few mins to perform 3 simulations at\n", + "for plasma_to_center_column_gap in [50, 100, 150]:\n", + "\n", + " # makes the 3d geometry\n", + " my_reactor = paramak.CenterColumnStudyReactor(\n", + " inner_bore_radial_thickness=20,\n", + " inboard_tf_leg_radial_thickness=50,\n", + " center_column_shield_radial_thickness_mid=50,\n", + " center_column_shield_radial_thickness_upper=100,\n", + " inboard_firstwall_radial_thickness=2,\n", + " divertor_radial_thickness=100,\n", + " inner_plasma_gap_radial_thickness=plasma_to_center_column_gap,\n", + " plasma_radial_thickness=200,\n", + " outer_plasma_gap_radial_thickness=90,\n", + " elongation=2.75,\n", + " triangularity=0.5,\n", + " plasma_gap_vertical_thickness=40,\n", + " center_column_arc_vertical_thickness=520,\n", + " rotation_angle=360\n", + " )\n", + " \n", + " # method is set to Trelis or Cubit by default to avoid overlaps in the geometry\n", + " # pymoab is used as it is open source and can be tested in the CI\n", + " # if you have Trelis or Cubit then this line can be deleted\n", + " my_reactor.method='pymoab'\n", + "\n", + " source = openmc.Source()\n", + " # sets the location of the source to x=0 y=0 z=0\n", + " source.space = openmc.stats.Point((my_reactor.major_radius, 0, 0))\n", + " # sets the direction to isotropic\n", + " source.angle = openmc.stats.Isotropic()\n", + " # sets the energy distribution to 100% 14MeV neutrons\n", + " source.energy = openmc.stats.Discrete([14e6], [1])\n", + "\n", + " # makes the neutronics model and assigns basic materials to each\n", + " # component\n", + " neutronics_model = paramak.NeutronicsModel(\n", + " geometry=my_reactor,\n", + " source=source,\n", + " materials={\n", + " 'DT_plasma': 'DT_plasma',\n", + " 'inboard_tf_coils_mat': 'eurofer',\n", + " 'center_column_shield_mat': 'eurofer',\n", + " 'divertor_mat': 'eurofer',\n", + " 'firstwall_mat': 'eurofer',\n", + " 'blanket_mat': 'Li4SiO4'},\n", + " cell_tallies=['heating'],\n", + " simulation_batches=2,\n", + " simulation_particles_per_batch=10, # this will need increasing to obtain accurate results\n", + " )\n", + "\n", + " # starts the neutronics simulation\n", + " neutronics_model.simulate()\n", + "\n", + " # converts the results to mega watts\n", + " total_heat_in_MW = neutronics_model.results['firstwall_mat_heating']['Watts']['result'] / 1e6\n", + "\n", + " # adds the results and inputs to a list\n", + " total_heats_in_MW.append(total_heat_in_MW)\n", + " plasma_to_center_column_gaps.append(plasma_to_center_column_gap)\n", + "\n", + "# plots the results\n", + "plt.scatter(plasma_to_center_column_gaps, total_heats_in_MW)\n", + "plt.xlabel('plasma_to_center_column_gap (cm)')\n", + "plt.ylabel('Heat on the inboard (MW)')\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/center_column_study_reactor_minimal.ipynb b/examples/center_column_study_reactor_minimal.ipynb new file mode 100644 index 0000000..ac5a14c --- /dev/null +++ b/examples/center_column_study_reactor_minimal.ipynb @@ -0,0 +1,98 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\"\"\"This is a minimal example that obtains the center column heating for a\n", + "parametric reactor.\"\"\"\n", + "\n", + "import openmc\n", + "import paramak\n", + "\n", + "\n", + "# makes the 3d geometry\n", + "my_reactor = paramak.CenterColumnStudyReactor(\n", + " inner_bore_radial_thickness=20,\n", + " inboard_tf_leg_radial_thickness=50,\n", + " center_column_shield_radial_thickness_mid=50,\n", + " center_column_shield_radial_thickness_upper=100,\n", + " inboard_firstwall_radial_thickness=20,\n", + " divertor_radial_thickness=100,\n", + " inner_plasma_gap_radial_thickness=80,\n", + " plasma_radial_thickness=200,\n", + " outer_plasma_gap_radial_thickness=90,\n", + " elongation=2.75,\n", + " triangularity=0.5,\n", + " plasma_gap_vertical_thickness=40,\n", + " center_column_arc_vertical_thickness=520,\n", + " rotation_angle=360\n", + ")\n", + "\n", + "# method is set to Trelis or Cubit by default to avoid overlaps in the geometry\n", + "# pymoab is used as it is open source and can be tested in the CI\n", + "# if you have Trelis or Cubit then this line can be deleted\n", + "my_reactor.method='pymoab'\n", + "\n", + "source = openmc.Source()\n", + "# sets the location of the source to x=0 y=0 z=0\n", + "source.space = openmc.stats.Point((my_reactor.major_radius, 0, 0))\n", + "# sets the direction to isotropic\n", + "source.angle = openmc.stats.Isotropic()\n", + "# sets the energy distribution to 100% 14MeV neutrons\n", + "source.energy = openmc.stats.Discrete([14e6], [1])\n", + "\n", + "# creates a neutronics model from the geometry and assigned materials\n", + "neutronics_model = paramak.NeutronicsModel(\n", + " geometry=my_reactor,\n", + " source=source,\n", + " materials={\n", + " 'inboard_tf_coils_mat': 'eurofer',\n", + " 'center_column_shield_mat': 'eurofer',\n", + " 'divertor_mat': 'eurofer',\n", + " 'firstwall_mat': 'eurofer',\n", + " 'blanket_mat': 'Li4SiO4'},\n", + " cell_tallies=['heating'],\n", + " simulation_batches=2,\n", + " simulation_particles_per_batch=10, # this will need increasing to obtain accurate results\n", + ")\n", + "\n", + "# starts the neutronics simulation\n", + "neutronics_model.simulate()\n", + "\n", + "# prints the results\n", + "print(neutronics_model.results)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/component_based_mesh_simulation.py b/examples/component_based_mesh_simulation.py new file mode 100644 index 0000000..b90ea16 --- /dev/null +++ b/examples/component_based_mesh_simulation.py @@ -0,0 +1,47 @@ +"""This example creates a curved center column with a few different sizes. +The shape is then converted into a neutronics geometry and the heat deposited +is simulated for a few different sizes of ceter column""" + +import matplotlib.pyplot as plt +import openmc +import paramak + + +def main(): + + # makes the component with a few different size mid radius values + my_shape = paramak.CenterColumnShieldHyperbola( + height=500, + inner_radius=50, + mid_radius=70, + outer_radius=100, + material_tag='center_column_shield_mat', + method='pymoab', + ) + + my_shape.export_stp('my_shape.stp') + + # makes the openmc neutron source at x,y,z 0, 0, 0 with isotropic + # diections + source = openmc.Source() + source.space = openmc.stats.Point((0, 0, 0)) + source.angle = openmc.stats.Isotropic() + + # converts the geometry into a neutronics geometry + my_model = paramak.NeutronicsModel( + geometry=my_shape, + source=source, + materials={'center_column_shield_mat': 'WB'}, # WB is tungsten boride + cell_tallies=['heating', 'TBR'], + mesh_tally_2d=['heating'], + mesh_tally_3d=['heating'], + simulation_batches=10, # should be increased for more accurate result + simulation_particles_per_batch=10, # settings are low to reduce time required + ) + + # performs an openmc simulation on the model + my_model.simulate() + + +if __name__ == "__main__": + main() diff --git a/examples/component_based_parameter_study.ipynb b/examples/component_based_parameter_study.ipynb new file mode 100644 index 0000000..0e35c0b --- /dev/null +++ b/examples/component_based_parameter_study.ipynb @@ -0,0 +1,110 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This example creates a curved center column with a few different sizes.\n", + "The shape is then converted into a neutronics geometry and the heat deposited\n", + "is simulated for a few different sizes of ceter column." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import openmc\n", + "import paramak\n", + "\n", + "\n", + "simulation_values = []\n", + "for mid_radius in [60, 70, 80]:\n", + "\n", + " # makes the component with a few different size mid radius values\n", + " my_shape = paramak.CenterColumnShieldHyperbola(\n", + " height=500,\n", + " inner_radius=50,\n", + " mid_radius=mid_radius,\n", + " outer_radius=100,\n", + " material_tag='center_column_shield_mat',\n", + " method='pymoab',\n", + " )\n", + "\n", + " my_shape.export_stp('my_shape' + str(mid_radius) + '.stp')\n", + "\n", + " # makes the openmc neutron source at x,y,z 0, 0, 0 with isotropic\n", + " # diections\n", + " source = openmc.Source()\n", + " source.space = openmc.stats.Point((0, 0, 0))\n", + " source.angle = openmc.stats.Isotropic()\n", + "\n", + " # converts the geometry into a neutronics geometry\n", + " my_model = paramak.NeutronicsModel(\n", + " geometry=my_shape,\n", + " source=source,\n", + " materials={\n", + " 'center_column_shield_mat': 'WB'},\n", + " # WB is tungsten boride\n", + " cell_tallies=['heating'],\n", + " simulation_batches=10, # should be increased for more accurate result\n", + " simulation_particles_per_batch=10, # settings are low to reduce time required\n", + " )\n", + "\n", + " # performs an openmc simulation on the model\n", + " my_model.simulate()\n", + "\n", + " # extracts the heat from the results dictionary\n", + " heat = my_model.results['center_column_shield_mat_heating']['Watts']['result']\n", + "\n", + " # adds the heat and the mid radius value to a list\n", + " simulation_values.append((mid_radius, heat))\n", + "\n", + "# plots the simualtion results vs the mid_radius used for the simulation\n", + "plt.plot(\n", + " [i[0] for i in simulation_values],\n", + " [i[1] for i in simulation_values],\n", + " '-p'\n", + ")\n", + "\n", + "# adds labels to the graph\n", + "plt.title(\"heating vs thickness\")\n", + "plt.xlabel(\"thickness (cm)\")\n", + "plt.ylabel(\"heating (watts)\")\n", + "\n", + "plt.savefig('heating_vs_thickness.svg')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/dagmc_logo.ipynb b/examples/dagmc_logo.ipynb new file mode 100644 index 0000000..ba56000 --- /dev/null +++ b/examples/dagmc_logo.ipynb @@ -0,0 +1,407 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import paramak" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Overwriting auto display for cadquery Workplane and Shape\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "ef616b44bad94747a1dc335b245c0e41", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "HBox(children=(VBox(children=(HBox(children=(Checkbox(value=False, description='Axes', indent=False, _dom_clas…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "white_part_1 = paramak.CenterColumnShieldCylinder(\n", + " inner_radius=0,\n", + " outer_radius=50,\n", + " height=10,\n", + " color=(1,1,1)\n", + ")\n", + "\n", + "white_part_1.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "ef885bb494b24c93968c1094df163c73", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "HBox(children=(VBox(children=(HBox(children=(Checkbox(value=False, description='Axes', indent=False, _dom_clas…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "blue_part_1 = paramak.CenterColumnShieldCylinder(\n", + " inner_radius=50,\n", + " outer_radius=50+82,\n", + " height=10,\n", + " rotation_angle=360*(2/3),\n", + " color=(0.09,0.212,0.365)\n", + ")\n", + "\n", + "blue_part_1.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "ec067ff0f2f74ff996843396ee012fb7", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "HBox(children=(VBox(children=(HBox(children=(Checkbox(value=False, description='Axes', indent=False, _dom_clas…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "red_part_1 = paramak.CenterColumnShieldCylinder(\n", + " inner_radius=50+82,\n", + " outer_radius=50+82+140,\n", + " height=10,\n", + " rotation_angle=360*(2/3),\n", + " color=(0.753, 0, 0)\n", + ")\n", + "\n", + "red_part_1.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "76e1b4207b6a4067bcc7329fcdd2bcc9", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "HBox(children=(VBox(children=(HBox(children=(Checkbox(value=False, description='Axes', indent=False, _dom_clas…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "yellow_part_1 = paramak.CenterColumnShieldCylinder(\n", + " inner_radius=50,\n", + " outer_radius=50+82+140,\n", + " height=10,\n", + " rotation_angle=360*(1/3),\n", + " azimuth_placement_angle=360-(360*(1/3)),\n", + " color=(0.463, 0.573, 0.235)\n", + ")\n", + "\n", + "yellow_part_1.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "c5db0af75cb14018a31d3dc2f0e4a20a", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "HBox(children=(VBox(children=(HBox(children=(Checkbox(value=False, description='Axes', indent=False, _dom_clas…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "white_part_2 = paramak.CenterColumnShieldCylinder(\n", + " inner_radius=50+82+140,\n", + " outer_radius=50+82+140+64,\n", + " height=10,\n", + " color=(1,1,1)\n", + ")\n", + "\n", + "white_part_2.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "f9fe455134504df4a3534d51e3c844e2", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "HBox(children=(VBox(children=(HBox(children=(Checkbox(value=False, description='Axes', indent=False, _dom_clas…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "red_part_2 = paramak.CenterColumnShieldCylinder(\n", + " inner_radius=50+82+140+64,\n", + " outer_radius=50+82+140+64+111,\n", + " height=10,\n", + " rotation_angle=360*(1/4),\n", + " azimuth_placement_angle=170,\n", + " color=(1,0,0)\n", + ")\n", + "\n", + "red_part_2.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "be38cb5d079544ba982ae1006ac388c7", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "HBox(children=(VBox(children=(HBox(children=(Checkbox(value=False, description='Axes', indent=False, _dom_clas…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "yellow_part_2 = paramak.CenterColumnShieldCylinder(\n", + " inner_radius=50+82+140+64,\n", + " outer_radius=50+82+140+64+111,\n", + " height=10,\n", + " rotation_angle=360*(1/4),\n", + " azimuth_placement_angle=0,\n", + " color=(0.463, 0.573, 0.235)\n", + ")\n", + "\n", + "yellow_part_2.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "38ee6893f28d4a69a04bfbf47e427bee", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "HBox(children=(VBox(children=(HBox(children=(Checkbox(value=False, description='Axes', indent=False, _dom_clas…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "white_part_3 = paramak.CenterColumnShieldCylinder(\n", + " inner_radius=50+82+140+64,\n", + " outer_radius=50+82+140+64+111,\n", + " height=10,\n", + " rotation_angle=360*(1/4),\n", + " azimuth_placement_angle=80,\n", + " color=(1,1,1),\n", + " name='white_part_3'\n", + ")\n", + "\n", + "white_part_3.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "2fd1caf0e87a4025bd8a17d31754b73b", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "HBox(children=(VBox(children=(HBox(children=(Checkbox(value=False, description='Axes', indent=False, _dom_clas…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "blue_part_2 = paramak.CenterColumnShieldCylinder(\n", + " inner_radius=50+82+140+64,\n", + " outer_radius=50+82+140+64+111,\n", + " height=10,\n", + " rotation_angle=360*(1/4),\n", + " azimuth_placement_angle=-10,\n", + " color=(0.09,0.212,0.365),\n", + " name='blue_part_2'\n", + ")\n", + "\n", + "blue_part_2.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "5bd23ac314e0417baee8228d4c5322e0", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "HBox(children=(VBox(children=(HBox(children=(Checkbox(value=False, description='Axes', indent=False, _dom_clas…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "yellow_part_2 = paramak.CenterColumnShieldCylinder(\n", + " inner_radius=50+82+140+64,\n", + " outer_radius=50+82+140+64+111,\n", + " height=10,\n", + " rotation_angle=360*(1/4),\n", + " azimuth_placement_angle=260,\n", + " color=(0.463, 0.573, 0.235),\n", + " name='yellow_part_2'\n", + ")\n", + "\n", + "yellow_part_2.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "8363e794e6204fb39f95c0ffa32bb69b", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "HBox(children=(VBox(children=(HBox(children=(Checkbox(value=False, description='Axes', indent=False, _dom_clas…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "all_parts = paramak.Reactor([yellow_part_1, yellow_part_2, white_part_1, white_part_2, white_part_3, red_part_1, red_part_2, blue_part_1, blue_part_2])\n", + "all_parts.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/heating_on_unstructured_mesh_shape.py b/examples/heating_on_unstructured_mesh_shape.py new file mode 100644 index 0000000..85fe8a3 --- /dev/null +++ b/examples/heating_on_unstructured_mesh_shape.py @@ -0,0 +1,23 @@ + +import paramak + + +rotated_spline = paramak.RotateSplineShape( + rotation_angle=180, + method='trelis', + points=[ + (500, 0), + (500, -20), + (400, -300), + (300, -300), + (400, 0), + (300, 300), + (400, 300), + (500, 20), + ], + tet_mesh=) + +rotated_spline.export_h5m( + # merge_tolerance, + # faceting_tolerance=, +) diff --git a/examples/make_unstructured_mesh.py b/examples/make_unstructured_mesh.py new file mode 100644 index 0000000..2f72faa --- /dev/null +++ b/examples/make_unstructured_mesh.py @@ -0,0 +1,139 @@ +#!/usr/env/python3 +import json +import os + + +def byteify(input): + if isinstance(input, dict): + return {byteify(key): byteify(value) + for key, value in input.iteritems()} + elif isinstance(input, list): + return [byteify(element) for element in input] + elif isinstance(input, unicode): + return input.encode("utf-8") + else: + return input + + +def save_tet_details_to_json_file( + geometry_details, + filename="mesh_details.json"): + for entry in geometry_details: + material = entry["material"] + tets_in_volumes = cubit.parse_cubit_list( + "tet", " in volume " + " ".join(entry["volumes"]) + ) + print("material ", material, " has ", len(tets_in_volumes), " tets") + entry["tet_ids"] = tets_in_volumes + with open(filename, "w") as outfile: + json.dump(geometry_details, outfile, indent=4) + + +def find_number_of_volumes_in_each_step_file(input_locations, basefolder): + body_ids = "" + volumes_in_each_step_file = [] + # all_groups=cubit.parse_cubit_list("group","all") + # starting_group_id = len(all_groups) + for entry in input_locations: + # starting_group_id = starting_group_id +1 + current_vols = cubit.parse_cubit_list("volume", "all") + print(os.path.join(basefolder, entry["stp_filename"])) + if entry["stp_filename"].endswith(".sat"): + import_type = "acis" + if entry["stp_filename"].endswith( + ".stp") or entry["stp_filename"].endswith(".step"): + import_type = "step" + short_file_name = os.path.split(entry["stp_filename"])[-1] + # print('short_file_name',short_file_name) + # cubit.cmd('import '+import_type+' "' + entry['stp_filename'] + '" separate_bodies no_surfaces no_curves no_vertices group "'+str(short_file_name)+'"') + cubit.cmd( + "import " + + import_type + + ' "' + + os.path.join(basefolder, entry["stp_filename"]) + + '" separate_bodies no_surfaces no_curves no_vertices ' + ) + all_vols = cubit.parse_cubit_list("volume", "all") + new_vols = set(current_vols).symmetric_difference(set(all_vols)) + new_vols = map(str, new_vols) + print("new_vols", new_vols, type(new_vols)) + current_bodies = cubit.parse_cubit_list("body", "all") + print("current_bodies", current_bodies) + # volumes_in_group = cubit.cmd('volume in group '+str(starting_group_id)) + # print('volumes_in_group',volumes_in_group,type(volumes_in_group)) + if len(new_vols) > 1: + cubit.cmd( + "unite vol " + + " ".join(new_vols) + + " with vol " + + " ".join(new_vols)) + all_vols = cubit.parse_cubit_list("volume", "all") + new_vols_after_unite = set( + current_vols).symmetric_difference(set(all_vols)) + new_vols_after_unite = map(str, new_vols_after_unite) + # cubit.cmd('group '+str(starting_group_id)+' copy rotate 45 about z repeat 7') + entry["volumes"] = new_vols_after_unite + cubit.cmd( + 'group "' + + short_file_name + + '" add volume ' + + " ".join( + entry["volumes"])) + # cubit.cmd('volume in group '+str(starting_group_id)+' copy rotate 45 about z repeat 7') + cubit.cmd("separate body all") + return input_locations + + +def imprint_and_merge_geometry(tolerance="1e-4"): + cubit.cmd("imprint body all") + cubit.cmd("merge tolerance " + tolerance) # optional as there is a default + cubit.cmd("merge vol all group_results") + cubit.cmd("graphics tol angle 3") + + +cubit.cmd("reset") + +cubit.cmd("set attribute on") + +with open("manifest.json") as json_file: + data = byteify(json.load(json_file)) + +input_locations = [] +for entry in data: + if "tet_mesh" in entry.keys(): + input_locations.append(entry) +geometry_details = find_number_of_volumes_in_each_step_file( + input_locations, str(os.path.abspath(".")) +) + +imprint_and_merge_geometry() + +current_vols = cubit.parse_cubit_list("volume", "all") + +cubit.cmd("Trimesher volume gradation 1.3") + +cubit.cmd("volume all size auto factor 5") +print(geometry_details) +for entry in geometry_details: + for volume in entry["volumes"]: + cubit.cmd( + "volume " + str(volume) + " size auto factor 6" + ) # this number is the size of the mesh 1 is small 10 is large + cubit.cmd( + "volume all scheme tetmesh proximity layers off geometric sizing on") + if "size" in entry["mesh"]: + cubit.cmd("volume " + str(volume) + " " + + entry["tet_mesh"]) # ' size 0.5' + else: + cubit.cmd("volume " + str(volume)) + cubit.cmd("mesh volume " + str(volume)) + + +cubit.cmd('export mesh "tet_mesh.exo" overwrite') +# cubit.cmd('export abaqus "tet_mesh.inp" overwrite') # asci format, not goood for large meshes +# cubit.cmd('save as "tet_mesh.cub" overwrite') # mbconvert code is older +# than the exo equivilent + +print("unstrutured mesh saved as tet_mesh.exo") + +save_tet_details_to_json_file(geometry_details) diff --git a/examples/openmc_logo_example.ipynb b/examples/openmc_logo_example.ipynb new file mode 100644 index 0000000..ab09243 --- /dev/null +++ b/examples/openmc_logo_example.ipynb @@ -0,0 +1,181 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This example creates a couple of CAD volumes, combines them into a reactor and then performs an OpenMC / DAGMC simulation on the geometry. The geometry made is the OpenMC logo. \n", + "\n", + "There is a corresponding video fro this example\n", + "https://youtu.be/40VARwD44FA" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import paramak" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "These two code blocks create a geometry and adds relevent meta data tags that are used in the neutronics simulation." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "red_part = paramak.ExtrudeMixedShape(\n", + " points=[\n", + " (-89.0, -50.0, 'straight'),\n", + " (-18.0, -50.0, 'straight'),\n", + " (30.0, 30.0, 'straight'),\n", + " (-12.0, 100.0, 'circle'),\n", + " (88.0, -50.0, 'circle'),\n", + " ],\n", + " distance=10,\n", + " color=(1,0,0),\n", + " name='red_part',\n", + " material_tag='red_part_material',\n", + " stp_filename='red_part.stp',\n", + " stl_filename='red_part.stl',\n", + ")\n", + "red_part.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "grey_part = paramak.ExtrudeMixedShape(\n", + " points=[\n", + " (-96.5, -27.0, 'straight'),\n", + " (-30.0, -27.0, 'straight'),\n", + " (3.0, 30.0, 'straight'),\n", + " (-30.0, 87.0, 'straight'),\n", + " (-50.0, 87.0, 'straight'),\n", + " (-30.0, 87.0, 'circle'),\n", + " (-97.0, 30.0, 'circle'),\n", + " ],\n", + " distance=10,\n", + " color=(0.5,0.5,0.5),\n", + " material_tag='grey_part_material',\n", + " stp_filename='grey_part.stp',\n", + " stl_filename='grey_part.stl',\n", + ")\n", + "\n", + "grey_part.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This combines the two volumes into a single Reactor object" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "both_parts = paramak.Reactor([red_part, grey_part])\n", + "both_parts.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This defines the neutron source used in the simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import openmc\n", + "source = openmc.Source()\n", + "source.space = openmc.stats.Point((0, 50, 0))\n", + "source.energy = openmc.stats.Discrete([14e6], [1])\n", + "source.angle = openmc.stats.Isotropic()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This combines the geometry withthe neutron source and materials to make a model. The specifies three types of tallies (cell, mesh 2d, mesh 3d) and the reaction to tally (tritium production and heating). Output png and vtk files will be produced for the mesh tallies." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "my_model = paramak.NeutronicsModel(\n", + " geometry=both_parts,\n", + " source=source,\n", + " simulation_batches=2, # this should be increased to get a better mesh tally result\n", + " simulation_particles_per_batch=10, # this should be increased to get a better mesh tally result\n", + " materials = {'red_part_material':'P91', 'grey_part_material':'Li4SiO4'},\n", + " mesh_tally_3d=['(n,Xt)','heating'],\n", + " mesh_tally_2d=['(n,Xt)','heating'],\n", + " cell_tallies=['(n,Xt)','heating', 'spectra'],\n", + ")\n", + "\n", + "my_model.simulate()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This shows the cell tally results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "my_model.results" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/segmented_blanket_ball_reactor.py b/examples/segmented_blanket_ball_reactor.py new file mode 100644 index 0000000..f804050 --- /dev/null +++ b/examples/segmented_blanket_ball_reactor.py @@ -0,0 +1,250 @@ +"""This example makes a reactor geometry, neutronics model and performs a TBR +simulation. A selection of materials are used from refrenced sources to +complete the neutronics model.""" + +import neutronics_material_maker as nmm +import openmc +import paramak + + +def make_model_and_simulate(): + """Makes a neutronics Reactor model and simulates the TBR with specified materials""" + + # based on + # http://www.euro-fusionscipub.org/wp-content/uploads/WPBBCP16_15535_submitted.pdf + firstwall_radial_thickness = 3.0 + firstwall_armour_material = "tungsten" + firstwall_coolant_material = "He" + firstwall_structural_material = "eurofer" + firstwall_armour_fraction = 0.106305 + firstwall_coolant_fraction = 0.333507 + firstwall_coolant_temperature_k = 400 + firstwall_coolant_pressure_Pa = 8e6 + firstwall_structural_fraction = 0.560188 + + firstwall_material = nmm.Material.from_mixture( + name="firstwall_mat", + materials=[ + nmm.Material.from_library( + name=firstwall_coolant_material, + temperature=firstwall_coolant_temperature_k, + pressure=firstwall_coolant_pressure_Pa, + ), + nmm.Material.from_library(name=firstwall_structural_material), + nmm.Material.from_library(name=firstwall_armour_material), + ], + fracs=[ + firstwall_coolant_fraction, + firstwall_structural_fraction, + firstwall_armour_fraction, + ], + percent_type="vo" + ) + + # based on + # https://www.sciencedirect.com/science/article/pii/S2352179118300437 + blanket_rear_wall_coolant_material = "H2O" + blanket_rear_wall_structural_material = "eurofer" + blanket_rear_wall_coolant_fraction = 0.3 + blanket_rear_wall_structural_fraction = 0.7 + # units of Kelvin, equivalent 200 degrees C + blanket_rear_wall_coolant_temperature = 473.15 + blanket_rear_wall_coolant_pressure = 1e6 # units of Pa + + blanket_rear_wall_material = nmm.Material.from_mixture( + name="blanket_rear_wall_mat", + materials=[ + nmm.Material.from_library( + name=blanket_rear_wall_coolant_material, + temperature=blanket_rear_wall_coolant_temperature, + pressure=blanket_rear_wall_coolant_pressure, + ), + nmm.Material.from_library( + name=blanket_rear_wall_structural_material), + ], + fracs=[ + blanket_rear_wall_coolant_fraction, + blanket_rear_wall_structural_fraction, + ], + percent_type="vo") + + # based on + # https://www.sciencedirect.com/science/article/pii/S2352179118300437 + blanket_lithium6_enrichment_percent = 60 + blanket_breeder_material = "Li4SiO4" + blanket_coolant_material = "He" + blanket_multiplier_material = "Be" + blanket_structural_material = "eurofer" + blanket_breeder_fraction = 0.15 + blanket_coolant_fraction = 0.05 + blanket_multiplier_fraction = 0.6 + blanket_structural_fraction = 0.2 + blanket_breeder_packing_fraction = 0.64 + blanket_multiplier_packing_fraction = 0.64 + blanket_coolant_temperature_k = 773.15 + blanket_coolant_pressure_Pa = 1e6 + blanket_breeder_temperature_k = 873.15 + blanket_breeder_pressure_Pa = 8e6 + + blanket_material = nmm.Material.from_mixture( + name="blanket_mat", + materials=[ + nmm.Material.from_library( + name=blanket_coolant_material, + temperature=blanket_coolant_temperature_k, + pressure=blanket_coolant_pressure_Pa, + ), + nmm.Material.from_library(name=blanket_structural_material), + nmm.Material.from_library( + name=blanket_multiplier_material, + packing_fraction=blanket_multiplier_packing_fraction, + ), + nmm.Material.from_library( + name=blanket_breeder_material, + enrichment=blanket_lithium6_enrichment_percent, + packing_fraction=blanket_breeder_packing_fraction, + temperature=blanket_breeder_temperature_k, + pressure=blanket_breeder_pressure_Pa, + ), + ], + fracs=[ + blanket_coolant_fraction, + blanket_structural_fraction, + blanket_multiplier_fraction, + blanket_breeder_fraction, + ], + percent_type="vo" + ) + + # based on + # https://www.sciencedirect.com/science/article/pii/S2352179118300437 + divertor_coolant_fraction = 0.57195798876 + divertor_structural_fraction = 0.42804201123 + divertor_coolant_material = "H2O" + divertor_structural_material = "tungsten" + divertor_coolant_temperature_k = 423.15 # equivalent to 150 degrees C + divertor_coolant_pressure_Pa = 5e6 + + divertor_material = nmm.Material.from_mixture( + name="divertor_mat", + materials=[ + nmm.Material.from_library( + name=divertor_coolant_material, + temperature=divertor_coolant_temperature_k, + pressure=divertor_coolant_pressure_Pa, + ), + nmm.Material.from_library(name=divertor_structural_material), + ], + fracs=[divertor_coolant_fraction, divertor_structural_fraction], + percent_type="vo" + ) + + # based on + # https://pdfs.semanticscholar.org/95fa/4dae7d82af89adf711b97e75a241051c7129.pdf + center_column_shield_coolant_fraction = 0.13 + center_column_shield_structural_fraction = 0.57 + center_column_shield_coolant_material = "H2O" + center_column_shield_structural_material = "tungsten" + center_column_shield_coolant_temperature_k = 423.15 # equivalent to 150 degrees C + center_column_shield_coolant_pressure_Pa = 5e6 + + center_column_shield_material = nmm.Material.from_mixture( + name="center_column_shield_mat", + materials=[ + nmm.Material.from_library( + name=center_column_shield_coolant_material, + temperature=center_column_shield_coolant_temperature_k, + pressure=center_column_shield_coolant_pressure_Pa, + ), + nmm.Material.from_library( + name=center_column_shield_structural_material), + ], + fracs=[ + center_column_shield_coolant_fraction, + center_column_shield_structural_fraction, + ], + percent_type="vo") + + # based on + # https://pdfs.semanticscholar.org/95fa/4dae7d82af89adf711b97e75a241051c7129.pdf + inboard_tf_coils_conductor_fraction = 0.57 + inboard_tf_coils_coolant_fraction = 0.05 + inboard_tf_coils_structure_fraction = 0.38 + inboard_tf_coils_conductor_material = "copper" + inboard_tf_coils_coolant_material = "He" + inboard_tf_coils_structure_material = "SS_316L_N_IG" + inboard_tf_coils_coolant_temperature_k = 303.15 # equivalent to 30 degrees C + inboard_tf_coils_coolant_pressure_Pa = 8e6 + + inboard_tf_coils_material = nmm.Material.from_mixture( + name="inboard_tf_coils_mat", + materials=[ + nmm.Material.from_library( + name=inboard_tf_coils_coolant_material, + temperature=inboard_tf_coils_coolant_temperature_k, + pressure=inboard_tf_coils_coolant_pressure_Pa, + ), + nmm.Material.from_library( + name=inboard_tf_coils_conductor_material), + nmm.Material.from_library( + name=inboard_tf_coils_structure_material), + ], + fracs=[ + inboard_tf_coils_coolant_fraction, + inboard_tf_coils_conductor_fraction, + inboard_tf_coils_structure_fraction, + ], + percent_type="vo") + + # makes the 3d geometry + my_reactor = paramak.BallReactor( + inner_bore_radial_thickness=1, + inboard_tf_leg_radial_thickness=30, + center_column_shield_radial_thickness=60, + divertor_radial_thickness=50, + inner_plasma_gap_radial_thickness=30, + plasma_radial_thickness=300, + outer_plasma_gap_radial_thickness=30, + firstwall_radial_thickness=firstwall_radial_thickness, + # http://www.euro-fusionscipub.org/wp-content/uploads/WPBBCP16_15535_submitted.pdf + blanket_radial_thickness=100, + blanket_rear_wall_radial_thickness=3, + elongation=2.75, + triangularity=0.5, + number_of_tf_coils=16, + rotation_angle=360, + ) + + source = openmc.Source() + # sets the location of the source to x=0 y=0 z=0 + source.space = openmc.stats.Point((my_reactor.major_radius, 0, 0)) + # sets the direction to isotropic + source.angle = openmc.stats.Isotropic() + # sets the energy distribution to 100% 14MeV neutrons + source.energy = openmc.stats.Discrete([14e6], [1]) + + # makes the neutronics material + neutronics_model = paramak.NeutronicsModel( + geometry=my_reactor, + source=source, + materials={ + 'inboard_tf_coils_mat': inboard_tf_coils_material, + 'center_column_shield_mat': center_column_shield_material, + 'divertor_mat': divertor_material, + 'firstwall_mat': firstwall_material, + 'blanket_mat': blanket_material, + 'blanket_rear_wall_mat': blanket_rear_wall_material}, + cell_tallies=['TBR'], + simulation_batches=5, + simulation_particles_per_batch=1e4, + ) + + # starts the neutronics simulation + neutronics_model.simulate() + + # prints the simulation results to screen + print('TBR', neutronics_model.results['TBR']) + + +if __name__ == "__main__": + make_model_and_simulate() diff --git a/examples/shape_with_gas_production.py b/examples/shape_with_gas_production.py new file mode 100644 index 0000000..1894172 --- /dev/null +++ b/examples/shape_with_gas_production.py @@ -0,0 +1,50 @@ +"""Demonstrates the use of reaction rates in the cell tally. +(n,Xp) is MT number 203 and scores all proton (hydrogen) production +(n,Xt) is MT number 205 and scores all tritium production +(n,Xa) is MT number 207 and scores all alpha paticle (helium) production +https://docs.openmc.org/en/latest/usersguide/tallies.html#scores +""" + + +import openmc +import paramak + + +def main(): + my_shape = paramak.CenterColumnShieldHyperbola( + height=500, + inner_radius=50, + mid_radius=60, + outer_radius=100, + material_tag='center_column_shield_mat', + method='pymoab' + ) + + # makes the openmc neutron source at x,y,z 0, 0, 0 with isotropic + # directions + source = openmc.Source() + source.space = openmc.stats.Point((0, 0, 0)) + source.energy = openmc.stats.Discrete([14e6], [1]) + source.angle = openmc.stats.Isotropic() + + # converts the geometry into a neutronics geometry + my_model = paramak.NeutronicsModel( + geometry=my_shape, + source=source, + materials={'center_column_shield_mat': 'Be'}, + cell_tallies=['(n,Xa)', '(n,Xt)', '(n,Xp)'], + mesh_tally_3d=['(n,Xa)', '(n,Xt)', '(n,Xp)'], + mesh_tally_2d=['(n,Xa)', '(n,Xt)', '(n,Xp)'], + simulation_batches=2, + simulation_particles_per_batch=10, + ) + + # performs an openmc simulation on the model + my_model.simulate() + + # this extracts the values from the results dictionary + print(my_model.results) + + +if __name__ == "__main__": + main() diff --git a/examples/shape_with_spectra_cell_tally.py b/examples/shape_with_spectra_cell_tally.py new file mode 100644 index 0000000..e9e3772 --- /dev/null +++ b/examples/shape_with_spectra_cell_tally.py @@ -0,0 +1,114 @@ + +import openmc +import paramak +import plotly.graph_objects as go + + +def main(): + my_shape = paramak.CenterColumnShieldHyperbola( + height=500, + inner_radius=50, + mid_radius=60, + outer_radius=100, + material_tag='center_column_shield_mat', + method='pymoab' + ) + + # makes the openmc neutron source at x,y,z 0, 0, 0 with isotropic + # directions + source = openmc.Source() + source.space = openmc.stats.Point((0, 0, 0)) + source.energy = openmc.stats.Discrete([14e6], [1]) + source.angle = openmc.stats.Isotropic() + + # converts the geometry into a neutronics geometry + my_model = paramak.NeutronicsModel( + geometry=my_shape, + source=source, + materials={'center_column_shield_mat': 'Be'}, + cell_tallies=['heating', 'flux', 'TBR', 'spectra'], + simulation_batches=2, + simulation_particles_per_batch=20, + ) + + # performs an openmc simulation on the model + output_filename = my_model.simulate() + + # this extracts the values from the results dictionary + energy_bins = my_model.results['center_column_shield_mat_photon_spectra']['Flux per source particle']['energy'] + neutron_spectra = my_model.results['center_column_shield_mat_neutron_spectra']['Flux per source particle']['result'] + photon_spectra = my_model.results['center_column_shield_mat_photon_spectra']['Flux per source particle']['result'] + + fig = go.Figure() + + # this sets the axis titles and range + fig.update_layout( + xaxis={'title': 'Energy (eV)', + 'range': (0, 14.1e6)}, + yaxis={'title': 'Neutrons per cm2 per source neutron'} + ) + + # this adds the neutron spectra line to the plot + fig.add_trace(go.Scatter( + x=energy_bins[:-85], # trims off the high energy range + y=neutron_spectra[:-85], # trims off the high energy range + name='neutron spectra', + line=dict(shape='hv') + ) + ) + + # this adds the photon spectra line to the plot + fig.add_trace(go.Scatter( + x=energy_bins[:-85], # trims off the high energy range + y=photon_spectra[:-85], # trims off the high energy range + name='photon spectra', + line=dict(shape='hv') + ) + ) + + # this adds the drop down menu fo log and linear scales + fig.update_layout( + updatemenus=[ + go.layout.Updatemenu( + buttons=list([ + dict( + args=[{ + "xaxis.type": 'lin', "yaxis.type": 'lin', + 'xaxis.range': (0, 14.1e6) + }], + label="linear(x) , linear(y)", + method="relayout" + ), + dict( + args=[{"xaxis.type": 'log', "yaxis.type": 'log'}], + label="log(x) , log(y)", + method="relayout" + ), + dict( + args=[{"xaxis.type": 'log', "yaxis.type": 'lin'}], + label="log(x) , linear(y)", + method="relayout" + ), + dict( + args=[{"xaxis.type": 'lin', "yaxis.type": 'log', + 'xaxis.range': (0, 14.1e6) + }], + label="linear(x) , log(y)", + method="relayout" + ) + ]), + pad={"r": 10, "t": 10}, + showactive=True, + x=0.5, + xanchor="left", + y=1.1, + yanchor="top" + ), + ] + ) + + fig.show() + + +if __name__ == "__main__": + main() diff --git a/examples/submersion_reactor.ipynb b/examples/submersion_reactor.ipynb new file mode 100644 index 0000000..09f280c --- /dev/null +++ b/examples/submersion_reactor.ipynb @@ -0,0 +1,152 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This is a example that obtains the tritium breeding ratio (TBR)\n", + "for a parametric submersion reactor and specified the faceting and merge\n", + "tolerance when creating the dagmc model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import neutronics_material_maker as nmm\n", + "import openmc\n", + "import paramak\n", + "\n", + "\n", + "# makes the 3d geometry from input parameters\n", + "my_reactor = paramak.SubmersionTokamak(\n", + " inner_bore_radial_thickness=30,\n", + " inboard_tf_leg_radial_thickness=30,\n", + " center_column_shield_radial_thickness=30,\n", + " divertor_radial_thickness=80,\n", + " inner_plasma_gap_radial_thickness=50,\n", + " plasma_radial_thickness=200,\n", + " outer_plasma_gap_radial_thickness=50,\n", + " firstwall_radial_thickness=30,\n", + " blanket_rear_wall_radial_thickness=30,\n", + " rotation_angle=359.9, # if method is changed to Trelis then this can be changed to 360\n", + " support_radial_thickness=50,\n", + " inboard_blanket_radial_thickness=30,\n", + " outboard_blanket_radial_thickness=30,\n", + " elongation=2.75,\n", + " triangularity=0.5,\n", + ")\n", + "\n", + "# pymoab is used as it is open source and can be tested in the CI\n", + "# if you have Trelis or Cubit then this line can be changed\n", + "my_reactor.export_h5m(method='pymoab')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This section remakes the blankt material (FLiBe) for a range of temperatures and obtains the TBR for each one." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def make_model_and_simulate(temperature):\n", + " \"\"\"Makes a neutronics Reactor model and simulates the flux\n", + "\n", + " Arguments:\n", + " temperature: the temperature of the submersion blanket in Kelivin\n", + " \"\"\"\n", + "\n", + " # this can just be set as a string as temperature is needed for this\n", + " # material\n", + " flibe = nmm.Material.from_library(name='FLiBe', temperature=temperature, temperature_to_neutronics_code=False)\n", + "\n", + " source = openmc.Source()\n", + " # sets the location of the source to x=0 y=0 z=0\n", + " source.space = openmc.stats.Point((my_reactor.major_radius, 0, 0))\n", + " # sets the direction to isotropic\n", + " source.angle = openmc.stats.Isotropic()\n", + " # sets the energy distribution to 100% 14MeV neutrons\n", + " source.energy = openmc.stats.Discrete([14e6], [1])\n", + "\n", + " # makes the neutronics model from the geometry and material allocations\n", + " neutronics_model = paramak.NeutronicsModel(\n", + " geometry=my_reactor,\n", + " source=source,\n", + " materials={\n", + " 'inboard_tf_coils_mat': 'eurofer',\n", + " 'center_column_shield_mat': 'eurofer',\n", + " 'divertor_mat': 'eurofer',\n", + " 'firstwall_mat': 'eurofer',\n", + " 'blanket_rear_wall_mat': 'eurofer',\n", + " 'blanket_mat': flibe,\n", + " 'supports_mat': 'eurofer'},\n", + " cell_tallies=['TBR'],\n", + " simulation_batches=2,\n", + " simulation_particles_per_batch=10, # to get more accurate results this will need increasing\n", + " )\n", + "\n", + " # simulate the neutronics model\n", + " neutronics_model.simulate(export_h5m=False) # uising dagmc.h5m previous created\n", + " return neutronics_model.results['TBR']['result']\n", + "\n", + "\n", + "# loops through different temperatures finding the TBR value at each one\n", + "tbr_values = []\n", + "temperature_values = [32, 100, 300, 500]\n", + "for temperature in temperature_values:\n", + " tbr = make_model_and_simulate(temperature)\n", + " tbr_values.append(tbr)\n", + "\n", + "# plots the results\n", + "plt.scatter(temperature_values, tbr_values)\n", + "plt.xlabel('FLiBe Temperature (degrees C)')\n", + "plt.ylabel('TBR')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/submersion_reactor_minimal.py b/examples/submersion_reactor_minimal.py new file mode 100644 index 0000000..0644ffe --- /dev/null +++ b/examples/submersion_reactor_minimal.py @@ -0,0 +1,66 @@ +"""This is a minimal example that obtains the particle flux in each component +for a parametric submersion reactor""" + +import neutronics_material_maker as nmm +import openmc +import paramak + + +def make_model_and_simulate(): + """Makes a neutronics Reactor model and simulates the flux""" + + # makes the 3d geometry from input parameters + my_reactor = paramak.SubmersionTokamak( + inner_bore_radial_thickness=30, + inboard_tf_leg_radial_thickness=30, + center_column_shield_radial_thickness=30, + divertor_radial_thickness=80, + inner_plasma_gap_radial_thickness=50, + plasma_radial_thickness=200, + outer_plasma_gap_radial_thickness=50, + firstwall_radial_thickness=30, + blanket_rear_wall_radial_thickness=30, + rotation_angle=180, + support_radial_thickness=50, + inboard_blanket_radial_thickness=30, + outboard_blanket_radial_thickness=30, + elongation=2.75, + triangularity=0.5, + ) + + # this can just be set as a string as temperature is needed for this + # material + flibe = nmm.Material.from_library(name='FLiBe', temperature=773.15) + + source = openmc.Source() + # sets the location of the source to x=0 y=0 z=0 + source.space = openmc.stats.Point((my_reactor.major_radius, 0, 0)) + # sets the direction to isotropic + source.angle = openmc.stats.Isotropic() + # sets the energy distribution to 100% 14MeV neutrons + source.energy = openmc.stats.Discrete([14e6], [1]) + + # makes the neutronics model from the geometry and material allocations + neutronics_model = paramak.NeutronicsModel( + geometry=my_reactor, + source=source, + materials={ + 'inboard_tf_coils_mat': 'eurofer', + 'center_column_shield_mat': 'eurofer', + 'divertor_mat': 'eurofer', + 'firstwall_mat': 'eurofer', + 'blanket_rear_wall_mat': 'eurofer', + 'blanket_mat': flibe, + 'supports_mat': 'eurofer'}, + cell_tallies=['flux'], + simulation_batches=5, + simulation_particles_per_batch=1e4, + ) + + # simulate the neutronics model + neutronics_model.simulate() + print(neutronics_model.results) + + +if __name__ == "__main__": + make_model_and_simulate() diff --git a/examples/text_example.ipynb b/examples/text_example.ipynb new file mode 100644 index 0000000..629801a --- /dev/null +++ b/examples/text_example.ipynb @@ -0,0 +1,84 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A simple example of text being " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import cadquery as cq\n", + "import openmc\n", + "import paramak\n", + "\n", + "text = cq.Workplane().text( \n", + " txt=\"The Paramak can easily convert CAD\\nto neutronics models for simulations.\",\n", + " fontsize=0.5,\n", + " distance=-1,\n", + " cut=True,\n", + " halign=\"left\",\n", + " valign=\"bottom\",\n", + " font=\"Sans\"\n", + ")\n", + "\n", + "# this creates an empty Shape and populates it with the CQ object\n", + "my_shape = paramak.Shape()\n", + "my_shape.solid = text\n", + "my_shape.method = 'pymoab'\n", + "my_shape.stl_filename = 'text.stl'\n", + "my_shape.material_tag = 'my_material'\n", + "my_shape.export_stp('text_cad.stp')\n", + "\n", + "coords = openmc.stats.Point((4, 1, 0.5))\n", + "energy = openmc.stats.Discrete([1.41E7], [1.0])\n", + "source = openmc.Source(space=coords, energy=energy)\n", + "\n", + "my_model = paramak.NeutronicsModel(\n", + " geometry=my_shape,\n", + " source=source,\n", + " materials={'my_material': 'Li4SiO4'},\n", + " mesh_tally_3d=['heating', '(n,Xt)'],\n", + " simulation_batches=2, # to get a good mesh tally result with will need increasing\n", + " mesh_3d_resolution=(850, 350, 10),\n", + " mesh_3d_corners=[(0, -1, -1), (8.5, 2.5, 0.5)],\n", + ")\n", + "\n", + "my_model.simulate()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}