From 0556db0c7f4a89a31da88aa30db684e673bce0d8 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 2 Apr 2024 02:12:56 +0100 Subject: [PATCH] draft solution --- ...strucutred_mesh_R2S_shut_down_dose_rate.py | 193 ++++++++++++++++++ 1 file changed, 193 insertions(+) create mode 100644 tasks/task_18_CAD_shut_down_dose_rate/1_unstrucutred_mesh_R2S_shut_down_dose_rate.py diff --git a/tasks/task_18_CAD_shut_down_dose_rate/1_unstrucutred_mesh_R2S_shut_down_dose_rate.py b/tasks/task_18_CAD_shut_down_dose_rate/1_unstrucutred_mesh_R2S_shut_down_dose_rate.py new file mode 100644 index 0000000..240c717 --- /dev/null +++ b/tasks/task_18_CAD_shut_down_dose_rate/1_unstrucutred_mesh_R2S_shut_down_dose_rate.py @@ -0,0 +1,193 @@ +from cad_to_dagmc import CadToDagmc +import cadquery as cq +import openmc +from matplotlib.colors import LogNorm + + +# makes a CAD geometry to use for the neutronics geometry +s = cq.Workplane("XY") +sPnts = [ + (2.75, 1.5), + (2.5, 1.75), + (2.0, 1.5), + (1.5, 1.0), + (1.0, 1.25), + (0.5, 1.0), + (0, 1.0), +] +r = s.lineTo(3.0, 0).lineTo(3.0, 1.0).spline(sPnts, includeCurrent=True).close() +result = r.extrude(0.5) + +my_model = CadToDagmc() + +my_model.add_cadquery_object(result) + +# this makes the tet mesh used for the unstructured mesh tally which is overlaid on the geometry +my_model.export_unstructured_mesh_file(filename="umesh.h5m", max_mesh_size=1, min_mesh_size=0.1) + +# this makes the surface mesh used for the material volume +my_model.export_dagmc_h5m_file(filename="dagmc.h5m", max_mesh_size=1, min_mesh_size=0.1) + +# the unstructured mesh to overlay on the DAGMC geometry +umesh = openmc.UnstructuredMesh("umesh.h5m", library="moab") + +# filters for use inf the tally +mesh_filter = openmc.MeshFilter(umesh) +neutron_particle_filter = openmc.ParticleFilter(['neutron']) +energy_filter = openmc.EnergyFilter.from_group_structure('CCFE-709') + +# sets up a neutron spectra tally on the unstructured mesh +tally = openmc.Tally(name="unstructured_mesh_neutron_spectra_tally") +tally.filters = [mesh_filter, energy_filter, neutron_particle_filter] +tally.scores = ["flux"] +tally.estimator = "tracklength" +my_tallies = openmc.Tallies([tally]) + +my_material = openmc.Material() +my_material.add_nuclide("H1", 1, percent_type="ao") +my_material.set_density("g/cm3", 0.001) +my_materials = openmc.Materials([my_material]) + +universe = openmc.DAGMCUniverse("dagmc.h5m").bounded_universe() +my_geometry = openmc.Geometry(universe) + +my_settings = openmc.Settings() +my_settings.batches = 10 +my_settings.inactive = 0 +my_settings.particles = 5000 +my_settings.run_mode = "fixed source" + +# Create a DT point source +my_source = openmc.IndependentSource() +my_source.space = openmc.stats.Point((0, 0, 0)) +my_source.angle = openmc.stats.Isotropic() +my_source.energy = openmc.stats.Discrete([14e6], [1]) +my_settings.source = my_source + +model = openmc.model.Model(my_geometry, my_materials, my_settings, my_tallies) +sp_filename = model.run() + + +sp = openmc.StatePoint(sp_filename) + +tally_result = sp.get_tally(name="unstrucutred_mesh_tally") + +# normally with regular meshes I would get the mesh from the tally +# but with unstrucutred meshes the tally does not contain the mesh +# however we can get it from the statepoint file +# umesh = tally_result.find_filter(openmc.MeshFilter) +umesh_from_sp = sp.meshes[1] + +# these trigger internal code in the mesh object so that its centroids and volumes become known. +# centroids and volumes are needed for the get_values and write_data_to_vtk steps +centroids = umesh_from_sp.centroids +mesh_vols = umesh_from_sp.volumes + +flux_in_each_group_for_each_voxel = tally_result.get_values(scores=["flux"], value="mean") + +all_sources = [] + +for flux_in_each_group, mesh_vol in zip(flux_in_each_group_for_each_voxel, mesh_vols): + micro_xs = openmc.deplete.MicroXS.from_multigroup_flux( + energies='CCFE-709', + multigroup_flux=flux_in_each_group, + temperature=294, # endf 8.0 has ['1200K', '2500K', '250K', '294K', '600K', '900K'] + chain_file= openmc.config['chain_file'], + nuclides=my_material.get_nuclides() + ) + + # constructing the operator, note we pass in the flux and micro xs + operator = openmc.deplete.IndependentOperator( + materials=my_materials, + fluxes=[sum(flux_in_each_group)*my_material.volume], # Flux in each group in [n-cm/src] for each domain + micros=[micro_xs], + reduce_chain=True, # reduced to only the isotopes present in depletable materials and their possible progeny + reduce_chain_level=5, + normalization_mode="source-rate" + ) + + integrator = openmc.deplete.PredictorIntegrator( + operator=operator, + timesteps=[5, 60, 60], + source_rates=[1e20, 0 , 0], # a 5 second pulse of neutrons followed by 120 seconds of decay + timestep_units='s' + ) + + integrator.integrate() + + results = openmc.deplete.Results.from_hdf5("depletion_results.h5") + last_time_step=result[-1] + activated_material = last_time_step.get_material(my_material.id) + + activated_material.volume = mesh_vol + + energy = activated_material.get_decay_photon_energy( + clip_tolerance = 1e-6, # cuts out a small fraction of the very low energy (and hence negligible dose contribution) photons + units = 'Bq', + ) + strength = energy.integral() + + +mesh_source = openmc.MeshSource( + mesh=umesh_from_sp, + sources=all_sources, +) + +my_gamma_settings = openmc.Settings() +my_gamma_settings.run_mode = "fixed source" +my_gamma_settings.batches = 100 +my_gamma_settings.particles = 100000 +my_gamma_settings.source = mesh_source + +energies, pSv_cm2 = openmc.data.dose_coefficients(particle="photon", geometry="AP") +dose_filter = openmc.EnergyFunctionFilter( + energies, pSv_cm2, interpolation="cubic" # interpolation method recommended by ICRP +) + +regularmesh = openmc.RegularMesh().from_domain( + my_geometry, + dimension=[10, 10, 10], # 10 voxels in each axis direction (x, y, z) +) + +particle_filter = openmc.ParticleFilter(["photon"]) +mesh_filter = openmc.MeshFilter(regularmesh) +dose_tally = openmc.Tally() +dose_tally.filters = [mesh_filter, dose_filter, particle_filter] +dose_tally.scores = ["flux"] +dose_tally.name = "photon_dose_on_mesh" +tallies = openmc.Tallies([dose_tally]) + + +model_gamma = openmc.Model(my_geometry, my_materials, my_gamma_settings, tallies) + +model_gamma.run(cwd="photons") + +# You may wish to plot the dose tally on a mesh, this package makes it easy to include the geometry with the mesh tally +from openmc_regular_mesh_plotter import plot_mesh_tally +with openmc.StatePoint('photons/statepoint.100.h5') as statepoint: + photon_tally = statepoint.get_tally(name="photon_dose_on_mesh") + + # normalising this tally is a little different to other examples as the source strength has been using units of photons per second. + # tally.mean is in units of pSv-cm3/source photon. + # as source strength is in photons_per_second this changes units to pSv-/second + + # multiplication by pico_to_micro converts from (pico) pSv/s to (micro) uSv/s + # dividing by mesh voxel volume cancles out the cm3 units + # could do the mesh volume scaling on the plot and vtk functions but doing it here instead + pico_to_micro = 1e-6 + seconds_to_hours = 60*60 + scaling_factor = (seconds_to_hours * pico_to_micro) / regularmesh.volumes[0][0][0] + + plot = plot_mesh_tally( + tally=photon_tally, + basis="xz", + # score='flux', # only one tally so can make use of default here + value="mean", + colorbar_kwargs={ + 'label': "Decay photon dose [µSv/h]", + }, + norm=LogNorm(), + volume_normalization=False, # this is done in the scaling_factor + scaling_factor=scaling_factor, + ) + plot.figure.savefig(f'shut_down_dose_map_timestep.png') \ No newline at end of file