Line data Source code
1 : /********************************************************************/
2 : /* SOFTWARE COPYRIGHT NOTIFICATION */
3 : /* Cardinal */
4 : /* */
5 : /* (c) 2021 UChicago Argonne, LLC */
6 : /* ALL RIGHTS RESERVED */
7 : /* */
8 : /* Prepared by UChicago Argonne, LLC */
9 : /* Under Contract No. DE-AC02-06CH11357 */
10 : /* With the U. S. Department of Energy */
11 : /* */
12 : /* Prepared by Battelle Energy Alliance, LLC */
13 : /* Under Contract No. DE-AC07-05ID14517 */
14 : /* With the U. S. Department of Energy */
15 : /* */
16 : /* See LICENSE for full restrictions */
17 : /********************************************************************/
18 :
19 : #ifdef ENABLE_OPENMC_COUPLING
20 : #include "MeshTally.h"
21 : #include "DisplacedProblem.h"
22 :
23 : #include "libmesh/replicated_mesh.h"
24 :
25 : registerMooseObject("CardinalApp", MeshTally);
26 :
27 : InputParameters
28 1089 : MeshTally::validParams()
29 : {
30 1089 : auto params = TallyBase::validParams();
31 1089 : params.addClassDescription("A class which implements unstructured mesh tallies.");
32 2178 : params.addParam<std::string>("mesh_template",
33 : "Mesh tally template for OpenMC when using mesh tallies; "
34 : "at present, this mesh must exactly match the mesh used in the "
35 : "[Mesh] block because a one-to-one copy "
36 : "is used to get OpenMC's tally results on the [Mesh].");
37 2178 : params.addParam<Point>("mesh_translation",
38 : "Coordinate to which this mesh should be "
39 : "translated. Units must match those used to define the [Mesh].");
40 :
41 : // The index of this tally into an array of mesh translations. Defaults to zero.
42 1089 : params.addPrivateParam<unsigned int>("instance", 0);
43 :
44 1089 : return params;
45 0 : }
46 :
47 726 : MeshTally::MeshTally(const InputParameters & parameters)
48 : : TallyBase(parameters),
49 1982 : _mesh_translation(isParamValid("mesh_translation") ? getParam<Point>("mesh_translation")
50 : : Point(0.0, 0.0, 0.0)),
51 1448 : _instance(getParam<unsigned int>("instance")),
52 2830 : _use_dof_map(_is_adaptive || isParamValid("block"))
53 : {
54 : bool nu_scatter =
55 724 : std::find(_tally_score.begin(), _tally_score.end(), "nu-scatter") != _tally_score.end();
56 :
57 : // Error check the estimators.
58 1448 : if (isParamValid("estimator"))
59 : {
60 14 : if (_estimator == openmc::TallyEstimator::TRACKLENGTH)
61 2 : paramError("estimator",
62 : "Tracklength estimators are currently incompatible with mesh tallies!");
63 : }
64 : else
65 1420 : _estimator = nu_scatter ? openmc::TallyEstimator::ANALOG : openmc::TallyEstimator::COLLISION;
66 :
67 : // Error check the mesh template.
68 724 : if (_openmc_problem.getMooseMesh().getMesh().allow_renumbering() &&
69 2 : !_openmc_problem.getMooseMesh().getMesh().is_replicated())
70 2 : mooseError(
71 : "Mesh tallies currently require 'allow_renumbering = false' to be set in the [Mesh]!");
72 :
73 1440 : if (isParamValid("mesh_template"))
74 : {
75 641 : if (_is_adaptive)
76 2 : paramError("mesh_template",
77 : "Adaptivity is not supported when loading a mesh from 'mesh_template'!");
78 :
79 1278 : if (isParamValid("block"))
80 2 : paramError("block",
81 : "Block restriction is currently not supported for mesh tallies which load a "
82 : "mesh from a file!");
83 :
84 637 : if (_openmc_problem.useDisplaced())
85 1 : paramError("mesh_template",
86 : "Tallying on a file-based mesh is not supported for moving-mesh cases as there is "
87 : "not a mechanism to update the mesh geometry. You can use a mesh tally for moving "
88 : "mesh cases by instead tallying directly on the [Mesh]. Simply do not provide the "
89 : "'mesh_template' parameter.");
90 :
91 1272 : _mesh_template_filename = &getParam<std::string>("mesh_template");
92 : }
93 : else
94 : {
95 79 : if (std::abs(_openmc_problem.scaling() - 1.0) > 1e-6)
96 2 : mooseError("Directly tallying on the [Mesh] is only supported for 'scaling' of unity. "
97 : "Instead, please make a file containing your tally mesh and set it with "
98 : "'mesh_template'. You can generate a mesh file corresponding to the [Mesh] "
99 2 : "by running:\n\ncardinal-opt -i " +
100 2 : _app.getFileName() + " --mesh-only");
101 :
102 : // for distributed meshes, each rank only owns a portion of the mesh information, but
103 : // OpenMC wants the entire mesh to be available on every rank. We might be able to add
104 : // this feature in the future, but will need to investigate
105 77 : if (!_openmc_problem.getMooseMesh().getMesh().is_replicated())
106 0 : mooseError("Directly tallying on the [Mesh] block by OpenMC is not yet supported "
107 : "for distributed meshes!");
108 :
109 154 : if (isParamValid("mesh_translation"))
110 0 : paramError("mesh_translation",
111 : "The mesh filter cannot be translated if directly tallying on the mesh "
112 : "provided in the [Mesh] block!");
113 : }
114 :
115 : /**
116 : * If the instance isn't zero this variable is a translated mesh tally. It will accumulate it's
117 : * scores in a different set of variables (the auxvars which are added by the first tally in a
118 : * sequence of mesh tallies), and so it doesn't need to create any auxvars.
119 : */
120 713 : if (_instance != 0)
121 356 : _tally_name = std::vector<std::string>();
122 713 : }
123 :
124 : std::pair<unsigned int, openmc::Filter *>
125 722 : MeshTally::spatialFilter()
126 : {
127 : // Create the OpenMC mesh which will be tallied on.
128 722 : if (!_mesh_template_filename)
129 : {
130 : auto msh =
131 85 : dynamic_cast<const libMesh::ReplicatedMesh *>(_openmc_problem.getMooseMesh().getMeshPtr());
132 85 : if (!msh)
133 0 : mooseError("Internal error: The mesh is not a replicated mesh.");
134 :
135 : // Adaptivity and block restriction both require a map from the element subsets we want to
136 : // tally on to the full mesh.
137 85 : if (_use_dof_map)
138 : {
139 50 : _bin_to_element_mapping.clear();
140 :
141 : auto begin = _tally_blocks.size() > 0
142 50 : ? msh->active_subdomain_set_elements_begin(_tally_blocks)
143 100 : : msh->active_elements_begin();
144 50 : auto end = _tally_blocks.size() > 0 ? msh->active_subdomain_set_elements_end(_tally_blocks)
145 50 : : msh->active_elements_end();
146 182692 : for (const auto & old_elem : libMesh::as_range(begin, end))
147 182642 : _bin_to_element_mapping.push_back(old_elem->id());
148 :
149 : _bin_to_element_mapping.shrink_to_fit();
150 : }
151 :
152 : // When block restriction is active we need to create a copy of the mesh which only contains
153 : // elements in the desired blocks.
154 85 : if (_tally_blocks.size() > 0)
155 : {
156 170 : _libmesh_mesh_copy = std::make_unique<libMesh::ReplicatedMesh>(
157 85 : _openmc_problem.comm(), _openmc_problem.getMooseMesh().dimension());
158 :
159 170 : msh->create_submesh(*_libmesh_mesh_copy.get(),
160 170 : msh->active_subdomain_set_elements_begin(_tally_blocks),
161 170 : msh->active_subdomain_set_elements_end(_tally_blocks));
162 : _libmesh_mesh_copy->allow_find_neighbors(true);
163 : _libmesh_mesh_copy->allow_renumbering(false);
164 85 : _libmesh_mesh_copy->prepare_for_use();
165 :
166 85 : openmc::model::meshes.emplace_back(
167 170 : std::make_unique<openmc::LibMesh>(*_libmesh_mesh_copy.get(), _openmc_problem.scaling()));
168 : }
169 : else
170 : {
171 0 : if (_is_adaptive)
172 0 : openmc::model::meshes.emplace_back(std::make_unique<openmc::AdaptiveLibMesh>(
173 0 : _openmc_problem.getMooseMesh().getMesh(), _openmc_problem.scaling()));
174 : else
175 0 : openmc::model::meshes.emplace_back(std::make_unique<openmc::LibMesh>(
176 0 : _openmc_problem.getMooseMesh().getMesh(), _openmc_problem.scaling()));
177 : }
178 : }
179 : else
180 637 : openmc::model::meshes.emplace_back(
181 1274 : std::make_unique<openmc::LibMesh>(*_mesh_template_filename, _openmc_problem.scaling()));
182 :
183 722 : _mesh_index = openmc::model::meshes.size() - 1;
184 722 : _mesh_template =
185 722 : dynamic_cast<openmc::UnstructuredMesh *>(openmc::model::meshes[_mesh_index].get());
186 :
187 : // by setting the ID to -1, OpenMC will automatically detect the next available ID
188 722 : _mesh_template->set_id(-1);
189 722 : _mesh_template->output_ = false;
190 :
191 722 : _mesh_filter = dynamic_cast<openmc::MeshFilter *>(openmc::Filter::create("mesh"));
192 722 : _mesh_filter->set_mesh(_mesh_index);
193 722 : _mesh_filter->set_translation({_mesh_translation(0), _mesh_translation(1), _mesh_translation(2)});
194 :
195 : // Validate the mesh filters to make sure we can run a copy transfer to the [Mesh].
196 722 : checkMeshTemplateAndTranslations();
197 :
198 716 : return std::make_pair(openmc::model::tally_filters.size() - 1, _mesh_filter);
199 : }
200 :
201 : void
202 33 : MeshTally::resetTally()
203 : {
204 33 : TallyBase::resetTally();
205 :
206 : // Erase the OpenMC mesh.
207 33 : openmc::model::meshes.erase(openmc::model::meshes.begin() + _mesh_index);
208 33 : }
209 :
210 : Real
211 1263 : MeshTally::storeResultsInner(const std::vector<unsigned int> & var_numbers,
212 : unsigned int local_score,
213 : unsigned int global_score,
214 : std::vector<xt::xtensor<double, 1>> tally_vals,
215 : bool norm_by_src_rate)
216 : {
217 : Real total = 0.0;
218 :
219 1263 : unsigned int mesh_offset = _instance * _mesh_filter->n_bins();
220 2772 : for (unsigned int ext_bin = 0; ext_bin < _num_ext_filter_bins; ++ext_bin)
221 : {
222 1214821 : for (decltype(_mesh_filter->n_bins()) e = 0; e < _mesh_filter->n_bins(); ++e)
223 : {
224 1213312 : Real unnormalized_tally = tally_vals[local_score](ext_bin * _mesh_filter->n_bins() + e);
225 :
226 : // divide each tally by the volume that it corresponds to in MOOSE
227 : // because we will apply it as a volumetric tally (per unit volume).
228 : // Because we require that the mesh template has units of cm based on the
229 : // mesh constructors in OpenMC, we need to adjust the division
230 1213312 : Real volumetric_tally = unnormalized_tally;
231 1213312 : volumetric_tally *= norm_by_src_rate
232 1213312 : ? _openmc_problem.tallyMultiplier(global_score) /
233 949860 : _mesh_template->volume(e) * _openmc_problem.scaling() *
234 : _openmc_problem.scaling() * _openmc_problem.scaling()
235 : : 1.0;
236 1213312 : total += _ext_bins_to_skip[ext_bin] ? 0.0 : unnormalized_tally;
237 :
238 1213312 : auto var = var_numbers[_num_ext_filter_bins * local_score + ext_bin];
239 1213312 : auto elem_id = _use_dof_map ? _bin_to_element_mapping[e] : mesh_offset + e;
240 1213312 : fillElementalAuxVariable(var, {elem_id}, volumetric_tally);
241 : }
242 : }
243 :
244 1263 : return total;
245 : }
246 :
247 : void
248 722 : MeshTally::checkMeshTemplateAndTranslations()
249 : {
250 : // we can do some rudimentary checking on the mesh template by comparing the centroid
251 : // coordinates compared to centroids in the [Mesh] (because right now, we just doing a simple
252 : // copy transfer that necessitates the meshes to have the same elements in the same order). In
253 : // other words, you might have two meshes that represent the same geometry, the element ordering
254 : // could be different.
255 722 : unsigned int mesh_offset = _instance * _mesh_filter->n_bins();
256 401370 : for (int e = 0; e < _mesh_filter->n_bins(); ++e)
257 : {
258 400654 : auto elem_id = _use_dof_map ? _bin_to_element_mapping[e] : mesh_offset + e;
259 400654 : auto elem_ptr = _openmc_problem.getMooseMesh().queryElemPtr(elem_id);
260 :
261 : // if element is not on this part of the distributed mesh, skip it
262 400654 : if (!elem_ptr)
263 60654 : continue;
264 :
265 340000 : const auto pt = _mesh_template->centroid(e);
266 340000 : Point centroid_template = {pt[0], pt[1], pt[2]};
267 :
268 : // The translation applied in OpenMC isn't actually registered in the mesh itself;
269 : // it is always added on to the point, so we need to do the same here
270 : centroid_template += _mesh_translation;
271 :
272 : // because the mesh template and [Mesh] may be in different units, we need
273 : // to adjust the [Mesh] by the scaling factor before doing a comparison.
274 340000 : Point centroid_mesh = elem_ptr->vertex_average() * _openmc_problem.scaling();
275 :
276 : // if the centroids are the same except for a factor of 'scaling', then we can
277 : // guess that the mesh_template is probably not in units of centimeters
278 340000 : if (_openmc_problem.hasScaling())
279 : {
280 : // if scaling was applied correctly, then each calculation of 'scaling' here should equal 1.
281 : // Otherwise, if they're all the same, then 'scaling_x' is probably the factor by which the
282 : // mesh_template needs to be multiplied, so we can print a helpful error message
283 : bool incorrect_scaling = true;
284 147240 : for (unsigned int j = 0; j < OpenMCCellAverageProblem::DIMENSION; ++j)
285 : {
286 110430 : Real scaling = centroid_mesh(j) / centroid_template(j);
287 110430 : incorrect_scaling = incorrect_scaling && !MooseUtils::absoluteFuzzyEqual(scaling, 1.0);
288 : }
289 :
290 36810 : if (incorrect_scaling)
291 2 : paramError("mesh_template",
292 : "The centroids of the 'mesh_template' differ from the "
293 2 : "centroids of the [Mesh] by a factor of " +
294 2 : Moose::stringify(centroid_mesh(0) / centroid_template(0)) +
295 : ".\nDid you forget that the 'mesh_template' must be in "
296 : "the same units as the [Mesh]?");
297 : }
298 :
299 : // check if centroids are the same
300 : bool different_centroids = false;
301 1359992 : for (unsigned int j = 0; j < OpenMCCellAverageProblem::DIMENSION; ++j)
302 1019994 : different_centroids = different_centroids ||
303 : !MooseUtils::absoluteFuzzyEqual(centroid_mesh(j), centroid_template(j));
304 :
305 339998 : if (different_centroids)
306 4 : paramError(
307 : "mesh_template",
308 4 : "Centroid for element " + Moose::stringify(elem_id) +
309 8 : " in the [Mesh] (cm): " + _openmc_problem.printPoint(centroid_mesh) +
310 4 : "\ndoes not match centroid for element " + Moose::stringify(e) +
311 8 : " in the 'mesh_template' with instance " + Moose::stringify(_instance) +
312 8 : " (cm): " + _openmc_problem.printPoint(centroid_template) +
313 : "!\n\nThe copy transfer requires that the [Mesh] and 'mesh_template' be identical.");
314 : }
315 716 : }
316 : #endif
|
