Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : // MOOSE includes
11 : #include "SampledOutput.h"
12 : #include "FEProblem.h"
13 : #include "DisplacedProblem.h"
14 : #include "MooseApp.h"
15 : #include "MoosePartitioner.h"
16 :
17 : #include "libmesh/distributed_mesh.h"
18 : #include "libmesh/equation_systems.h"
19 : #include "libmesh/mesh_function.h"
20 : #include "libmesh/explicit_system.h"
21 :
22 : using namespace libMesh;
23 :
24 : InputParameters
25 163904 : SampledOutput::validParams()
26 : {
27 :
28 : // Get the parameters from the parent object
29 163904 : InputParameters params = AdvancedOutput::validParams();
30 491712 : params.addParam<unsigned int>("refinements",
31 327808 : 0,
32 : "Number of uniform refinements for oversampling "
33 : "(refinement levels beyond any level of "
34 : "refinements already applied on the regular mesh)");
35 163904 : params.addParam<Point>("position",
36 : "Set a positional offset, this vector will get added to the "
37 : "nodal coordinates to move the domain.");
38 163904 : params.addParam<MeshFileName>("file", "The name of the mesh file to read, for oversampling");
39 163904 : params.addParam<std::vector<SubdomainName>>(
40 : "sampling_blocks", "The list of blocks to restrict the mesh sampling to");
41 491712 : params.addParam<bool>(
42 : "serialize_sampling",
43 327808 : true,
44 : "If set to true, all sampled output (see sampling parameters) will be done "
45 : "on rank 0. This option is useful to debug suspected parallel output issues");
46 :
47 : // **** DEPRECATED PARAMETERS ****
48 491712 : params.addDeprecatedParam<bool>("append_oversample",
49 327808 : false,
50 : "Append '_oversample' to the output file base",
51 : "This parameter is deprecated. To append '_oversample' utilize "
52 : "the output block name or the 'file_base'");
53 :
54 : // 'Oversampling' Group
55 163904 : params.addParamNamesToGroup("refinements position file sampling_blocks serialize_sampling",
56 : "Modified Mesh Sampling");
57 :
58 163904 : return params;
59 0 : }
60 :
61 37913 : SampledOutput::SampledOutput(const InputParameters & parameters)
62 : : AdvancedOutput(parameters),
63 37901 : _refinements(getParam<unsigned int>("refinements")),
64 37901 : _using_external_sampling_file(isParamValid("file")),
65 37901 : _change_position(isParamValid("position")),
66 37769 : _use_sampled_output(_refinements > 0 || _using_external_sampling_file ||
67 75670 : isParamValid("sampling_blocks") || _change_position),
68 37901 : _position(_change_position ? getParam<Point>("position") : Point()),
69 37901 : _sampling_mesh_changed(true),
70 37901 : _mesh_subdomains_match(true),
71 113715 : _serialize(getParam<bool>("serialize_sampling"))
72 : {
73 37901 : }
74 :
75 : void
76 37438 : SampledOutput::initialSetup()
77 : {
78 37438 : AdvancedOutput::initialSetup();
79 :
80 : // Creates and initializes the sampling mesh
81 37434 : initSample();
82 37434 : }
83 :
84 : void
85 2609699 : SampledOutput::outputStep(const ExecFlagType & type)
86 : {
87 : // Output is not allowed
88 2609699 : if (!_allow_output && type != EXEC_FORCED)
89 51111 : return;
90 :
91 : // If recovering disable output of initial condition, it was already output
92 2558588 : if (type == EXEC_INITIAL && _app.isRecovering())
93 0 : return;
94 :
95 : // Return if the current output is not on the desired interval
96 2558588 : if (type != EXEC_FINAL && !onInterval())
97 47893 : return;
98 :
99 : // store current simulation time
100 2510695 : _last_output_simulation_time = _time;
101 :
102 : // set current type
103 2510695 : _current_execute_flag = type;
104 :
105 : // Call the output method
106 2510695 : if (shouldOutput())
107 : {
108 172199 : TIME_SECTION("outputStep", 2, "Outputting Step");
109 172199 : updateSample();
110 172199 : output();
111 172199 : }
112 :
113 2510695 : _current_execute_flag = EXEC_NONE;
114 : }
115 :
116 36481 : SampledOutput::~SampledOutput()
117 : {
118 : // TODO: Remove once libmesh Issue #1184 is fixed
119 36481 : _sampling_es.reset();
120 36481 : _sampling_mesh_ptr.reset();
121 36481 : }
122 :
123 : void
124 5194 : SampledOutput::meshChanged()
125 : {
126 5194 : _sampling_mesh_changed = true;
127 5194 : }
128 :
129 : void
130 37434 : SampledOutput::initSample()
131 : {
132 : // Perform the mesh cloning, if needed
133 37434 : if (!_use_sampled_output)
134 37241 : return;
135 :
136 193 : cloneMesh();
137 :
138 : // Re-position the sampling mesh
139 193 : if (_change_position)
140 2508 : for (auto & node : _mesh_ptr->getMesh().node_ptr_range())
141 2508 : *node += _position;
142 :
143 : // Perform the mesh refinement
144 193 : if (_refinements > 0)
145 : {
146 132 : MeshRefinement mesh_refinement(_mesh_ptr->getMesh());
147 :
148 : // We want original and refined partitioning to match so we can
149 : // query from one to the other safely on distributed meshes.
150 132 : _mesh_ptr->getMesh().skip_partitioning(true);
151 132 : mesh_refinement.uniformly_refine(_refinements);
152 :
153 : // Note that nodesets are not propagated with mesh refinement, unless you built the nodesets
154 : // from the sidesets again, which is what happens for the regular mesh with initial refinement
155 132 : }
156 :
157 : // We can't allow renumbering if we want to output multiple time
158 : // steps to the same Exodus file
159 193 : _mesh_ptr->getMesh().allow_renumbering(false);
160 :
161 : // This should be called after changing the mesh (block restriction for example)
162 193 : if (_change_position || (_refinements > 0) || isParamValid("sampling_blocks"))
163 182 : _sampling_mesh_ptr->meshChanged();
164 :
165 : // Create the new EquationSystems
166 193 : _sampling_es = std::make_unique<EquationSystems>(_mesh_ptr->getMesh());
167 193 : _es_ptr = _sampling_es.get();
168 :
169 : // Reference the system from which we are copying
170 193 : EquationSystems & source_es = _problem_ptr->es();
171 :
172 : // If we're going to be copying from that system later, we need to keep its
173 : // original elements as ghost elements even if it gets grossly
174 : // repartitioned, since we can't repartition the sample mesh to
175 : // match.
176 : // FIXME: this is not enough. It assumes our initial partition of the sampling mesh
177 : // and the source mesh match. But that's usually only true in the 'refinement' case,
178 : // not with an arbitrary sampling mesh file
179 193 : DistributedMesh * dist_mesh = dynamic_cast<DistributedMesh *>(&source_es.get_mesh());
180 193 : if (dist_mesh)
181 : {
182 650 : for (auto & elem : dist_mesh->active_local_element_ptr_range())
183 650 : dist_mesh->add_extra_ghost_elem(elem);
184 : }
185 :
186 : // Initialize the _mesh_functions vector
187 193 : const auto num_systems = source_es.n_systems();
188 193 : _mesh_functions.resize(num_systems);
189 :
190 : // Keep track of the variable numbering in both regular and sampled system
191 193 : _variable_numbers_in_system.resize(num_systems);
192 :
193 : // Get the list of nodal and elemental output data
194 193 : const auto & nodal_data = getNodalVariableOutput();
195 193 : const auto & elemental_data = getElementalVariableOutput();
196 :
197 : // Loop over the number of systems
198 579 : for (const auto sys_num : make_range(num_systems))
199 : {
200 : // Reference to the current system
201 386 : const auto & source_sys = source_es.get_system(sys_num);
202 :
203 : // Add the system to the new EquationsSystems
204 386 : ExplicitSystem & dest_sys = _sampling_es->add_system<ExplicitSystem>(source_sys.name());
205 :
206 : // Loop through the variables in the System
207 386 : const auto num_vars = source_sys.n_vars();
208 386 : unsigned int num_actual_vars = 0;
209 386 : if (num_vars > 0)
210 : {
211 265 : if (_serialize)
212 239 : _serialized_solution = NumericVector<Number>::build(_communicator);
213 :
214 : // Add the variables to the system... simultaneously creating MeshFunctions for them.
215 669 : for (const auto var_num : make_range(num_vars))
216 : {
217 : // Is the variable supposed to be output?
218 404 : const auto & var_name = source_sys.variable_name(var_num);
219 404 : if (!nodal_data.count(var_name) && !elemental_data.count(var_name))
220 159 : continue;
221 :
222 : // We do what we can to preserve the block restriction
223 : const std::set<SubdomainID> * subdomains;
224 271 : std::set<SubdomainID> restricted_subdomains;
225 271 : if (_using_external_sampling_file && !_mesh_subdomains_match)
226 0 : subdomains = nullptr;
227 : else
228 : {
229 271 : subdomains = &source_sys.variable(var_num).active_subdomains();
230 : // Reduce the block restriction if the output is block restricted
231 271 : if (isParamValid("sampling_blocks") && !subdomains->empty())
232 : {
233 : const auto & sampling_blocks = _sampling_mesh_ptr->getSubdomainIDs(
234 52 : getParam<std::vector<SubdomainName>>("sampling_blocks"));
235 52 : set_intersection(subdomains->begin(),
236 : subdomains->end(),
237 : sampling_blocks.begin(),
238 : sampling_blocks.end(),
239 : std::inserter(restricted_subdomains, restricted_subdomains.begin()));
240 52 : subdomains = &restricted_subdomains;
241 :
242 : // None of the subdomains are included in the sampling, might as well skip
243 52 : if (subdomains->empty())
244 : {
245 26 : hideAdditionalVariable(nodal_data.count(var_name) ? "nodal" : "elemental", var_name);
246 26 : continue;
247 : }
248 52 : }
249 : }
250 :
251 : // We are going to add the variable, let's count it
252 245 : _variable_numbers_in_system[sys_num].push_back(var_num);
253 245 : num_actual_vars++;
254 :
255 : // Add the variable. We essentially support nodal variables and constant monomials
256 245 : const FEType & fe_type = source_sys.variable_type(var_num);
257 245 : if (isSampledAtNodes(fe_type))
258 : {
259 180 : dest_sys.add_variable(source_sys.variable_name(var_num), fe_type, subdomains);
260 180 : if (dist_mesh && !_serialize)
261 0 : paramError("serialize_sampling",
262 : "Variables sampled as nodal currently require serialization with a "
263 : "distributed mesh.");
264 : }
265 : else
266 : {
267 65 : const auto & var_name = source_sys.variable_name(var_num);
268 65 : if (fe_type != FEType(CONSTANT, MONOMIAL))
269 : {
270 0 : mooseInfoRepeated("Sampled output projects variable '" + var_name +
271 : "' onto a constant monomial");
272 0 : if (!_serialize)
273 0 : paramWarning("serialize_sampling",
274 : "Projection without serialization may fail with insufficient ghosting. "
275 : "Consider setting 'serialize_sampling' to true.");
276 : }
277 65 : dest_sys.add_variable(var_name, FEType(CONSTANT, MONOMIAL), subdomains);
278 : }
279 : // Note: we could do more, using the generic projector. But exodus output of higher order
280 : // or more exotic variables is limited anyway
281 271 : }
282 :
283 : // Size for the actual number of variables output
284 265 : _mesh_functions[sys_num].resize(num_actual_vars);
285 : }
286 : }
287 :
288 : // Initialize the newly created EquationSystem
289 193 : _sampling_es->init();
290 : }
291 :
292 : void
293 172199 : SampledOutput::updateSample()
294 : {
295 : // Do nothing if oversampling and changing position are not enabled
296 172199 : if (!_use_sampled_output)
297 171632 : return;
298 :
299 : // We need the mesh functions to extend the whole domain so we serialize both the mesh and the
300 : // solution. We need this because the partitioning of the sampling mesh may not match the
301 : // partitioning of the source mesh
302 567 : if (_serialize)
303 : {
304 543 : _problem_ptr->mesh().getMesh().gather_to_zero();
305 543 : _mesh_ptr->getMesh().gather_to_zero();
306 : }
307 :
308 : // Get a reference to actual equation system
309 567 : EquationSystems & source_es = _problem_ptr->es();
310 567 : const auto num_systems = source_es.n_systems();
311 :
312 : // Loop through each system
313 1701 : for (const auto sys_num : make_range(num_systems))
314 : {
315 1134 : if (!_mesh_functions[sys_num].empty())
316 : {
317 : // Get references to the source and destination systems
318 615 : System & source_sys = source_es.get_system(sys_num);
319 615 : System & dest_sys = _sampling_es->get_system(sys_num);
320 :
321 : // Update the solution for the sampling mesh
322 615 : if (_serialize)
323 : {
324 591 : _serialized_solution->clear();
325 591 : _serialized_solution->init(source_sys.n_dofs(), false, SERIAL);
326 : // Pull down a full copy of this vector on every processor so we can get values in
327 : // parallel
328 591 : source_sys.solution->localize(*_serialized_solution);
329 : }
330 :
331 : // Update the mesh functions
332 1278 : for (const auto var_num : index_range(_mesh_functions[sys_num]))
333 : {
334 663 : const auto original_var_num = _variable_numbers_in_system[sys_num][var_num];
335 :
336 : // If the mesh has changed, the MeshFunctions need to be re-built, otherwise simply clear
337 : // it for re-initialization
338 : // TODO: inherit from MeshChangedInterface and rebuild mesh functions on meshChanged()
339 663 : if (!_mesh_functions[sys_num][var_num] || _sampling_mesh_changed)
340 912 : _mesh_functions[sys_num][var_num] = std::make_unique<MeshFunction>(
341 : source_es,
342 304 : _serialize ? *_serialized_solution : *source_sys.solution,
343 : source_sys.get_dof_map(),
344 304 : original_var_num);
345 : else
346 359 : _mesh_functions[sys_num][var_num]->clear();
347 :
348 : // Initialize the MeshFunctions for application to the sampled solution
349 663 : _mesh_functions[sys_num][var_num]->init();
350 :
351 : // Mesh functions are still defined on the original mesh, which might not fully overlap
352 : // with the sampling mesh. We don't want to error with a libMesh assert on the out of mesh
353 : // mode
354 663 : _mesh_functions[sys_num][var_num]->enable_out_of_mesh_mode(-1e6);
355 : }
356 :
357 : // Fill solution vectors by evaluating mesh functions on sampling mesh
358 1107 : for (const auto var_num : index_range(_mesh_functions[sys_num]))
359 : {
360 : // we serialized the mesh and the solution vector, we might as well just do this only on
361 : // processor 0.
362 651 : if (_serialize && processor_id() > 0)
363 159 : break;
364 :
365 492 : const auto original_var_num = _variable_numbers_in_system[sys_num][var_num];
366 492 : const FEType & fe_type = source_sys.variable_type(original_var_num);
367 : // we use the original variable block restriction for sampling
368 492 : const auto * var_blocks = &source_sys.variable(original_var_num).active_subdomains();
369 : // NOTE: if we have overlapping domains between the sampling mesh and the source mesh
370 : // we would get a value from the source mesh domain. We could further restrict this
371 : // block restriction with the sampling mesh block restriction to prevent this.
372 :
373 : // Loop over the mesh, nodes for nodal data, elements for element data
374 492 : if (isSampledAtNodes(fe_type))
375 : {
376 792 : for (const auto & node : (_serialize ? _mesh_ptr->getMesh().node_ptr_range()
377 1131336 : : _mesh_ptr->getMesh().local_node_ptr_range()))
378 : {
379 : // Avoid working on ghosted dofs
380 1130544 : if (node->n_dofs(sys_num, var_num) &&
381 565272 : (_serialize || processor_id() == node->processor_id()))
382 : {
383 : // the node has to be within the domain of the mesh function
384 565272 : DenseVector<Real> value(1);
385 565272 : if (var_blocks->size())
386 720 : (*_mesh_functions[sys_num][var_num])(
387 720 : *node - _position, /*time*/ 0., value, var_blocks);
388 : else
389 564912 : value[0] = (*_mesh_functions[sys_num][var_num])(*node - _position);
390 :
391 565272 : if (value[0] != -1e6)
392 565272 : dest_sys.solution->set(node->dof_number(sys_num, var_num, /*comp=*/0), value[0]);
393 : else
394 0 : mooseDoOnce(mooseWarning(
395 : "Sampling at location ",
396 : *node - _position,
397 : " by process ",
398 : std::to_string(processor_id()),
399 : " was outside the problem mesh.\nThis message will not be repeated"));
400 565272 : }
401 396 : }
402 : }
403 : else
404 : {
405 96 : const auto elem_range = _serialize
406 96 : ? _mesh_ptr->getMesh().active_element_ptr_range()
407 96 : : _mesh_ptr->getMesh().active_local_element_ptr_range();
408 1716 : for (const auto & elem : elem_range)
409 : {
410 1620 : if (elem->n_dofs(sys_num, var_num) &&
411 810 : (_serialize || processor_id() == elem->processor_id()))
412 : {
413 810 : DenseVector<Real> value(1);
414 810 : if (var_blocks->size())
415 324 : (*_mesh_functions[sys_num][var_num])(
416 324 : elem->true_centroid() - _position, /*time*/ 0., value, var_blocks);
417 : else
418 648 : value[0] = (*_mesh_functions[sys_num][var_num])(elem->true_centroid() - _position);
419 :
420 810 : if (value[0] != -1e6)
421 810 : dest_sys.solution->set(elem->dof_number(sys_num, var_num, /*comp=*/0), value[0]);
422 : else
423 0 : mooseDoOnce(mooseWarning(
424 : "Sampling at location ",
425 : elem->true_centroid() - _position,
426 : " was outside the problem mesh.\nThis message will not be repeated."));
427 810 : }
428 96 : }
429 96 : }
430 : }
431 :
432 : // We modified the solution vector directly, we have to close it
433 615 : dest_sys.solution->close();
434 : }
435 : }
436 :
437 : // Set this to false so that new output files are not created, since the sampling mesh
438 : // doesn't actually change
439 567 : _sampling_mesh_changed = false;
440 : }
441 :
442 : void
443 193 : SampledOutput::cloneMesh()
444 : {
445 : // Create the new mesh from a file
446 193 : if (isParamValid("file"))
447 : {
448 11 : InputParameters mesh_params = _app.getFactory().getValidParams("FileMesh");
449 11 : mesh_params.applyParameters(parameters(), {}, true);
450 11 : mesh_params.set<bool>("nemesis") = false;
451 : _sampling_mesh_ptr =
452 11 : _app.getFactory().createUnique<MooseMesh>("FileMesh", "output_problem_mesh", mesh_params);
453 11 : _sampling_mesh_ptr->allowRecovery(false); // We actually want to reread the initial mesh
454 11 : _sampling_mesh_ptr->init();
455 11 : }
456 : // Clone the existing mesh
457 : else
458 : {
459 182 : if (_app.isRecovering())
460 3 : mooseWarning("Recovering or Restarting with oversampling may not work (especially with "
461 : "adapted meshes)!! Refs #2295");
462 182 : _sampling_mesh_ptr = _mesh_ptr->safeClone();
463 : }
464 :
465 : // Remove unspecified blocks
466 193 : if (isParamValid("sampling_blocks"))
467 : {
468 : // Remove all elements not in the blocks
469 39 : const auto & blocks_to_keep_names = getParam<std::vector<SubdomainName>>("sampling_blocks");
470 39 : const auto & blocks_to_keep = _sampling_mesh_ptr->getSubdomainIDs(blocks_to_keep_names);
471 1239 : for (const auto & elem_ptr : _sampling_mesh_ptr->getMesh().element_ptr_range())
472 600 : if (std::find(blocks_to_keep.begin(), blocks_to_keep.end(), elem_ptr->subdomain_id()) ==
473 1200 : blocks_to_keep.end())
474 300 : _sampling_mesh_ptr->getMesh().delete_elem(elem_ptr);
475 :
476 : // Deleting elements and isolated nodes would cause renumbering. Not renumbering might help
477 : // user examining the sampling mesh and the regular mesh. Also if we end up partitioning the
478 : // elements, the node partitioning is unlikely to match if the element numbering is different.
479 : // Still not enough of a guarantee, because of deleted elements the node partitioning could be
480 : // different. We will rely on ghosting to make it work
481 39 : _sampling_mesh_ptr->getMesh().allow_renumbering(false);
482 39 : }
483 :
484 : // Set a partitioner
485 193 : if (!_serialize)
486 : {
487 13 : _sampling_mesh_ptr->setIsCustomPartitionerRequested(true);
488 13 : InputParameters partition_params = _app.getFactory().getValidParams("CopyMeshPartitioner");
489 13 : partition_params.set<MooseMesh *>("mesh") = _sampling_mesh_ptr.get();
490 13 : partition_params.set<MooseMesh *>("source_mesh") = _mesh_ptr;
491 13 : std::shared_ptr<MoosePartitioner> mp = _factory.create<MoosePartitioner>(
492 13 : "CopyMeshPartitioner", "sampled_output_part", partition_params);
493 13 : _sampling_mesh_ptr->setCustomPartitioner(mp.get());
494 :
495 13 : _sampling_mesh_ptr->getMesh().prepare_for_use();
496 : // this should be called by prepare_for_use, but is not.
497 : // it also requires a prior call to prepare_for_use()
498 13 : mp->partition(_sampling_mesh_ptr->getMesh(), comm().size());
499 13 : }
500 :
501 : // Prepare mesh, needed for the mesh functions
502 193 : if (_using_external_sampling_file)
503 11 : _sampling_mesh_ptr->prepare(/*mesh to clone*/ nullptr);
504 182 : else if (_serialize && isParamValid("sampling_blocks"))
505 : // TODO: constraints have not been initialized?
506 26 : _sampling_mesh_ptr->getMesh().prepare_for_use();
507 :
508 193 : if (_serialize)
509 : // we want to avoid re-partitioning, as we will serialize anyway
510 180 : _sampling_mesh_ptr->getMesh().skip_partitioning(true);
511 :
512 : // Make sure that the mesh pointer points to the newly cloned mesh
513 193 : _mesh_ptr = _sampling_mesh_ptr.get();
514 :
515 : // Check the source and target mesh in case their subdomains match
516 193 : const std::vector<SubdomainID> mesh_subdomain_ids_vec(_mesh_ptr->meshSubdomains().begin(),
517 386 : _mesh_ptr->meshSubdomains().end());
518 : const std::vector<SubdomainID> initial_mesh_subdomain_ids_vec(
519 193 : _problem_ptr->mesh().meshSubdomains().begin(), _problem_ptr->mesh().meshSubdomains().end());
520 193 : _mesh_subdomains_match =
521 204 : (_mesh_ptr->meshSubdomains() == _problem_ptr->mesh().meshSubdomains() &&
522 204 : _mesh_ptr->getSubdomainNames(mesh_subdomain_ids_vec) ==
523 204 : _problem_ptr->mesh().getSubdomainNames(initial_mesh_subdomain_ids_vec));
524 193 : if (_using_external_sampling_file && !_mesh_subdomains_match)
525 0 : mooseInfoRepeated("Variable block restriction disabled in sampled output due to non-matching "
526 : "subdomain names and ids");
527 193 : }
528 :
529 : void
530 28429 : SampledOutput::setFileBaseInternal(const std::string & file_base)
531 : {
532 28429 : AdvancedOutput::setFileBaseInternal(file_base);
533 : // ** DEPRECATED SUPPORT **
534 28429 : if (getParam<bool>("append_oversample"))
535 0 : _file_base += "_oversample";
536 28429 : }
537 :
538 : bool
539 737 : SampledOutput::isSampledAtNodes(const FEType & fe_type) const
540 : {
541 : // This is the same criterion as in MooseVariableData
542 737 : const auto continuity = FEInterface::get_continuity(fe_type);
543 737 : return (continuity == C_ZERO || continuity == C_ONE);
544 : }
|
