diff --git a/examples/tokamak_structured_mesh_simulation.py b/examples/tokamak_structured_mesh_simulation.py new file mode 100644 index 00000000..a709b2a4 --- /dev/null +++ b/examples/tokamak_structured_mesh_simulation.py @@ -0,0 +1,132 @@ +import openmc + +import paramak + +openmc.config['cross_sections'] = '/nuclear_data/cross_sections.xml' # change this to the dir of your cross_sections.xml + +my_reactor = paramak.tokamak_from_plasma( + radial_build=[ + (paramak.LayerType.GAP, 10), + (paramak.LayerType.SOLID, 30), # layer_1 + (paramak.LayerType.SOLID, 50), # layer_2 + (paramak.LayerType.SOLID, 10), # layer_3 + (paramak.LayerType.SOLID, 120), # layer_4 + (paramak.LayerType.SOLID, 20), # layer_5 + (paramak.LayerType.GAP, 60), + (paramak.LayerType.PLASMA, 300), # plasma + (paramak.LayerType.GAP, 60), + (paramak.LayerType.SOLID, 20), # layer_3 + (paramak.LayerType.SOLID, 120), # layer_4 + (paramak.LayerType.SOLID, 10), # layer_5 + ], + elongation=2, + triangularity=0.55, + rotation_angle=180, +) +my_reactor.save(f"tokamak_minimal.step") +print(f"Saved as tokamak_minimal.step") + +# this prints all the names of the parts in the reactor +for part in my_reactor: + print(part[1]) + +from cad_to_dagmc import CadToDagmc +my_model = CadToDagmc() +material_tags = ["layer_1", "layer_2", "layer_3", "layer_4", "layer_5", "plasma"] + +# 6 material tags as inner and outer layers are one solid there are only 6 solids in model +my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) + +h5m_filename = "dagmc.h5m" +# my_model.export_dagmc_h5m_file(filename=h5m_filename, min_mesh_size=10.0, max_mesh_size=20.0) + +# simplified material definitions have been used to keen this example minimal +mat_layer_1 = openmc.Material(name='layer_1') +mat_layer_1.add_element('Cu', 1, "ao") +mat_layer_1.set_density("g/cm3", 7) + +mat_layer_2 = openmc.Material(name='layer_2') +mat_layer_2.add_nuclide('W186', 1, "ao") +mat_layer_2.set_density("g/cm3", 0.01) + +mat_layer_3 = openmc.Material(name='layer_3') +mat_layer_3.add_nuclide('Fe56', 1, "ao") +mat_layer_3.set_density("g/cm3", 7) + +mat_layer_4 = openmc.Material(name='layer_4') +mat_layer_4.add_element('Li', 1, "ao") +mat_layer_4.set_density("g/cm3", 0.5) + +mat_layer_5 = openmc.Material(name='layer_5') +mat_layer_5.add_nuclide('Fe56', 1, "ao") +mat_layer_5.set_density("g/cm3", 7) + +mat_plasma = openmc.Material(name='plasma') +mat_plasma.add_nuclide('Fe56', 1, "ao") +mat_plasma.set_density("g/cm3", 7) + + +materials = openmc.Materials([mat_layer_1, mat_layer_2, mat_layer_3, mat_layer_4, mat_layer_5, mat_plasma]) + + + +dag_univ = openmc.DAGMCUniverse(filename=h5m_filename) +bound_dag_univ = dag_univ.bounded_universe() +geometry = openmc.Geometry(root=bound_dag_univ) + +# initializes a new source object +my_source = openmc.IndependentSource() + +# initialises a new source object +my_source = openmc.IndependentSource() + +# the distribution of radius is just a single value +radius = openmc.stats.Discrete([10], [1]) + +# the distribution of source z values is just a single value +z_values = openmc.stats.Discrete([0], [1]) + +# the distribution of source azimuthal angles values is a uniform distribution between 0 and 2 Pi +angle = openmc.stats.Uniform(a=0., b=2* 3.14159265359) + +# this makes the ring source using the three distributions and a radius +my_source.space = openmc.stats.CylindricalIndependent(r=radius, phi=angle, z=z_values, origin=(0.0, 0.0, 0.0)) + +# sets the direction to isotropic +my_source.angle = openmc.stats.Isotropic() + +# sets the energy distribution to a Muir distribution neutrons +my_source.energy = openmc.stats.muir(e0=14080000.0, m_rat=5.0, kt=20000.0) + + +# specifies the simulation computational intensity +settings = openmc.Settings() +settings.batches = 10 +settings.particles = 10000 +settings.run_mode = "fixed source" +settings.source = my_source + +# adds a tally to record the heat deposited in a mesh overlaid over entire geometry + +mesh = openmc.RegularMesh().from_domain(bound_dag_univ) +mesh.dimension=[100, 100, 100] +mesh_tally = openmc.Tally(name="flux") +mesh_tally.filters=[openmc.MeshFilter(mesh)] +mesh_tally.scores = ["flux"] +tallies = openmc.Tallies([mesh_tally]) + +# # builds the openmc model +my_model = openmc.Model(materials=materials, geometry=geometry, settings=settings, tallies=tallies) + +# # starts the simulation +output_file = my_model.run() + +# # loads up the output file from the simulation +statepoint = openmc.StatePoint(output_file) + +my_flux_mesh_tally = statepoint.get_tally(name="flux") + +mesh.write_data_to_vtk( + datasets={'mean': my_flux_mesh_tally.mean.flatten()}, + filename='mesh_tally_results.vtk', +) \ No newline at end of file diff --git a/examples/tokamak_unstructured_mesh_simulation.py b/examples/tokamak_unstructured_mesh_simulation.py index 29ba178d..56d9b957 100644 --- a/examples/tokamak_unstructured_mesh_simulation.py +++ b/examples/tokamak_unstructured_mesh_simulation.py @@ -39,8 +39,8 @@ my_model = CadToDagmc() material_tags = ["mat1"] * 6 # as inner and outer layers are one solid there are only 6 solids in model my_model.add_cadquery_object(cadquery_object=my_reactor, material_tags=material_tags) -my_model.export_dagmc_h5m_file(min_mesh_size=3.0, max_mesh_size=20.0, filename="dagmc.h5m") -my_model.export_unstructured_mesh_file(min_mesh_size=3.0, max_mesh_size=20.0, filename="tokamak_minimal.vtk") +my_model.export_dagmc_h5m_file(min_mesh_size=10.0, max_mesh_size=20.0, filename="dagmc.h5m") +my_model.export_unstructured_mesh_file(min_mesh_size=10.0, max_mesh_size=20.0, filename="tokamak_minimal.vtk") h5m_filename = "dagmc.h5m" flux = transport_particles_on_h5m_geometry(