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