diff --git a/.circleci/config.yml b/.circleci/config.yml index 79898b9..4275b0c 100755 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -18,17 +18,17 @@ jobs: - run: name: run test_NeutronicModel command: - pytest tests/test_NeutronicModel.py -v --cov=paramak_neutronics --cov-append --cov-report term --cov-report xml + pytest tests/test_neutronic_model.py -v --cov=paramak_neutronics --cov-append --cov-report term --cov-report xml - run: - name: run test_Reactor_neutronics + name: run test_shape_neutronics command: - pytest tests/test_Reactor_neutronics.py -v --cov=paramak_neutronics --cov-append --cov-report term --cov-report xml + pytest tests/test_shape_neutronics.py -v --cov=paramak_neutronics --cov-append --cov-report term --cov-report xml - run: - name: run test_Shape_neutronics + name: run test_reactor_neutronics command: - pytest tests/test_Shape_neutronics.py -v --cov=paramak_neutronics --cov-append --cov-report term --cov-report xml + pytest tests/test_reactor_neutronics.py -v --cov=paramak_neutronics --cov-append --cov-report term --cov-report xml - run: name: run test_neutronics_utils diff --git a/Dockerfile b/Dockerfile index 368fd40..0a40ada 100755 --- a/Dockerfile +++ b/Dockerfile @@ -29,6 +29,7 @@ ENV LANG=C.UTF-8 LC_ALL=C.UTF-8 \ CC=/usr/bin/mpicc CXX=/usr/bin/mpicxx \ DEBIAN_FRONTEND=noninteractive +RUN apt-get --allow-releaseinfo-change update RUN apt-get update -y && \ apt-get upgrade -y @@ -126,7 +127,8 @@ RUN git clone --shallow-submodules --single-branch --branch main --depth 1 https RUN mkdir DAGMC && \ cd DAGMC && \ # git clone --single-branch --branch 3.2.0 --depth 1 https://github.com/svalinn/DAGMC.git && \ - git clone --shallow-submodules --single-branch --branch develop --depth 1 https://github.com/svalinn/DAGMC.git && \ + # git clone --shallow-submodules --single-branch --branch develop --depth 1 https://github.com/svalinn/DAGMC.git && \ + git clone --shallow-submodules --single-branch --branch dd_bbox --depth 1 https://github.com/pshriwise/DAGMC.git && \ mkdir build && \ cd build && \ cmake ../DAGMC -DBUILD_TALLY=ON \ @@ -165,15 +167,14 @@ ENV OPENMC_CROSS_SECTIONS=/nuclear_data/cross_sections.xml ENV PATH="/MOAB/build/bin:${PATH}" ENV PATH="/DAGMC/bin:${PATH}" -RUN pip install paramak FROM dependencies as final COPY run_tests.sh run_tests.sh COPY paramak_neutronics paramak_neutronics/ COPY setup.py setup.py -COPY examples examples/ -COPY tests tests/ COPY README.md README.md +COPY tests tests/ +COPY examples examples/ RUN python setup.py install diff --git a/examples/ball_reactor.ipynb b/examples/ball_reactor.ipynb index 6e478cb..6097bec 100644 --- a/examples/ball_reactor.ipynb +++ b/examples/ball_reactor.ipynb @@ -2,14 +2,16 @@ "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." - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "import neutronics_material_maker as nmm\n", "import openmc\n", @@ -61,7 +63,7 @@ "\n", "# makes the neutronics material\n", "neutronics_model = paramak_neutronics.NeutronicsModel(\n", - " geometry=my_reactor,\n", + " h5m_filename=my_reactor.export_h5m(),\n", " source=source,\n", " materials={\n", " 'inboard_tf_coils_mat': 'copper',\n", @@ -80,28 +82,20 @@ "\n", "# prints the results to screen\n", "print('TBR', neutronics_model.results['TBR'])" - ], - "outputs": [], - "metadata": {} + ] } ], "metadata": { + "interpreter": { + "hash": "71893251c0b0cd78b546a9c905312457257aecf63dc5b51c6968316964f80317" + }, "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", + "display_name": "Python 3.8.10 64-bit ('cq': conda)", "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.10" + "version": "" } }, "nbformat": 4, diff --git a/examples/ball_reactor_minimal.ipynb b/examples/ball_reactor_minimal.ipynb index 0187b64..afad992 100644 --- a/examples/ball_reactor_minimal.ipynb +++ b/examples/ball_reactor_minimal.ipynb @@ -2,14 +2,16 @@ "cells": [ { "cell_type": "markdown", + "metadata": {}, "source": [ "This is a minimal example that obtains the TBR (Tritium Breeding Ratio) for a parametric ball reactor" - ], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "import openmc\n", "import paramak\n", @@ -49,7 +51,7 @@ "\n", "# makes the neutronics model from the geometry and material allocations\n", "neutronics_model = paramak_neutronics.NeutronicsModel(\n", - " geometry=my_reactor,\n", + " h5m_filename=my_reactor.export_h5m(),\n", " source=source,\n", " materials={\n", " 'inboard_tf_coils_mat': 'eurofer',\n", @@ -66,35 +68,27 @@ "# simulate the neutronics model\n", "neutronics_model.simulate()\n", "print(neutronics_model.results)\n" - ], - "outputs": [], - "metadata": {} + ] }, { "cell_type": "code", "execution_count": null, - "source": [], + "metadata": {}, "outputs": [], - "metadata": {} + "source": [] } ], "metadata": { + "interpreter": { + "hash": "71893251c0b0cd78b546a9c905312457257aecf63dc5b51c6968316964f80317" + }, "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", + "display_name": "Python 3.8.10 64-bit ('cq': conda)", "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.10" + "version": "" } }, "nbformat": 4, diff --git a/examples/ball_reactor_source_plot.ipynb b/examples/ball_reactor_source_plot.ipynb index 8c05f3f..90bd93e 100644 --- a/examples/ball_reactor_source_plot.ipynb +++ b/examples/ball_reactor_source_plot.ipynb @@ -82,7 +82,7 @@ "\n", "# makes the neutronics material\n", "neutronics_model = paramak_neutronics.NeutronicsModel(\n", - " geometry=my_reactor,\n", + " h5m_filename=my_reactor.export_h5m(),\n", " source=source,\n", " materials={\n", " 'inboard_tf_coils_mat': 'copper',\n", @@ -132,24 +132,18 @@ } ], "metadata": { + "interpreter": { + "hash": "71893251c0b0cd78b546a9c905312457257aecf63dc5b51c6968316964f80317" + }, "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", + "display_name": "Python 3.8.10 64-bit ('cq': conda)", "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.10" + "version": "" } }, "nbformat": 4, "nbformat_minor": 4 -} +} \ No newline at end of file diff --git a/examples/center_column_study_reactor.ipynb b/examples/center_column_study_reactor.ipynb index 3c2b4a2..2fabdb5 100644 --- a/examples/center_column_study_reactor.ipynb +++ b/examples/center_column_study_reactor.ipynb @@ -58,7 +58,7 @@ " # makes the neutronics model and assigns basic materials to each\n", " # component\n", " neutronics_model = paramak_neutronics.NeutronicsModel(\n", - " geometry=my_reactor,\n", + " h5m_filename=my_reactor.export_h5m(),\n", " source=source,\n", " materials={\n", " 'DT_plasma': 'DT_plasma',\n", @@ -91,24 +91,18 @@ } ], "metadata": { + "interpreter": { + "hash": "71893251c0b0cd78b546a9c905312457257aecf63dc5b51c6968316964f80317" + }, "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", + "display_name": "Python 3.8.10 64-bit ('cq': conda)", "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.10" + "version": "" } }, "nbformat": 4, "nbformat_minor": 4 -} +} \ No newline at end of file diff --git a/examples/center_column_study_reactor_minimal.ipynb b/examples/center_column_study_reactor_minimal.ipynb index 8e1ef86..a57b40d 100644 --- a/examples/center_column_study_reactor_minimal.ipynb +++ b/examples/center_column_study_reactor_minimal.ipynb @@ -46,7 +46,7 @@ "\n", "# creates a neutronics model from the geometry and assigned materials\n", "neutronics_model = paramak_neutronics.NeutronicsModel(\n", - " geometry=my_reactor,\n", + " h5m_filename=my_reactor.export_h5m(),\n", " source=source,\n", " materials={\n", " 'inboard_tf_coils_mat': 'eurofer',\n", @@ -75,24 +75,18 @@ } ], "metadata": { + "interpreter": { + "hash": "71893251c0b0cd78b546a9c905312457257aecf63dc5b51c6968316964f80317" + }, "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", + "display_name": "Python 3.8.10 64-bit ('cq': conda)", "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.10" + "version": "" } }, "nbformat": 4, "nbformat_minor": 4 -} +} \ No newline at end of file diff --git a/examples/component_based_mesh_simulation.py b/examples/component_based_mesh_simulation.py index b90ea16..ec40656 100644 --- a/examples/component_based_mesh_simulation.py +++ b/examples/component_based_mesh_simulation.py @@ -29,7 +29,7 @@ def main(): # converts the geometry into a neutronics geometry my_model = paramak.NeutronicsModel( - geometry=my_shape, + h5m_filename=my_shape.export_h5m(), source=source, materials={'center_column_shield_mat': 'WB'}, # WB is tungsten boride cell_tallies=['heating', 'TBR'], diff --git a/examples/component_based_parameter_study.ipynb b/examples/component_based_parameter_study.ipynb index ad6f247..d824435 100644 --- a/examples/component_based_parameter_study.ipynb +++ b/examples/component_based_parameter_study.ipynb @@ -43,7 +43,7 @@ "\n", " # converts the geometry into a neutronics geometry\n", " my_model = paramak_neutronics.NeutronicsModel(\n", - " geometry=my_shape,\n", + " h5m_filename=my_shape.export_h5m(),\n", " source=source,\n", " materials={\n", " 'center_column_shield_mat': 'WB'},\n", @@ -87,24 +87,18 @@ } ], "metadata": { + "interpreter": { + "hash": "71893251c0b0cd78b546a9c905312457257aecf63dc5b51c6968316964f80317" + }, "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", + "display_name": "Python 3.8.10 64-bit ('cq': conda)", "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.10" + "version": "" } }, "nbformat": 4, "nbformat_minor": 4 -} +} \ No newline at end of file diff --git a/examples/make_unstructured_mesh.py b/examples/make_unstructured_mesh.py deleted file mode 100644 index 2f72faa..0000000 --- a/examples/make_unstructured_mesh.py +++ /dev/null @@ -1,139 +0,0 @@ -#!/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 index bf51b61..bd14bc1 100644 --- a/examples/openmc_logo_example.ipynb +++ b/examples/openmc_logo_example.ipynb @@ -128,7 +128,7 @@ "outputs": [], "source": [ "my_model = paramak_neutronics.NeutronicsModel(\n", - " geometry=both_parts,\n", + " h5m_filename=both_parts.export_h5m(),\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", @@ -179,4 +179,4 @@ }, "nbformat": 4, "nbformat_minor": 4 -} +} \ No newline at end of file diff --git a/examples/segmented_blanket_ball_reactor.py b/examples/segmented_blanket_ball_reactor.py index f804050..93adc45 100644 --- a/examples/segmented_blanket_ball_reactor.py +++ b/examples/segmented_blanket_ball_reactor.py @@ -225,7 +225,7 @@ def make_model_and_simulate(): # makes the neutronics material neutronics_model = paramak.NeutronicsModel( - geometry=my_reactor, + h5m_filename=my_reactor.export_h5m(), source=source, materials={ 'inboard_tf_coils_mat': inboard_tf_coils_material, diff --git a/examples/shape_with_gas_production.py b/examples/shape_with_gas_production.py index 1894172..623926e 100644 --- a/examples/shape_with_gas_production.py +++ b/examples/shape_with_gas_production.py @@ -29,7 +29,7 @@ def main(): # converts the geometry into a neutronics geometry my_model = paramak.NeutronicsModel( - geometry=my_shape, + h5m_filename=my_shape.export_h5m(), source=source, materials={'center_column_shield_mat': 'Be'}, cell_tallies=['(n,Xa)', '(n,Xt)', '(n,Xp)'], diff --git a/examples/shape_with_spectra_cell_tally.py b/examples/shape_with_spectra_cell_tally.py index e9e3772..7adbdc2 100644 --- a/examples/shape_with_spectra_cell_tally.py +++ b/examples/shape_with_spectra_cell_tally.py @@ -23,7 +23,7 @@ def main(): # converts the geometry into a neutronics geometry my_model = paramak.NeutronicsModel( - geometry=my_shape, + h5m_filename=my_shape.export_h5m(), source=source, materials={'center_column_shield_mat': 'Be'}, cell_tallies=['heating', 'flux', 'TBR', 'spectra'], diff --git a/examples/submersion_reactor.ipynb b/examples/submersion_reactor.ipynb index 7e1ad58..b5e1309 100644 --- a/examples/submersion_reactor.ipynb +++ b/examples/submersion_reactor.ipynb @@ -79,7 +79,7 @@ "\n", " # makes the neutronics model from the geometry and material allocations\n", " neutronics_model = paramak_neutronics.NeutronicsModel(\n", - " geometry=my_reactor,\n", + " h5m_filename=my_reactor.export_h5m(),\n", " source=source,\n", " materials={\n", " 'inboard_tf_coils_mat': 'eurofer',\n", @@ -129,24 +129,18 @@ } ], "metadata": { + "interpreter": { + "hash": "71893251c0b0cd78b546a9c905312457257aecf63dc5b51c6968316964f80317" + }, "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", + "display_name": "Python 3.8.10 64-bit ('cq': conda)", "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.10" + "version": "" } }, "nbformat": 4, "nbformat_minor": 4 -} +} \ No newline at end of file diff --git a/examples/submersion_reactor_minimal.py b/examples/submersion_reactor_minimal.py index 0644ffe..89a6482 100644 --- a/examples/submersion_reactor_minimal.py +++ b/examples/submersion_reactor_minimal.py @@ -42,7 +42,7 @@ def make_model_and_simulate(): # makes the neutronics model from the geometry and material allocations neutronics_model = paramak.NeutronicsModel( - geometry=my_reactor, + h5m_filename=my_reactor.export_h5m(), source=source, materials={ 'inboard_tf_coils_mat': 'eurofer', diff --git a/examples/text_example.ipynb b/examples/text_example.ipynb index 4a1249d..61a1f05 100644 --- a/examples/text_example.ipynb +++ b/examples/text_example.ipynb @@ -41,7 +41,7 @@ "source = openmc.Source(space=coords, energy=energy)\n", "\n", "my_model = paramak_neutronics.NeutronicsModel(\n", - " geometry=my_shape,\n", + " h5m_filename=my_shape.export_h5m(),\n", " source=source,\n", " materials={'my_material': 'Li4SiO4'},\n", " mesh_tally_3d=['heating', '(n,Xt)'],\n", @@ -62,24 +62,18 @@ } ], "metadata": { + "interpreter": { + "hash": "71893251c0b0cd78b546a9c905312457257aecf63dc5b51c6968316964f80317" + }, "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", + "display_name": "Python 3.8.10 64-bit ('cq': conda)", "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.10" + "version": "" } }, "nbformat": 4, "nbformat_minor": 4 -} +} \ No newline at end of file diff --git a/install_scripts/ubuntu-20.04-install.sh b/install_scripts/ubuntu-20.04-install.sh index e5793e0..aa631d6 100644 --- a/install_scripts/ubuntu-20.04-install.sh +++ b/install_scripts/ubuntu-20.04-install.sh @@ -12,6 +12,7 @@ printf '\nexport PATH="/opt/openmc/bin:$PATH"' >> ~/.bashrc printf '\nexport PATH="/opt/DAGMC/bin:$PATH"' >> ~/.bashrc printf '\nexport PATH="/opt/MOAB/bin:$PATH"' >> ~/.bashrc printf '\nexport LD_LIBRARY_PATH="/opt/openmc/lib:$LD_LIBRARY_PATH"' >> ~/.bashrc +printf '\nexport LD_LIBRARY_PATH="/opt/MOAB/lib:$LD_LIBRARY_PATH"' >> ~/.bashrc # CC=/usr/bin/mpicc # CXX=/usr/bin/mpicxx @@ -90,19 +91,20 @@ make -j"$compile_cores" install # Clone and install MOAB pip install --upgrade numpy cython -mkdir /opt/MOAB +git clone --single-branch --branch 5.3.0 --depth 1 https://bitbucket.org/fathomteam/moab.git /opt/MOAB/moab cd /opt/MOAB -git clone --single-branch --branch 5.2.1 --depth 1 https://bitbucket.org/fathomteam/moab.git mkdir build -cd build -cmake ../moab -DENABLE_HDF5=ON -DENABLE_NETCDF=ON -DENABLE_FORTRAN=OFF -DENABLE_BLASLAPACK=OFF -DBUILD_SHARED_LIBS=OFF -DCMAKE_INSTALL_PREFIX=/opt/MOAB -make -j"$compile_cores" -make -j"$compile_cores" install -rm -rf * -cmake ../moab -DENABLE_HDF5=ON -DENABLE_PYMOAB=ON -DENABLE_FORTRAN=OFF -DBUILD_SHARED_LIBS=ON -DENABLE_BLASLAPACK=OFF -DCMAKE_INSTALL_PREFIX=/opt/MOAB +cd /opt/MOAB/build +# this double build was needed in earlier versions of moab +# cmake ../moab -DENABLE_HDF5=ON -DHDF5_ROOT=/usr/lib/x86_64-linux-gnu/hdf5/serial -DENABLE_NETCDF=ON -DENABLE_FORTRAN=OFF -DENABLE_BLASLAPACK=OFF -DBUILD_SHARED_LIBS=OFF -DCMAKE_INSTALL_PREFIX=/opt/MOAB +# make -j"$compile_cores" +# make -j"$compile_cores" install +# rm -rf * +# -DHDF5_ROOT=$HDF5_DIR might be needed to to avoid conflicts with Cubit +cmake ../moab -DENABLE_HDF5=ON -DHDF5_ROOT=/usr/lib/x86_64-linux-gnu/hdf5/serial -DENABLE_PYMOAB=ON -DENABLE_FORTRAN=OFF -DBUILD_SHARED_LIBS=ON -DENABLE_BLASLAPACK=OFF -DCMAKE_INSTALL_PREFIX=/opt/MOAB make -j"$compile_cores" make -j"$compile_cores" install -cd pymoab +cd /opt/MOAB/build/pymoab bash install.sh python setup.py install @@ -111,7 +113,7 @@ python setup.py install git clone --single-branch --branch main https://github.com/pshriwise/double-down.git /opt/double-down cd /opt/double-down mkdir build -cd build +cd /opt/double-down/build cmake .. -DMOAB_DIR=/opt/MOAB -DCMAKE_INSTALL_PREFIX=.. -DEMBREE_DIR=/opt/embree make -j"$compile_cores" make -j"$compile_cores" install @@ -119,12 +121,13 @@ make -j"$compile_cores" install # Clone and install DAGMC -cd /opt -mkdir DAGMC -cd DAGMC +# cd /opt +# mkdir DAGMC +# cd DAGMC # TODO change to tagged release # git clone --single-branch --branch 3.2.0 --depth 1 https://github.com/svalinn/DAGMC.git -git clone --single-branch --branch develop --depth 1 https://github.com/svalinn/DAGMC.git +git clone --single-branch --branch develop --depth 1 https://github.com/svalinn/DAGMC.git /opt/DAGMC/DAGMC +cd /opt/DAGMC mkdir build cd build cmake ../DAGMC -DBUILD_TALLY=ON \ @@ -135,6 +138,7 @@ cmake ../DAGMC -DBUILD_TALLY=ON \ -DCMAKE_INSTALL_PREFIX=/opt/DAGMC/ \ -DDOUBLE_DOWN_DIR=/opt/double-down make -j"$compile_cores" install +# optional space saving to delete files rm -rf /opt/DAGMC/DAGMC /opt/DAGMC/build # Clone and install OpenMC with DAGMC @@ -152,8 +156,9 @@ sudo make -j"$compile_cores" install cd /opt/openmc pip install -e . +# vtk is an optional dependency pip install vtk -pip install neutronics_material_maker + # installs python packages for nuclear data pip install openmc_data_downloader @@ -188,7 +193,7 @@ sudo dpkg -i coreform-cubit-2021.5.deb wget https://github.com/svalinn/Cubit-plugin/releases/download/0.1.0/svalinn-plugin_ubuntu-20.04_cubit_2021.5.tgz sudo tar -xzvf svalinn-plugin_ubuntu-20.04_cubit_2021.5.tgz -C /opt/Coreform-Cubit-2021.5 -# check all packages are installed +# check all following packages are installed # sudo apt-get install -y libx11-6 # sudo apt-get install -y libxt6 # sudo apt-get install -y libgl1 diff --git a/paramak_neutronics/__init__.py b/paramak_neutronics/__init__.py index 1e4864d..43d6383 100644 --- a/paramak_neutronics/__init__.py +++ b/paramak_neutronics/__init__.py @@ -1,4 +1,3 @@ - from .neutronics_utils import get_neutronics_results_from_statepoint_file from .neutronics_utils import find_material_groups_in_h5m from .neutronics_utils import find_volume_ids_in_h5m diff --git a/paramak_neutronics/neutronics_model.py b/paramak_neutronics/neutronics_model.py index 883d221..c8b3876 100644 --- a/paramak_neutronics/neutronics_model.py +++ b/paramak_neutronics/neutronics_model.py @@ -1,33 +1,34 @@ - import json -import os import warnings from pathlib import Path from typing import List, Optional, Tuple, Union -import paramak -from .neutronics_utils import get_neutronics_results_from_statepoint_file -from .neutronics_utils import (create_inital_particles, - silently_remove_file, - extract_points_from_initial_source) +import openmc +import openmc.lib # needed to find bounding box of h5m file +import plotly.graph_objects as go +from openmc.data import REACTION_MT, REACTION_NAME from paramak.utils import plotly_trace -try: - import openmc - from openmc.data import REACTION_MT, REACTION_NAME -except ImportError: - warnings.warn('OpenMC not found, NeutronicsModelFromReactor.simulate \ - method not available', UserWarning) +from .neutronics_utils import ( + create_inital_particles, + extract_points_from_initial_source, + get_neutronics_results_from_statepoint_file, + silently_remove_file, + load_moab_file, +) try: import neutronics_material_maker as nmm except ImportError: - warnings.warn("neutronics_material_maker not found, \ + warnings.warn( + "neutronics_material_maker not found, \ NeutronicsModelFromReactor.materials can't accept strings or \ - neutronics_material_maker objects", UserWarning) + neutronics_material_maker objects", + UserWarning, + ) -class NeutronicsModel(): +class NeutronicsModel: """Creates a neutronics model of the provided shape geometry with assigned materials, source and neutronics tallies. @@ -67,11 +68,11 @@ class NeutronicsModel(): computational intensity is required to converge each mesh element. mesh_2d_corners: The upper and lower corner locations for the 2d mesh. This sets the location of the mesh. Defaults to None which - uses the NeutronicsModel.geometry.largest_dimension property to set + uses the NeutronicsModel.largest_dimension property to set the corners. mesh_3d_corners: The upper and lower corner locations for the 3d mesh. This sets the location of the mesh. Defaults to None which - uses the NeutronicsModel.geometry.largest_dimension property to set + uses the NeutronicsModel.largest_dimension property to set the corners. fusion_power: the power in watts emitted by the fusion reaction recalling that each DT fusion reaction emitts 17.6 MeV or @@ -80,7 +81,7 @@ class NeutronicsModel(): def __init__( self, - geometry: Union[paramak.Reactor, paramak.Shape], + h5m_filename: str, source, materials: dict, simulation_batches: Optional[int] = 100, @@ -90,18 +91,21 @@ def __init__( mesh_tally_3d: Optional[List[str]] = None, mesh_2d_resolution: Optional[Tuple[int, int, int]] = (400, 400), mesh_3d_resolution: Optional[Tuple[int, int, int]] = (100, 100, 100), - mesh_2d_corners: Optional[Tuple[Tuple[float, float, - float], Tuple[float, float, float]]] = None, - mesh_3d_corners: Optional[Tuple[Tuple[float, float, - float], Tuple[float, float, float]]] = None, + mesh_2d_corners: Optional[ + Tuple[Tuple[float, float, float], Tuple[float, float, float]] + ] = None, + mesh_3d_corners: Optional[ + Tuple[Tuple[float, float, float], Tuple[float, float, float]] + ] = None, fusion_power: Optional[float] = 1e9, photon_transport: Optional[bool] = True, # convert from watts to activity source_activity max_lost_particles: Optional[int] = 10, + largest_dimension: Optional[float] = 1000, ): - + self.largest_dimension = largest_dimension self.materials = materials - self.geometry = geometry + self.h5m_filename = h5m_filename self.source = source self.cell_tallies = cell_tallies self.mesh_tally_2d = mesh_tally_2d @@ -124,17 +128,15 @@ def __init__( self.statepoint_filename = None @property - def geometry(self): - return self._geometry + def h5m_filename(self): + return self._h5m_filename - @geometry.setter - def geometry(self, value): - if isinstance(value, (paramak.Shape, paramak.Reactor, type(None))): - self._geometry = value + @h5m_filename.setter + def h5m_filename(self, value): + if isinstance(value, str): + self._h5m_filename = value else: - raise TypeError( - "NeutronicsModelFromReactor.geometry should be a \ - paramak.Shape(), paramak.Reactor()") + raise TypeError("NeutronicsModelFromReactor.geometry should be a string") @property def source(self): @@ -145,7 +147,8 @@ def source(self, value): if not isinstance(value, (openmc.Source, type(None))): raise TypeError( "NeutronicsModelFromReactor.source should be an \ - openmc.Source() object") + openmc.Source() object" + ) self._source = value @property @@ -158,16 +161,21 @@ def cell_tallies(self, value): if not isinstance(value, list): raise TypeError( "NeutronicsModelFromReactor.cell_tallies should be a\ - list") - output_options = ['TBR', 'heating', 'flux', 'spectra', 'absorption'] + \ - list(REACTION_MT.keys()) + list(REACTION_NAME.keys()) + list" + ) + output_options = ( + ["TBR", "heating", "flux", "spectra", "absorption"] + + list(REACTION_MT.keys()) + + list(REACTION_NAME.keys()) + ) for entry in value: if entry not in output_options: raise ValueError( "NeutronicsModelFromReactor.cell_tallies argument", entry, "not allowed, the following options are supported", - output_options) + output_options, + ) self._cell_tallies = value @property @@ -180,16 +188,21 @@ def mesh_tally_2d(self, value): if not isinstance(value, list): raise TypeError( "NeutronicsModelFromReactor.mesh_tally_2d should be a\ - list") - output_options = ['heating', 'flux', 'absorption'] + \ - list(REACTION_MT.keys()) + list(REACTION_NAME.keys()) + list" + ) + output_options = ( + ["heating", "flux", "absorption"] + + list(REACTION_MT.keys()) + + list(REACTION_NAME.keys()) + ) for entry in value: if entry not in output_options: raise ValueError( "NeutronicsModelFromReactor.mesh_tally_2d argument", entry, "not allowed, the following options are supported", - output_options) + output_options, + ) self._mesh_tally_2d = value @property @@ -202,16 +215,21 @@ def mesh_tally_3d(self, value): if not isinstance(value, list): raise TypeError( "NeutronicsModelFromReactor.mesh_tally_3d should be a\ - list") - output_options = ['heating', 'flux', 'absorption'] + \ - list(REACTION_MT.keys()) + list(REACTION_NAME.keys()) + list" + ) + output_options = ( + ["heating", "flux", "absorption"] + + list(REACTION_MT.keys()) + + list(REACTION_NAME.keys()) + ) for entry in value: if entry not in output_options: raise ValueError( "NeutronicsModelFromReactor.mesh_tally_3d argument", entry, "not allowed, the following options are supported", - output_options) + output_options, + ) self._mesh_tally_3d = value @property @@ -221,8 +239,10 @@ def materials(self): @materials.setter def materials(self, value): if not isinstance(value, dict): - raise TypeError("NeutronicsModelFromReactor.materials should be a\ - dictionary") + raise TypeError( + "NeutronicsModelFromReactor.materials should be a\ + dictionary" + ) self._materials = value @property @@ -235,11 +255,10 @@ def simulation_batches(self, value): value = int(value) if not isinstance(value, int): raise TypeError( - "NeutronicsModelFromReactor.simulation_batches should be an int") - if value < 2: - raise ValueError( - "The minimum of setting for simulation_batches is 2" + "NeutronicsModelFromReactor.simulation_batches should be an int" ) + if value < 2: + raise ValueError("The minimum of setting for simulation_batches is 2") self._simulation_batches = value @property @@ -253,13 +272,15 @@ def simulation_particles_per_batch(self, value): if not isinstance(value, int): raise TypeError( "NeutronicsModelFromReactor.simulation_particles_per_batch\ - should be an int") + should be an int" + ) self._simulation_particles_per_batch = value def create_material(self, material_tag: str, material_entry): if isinstance(material_entry, str): openmc_material = nmm.Material.from_library( - name=material_entry, material_id=None).openmc_material + name=material_entry, material_id=None + ).openmc_material elif isinstance(material_entry, openmc.Material): # sets the material name in the event that it had not been set openmc_material = material_entry @@ -267,9 +288,13 @@ def create_material(self, material_tag: str, material_entry): # sets the material tag in the event that it had not been set openmc_material = material_entry.openmc_material else: - raise TypeError("materials must be either a str, \ + raise TypeError( + "materials must be either a str, \ openmc.Material, nmm.MultiMaterial or nmm.Material object \ - not a ", type(material_entry), material_entry) + not a ", + type(material_entry), + material_entry, + ) openmc_material.name = material_tag return openmc_material @@ -288,12 +313,11 @@ def create_openmc_materials(self): # "material has been added that is not needed for this \ # reactor model", reactor_material) - silently_remove_file('materials.xml') + silently_remove_file("materials.xml") openmc_materials = {} for material_tag, material_entry in self.materials.items(): - openmc_material = self.create_material( - material_tag, material_entry) + openmc_material = self.create_material(material_tag, material_entry) openmc_materials[material_tag] = openmc_material self.openmc_materials = openmc_materials @@ -301,22 +325,67 @@ def create_openmc_materials(self): self.mats = openmc.Materials(list(self.openmc_materials.values())) self.mats.export_to_xml() - return self.mats + def find_bounding_box(self): + """Computes the bounding box of the DAGMC geometry""" + + if not Path(self.h5m_filename).is_file: + msg = f'h5m file with filename {self.h5m_filename} not found' + raise FileNotFoundError(msg) + dag_univ = openmc.DAGMCUniverse(self.h5m_filename, auto_geom_ids=False) + + geometry = openmc.Geometry(root=dag_univ) + geometry.root_universe = dag_univ + geometry.export_to_xml() + + # exports materials.xml + # replace this with a empty materisl with the correct names + self.create_openmc_materials() + # openmc.Materials().export_to_xml() + + openmc.Plots().export_to_xml() + + # a minimal settings .xml to allow openmc to init + settings = openmc.Settings() + settings.verbosity = 1 + settings.batches = 1 + settings.particles = 1 + settings.export_to_xml() + + # The -p runs in plotting mode which avoids the check that OpenMC does + # when looking for boundary surfaces and therefore avoids this error + # ERROR: No boundary conditions were applied to any surfaces! + openmc.lib.init(["-p"]) + + bbox = openmc.lib.global_bounding_box() + openmc.lib.finalize() + + silently_remove_file("settings.xml") + silently_remove_file("plots.xml") + silently_remove_file("geometry.xml") + silently_remove_file("materials.xml") + return bbox + + # def build_csg_graveyard(self): + def export_xml( - self, - simulation_batches: Optional[int] = None, - source=None, - max_lost_particles: Optional[int] = None, - simulation_particles_per_batch: Optional[int] = None, - mesh_tally_3d: Optional[float] = None, - mesh_tally_2d: Optional[float] = None, - cell_tallies: Optional[float] = None, - mesh_2d_resolution: Optional[Tuple[int, int, int]] = None, - mesh_3d_resolution: Optional[Tuple[int, int, int]] = None, - mesh_2d_corners: Optional[Tuple[Tuple[float, float, float], Tuple[float, float, float]]] = None, - mesh_3d_corners: Optional[Tuple[Tuple[float, float, float], Tuple[float, float, float]]] = None, + self, + simulation_batches: Optional[int] = None, + source=None, + max_lost_particles: Optional[int] = None, + simulation_particles_per_batch: Optional[int] = None, + mesh_tally_3d: Optional[float] = None, + mesh_tally_2d: Optional[float] = None, + cell_tallies: Optional[float] = None, + mesh_2d_resolution: Optional[Tuple[int, int, int]] = None, + mesh_3d_resolution: Optional[Tuple[int, int, int]] = None, + mesh_2d_corners: Optional[ + Tuple[Tuple[float, float, float], Tuple[float, float, float]] + ] = None, + mesh_3d_corners: Optional[ + Tuple[Tuple[float, float, float], Tuple[float, float, float]] + ] = None, ): """Uses OpenMC python API to make a neutronics model, including tallies (cell_tallies and mesh_tally_2d), simulation settings (batches, @@ -396,19 +465,17 @@ def export_xml( mesh_3d_corners = self.mesh_3d_corners # this removes any old file from previous simulations - silently_remove_file('geometry.xml') - silently_remove_file('settings.xml') - silently_remove_file('tallies.xml') + silently_remove_file("geometry.xml") + silently_remove_file("settings.xml") + silently_remove_file("tallies.xml") # materials.xml is removed in this function self.create_openmc_materials() # this is the underlying geometry container that is filled with the - # faceteted DGAMC CAD model - dag_univ = openmc.DAGMCUniverse("dagmc.h5m") + # faceteted DAGMC CAD model + dag_univ = openmc.DAGMCUniverse(self.h5m_filename) geom = openmc.Geometry(root=dag_univ) - # self.universe = openmc.Universe() - # geom = openmc.Geometry(self.universe) # settings for the number of neutrons to simulate settings = openmc.Settings() @@ -416,7 +483,6 @@ def export_xml( settings.inactive = 0 settings.particles = self.simulation_particles_per_batch settings.run_mode = "fixed source" - # settings.dagmc = True settings.photon_transport = self.photon_transport settings.source = self.source @@ -426,19 +492,19 @@ def export_xml( self.tallies = openmc.Tallies() if self.mesh_tally_3d is not None: - mesh_xyz = openmc.RegularMesh(mesh_id=1, name='3d_mesh') + mesh_xyz = openmc.RegularMesh(mesh_id=1, name="3d_mesh") mesh_xyz.dimension = self.mesh_3d_resolution if self.mesh_3d_corners is None: mesh_xyz.lower_left = [ - -self.geometry.largest_dimension, - -self.geometry.largest_dimension, - -self.geometry.largest_dimension + -self.largest_dimension, + -self.largest_dimension, + -self.largest_dimension, ] mesh_xyz.upper_right = [ - self.geometry.largest_dimension, - self.geometry.largest_dimension, - self.geometry.largest_dimension + self.largest_dimension, + self.largest_dimension, + self.largest_dimension, ] else: mesh_xyz.lower_left = self.mesh_3d_corners[0] @@ -448,7 +514,7 @@ def export_xml( score = standard_tally prefix = standard_tally mesh_filter = openmc.MeshFilter(mesh_xyz) - tally = openmc.Tally(name=prefix + '_on_3D_mesh') + tally = openmc.Tally(name=prefix + "_on_3D_mesh") tally.filters = [mesh_filter] tally.scores = [score] self.tallies.append(tally) @@ -456,71 +522,71 @@ def export_xml( if self.mesh_tally_2d is not None: # Create mesh which will be used for tally - mesh_xz = openmc.RegularMesh(mesh_id=2, name='2d_mesh_xz') + mesh_xz = openmc.RegularMesh(mesh_id=2, name="2d_mesh_xz") mesh_xz.dimension = [ self.mesh_2d_resolution[1], 1, - self.mesh_2d_resolution[0] + self.mesh_2d_resolution[0], ] if self.mesh_2d_corners is None: mesh_xz.lower_left = [ - -self.geometry.largest_dimension, + -self.largest_dimension, -1, - -self.geometry.largest_dimension + -self.largest_dimension, ] mesh_xz.upper_right = [ - self.geometry.largest_dimension, + self.largest_dimension, 1, - self.geometry.largest_dimension + self.largest_dimension, ] else: mesh_xz.lower_left = self.mesh_2d_corners[0] mesh_xz.upper_right = self.mesh_2d_corners[1] - mesh_xy = openmc.RegularMesh(mesh_id=3, name='2d_mesh_xy') + mesh_xy = openmc.RegularMesh(mesh_id=3, name="2d_mesh_xy") mesh_xy.dimension = [ self.mesh_2d_resolution[1], self.mesh_2d_resolution[0], - 1 + 1, ] if self.mesh_2d_corners is None: mesh_xy.lower_left = [ - -self.geometry.largest_dimension, - -self.geometry.largest_dimension, - -1 + -self.largest_dimension, + -self.largest_dimension, + -1, ] mesh_xy.upper_right = [ - self.geometry.largest_dimension, - self.geometry.largest_dimension, - 1 + self.largest_dimension, + self.largest_dimension, + 1, ] else: mesh_xy.lower_left = self.mesh_2d_corners[0] mesh_xy.upper_right = self.mesh_2d_corners[1] - mesh_yz = openmc.RegularMesh(mesh_id=4, name='2d_mesh_yz') + mesh_yz = openmc.RegularMesh(mesh_id=4, name="2d_mesh_yz") mesh_yz.dimension = [ 1, self.mesh_2d_resolution[1], - self.mesh_2d_resolution[0] + self.mesh_2d_resolution[0], ] if self.mesh_2d_corners is None: mesh_yz.lower_left = [ -1, - -self.geometry.largest_dimension, - -self.geometry.largest_dimension + -self.largest_dimension, + -self.largest_dimension, ] mesh_yz.upper_right = [ 1, - self.geometry.largest_dimension, - self.geometry.largest_dimension + self.largest_dimension, + self.largest_dimension, ] else: mesh_yz.lower_left = self.mesh_2d_corners[0] @@ -531,9 +597,10 @@ def export_xml( prefix = standard_tally for mesh_filter, plane in zip( - [mesh_xz, mesh_xy, mesh_yz], ['xz', 'xy', 'yz']): + [mesh_xz, mesh_xy, mesh_yz], ["xz", "xy", "yz"] + ): mesh_filter = openmc.MeshFilter(mesh_filter) - tally = openmc.Tally(name=prefix + '_on_2D_mesh_' + plane) + tally = openmc.Tally(name=prefix + "_on_2D_mesh_" + plane) tally.filters = [mesh_filter] tally.scores = [score] self.tallies.append(tally) @@ -541,33 +608,31 @@ def export_xml( if self.cell_tallies is not None: for standard_tally in self.cell_tallies: - if standard_tally == 'TBR': - score = '(n,Xt)' # where X is a wild card - sufix = 'TBR' - tally = openmc.Tally(name='TBR') + if standard_tally == "TBR": + score = "(n,Xt)" # where X is a wild card + sufix = "TBR" + tally = openmc.Tally(name="TBR") tally.scores = [score] self.tallies.append(tally) self._add_tally_for_every_material(sufix, score) - elif standard_tally == 'spectra': + elif standard_tally == "spectra": - energy_bins = openmc.mgxs.GROUP_STRUCTURES['CCFE-709'] + energy_bins = openmc.mgxs.GROUP_STRUCTURES["CCFE-709"] energy_filter = openmc.EnergyFilter(energy_bins) - neutron_particle_filter = openmc.ParticleFilter([ - 'neutron']) + neutron_particle_filter = openmc.ParticleFilter(["neutron"]) self._add_tally_for_every_material( - 'neutron_spectra', - 'flux', - [neutron_particle_filter, energy_filter] + "neutron_spectra", + "flux", + [neutron_particle_filter, energy_filter], ) if self.photon_transport is True: - photon_particle_filter = openmc.ParticleFilter([ - 'photon']) + photon_particle_filter = openmc.ParticleFilter(["photon"]) self._add_tally_for_every_material( - 'photon_spectra', - 'flux', - [photon_particle_filter, energy_filter] + "photon_spectra", + "flux", + [photon_particle_filter, energy_filter], ) else: score = standard_tally @@ -575,8 +640,7 @@ def export_xml( self._add_tally_for_every_material(sufix, score) # make the model from geometry, materials, settings and tallies - model = openmc.model.Model( - geom, self.mats, settings, self.tallies) + model = openmc.model.Model(geom, self.mats, settings, self.tallies) geom.export_to_xml() settings.export_to_xml() @@ -585,8 +649,9 @@ def export_xml( self.model = model return model - def _add_tally_for_every_material(self, sufix: str, score: str, - additional_filters: List = None) -> None: + def _add_tally_for_every_material( + self, sufix: str, score: str, additional_filters: List = None + ) -> None: """Adds a tally to self.tallies for every material. Arguments: @@ -597,20 +662,20 @@ def _add_tally_for_every_material(self, sufix: str, score: str, if additional_filters is None: additional_filters = [] for key, value in self.openmc_materials.items(): - if key != 'DT_plasma': + if key != "DT_plasma": material_filter = openmc.MaterialFilter(value) - tally = openmc.Tally(name=key + '_' + sufix) + tally = openmc.Tally(name=key + "_" + sufix) tally.filters = [material_filter] + additional_filters tally.scores = [score] self.tallies.append(tally) def simulate( - self, - verbose: Optional[bool] = True, - cell_tally_results_filename: Optional[str] = 'results.json', - threads: Optional[int] = None, - export_h5m: Optional[bool] = True, - export_xml: Optional[bool] = True, + self, + verbose: Optional[bool] = True, + cell_tally_results_filename: Optional[str] = "results.json", + threads: Optional[int] = None, + export_h5m: Optional[bool] = True, + export_xml: Optional[bool] = True, ) -> str: """Run the OpenMC simulation. Deletes existing simulation output (summary.h5) if files exists. @@ -641,51 +706,52 @@ def simulate( if export_xml is True: self.export_xml() - if export_h5m is True: - silently_remove_file('dagmc.h5') - self.geometry.export_h5m() - # checks all the nessecary files are found - for required_file in ['geometry.xml', 'materials.xml', 'settings.xml', - 'tallies.xml']: + for required_file in [ + "geometry.xml", + "materials.xml", + "settings.xml", + "tallies.xml", + ]: if not Path(required_file).is_file(): msg = "{} file was not found. Please set export_xml \ to True or use the export_xml() \ - method to create the xml files".format(required_file) + method to create the xml files".format( + required_file + ) raise FileNotFoundError(msg) - if not Path('dagmc.h5m').is_file(): - msg = """dagmc.h5m file was not found. Please set export_h5m to - True or use the export_h5m() methods to create the dagmc.h5m - file""" + if not Path(self.h5m_filename).is_file(): + msg = f"""{self.h5m_filename} file was not found. Please set + export_h5m to True or use the export_h5m() methods to create + the dagmc.h5m file""" raise FileNotFoundError(msg) # Deletes summary.h5m if it already exists. # This avoids permission problems when trying to overwrite the file - silently_remove_file('summary.h5') - silently_remove_file('statepoint.' + str(self.simulation_batches) + '.h5') + silently_remove_file("summary.h5") + silently_remove_file("statepoint." + str(self.simulation_batches) + ".h5") - self.statepoint_filename = self.model.run( - output=verbose, threads=threads - ) + self.statepoint_filename = self.model.run(output=verbose, threads=threads) self.results = get_neutronics_results_from_statepoint_file( - statepoint_filename=self.statepoint_filename, - fusion_power=self.fusion_power + statepoint_filename=self.statepoint_filename, fusion_power=self.fusion_power ) - with open(cell_tally_results_filename, 'w') as outfile: + with open(cell_tally_results_filename, "w") as outfile: json.dump(self.results, outfile, indent=4, sort_keys=True) return self.statepoint_filename def export_html( - self, - filename: Optional[str] = "neutronics_model.html", - facet_splines: Optional[bool] = True, - facet_circles: Optional[bool] = True, - tolerance: Optional[float] = 1., - view_plane: Optional[str] = 'RZ', - number_of_source_particles: Optional[int] = 1000): + self, + figure: go.Figure() = go.Figure(), + filename: Optional[str] = "neutronics_model.html", + # facet_splines: Optional[bool] = True, + # facet_circles: Optional[bool] = True, + # tolerance: Optional[float] = 1., + view_plane: Optional[str] = "RZ", + number_of_source_particles: Optional[int] = 1000, + ): """Creates a html graph representation of the points for the Shape objects that make up the reactor and optionally the source. Shapes are colored by their .color property. Shapes are also labelled by their @@ -693,14 +759,13 @@ def export_html( added. Args: + figure: The Plotly figure to add the source points to. + Paramak.export_html() returns a go.Figure() objct that can be + passed in here and have source points added to it. Otherwise + this defaults to plotly.graph_objects.Figure() which provides + an empty figure for source points. filename: the filename used to save the html graph. Defaults to neutronics_model.html - facet_splines: If True then spline edges will be faceted. Defaults - to True. - facet_circles: If True then circle edges will be faceted. Defaults - to True. - tolerance: faceting toleranceto use when faceting cirles and - splines. Defaults to 1e-3. view_plane: The plane to project. Options are 'XZ', 'XY', 'YZ', 'YX', 'ZY', 'ZX', 'RZ' and 'XYZ'. Defaults to 'RZ'. Defaults to 'RZ'. @@ -709,27 +774,13 @@ def export_html( plotly.Figure(): figure object """ - fig = self.geometry.export_html( - filename=None, - facet_splines=facet_splines, - facet_circles=facet_circles, - tolerance=tolerance, - view_plane=view_plane, - ) - if number_of_source_particles != 0: source_filename = create_inital_particles( - self.source, number_of_source_particles) - points = extract_points_from_initial_source( - source_filename, view_plane) - - fig.add_trace( - plotly_trace( - points=points, - mode="markers", - name='source' - ) + self.source, number_of_source_particles ) + points = extract_points_from_initial_source(source_filename, view_plane) + + figure.add_trace(plotly_trace(points=points, mode="markers", name="source")) if filename is not None: @@ -740,8 +791,8 @@ def export_html( if path_filename.suffix != ".html": path_filename = path_filename.with_suffix(".html") - fig.write_html(str(path_filename)) + figure.write_html(str(path_filename)) print("Exported html graph to ", path_filename) - return fig + return figure diff --git a/paramak_neutronics/neutronics_utils.py b/paramak_neutronics/neutronics_utils.py index 42b4a59..7326acd 100644 --- a/paramak_neutronics/neutronics_utils.py +++ b/paramak_neutronics/neutronics_utils.py @@ -1,7 +1,6 @@ - +import errno import math import os -import errno import subprocess import warnings from collections import defaultdict @@ -11,17 +10,19 @@ import defusedxml.ElementTree as ET import matplotlib.pyplot as plt import numpy as np +import openmc from pymoab import core, types -try: - import openmc -except ImportError: - warnings.warn('OpenMC not found, create_inital_particles \ - method not available', UserWarning) + +def load_moab_file(filename: str): + """Loads a h5m into a Moab Core object and returns the object""" + moab_core = core.Core() + moab_core.load_file(filename) + return moab_core def silently_remove_file(filename: str): - """Allows files to be deleted without printing warning messages int the + """Allows files to be deleted without printing warning messages int the terminal. input XML files for OpenMC are deleted prior to running simulations and many not exist.""" try: @@ -30,9 +31,7 @@ def silently_remove_file(filename: str): pass # in some cases the file will not exist -def find_volume_ids_in_h5m( - filename: Optional[str] = 'dagmc.h5m' -) -> List[str]: +def find_volume_ids_in_h5m(filename: Optional[str] = "dagmc.h5m") -> List[str]: """Reads in a DAGMC h5m file and uses PyMoab to find the volume ids of the volumes in the file @@ -44,33 +43,34 @@ def find_volume_ids_in_h5m( """ # create a new PyMOAB instance and load the specified DAGMC file - mb = core.Core() - mb.load_file(filename) + moab_core = load_moab_file(filename) # retrieve the category tag on the instance try: - cat_tag = mb.tag_get_handle(types.CATEGORY_TAG_NAME) + cat_tag = moab_core.tag_get_handle(types.CATEGORY_TAG_NAME) except types.MB_ENTITY_NOT_FOUND: - raise RuntimeError("The category tag could not be found in the PyMOAB instance." - "Please check that the DAGMC file has been loaded.") + msg = ( + "The category tag could not be found in the PyMOAB instance." + "Please check that the DAGMC file has been loaded." + ) + raise RuntimeError(msg) # get the id tag - gid_tag = mb.tag_get_handle(types.GLOBAL_ID_TAG_NAME) + gid_tag = moab_core.tag_get_handle(types.GLOBAL_ID_TAG_NAME) # get the set of entities using the provided category tag name # (0 means search on the instance's root set) - ents = mb.get_entities_by_type_and_tag(0, types.MBENTITYSET, [cat_tag], ["Volume"]) + ents = moab_core.get_entities_by_type_and_tag( + 0, types.MBENTITYSET, [cat_tag], ["Volume"] + ) # retrieve the IDs of the entities - ids = mb.tag_get_data(gid_tag, ents).flatten() + ids = moab_core.tag_get_data(gid_tag, ents).flatten() return sorted(list(ids)) - -def find_material_groups_in_h5m( - filename: Optional[str] = 'dagmc.h5m' -) -> List[str]: +def find_material_groups_in_h5m(filename: Optional[str] = "dagmc.h5m") -> List[str]: """Reads in a DAGMC h5m file and uses mbsize to find the names of the material groups in the file @@ -91,86 +91,19 @@ def find_material_groups_in_h5m( raise ValueError( "mbsize failed, check MOAB is install and the MOAB/build/bin " "folder is in the path directory (Linux and Mac) or set as an " - "enviromental varible (Windows)") + "enviromental varible (Windows)" + ) list_of_mats = terminal_output.split() - list_of_mats = list(filter(lambda a: a != '=', list_of_mats)) - list_of_mats = list(filter(lambda a: a != 'NAME', list_of_mats)) - list_of_mats = list(filter(lambda a: a != 'EXTRA_NAME0', list_of_mats)) + list_of_mats = list(filter(lambda a: a != "=", list_of_mats)) + list_of_mats = list(filter(lambda a: a != "NAME", list_of_mats)) + list_of_mats = list(filter(lambda a: a != "EXTRA_NAME0", list_of_mats)) list_of_mats = list(set(list_of_mats)) return list_of_mats -def remove_tag_from_h5m_file( - input_h5m_filename: Optional[str] = 'dagmc.h5m', - output_h5m_filename: Optional[str] = 'dagmc_removed_tag.h5m', - tag_to_remove: Optional[str] = 'graveyard', -) -> str: - """Removes a specific tag from a dagmc h5m file and saves the remaining - geometry as a new h5m file. Useful for visulising the geometry by removing - the graveyard tag and then the vtk file can be made without a bounding box - graveyard obstructing the view. Adapted from - https://github.com/svalinn/DAGMC-viz source code - - Arguments: - input_h5m_filename: The name of the h5m file to remove the graveyard from - output_h5m_filename: The name of the outfile h5m without a graveyard - - Returns: - filename of the new dagmc h5m file with the tags removed - """ - - try: - from pymoab import core, types - from pymoab.types import MBENTITYSET - except ImportError: - raise ImportError( - 'PyMoab not found, remove_tag_from_h5m_file method is not ' - 'available' - ) - - moab_core = core.Core() - moab_core.load_file(input_h5m_filename) - - tag_name = moab_core.tag_get_handle(str(types.NAME_TAG_NAME)) - - tag_category = moab_core.tag_get_handle(str(types.CATEGORY_TAG_NAME)) - root = moab_core.get_root_set() - - # An array of tag values to be matched for entities returned by the - # following call. - group_tag_values = np.array(["Group"]) - - # Retrieve all EntitySets with a category tag of the user input value. - group_categories = list(moab_core.get_entities_by_type_and_tag( - root, MBENTITYSET, tag_category, group_tag_values)) - - # Retrieve all EntitySets with a name tag. - group_names = moab_core.tag_get_data(tag_name, group_categories, flat=True) - - # Find the EntitySet whose name includes tag provided - sets_to_remove = [ - group_set for group_set, - name in zip( - group_categories, - group_names) if tag_to_remove in str( - name.lower())] - - # Remove the graveyard EntitySet from the data. - groups_to_write = [ - group_set for group_set in group_categories if group_set not in sets_to_remove] - - moab_core.write_file(output_h5m_filename, output_sets=groups_to_write) - - return output_h5m_filename - - -def _save_2d_mesh_tally_as_png( - score: str, - filename: str, - tally -) -> str: +def _save_2d_mesh_tally_as_png(score: str, filename: str, tally) -> str: """Extracts 2D mesh tally results from a tally and saves the result as a png image. @@ -182,13 +115,6 @@ def _save_2d_mesh_tally_as_png( resutls from. """ - try: - import openmc - except ImportError as err: - raise err( - 'openmc not found, _save_2d_mesh_tally_as_png method is not \ - available') - my_slice = tally.get_slice(scores=[score]) tally_filter = tally.find_filter(filter_type=openmc.MeshFilter) shape = tally_filter.mesh.dimension.tolist() @@ -203,8 +129,7 @@ def _save_2d_mesh_tally_as_png( def get_neutronics_results_from_statepoint_file( - statepoint_filename: str, - fusion_power: Optional[float] = None + statepoint_filename: str, fusion_power: Optional[float] = None ) -> dict: """Reads the statepoint file from the neutronics simulation and extracts the tally results. @@ -218,13 +143,6 @@ def get_neutronics_results_from_statepoint_file( dict: a dictionary of the simulation results """ - try: - import openmc - except ImportError as err: - raise err( - 'openmc not found, get_neutronics_results_from_statepoint_file \ - method is not available') - # open the results file statepoint = openmc.StatePoint(statepoint_filename) @@ -233,83 +151,76 @@ def get_neutronics_results_from_statepoint_file( # access the tallies for tally in statepoint.tallies.values(): - if tally.name.endswith('TBR'): + if tally.name.endswith("TBR"): data_frame = tally.get_pandas_dataframe() tally_result = data_frame["mean"].sum() - tally_std_dev = data_frame['std. dev.'].sum() + tally_std_dev = data_frame["std. dev."].sum() results[tally.name] = { - 'result': tally_result, - 'std. dev.': tally_std_dev, + "result": tally_result, + "std. dev.": tally_std_dev, } - elif tally.name.endswith('heating'): + elif tally.name.endswith("heating"): data_frame = tally.get_pandas_dataframe() tally_result = data_frame["mean"].sum() - tally_std_dev = data_frame['std. dev.'].sum() - results[tally.name]['MeV per source particle'] = { - 'result': tally_result / 1e6, - 'std. dev.': tally_std_dev / 1e6, + tally_std_dev = data_frame["std. dev."].sum() + results[tally.name]["MeV per source particle"] = { + "result": tally_result / 1e6, + "std. dev.": tally_std_dev / 1e6, } if fusion_power is not None: - results[tally.name]['Watts'] = { - 'result': tally_result * 1.602176487e-19 * (fusion_power / ((17.58 * 1e6) / 6.2415090744e18)), - 'std. dev.': tally_std_dev * 1.602176487e-19 * (fusion_power / ((17.58 * 1e6) / 6.2415090744e18)), + results[tally.name]["Watts"] = { + "result": tally_result + * 1.602176487e-19 + * (fusion_power / ((17.58 * 1e6) / 6.2415090744e18)), + "std. dev.": tally_std_dev + * 1.602176487e-19 + * (fusion_power / ((17.58 * 1e6) / 6.2415090744e18)), } - elif tally.name.endswith('flux'): + elif tally.name.endswith("flux"): data_frame = tally.get_pandas_dataframe() tally_result = data_frame["mean"].sum() - tally_std_dev = data_frame['std. dev.'].sum() - results[tally.name]['Flux per source particle'] = { - 'result': tally_result, - 'std. dev.': tally_std_dev, + tally_std_dev = data_frame["std. dev."].sum() + results[tally.name]["Flux per source particle"] = { + "result": tally_result, + "std. dev.": tally_std_dev, } - elif tally.name.endswith('spectra'): + elif tally.name.endswith("spectra"): data_frame = tally.get_pandas_dataframe() tally_result = data_frame["mean"] - tally_std_dev = data_frame['std. dev.'] - results[tally.name]['Flux per source particle'] = { - 'energy': openmc.mgxs.GROUP_STRUCTURES['CCFE-709'].tolist(), - 'result': tally_result.tolist(), - 'std. dev.': tally_std_dev.tolist(), + tally_std_dev = data_frame["std. dev."] + results[tally.name]["Flux per source particle"] = { + "energy": openmc.mgxs.GROUP_STRUCTURES["CCFE-709"].tolist(), + "result": tally_result.tolist(), + "std. dev.": tally_std_dev.tolist(), } - elif '_on_2D_mesh' in tally.name: - score = tally.name.split('_')[0] + elif "_on_2D_mesh" in tally.name: + score = tally.name.split("_")[0] _save_2d_mesh_tally_as_png( score=score, tally=tally, - filename=tally.name.replace( - '(', - '').replace( - ')', - '').replace( - ',', - '-')) - - elif '_on_3D_mesh' in tally.name: + filename=tally.name.replace("(", "").replace(")", "").replace(",", "-"), + ) + + elif "_on_3D_mesh" in tally.name: mesh_id = 1 mesh = statepoint.meshes[mesh_id] xs = np.linspace( - mesh.lower_left[0], - mesh.upper_right[0], - mesh.dimension[0] + 1 + mesh.lower_left[0], mesh.upper_right[0], mesh.dimension[0] + 1 ) ys = np.linspace( - mesh.lower_left[1], - mesh.upper_right[1], - mesh.dimension[1] + 1 + mesh.lower_left[1], mesh.upper_right[1], mesh.dimension[1] + 1 ) zs = np.linspace( - mesh.lower_left[2], - mesh.upper_right[2], - mesh.dimension[2] + 1 + mesh.lower_left[2], mesh.upper_right[2], mesh.dimension[2] + 1 ) tally = statepoint.get_tally(name=tally.name) @@ -322,7 +233,7 @@ def get_neutronics_results_from_statepoint_file( for content in [data, error]: for counter, i in enumerate(content): if math.isnan(i): - content[counter] = 0. + content[counter] = 0.0 write_3d_mesh_tally_to_vtk( xs=xs, @@ -331,36 +242,31 @@ def get_neutronics_results_from_statepoint_file( tally_label=tally.name, tally_data=data, error_data=error, - outfile=tally.name.replace( - '(', - '').replace( - ')', - '').replace( - ',', - '-') + - '.vtk') + outfile=tally.name.replace("(", "").replace(")", "").replace(",", "-") + + ".vtk", + ) else: # this must be a standard score cell tally data_frame = tally.get_pandas_dataframe() tally_result = data_frame["mean"].sum() - tally_std_dev = data_frame['std. dev.'].sum() - results[tally.name]['events per source particle'] = { - 'result': tally_result, - 'std. dev.': tally_std_dev, + tally_std_dev = data_frame["std. dev."].sum() + results[tally.name]["events per source particle"] = { + "result": tally_result, + "std. dev.": tally_std_dev, } return results def write_3d_mesh_tally_to_vtk( - xs: np.linspace, - ys: np.linspace, - zs: np.linspace, - tally_data: List[float], - error_data: Optional[List[float]] = None, - outfile: Optional[str] = '3d_mesh_tally_data.vtk', - tally_label: Optional[str] = '3d_mesh_tally_data', + xs: np.linspace, + ys: np.linspace, + zs: np.linspace, + tally_data: List[float], + error_data: Optional[List[float]] = None, + outfile: Optional[str] = "3d_mesh_tally_data.vtk", + tally_label: Optional[str] = "3d_mesh_tally_data", ) -> str: """Converts regular 3d data into a vtk file for visualising the data. Programs that can visualise vtk files include Paraview @@ -385,8 +291,10 @@ def write_3d_mesh_tally_to_vtk( try: import vtk except (ImportError, ModuleNotFoundError): - msg = "Conversion to VTK requested," \ + msg = ( + "Conversion to VTK requested," "but the Python VTK module is not installed. Try pip install pyvtk" + ) raise ImportError(msg) vtk_box = vtk.vtkRectilinearGrid() @@ -394,17 +302,17 @@ def write_3d_mesh_tally_to_vtk( vtk_box.SetDimensions(len(xs), len(ys), len(zs)) vtk_x_array = vtk.vtkDoubleArray() - vtk_x_array.SetName('x-coords') + vtk_x_array.SetName("x-coords") vtk_x_array.SetArray(xs, len(xs), True) vtk_box.SetXCoordinates(vtk_x_array) vtk_y_array = vtk.vtkDoubleArray() - vtk_y_array.SetName('y-coords') + vtk_y_array.SetName("y-coords") vtk_y_array.SetArray(ys, len(ys), True) vtk_box.SetYCoordinates(vtk_y_array) vtk_z_array = vtk.vtkDoubleArray() - vtk_z_array.SetName('z-coords') + vtk_z_array.SetName("z-coords") vtk_z_array.SetArray(zs, len(zs), True) vtk_box.SetZCoordinates(vtk_z_array) @@ -428,17 +336,14 @@ def write_3d_mesh_tally_to_vtk( writer.SetInputData(vtk_box) - print('Writing %s' % outfile) + print("Writing %s" % outfile) writer.Write() return outfile -def create_inital_particles( - source, - number_of_source_particles: int = 2000 -) -> str: +def create_inital_particles(source, number_of_source_particles: int = 2000) -> str: """Accepts an openmc source and creates an inital_source.h5 that can be used to find intial xyz, direction and energy of the partice source. @@ -459,7 +364,7 @@ def create_inital_particles( # GEOMETRY # just a minimal geometry - outer_surface = openmc.Sphere(r=100000, boundary_type='vacuum') + outer_surface = openmc.Sphere(r=100000, boundary_type="vacuum") cell = openmc.Cell(region=-outer_surface) universe = openmc.Universe(cells=[cell]) geom = openmc.Geometry(universe) @@ -479,11 +384,11 @@ def create_inital_particles( model = openmc.model.Model(geom, mats, sett) - silently_remove_file('settings.xml') - silently_remove_file('materials.xml') - silently_remove_file('geometry.xml') - silently_remove_file('settings.xml') - silently_remove_file('tallies.xml') + silently_remove_file("settings.xml") + silently_remove_file("materials.xml") + silently_remove_file("geometry.xml") + silently_remove_file("settings.xml") + silently_remove_file("tallies.xml") model.export_to_xml() # this just adds write_initial_source == True to the settings.xml @@ -504,8 +409,7 @@ def create_inital_particles( def extract_points_from_initial_source( - input_filename: str = 'initial_source.h5', - view_plane: str = 'RZ' + input_filename: str = "initial_source.h5", view_plane: str = "RZ" ) -> list: """Reads in an inital source h5 file (generated by OpenMC), extracts point and projects them onto a view plane. @@ -521,31 +425,30 @@ def extract_points_from_initial_source( list: list of points extracted """ import h5py - h5_file = h5py.File(input_filename, 'r') - dset = h5_file['source_bank'] + + h5_file = h5py.File(input_filename, "r") + dset = h5_file["source_bank"] points = [] for particle in dset: - if view_plane == 'XZ': + if view_plane == "XZ": points.append((particle[0][0], particle[0][2])) - elif view_plane == 'XY': + elif view_plane == "XY": points.append((particle[0][0], particle[0][1])) - elif view_plane == 'YZ': + elif view_plane == "YZ": points.append((particle[0][1], particle[0][2])) - elif view_plane == 'YX': + elif view_plane == "YX": points.append((particle[0][1], particle[0][0])) - elif view_plane == 'ZY': + elif view_plane == "ZY": points.append((particle[0][2], particle[0][1])) - elif view_plane == 'ZX': + elif view_plane == "ZX": points.append((particle[0][2], particle[0][0])) - elif view_plane == 'RZ': - xy_coord = math.pow(particle[0][0], 2) + \ - math.pow(particle[0][1], 2) + elif view_plane == "RZ": + xy_coord = math.pow(particle[0][0], 2) + math.pow(particle[0][1], 2) points.append((math.sqrt(xy_coord), particle[0][2])) - elif view_plane == 'XYZ': + elif view_plane == "XYZ": points.append((particle[0][0], particle[0][1], particle[0][2])) else: - raise ValueError('view_plane value of ', view_plane, - ' is not supported') + raise ValueError("view_plane value of ", view_plane, " is not supported") return points diff --git a/run_tests.sh b/run_tests.sh index 1a60ee0..2000950 100644 --- a/run_tests.sh +++ b/run_tests.sh @@ -1,7 +1,7 @@ pytest tests/notebook_testing.py -v --cov=paramak_neutronics --cov-append --cov-report term --cov-report xml -pytest tests/test_NeutronicModel.py -v --cov=paramak_neutronics --cov-append --cov-report term --cov-report xml -pytest tests/test_Reactor_neutronics.py -v --cov=paramak_neutronics --cov-append --cov-report term --cov-report xml -pytest tests/test_Shape_neutronics.py -v --cov=paramak_neutronics --cov-append --cov-report term --cov-report xml +pytest tests/test_neutronics_model.py -v --cov=paramak_neutronics --cov-append --cov-report term --cov-report xml +pytest tests/test_reactor_neutronics.py -v --cov=paramak_neutronics --cov-append --cov-report term --cov-report xml +pytest tests/test_shape_neutronics.py -v --cov=paramak_neutronics --cov-append --cov-report term --cov-report xml pytest tests/test_example_neutronics_simulations.py -v --cov=paramak_neutronics --cov-append --cov-report term --cov-report xml pytest tests/test_neutronics_utils.py -v --cov=paramak_neutronics --cov-append --cov-report term --cov-report xml diff --git a/setup.py b/setup.py index 05dfa55..426acae 100644 --- a/setup.py +++ b/setup.py @@ -27,6 +27,7 @@ "remove_dagmc_tags", "neutronics_material_maker", "vtk", + "neutronics_material_maker", ], # openmc, dagmc, pymoab are also needed but not avaible on PyPi ) diff --git a/tests/test_Reactor_neutronics.py b/tests/test_Reactor_neutronics.py deleted file mode 100644 index 9810ea7..0000000 --- a/tests/test_Reactor_neutronics.py +++ /dev/null @@ -1,88 +0,0 @@ - -import os -import unittest -from pathlib import Path - -import paramak - - -class TestReactorNeutronics(unittest.TestCase): - - def test_export_h5m(self): - """Creates a Reactor object consisting of two shapes and checks a h5m - file of the reactor can be exported using the export_h5m method.""" - - os.system('rm small_dagmc.h5m') - os.system('rm small_dagmc_without_graveyard.h5m') - os.system('rm small_dagmc_with_graveyard.h5m') - os.system('rm large_dagmc.h5m') - test_shape = paramak.RotateStraightShape( - points=[(0, 0), (0, 20), (20, 20)], - material_tag='mat1') - test_shape2 = paramak.RotateSplineShape( - points=[(0, 0), (0, 20), (20, 20)], - material_tag='mat2') - test_shape.rotation_angle = 360 - test_reactor = paramak.Reactor([test_shape, test_shape2]) - test_reactor.export_h5m( - filename='small_dagmc.h5m', - faceting_tolerance=0.01 - ) - test_reactor.export_h5m( - filename='small_dagmc_without_graveyard.h5m', - faceting_tolerance=0.01, - include_graveyard=False - ) - test_reactor.export_h5m( - filename='small_dagmc_with_graveyard.h5m', - faceting_tolerance=0.01, - include_graveyard=True - ) - test_reactor.export_h5m( - filename='large_dagmc.h5m', - faceting_tolerance=0.001 - ) - - assert Path("small_dagmc.h5m").exists() is True - assert Path("small_dagmc_with_graveyard.h5m").exists() is True - assert Path("large_dagmc.h5m").exists() is True - assert Path("large_dagmc.h5m").stat().st_size > Path( - "small_dagmc.h5m").stat().st_size - assert Path("small_dagmc_without_graveyard.h5m").stat( - ).st_size < Path("small_dagmc.h5m").stat().st_size - - def test_export_h5m_without_extension(self): - """Tests that the code appends .h5m to the end of the filename""" - - os.system('rm out.h5m') - test_shape = paramak.RotateStraightShape( - points=[(0, 0), (0, 20), (20, 20)], - material_tag='mat1') - test_shape2 = paramak.RotateSplineShape( - points=[(0, 0), (0, 20), (20, 20)], - material_tag='mat2') - test_shape.rotation_angle = 360 - test_reactor = paramak.Reactor([test_shape, test_shape2]) - test_reactor.export_h5m(filename='out', faceting_tolerance=0.01) - assert Path("out.h5m").exists() is True - os.system('rm out.h5m') - - def test_make_graveyard_accepts_offset_from_graveyard(self): - """Creates a graveyard for a reactor and sets the graveyard_offset. - Checks that the Reactor.graveyard_offset property is set""" - - test_shape = paramak.RotateStraightShape( - points=[(0, 0), (0, 20), (20, 20)], - material_tag='mat1') - test_shape2 = paramak.RotateSplineShape( - points=[(0, 0), (0, 20), (20, 20)], - material_tag='mat2') - test_shape.rotation_angle = 360 - test_reactor = paramak.Reactor([test_shape, test_shape2]) - test_reactor.graveyard_offset == 101 - graveyard = test_reactor.make_graveyard() - assert graveyard.volume > 0 - - -if __name__ == "__main__": - unittest.main() diff --git a/tests/test_NeutronicModel.py b/tests/test_neutronics_model.py similarity index 94% rename from tests/test_NeutronicModel.py rename to tests/test_neutronics_model.py index cb77490..c5837c7 100644 --- a/tests/test_NeutronicModel.py +++ b/tests/test_neutronics_model.py @@ -23,6 +23,8 @@ def setUp(self): method='pymoab' ) + self.h5m_filename = self.my_shape.export_h5m() + # makes the openmc neutron source at x,y,z 0, 0, 0 with isotropic # directions and 14MeV neutrons source = openmc.Source() @@ -37,14 +39,12 @@ def simulation_with_previous_h5m_file(self): os.system('rm *.h5m') my_model = paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 'WC'}, ) - self.my_shape.export_h5m() - - my_model.simulate(export_h5m=False) + my_model.simulate() my_model.results is not None @@ -57,7 +57,7 @@ def test_neutronics_component_simulation_with_openmc_mat(self): # converts the geometry into a neutronics geometry my_model = paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': test_mat}, cell_tallies=['heating'], @@ -84,7 +84,7 @@ def test_neutronics_component_simulation_with_nmm(self): # converts the geometry into a neutronics geometry my_model = paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': test_mat}, cell_tallies=['heating'], @@ -116,7 +116,7 @@ def test_cell_tally_output_file_creation(self): # converts the geometry into a neutronics geometry # this simulation has no tally to test this edge case my_model = paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': test_mat}, simulation_batches=2, @@ -140,7 +140,7 @@ def test_incorrect_faceting_tolerance(self): def incorrect_faceting_tolerance(): """Sets faceting_tolerance as a string which should raise an error""" paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 'eurofer'}, faceting_tolerance='coucou' @@ -156,7 +156,7 @@ def test_incorrect_merge_tolerance(self): def incorrect_merge_tolerance(): """Set merge_tolerance as a string which should raise an error""" paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 'eurofer'}, merge_tolerance='coucou' @@ -172,7 +172,7 @@ def test_incorrect_cell_tallies(self): def incorrect_cell_tallies(): """Set a cell tally that is not accepted which should raise an error""" paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 'eurofer'}, cell_tallies=['coucou'], @@ -188,7 +188,7 @@ def test_incorrect_cell_tally_type(self): def incorrect_cell_tally_type(): """Set a cell tally that is the wrong type which should raise an error""" paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 'eurofer'}, cell_tallies=1, @@ -204,7 +204,7 @@ def test_incorrect_mesh_tally_2d(self): def incorrect_mesh_tally_2d(): """Set a mesh_tally_2d that is not accepted which should raise an error""" paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 'eurofer'}, mesh_tally_2d=['coucou'], @@ -220,7 +220,7 @@ def test_incorrect_mesh_tally_2d_type(self): def incorrect_mesh_tally_2d_type(): """Set a mesh_tally_2d that is the wrong type which should raise an error""" paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 'eurofer'}, mesh_tally_2d=1, @@ -236,7 +236,7 @@ def test_incorrect_mesh_tally_3d(self): def incorrect_mesh_tally_3d(): """Set a mesh_tally_3d that is not accepted which should raise an error""" paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 'eurofer'}, mesh_tally_3d=['coucou'], @@ -252,7 +252,7 @@ def test_incorrect_mesh_tally_3d_type(self): def incorrect_mesh_tally_3d_type(): """Set a mesh_tally_3d that is the wrong type which should raise an error""" paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 'eurofer'}, mesh_tally_3d=1, @@ -268,7 +268,7 @@ def test_incorrect_materials(self): def incorrect_materials(): """Set a material as a string which should raise an error""" paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials='coucou', ) @@ -283,7 +283,7 @@ def test_incorrect_materials_type(self): def incorrect_materials_type(): """Sets a material as an int which should raise an error""" test_model = paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 23}, ) @@ -300,7 +300,7 @@ def test_incorrect_simulation_batches_to_small(self): def incorrect_simulation_batches_to_small(): """Sets simulation batch below 2 which should raise an error""" paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 'eurofer'}, simulation_batches=1 @@ -316,7 +316,7 @@ def test_incorrect_simulation_batches_wrong_type(self): def incorrect_simulation_batches_wrong_type(): """Sets simulation_batches as a string which should raise an error""" paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 'eurofer'}, simulation_batches='one' @@ -332,7 +332,7 @@ def test_incorrect_simulation_particles_per_batch_wrong_type(self): def incorrect_simulation_particles_per_batch_wrong_type(): """Sets simulation_particles_per_batch below 2 which should raise an error""" paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 'eurofer'}, simulation_particles_per_batch='one' @@ -353,7 +353,7 @@ def test_neutronics_component_cell_simulation_heating(self): # converts the geometry into a neutronics geometry my_model = paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': mat}, cell_tallies=['heating', 'flux', 'TBR', 'spectra'], @@ -396,7 +396,7 @@ def test_neutronics_component_2d_mesh_simulation(self): # converts the geometry into a neutronics geometry my_model = paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 'Be'}, mesh_tally_2d=['heating'], @@ -423,7 +423,7 @@ def test_neutronics_component_3d_mesh_simulation(self): # converts the geometry into a neutronics geometry my_model = paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 'Be'}, mesh_tally_3d=['heating', '(n,Xt)'], @@ -451,7 +451,7 @@ def test_batches_and_particles_convert_to_int(self): # converts the geometry into a neutronics geometry my_model = paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 'Be'}, simulation_batches=3.1, @@ -472,7 +472,7 @@ def test_neutronics_component_3d_and_2d_mesh_simulation(self): # converts the geometry into a neutronics geometry my_model = paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 'Be'}, mesh_tally_3d=['heating'], @@ -503,7 +503,7 @@ def test_neutronics_component_3d_and_2d_mesh_simulation_with_corner_points( # converts the geometry into a neutronics geometry my_model = paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 'Be'}, mesh_tally_3d=['heating'], @@ -545,7 +545,7 @@ def test_reactor_from_shapes_cell_tallies(self): test_reactor = paramak.Reactor([test_shape, test_shape2]) neutronics_model = paramak_neutronics.NeutronicsModel( - geometry=test_reactor, + h5m_filename=test_reactor.export_h5m(), source=self.source, materials={ 'mat1': 'copper', @@ -578,7 +578,7 @@ def test_reactor_from_shapes_2d_mesh_tallies(self): test_reactor = paramak.Reactor([test_shape, test_shape2]) neutronics_model = paramak_neutronics.NeutronicsModel( - geometry=test_reactor, + h5m_filename=test_reactor.export_h5m(), source=self.source, materials={ 'mat1': 'copper', @@ -612,7 +612,7 @@ def test_export_h5m_error_handling(): and tests the TBR simulation value""" self.my_shape.method = 'cubit' paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 'WC'}, ) @@ -630,7 +630,7 @@ def test_missing_xml_file_error_handling(): with a FileNotFoundError""" my_model = paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 'WC'}, ) @@ -654,7 +654,7 @@ def test_missing_h5m_file_error_handling(): with a FileNotFoundError""" my_model = paramak_neutronics.NeutronicsModel( - geometry=self.my_shape, + h5m_filename=self.h5m_filename, source=self.source, materials={'center_column_shield_mat': 'WC'}, ) @@ -717,8 +717,8 @@ def test_neutronics_model_attributes(self): # makes the neutronics material neutronics_model = paramak_neutronics.NeutronicsModel( - geometry=self.my_reactor, source=openmc.Source(), + h5m_filename='placeholder.h5m', materials={ 'inboard_tf_coils_mat': 'copper', 'center_column_shield_mat': 'WC', @@ -731,7 +731,7 @@ def test_neutronics_model_attributes(self): simulation_particles_per_batch=12, ) - assert neutronics_model.geometry == self.my_reactor + assert neutronics_model.h5m_filename == 'placeholder.h5m' assert neutronics_model.materials == { 'inboard_tf_coils_mat': 'copper', diff --git a/tests/test_reactor_neutronics.py b/tests/test_reactor_neutronics.py new file mode 100644 index 0000000..4245d7c --- /dev/null +++ b/tests/test_reactor_neutronics.py @@ -0,0 +1,74 @@ +import os +import unittest +from pathlib import Path + +import openmc +import paramak +import paramak_neutronics +import pytest + + +class TestNeutronicsModelWithReactor(unittest.TestCase): + """Tests Shape object arguments that involve neutronics usage""" + + def setUp(self): + self.test_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, + ) + + def test_bounding_box_size(self): + + h5m_filename = self.test_reactor.export_h5m_with_pymoab( + include_graveyard=False, + faceting_tolerance=1e-1 + ) + + # makes the openmc neutron source at x,y,z 0, 0, 0 with isotropic + # directions and 14MeV neutrons + source = openmc.Source() + source.space = openmc.stats.Point((0, 0, 0)) + source.angle = openmc.stats.Isotropic() + source.energy = openmc.stats.Discrete([14e6], [1]) + + h5m_filename='dagmc.h5m' + my_model = paramak_neutronics.NeutronicsModel( + h5m_filename=h5m_filename, + source=source, + materials={ + 'center_column_shield_mat': 'Be', + 'firstwall_mat': 'Be', + 'blanket_mat': 'Be', + 'divertor_mat': 'Be', + 'supports_mat': 'Be', + 'blanket_rear_wall_mat': 'Be', + 'inboard_tf_coils_mat': 'Be' + }, + simulation_batches=3, + simulation_particles_per_batch=2 + ) + + bounding_box=my_model.find_bounding_box() + + assert len(bounding_box) == 2 + assert len(bounding_box[0]) == 3 + assert len(bounding_box[1]) == 3 + assert bounding_box[0][0] == pytest.approx(-540, abs=0.2) + assert bounding_box[0][1] == pytest.approx(0, abs=0.2) + assert bounding_box[0][2] == pytest.approx(-415, abs=0.2) + assert bounding_box[1][0] == pytest.approx(540, abs=0.2) + assert bounding_box[1][1] == pytest.approx(540, abs=0.2) + assert bounding_box[1][2] == pytest.approx(415, abs=0.2) diff --git a/tests/test_Shape_neutronics.py b/tests/test_shape_neutronics.py similarity index 85% rename from tests/test_Shape_neutronics.py rename to tests/test_shape_neutronics.py index a6adc0e..48af2c5 100644 --- a/tests/test_Shape_neutronics.py +++ b/tests/test_shape_neutronics.py @@ -21,7 +21,8 @@ def setUp(self): (70, 50, "circle"), (60, 25, "circle"), (70, 0, "straight")], - distance=50 + distance=50, + material_tag='test_shape' ) def test_export_h5m_creates_file(self): @@ -76,6 +77,41 @@ def test_graveyard_offset_increases_voulme(self): large_offset = self.test_shape.graveyard.volume assert small_offset < large_offset + def test_bounding_box_size(self): + + h5m_filename = self.test_shape.export_h5m_with_pymoab( + include_graveyard=False, + faceting_tolerance=1e-1 + ) + + # makes the openmc neutron source at x,y,z 0, 0, 0 with isotropic + # directions and 14MeV neutrons + source = openmc.Source() + source.space = openmc.stats.Point((0, 0, 0)) + source.angle = openmc.stats.Isotropic() + source.energy = openmc.stats.Discrete([14e6], [1]) + + h5m_filename='dagmc.h5m' + my_model = paramak_neutronics.NeutronicsModel( + h5m_filename=h5m_filename, + source=source, + materials={'test_shape': 'Be'}, + simulation_batches=3, + simulation_particles_per_batch=2 + ) + + bounding_box=my_model.find_bounding_box() + + assert len(bounding_box) == 2 + assert len(bounding_box[0]) == 3 + assert len(bounding_box[1]) == 3 + assert bounding_box[0][0] == pytest.approx(50, abs=0.1) + assert bounding_box[0][1] == pytest.approx(-25, abs=0.1) + assert bounding_box[0][2] == pytest.approx(0, abs=0.1) + assert bounding_box[1][0] == pytest.approx(70, abs=0.1) + assert bounding_box[1][1] == pytest.approx(25, abs=0.1) + assert bounding_box[1][2] == pytest.approx(70, abs=0.1) + class TestSimulationResultsVsCsg(unittest.TestCase): """Makes a geometry in the paramak and in CSG geometry, simulates and @@ -216,7 +252,7 @@ def simulate_cylinder_cask_cad( ) my_model = paramak_neutronics.NeutronicsModel( - geometry=my_geometry, + h5m_filename=my_geometry.export_h5m(), source=source, simulation_batches=batches, simulation_particles_per_batch=particles,