|
| 1 | +import openmc |
| 2 | +import openmc.deplete |
| 3 | +from openmc_model import openmc_model |
| 4 | +import numpy as np |
| 5 | +from pathlib import Path |
| 6 | + |
| 7 | +# Fetch model |
| 8 | +model = openmc_model() |
| 9 | + |
| 10 | +model.settings.particles = int(1e1) |
| 11 | + |
| 12 | +# r2s depletion |
| 13 | +openmc.config['chain_file'] = '/home/segantin/openmc_models/CROSS_SECTIONS/chain_endfb80_sfr.xml' |
| 14 | + |
| 15 | +# deplete all materials & collect all nuclides |
| 16 | +all_nuclides = set() |
| 17 | +for i, m in enumerate(model.materials): |
| 18 | + model.materials[i].depletable = True |
| 19 | + all_nuclides.update(m.get_nuclides()) |
| 20 | + |
| 21 | +umesh = openmc.UnstructuredMesh("libra-one_nosrc.mesh.h5m", library="moab") |
| 22 | + |
| 23 | +flux_in_each_voxel, micro_xs = openmc.deplete.get_microxs_and_flux( |
| 24 | + model=model, |
| 25 | + domains=umesh, |
| 26 | + energies=[0, 30e6], # one energy bin from 0 to 30MeV |
| 27 | + chain_file=openmc.config['chain_file'], |
| 28 | + # needed otherwise the statepoint file is produced in an unknown temporary directory |
| 29 | + run_kwargs={'cwd': '.'}, |
| 30 | + nuclides=list(all_nuclides) # Convert set to list |
| 31 | +) |
| 32 | + |
| 33 | +print('\\DONE 1!\\') |
| 34 | +# Read in the unstructured from the statepoint, this contains additional information (centroids and volumes) compared to the umesh object |
| 35 | +sp_filename = f'statepoint.100.h5' |
| 36 | +sp = openmc.StatePoint(sp_filename) |
| 37 | +print('\\DONE 2!\\') |
| 38 | +# normally with regular meshes I would get the mesh from the tally |
| 39 | +# but with unstructured meshes the tally does not contain the mesh |
| 40 | +# however we can get it from the statepoint file |
| 41 | +umesh_from_sp = sp.meshes[umesh.id] |
| 42 | +# reading a unstructured mesh from the statepoint trigger internal code in the mesh |
| 43 | +# object so that its centroids and volumes become known. |
| 44 | +# centroids and volumes are needed for the get_values and write_data_to_vtk steps |
| 45 | +centroids = umesh_from_sp.centroids |
| 46 | +mesh_vols = umesh_from_sp.volumes |
| 47 | + |
| 48 | +# Calcualte the material volumes for each mesh element from the unstructured mesh |
| 49 | +mat_vols = umesh_from_sp.material_volumes(model=model, n_samples=1_000_000) |
| 50 | + |
| 51 | +print('\\DONE 3!\\') |
| 52 | + |
| 53 | +# Make a new fresh material for every tet in the unstructured mesh. |
| 54 | +# Assign the material volume for each tet as the volume is needed to deplete the material. |
| 55 | +# Get material IDs from my_materials object |
| 56 | +material_ids = [mat.id for mat in model.materials] |
| 57 | + |
| 58 | +materials_for_every_mesh_voxel = [] |
| 59 | +for i in range(len(mat_vols[material_ids[0]])): |
| 60 | + material_id = next( |
| 61 | + (mid for mid in material_ids if mat_vols[mid][i] > 0), None) |
| 62 | + if material_id is None: |
| 63 | + # handle this case if needed |
| 64 | + continue |
| 65 | + |
| 66 | + material = next(mat for mat in model.materials if mat.id == material_id) |
| 67 | + |
| 68 | + # Create a new material instance for this voxel |
| 69 | + new_mat = material.clone() |
| 70 | + new_mat.id = i |
| 71 | + # Use the volume of this material in this voxel from mat_vol |
| 72 | + new_mat.volume = mat_vols[material_id][i] |
| 73 | + materials_for_every_mesh_voxel.append(new_mat) |
| 74 | + |
| 75 | +print('\\DONE 4!\\') |
| 76 | +# Define irradiation and cooling time steps. |
| 77 | +# Set source rates to zero during decay-only steps. |
| 78 | +timesteps = np.diff(np.logspace(0.1, 7, num=10)) |
| 79 | +timesteps = np.insert(timesteps, 0, 10) # [s] add 10s of shot time |
| 80 | +source_rates = np.zeros(len(timesteps)) |
| 81 | +# n/s # [n/s] initial source rate // Full power H-mode DD pulse is 7E18 n/s |
| 82 | +source_rates[0] = 1e9 # [n/s] |
| 83 | + |
| 84 | +# Perform the activation / depletion / transmutation of all the materials |
| 85 | +# constructing the operator, note we pass in the flux and micro xs |
| 86 | +operator = openmc.deplete.IndependentOperator( |
| 87 | + materials=openmc.Materials(materials_for_every_mesh_voxel), |
| 88 | + # Flux in each group in [n-cm/src] for each domain |
| 89 | + fluxes=[flux[0] for flux in flux_in_each_voxel], |
| 90 | + micros=micro_xs, |
| 91 | + reduce_chain_level=5, |
| 92 | + normalization_mode="source-rate" |
| 93 | +) |
| 94 | + |
| 95 | +integrator = openmc.deplete.PredictorIntegrator( |
| 96 | + operator=operator, |
| 97 | + timesteps=timesteps, |
| 98 | + # a 5 second pulse of neutrons followed by 120 seconds of decay |
| 99 | + source_rates=source_rates, |
| 100 | + timestep_units='s' |
| 101 | +) |
| 102 | + |
| 103 | +integrator.integrate() |
| 104 | + |
| 105 | +print('\\DONE 5!\\') |
| 106 | + |
| 107 | +# Make a decay photon flux tally on a regular mesh for the photon / gamma dose. |
| 108 | +# The tallies will be used in each of the gamma simulations to see the shutdown dose |
| 109 | +# Include regular mesh photon flux tally |
| 110 | +mesh = openmc.RegularMesh().from_domain( |
| 111 | + model.geometry, |
| 112 | + dimension=[50, 50, 50], |
| 113 | + # 100 voxels in x and y axis directions and 1 voxel in z as we want a xy plot |
| 114 | +) |
| 115 | + |
| 116 | +rm_filter = openmc.MeshFilter(mesh) |
| 117 | +photon_filter = openmc.ParticleFilter(['photon']) |
| 118 | +tally = openmc.Tally(name='gflux_rm') |
| 119 | +tally.filters = [photon_filter, rm_filter] |
| 120 | +tally.scores = ['flux'] |
| 121 | +model.tallies.append(tally) |
| 122 | + |
| 123 | +print('\\DONE 6!\\') |
| 124 | +# We will collect the gamma source for all cooling time steps |
| 125 | +# Extract all the materials and get their gamma emission spectrum |
| 126 | +# Turn these gamma spectra into source terms for later use |
| 127 | +results = openmc.deplete.Results.from_hdf5("depletion_results.h5") |
| 128 | + |
| 129 | +all_mesh_sources = [] |
| 130 | +# skip the first time step as it is the irradiation step |
| 131 | +for i_cool in range(1, len(timesteps)): |
| 132 | + all_sources = [] |
| 133 | + for i, mesh_vol in enumerate(mesh_vols): |
| 134 | + material_id = str(i) |
| 135 | + |
| 136 | + activated_material = results[i_cool].get_material(material_id) |
| 137 | + activated_material.volume = mesh_vol |
| 138 | + energy = activated_material.get_decay_photon_energy( |
| 139 | + clip_tolerance=1e-6, |
| 140 | + units='Bq', |
| 141 | + ) |
| 142 | + |
| 143 | + if energy: |
| 144 | + strength = energy.integral() |
| 145 | + # for the strength == None case |
| 146 | + else: |
| 147 | + strength = 0 |
| 148 | + |
| 149 | + my_source = openmc.IndependentSource( |
| 150 | + energy=energy, |
| 151 | + particle="photon", |
| 152 | + strength=strength, |
| 153 | + # constraints={'domains':my_material} |
| 154 | + ) |
| 155 | + |
| 156 | + all_sources.append(my_source) |
| 157 | + |
| 158 | + # Make a mesh source out of the IndependentSource just made in the inner loop |
| 159 | + mesh_source = openmc.MeshSource( |
| 160 | + mesh=umesh_from_sp, |
| 161 | + sources=all_sources, |
| 162 | + ) |
| 163 | + |
| 164 | + all_mesh_sources.append(mesh_source) |
| 165 | + |
| 166 | +# Makes and runs a simulation model for each time meshsource that has been made. |
| 167 | +# Make simulation settings for the gamma transport simulation |
| 168 | +my_gamma_settings = openmc.Settings() |
| 169 | +my_gamma_settings.run_mode = "fixed source" |
| 170 | +my_gamma_settings.batches = 100 |
| 171 | +my_gamma_settings.particles = int(1e3) |
| 172 | +my_gamma_settings.output = {'summary': False} |
| 173 | +my_gamma_settings.photon_transport = True |
| 174 | + |
| 175 | +all_gamma_sp_filename = [] |
| 176 | + |
| 177 | +for mesh_source in all_mesh_sources: |
| 178 | + |
| 179 | + my_gamma_settings.source = mesh_source |
| 180 | + |
| 181 | + # here we use the same pristine materials from before neutron irradaiton as the burnup is low |
| 182 | + # and the materials have not changed much so they would not perterb the neutron spectrum significantly |
| 183 | + # you could also use the activated materials from the depletion results but this would significantly slow the simulation down |
| 184 | + model_gamma = openmc.Model( |
| 185 | + model.geometry, model.materials, my_gamma_settings, model.tallies) |
| 186 | + |
| 187 | + # Make the model for the gamma / photon transport and run the simulation |
| 188 | + # a folder will be made for each photon transposrt, XML files will be saves there as well as the plot |
| 189 | + # Create directories with parents=True to ensure parent directories are created |
| 190 | + output_dir = Path(f"photons/photon_at_time_{i_cool}") |
| 191 | + output_dir.mkdir(parents=True, exist_ok=True) |
| 192 | + |
| 193 | + gamma_sp_filename = model_gamma.run(cwd=str(output_dir)) |
| 194 | + all_gamma_sp_filename.append(gamma_sp_filename) |
| 195 | + |
| 196 | + |
| 197 | +print('\\DONE 7!\\') |
0 commit comments