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