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 : #include "Adaptivity.h"
11 :
12 : #include "AuxiliarySystem.h"
13 : #include "DisplacedProblem.h"
14 : #include "FEProblem.h"
15 : #include "FlagElementsThread.h"
16 : #include "MooseMesh.h"
17 : #include "NonlinearSystemBase.h"
18 : #include "UpdateErrorVectorsThread.h"
19 :
20 : // libMesh
21 : #include "libmesh/equation_systems.h"
22 : #include "libmesh/kelly_error_estimator.h"
23 : #include "libmesh/patch_recovery_error_estimator.h"
24 : #include "libmesh/fourth_error_estimators.h"
25 : #include "libmesh/parallel.h"
26 : #include "libmesh/error_vector.h"
27 : #include "libmesh/distributed_mesh.h"
28 :
29 : using namespace libMesh;
30 :
31 : #ifdef LIBMESH_ENABLE_AMR
32 :
33 61643 : Adaptivity::Adaptivity(FEProblemBase & fe_problem)
34 : : ConsoleStreamInterface(fe_problem.getMooseApp()),
35 61643 : PerfGraphInterface(fe_problem.getMooseApp().perfGraph(), "Adaptivity"),
36 61643 : ParallelObject(fe_problem.getMooseApp()),
37 61643 : _fe_problem(fe_problem),
38 123286 : _mesh(_fe_problem.mesh()),
39 61643 : _mesh_refinement_on(false),
40 61643 : _initialized(false),
41 61643 : _initial_steps(0),
42 61643 : _steps(0),
43 61643 : _print_mesh_changed(false),
44 61643 : _t(_fe_problem.time()),
45 61643 : _step(_fe_problem.timeStep()),
46 61643 : _interval(1),
47 61643 : _start_time(-std::numeric_limits<Real>::max()),
48 61643 : _stop_time(std::numeric_limits<Real>::max()),
49 61643 : _cycles_per_step(1),
50 61643 : _use_new_system(false),
51 61643 : _adaptivity_type(AdaptivityType::H),
52 61643 : _max_h_level(0),
53 308215 : _recompute_markers_during_cycles(false)
54 : {
55 61643 : }
56 :
57 58563 : Adaptivity::~Adaptivity() {}
58 :
59 : void
60 2247 : Adaptivity::init(const unsigned int steps,
61 : const unsigned int initial_steps,
62 : const AdaptivityType adaptivity_type)
63 : {
64 : // Get the pointer to the DisplacedProblem, this cannot be done at construction because
65 : // DisplacedProblem
66 : // does not exist at that point.
67 2247 : _displaced_problem = _fe_problem.getDisplacedProblem();
68 :
69 2247 : _mesh_refinement = std::make_unique<MeshRefinement>(_mesh);
70 2247 : _error = std::make_unique<ErrorVector>();
71 :
72 2247 : EquationSystems & es = _fe_problem.es();
73 4494 : es.parameters.set<bool>("adaptivity") = true;
74 :
75 2247 : _initial_steps = initial_steps;
76 2247 : _steps = steps;
77 :
78 2247 : _adaptivity_type = adaptivity_type;
79 :
80 2247 : _mesh_refinement_on = true;
81 :
82 2247 : if (_adaptivity_type == AdaptivityType::P || _adaptivity_type == AdaptivityType::HP)
83 228 : _mesh.doingPRefinement(true);
84 :
85 2247 : if (_adaptivity_type == AdaptivityType::HP)
86 : {
87 24 : if (_fe_problem.numSolverSystems() > 1)
88 0 : mooseError("HP adaptivity doesn't support multiple solver systems because currently only the "
89 : "zeroth solver system solution is analyzed for determining whether to toggle "
90 : "h-refinement flags to p-refinement flags");
91 :
92 24 : _sibling_coupling = std::make_unique<SiblingCoupling>();
93 48 : _fe_problem.getSolverSystem(0).system().get_dof_map().add_algebraic_ghosting_functor(
94 24 : *_sibling_coupling);
95 : }
96 :
97 4494 : _mesh_refinement->set_periodic_boundaries_ptr(
98 2247 : _fe_problem.getSolverSystem(/*nl_sys=*/0).dofMap().get_periodic_boundaries());
99 :
100 2247 : if (_displaced_problem)
101 : {
102 124 : EquationSystems & displaced_es = _displaced_problem->es();
103 248 : displaced_es.parameters.set<bool>("adaptivity") = true;
104 :
105 124 : if (!_displaced_mesh_refinement)
106 124 : _displaced_mesh_refinement = std::make_unique<MeshRefinement>(_displaced_problem->mesh());
107 :
108 : // The periodic boundaries pointer allows the MeshRefinement
109 : // object to determine elements which are "topological" neighbors,
110 : // i.e. neighbors across periodic boundaries, for the purposes of
111 : // refinement.
112 248 : _displaced_mesh_refinement->set_periodic_boundaries_ptr(
113 124 : _fe_problem.getSolverSystem(/*nl_sys=*/0).dofMap().get_periodic_boundaries());
114 :
115 : // TODO: This is currently an empty function on the DisplacedProblem... could it be removed?
116 124 : _displaced_problem->initAdaptivity();
117 : }
118 :
119 : // indicate the Adaptivity system has been initialized
120 2247 : _initialized = true;
121 2247 : }
122 :
123 : void
124 422 : Adaptivity::setErrorEstimator(const MooseEnum & error_estimator_name)
125 : {
126 422 : if (error_estimator_name == "KellyErrorEstimator")
127 416 : _error_estimator = std::make_unique<KellyErrorEstimator>();
128 6 : else if (error_estimator_name == "LaplacianErrorEstimator")
129 0 : _error_estimator = std::make_unique<LaplacianErrorEstimator>();
130 6 : else if (error_estimator_name == "PatchRecoveryErrorEstimator")
131 6 : _error_estimator = std::make_unique<PatchRecoveryErrorEstimator>();
132 : else
133 0 : mooseError(std::string("Unknown error_estimator selection: ") +
134 0 : std::string(error_estimator_name));
135 422 : }
136 :
137 : void
138 12 : Adaptivity::setErrorNorm(SystemNorm & sys_norm)
139 : {
140 : mooseAssert(_error_estimator, "error_estimator not initialized. Did you call init_adaptivity()?");
141 12 : _error_estimator->error_norm = sys_norm;
142 12 : }
143 :
144 : bool
145 5762 : Adaptivity::adaptMesh(std::string marker_name /*=std::string()*/)
146 : {
147 28810 : TIME_SECTION("adaptMesh", 3, "Adapting Mesh");
148 :
149 : // If the marker name is supplied, use it. Otherwise, use the one in _marker_variable_name
150 5762 : if (marker_name.empty())
151 4882 : marker_name = _marker_variable_name;
152 :
153 5762 : bool mesh_changed = false;
154 :
155 : // If mesh adaptivity is carried out in a distributed (scalable) way
156 5762 : bool distributed_adaptivity = false;
157 :
158 5762 : if (_use_new_system)
159 : {
160 4890 : if (!marker_name.empty()) // Only flag if a marker variable name has been set
161 : {
162 4857 : _mesh_refinement->clean_refinement_flags();
163 :
164 4857 : std::vector<Number> serialized_solution;
165 :
166 4857 : auto distributed_mesh = dynamic_cast<DistributedMesh *>(&_fe_problem.mesh().getMesh());
167 :
168 : // Element range
169 4857 : std::unique_ptr<ConstElemRange> all_elems;
170 : // If the mesh is distributed and we do not do "gather to zero" or "allgather".
171 : // Then it is safe to not serialize solution.
172 : // Some output idiom (Exodus) will do "gather to zero". That being said,
173 : // if you have exodus output on, mesh adaptivty is not scalable.
174 4857 : if (distributed_mesh && !distributed_mesh->is_serial_on_zero())
175 : {
176 : // We update here to make sure local solution is up-to-date
177 904 : _fe_problem.getAuxiliarySystem().update();
178 904 : distributed_adaptivity = true;
179 :
180 : // We can not assume that geometric and algebraic ghosting functors cover
181 : // the same set of elements/nodes. That being said, in general,
182 : // we would expect G(e) > A(e). Here G(e) is the set of elements reserved
183 : // by the geometric ghosting functors, and A(e) corresponds to
184 : // the one covered by the algebraic ghosting functors.
185 : // Therefore, we have to work only on local elements instead of
186 : // ghosted + local elements. The ghosted solution might not be enough
187 : // for ghosted+local elements. But it is always sufficient for local elements.
188 : // After we set markers for all local elements, we will do a global
189 : // communication to sync markers for ghosted elements from their owners.
190 1808 : all_elems = std::make_unique<ConstElemRange>(
191 1808 : _fe_problem.mesh().getMesh().active_local_elements_begin(),
192 2712 : _fe_problem.mesh().getMesh().active_local_elements_end());
193 : }
194 : else // This is not scalable but it might be useful for small-size problems
195 : {
196 3953 : _fe_problem.getAuxiliarySystem().solution().close();
197 3953 : _fe_problem.getAuxiliarySystem().solution().localize(serialized_solution);
198 3953 : distributed_adaptivity = false;
199 :
200 : // For a replicated mesh or a serialized distributed mesh, the solution
201 : // is serialized to everyone. Then we update markers for all active elements.
202 : // In this case, we can avoid a global communication to update mesh.
203 : // I do not know if it is a good idea, but it the old code behavior.
204 : // We might not care about much since a replicated mesh
205 : // or a serialized distributed mesh is not scalable anyway.
206 : all_elems =
207 7906 : std::make_unique<ConstElemRange>(_fe_problem.mesh().getMesh().active_elements_begin(),
208 11859 : _fe_problem.mesh().getMesh().active_elements_end());
209 : }
210 :
211 : FlagElementsThread fet(
212 4857 : _fe_problem, serialized_solution, _max_h_level, marker_name, !distributed_adaptivity);
213 4857 : Threads::parallel_reduce(*all_elems, fet);
214 4857 : _fe_problem.getAuxiliarySystem().solution().close();
215 4857 : }
216 : }
217 : else
218 : {
219 872 : if (_fe_problem.numSolverSystems() > 1)
220 0 : mooseError("The old adaptivity system based on libMesh error estimators does not currently "
221 : "support multiple solver systems because we currently only use the zeroth solver "
222 : "system solution for estimating errors");
223 :
224 : // Compute the error for each active element
225 872 : _error_estimator->estimate_error(_fe_problem.getSolverSystem(/*nl_sys=*/0).system(), *_error);
226 :
227 : // Flag elements to be refined and coarsened
228 872 : _mesh_refinement->flag_elements_by_error_fraction(*_error);
229 : }
230 :
231 : // Moving some of h flagged elements to p flagged based on the
232 : // local smoothness and prior h & p error estimates
233 5762 : if (_adaptivity_type == AdaptivityType::HP)
234 : {
235 66 : _hp_coarsen_test = std::make_unique<HPCoarsenTest>();
236 66 : _hp_coarsen_test->select_refinement(_fe_problem.getSolverSystem(/*nl_sys=*/0).system());
237 : }
238 :
239 : // If the DisplacedProblem is active, undisplace the DisplacedMesh
240 : // in preparation for refinement. We can't safely refine the
241 : // DisplacedMesh directly, since the Hilbert keys computed on the
242 : // inconsistenly-displaced Mesh are different on different
243 : // processors, leading to inconsistent Hilbert keys. We must do
244 : // this before the undisplaced Mesh is refined, so that the
245 : // element and node numbering is still consistent.
246 5762 : if (_displaced_problem)
247 510 : _displaced_problem->undisplaceMesh();
248 :
249 : // If markers are added to only local elements,
250 : // we sync them here.
251 5762 : if (distributed_adaptivity)
252 904 : _mesh_refinement->make_flags_parallel_consistent();
253 :
254 : // Sync flags from the reference mesh
255 5762 : if (_displaced_problem)
256 510 : for (auto * const displaced_elem :
257 337800 : _displaced_problem->mesh().getMesh().active_element_ptr_range())
258 : {
259 336780 : const auto * const reference_elem = _fe_problem.mesh().elemPtr(displaced_elem->id());
260 336780 : displaced_elem->set_refinement_flag(reference_elem->refinement_flag());
261 336780 : displaced_elem->set_p_refinement_flag(reference_elem->p_refinement_flag());
262 510 : }
263 :
264 5762 : if (_adaptivity_type == AdaptivityType::P)
265 253 : _mesh_refinement->switch_h_to_p_refinement();
266 :
267 : // Perform refinement and coarsening
268 5762 : mesh_changed = _mesh_refinement->refine_and_coarsen_elements();
269 :
270 5762 : if (_displaced_problem && mesh_changed)
271 : {
272 387 : if (_adaptivity_type == AdaptivityType::P)
273 0 : _displaced_mesh_refinement->switch_h_to_p_refinement();
274 :
275 : #ifndef NDEBUG
276 : bool displaced_mesh_changed =
277 : #endif
278 387 : _displaced_mesh_refinement->refine_and_coarsen_elements();
279 :
280 : // Since the undisplaced mesh changed, the displaced mesh better have changed!
281 : mooseAssert(displaced_mesh_changed, "Undisplaced mesh changed, but displaced mesh did not!");
282 : }
283 :
284 5762 : if (mesh_changed && _print_mesh_changed)
285 : {
286 0 : _console << "\nMesh Changed:\n";
287 0 : _mesh.printInfo();
288 0 : _console << std::flush;
289 : }
290 :
291 5762 : return mesh_changed;
292 5762 : }
293 :
294 : bool
295 1055 : Adaptivity::initialAdaptMesh()
296 : {
297 1055 : return adaptMesh(_initial_marker_variable_name);
298 : }
299 :
300 : void
301 3698 : Adaptivity::uniformRefine(MooseMesh * mesh, unsigned int level /*=libMesh::invalid_uint*/)
302 : {
303 : mooseAssert(mesh, "Mesh pointer must not be NULL");
304 :
305 : // NOTE: we are using a separate object here, since adaptivity may not be on, but we need to be
306 : // able to do refinements
307 3698 : MeshRefinement mesh_refinement(*mesh);
308 3698 : if (level == libMesh::invalid_uint)
309 3614 : level = mesh->uniformRefineLevel();
310 :
311 : // Skip deletion and repartition will make uniform refinements run more
312 : // efficiently, but at the same time, there might be extra ghosting elements.
313 : // The number of layers of additional ghosting elements depends on the number
314 : // of uniform refinement levels. This should happen only when you have a "fine enough"
315 : // coarse mesh and want to refine the mesh by a few levels. Otherwise, it might
316 : // introduce an unbalanced workload and too large ghosting domain.
317 3698 : if (mesh->skipDeletionRepartitionAfterRefine())
318 : {
319 24 : mesh->getMesh().skip_partitioning(true);
320 24 : mesh->getMesh().allow_remote_element_removal(false);
321 24 : mesh->needsRemoteElemDeletion(false);
322 : }
323 :
324 3698 : mesh_refinement.uniformly_refine(level);
325 3698 : }
326 :
327 : void
328 10 : Adaptivity::uniformRefineWithProjection()
329 : {
330 50 : TIME_SECTION("uniformRefineWithProjection", 2, "Uniformly Refining and Reprojecting");
331 :
332 : // NOTE: we are using a separate object here, since adaptivity may not be on, but we need to be
333 : // able to do refinements
334 10 : MeshRefinement mesh_refinement(_mesh);
335 10 : unsigned int level = _mesh.uniformRefineLevel();
336 10 : MeshRefinement displaced_mesh_refinement(_displaced_problem ? _displaced_problem->mesh() : _mesh);
337 :
338 : // we have to go step by step so EquationSystems::reinit() won't freak out
339 20 : for (unsigned int i = 0; i < level; i++)
340 : {
341 : // See comment above about why refining the displaced mesh is potentially unsafe.
342 10 : if (_displaced_problem)
343 0 : _displaced_problem->undisplaceMesh();
344 :
345 10 : mesh_refinement.uniformly_refine(1);
346 :
347 10 : if (_displaced_problem)
348 0 : displaced_mesh_refinement.uniformly_refine(1);
349 10 : _fe_problem.meshChanged(
350 : /*intermediate_step=*/false, /*contract_mesh=*/true, /*clean_refinement_flags=*/true);
351 : }
352 10 : }
353 :
354 : void
355 206 : Adaptivity::setAdaptivityOn(bool state)
356 : {
357 : // check if Adaptivity has been initialized before turning on
358 206 : if (state == true && !_initialized)
359 0 : mooseError("Mesh adaptivity system not available");
360 :
361 206 : _mesh_refinement_on = state;
362 206 : }
363 :
364 : void
365 2247 : Adaptivity::setTimeActive(Real start_time, Real stop_time)
366 : {
367 2247 : _start_time = start_time;
368 2247 : _stop_time = stop_time;
369 2247 : }
370 :
371 : void
372 1825 : Adaptivity::setUseNewSystem()
373 : {
374 1825 : _use_new_system = true;
375 1825 : }
376 :
377 : void
378 1142 : Adaptivity::setMarkerVariableName(std::string marker_field)
379 : {
380 1142 : _marker_variable_name = marker_field;
381 1142 : }
382 :
383 : void
384 670 : Adaptivity::setInitialMarkerVariableName(std::string marker_field)
385 : {
386 670 : _initial_marker_variable_name = marker_field;
387 670 : }
388 :
389 : ErrorVector &
390 530 : Adaptivity::getErrorVector(const std::string & indicator_field)
391 : {
392 : // Insert or retrieve error vector
393 530 : auto insert_pair = moose_try_emplace(
394 530 : _indicator_field_to_error_vector, indicator_field, std::make_unique<ErrorVector>());
395 1060 : return *insert_pair.first->second;
396 : }
397 :
398 : void
399 6477 : Adaptivity::updateErrorVectors()
400 : {
401 32385 : TIME_SECTION("updateErrorVectors", 5, "Updating Error Vectors");
402 :
403 : // Resize all of the ErrorVectors in case the mesh has changed
404 8731 : for (const auto & it : _indicator_field_to_error_vector)
405 : {
406 2254 : ErrorVector & vec = *(it.second);
407 2254 : vec.assign(_mesh.getMesh().max_elem_id(), 0);
408 : }
409 :
410 : // Fill the vectors with the local contributions
411 6477 : UpdateErrorVectorsThread uevt(_fe_problem, _indicator_field_to_error_vector);
412 6477 : Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), uevt);
413 :
414 : // Now sum across all processors
415 8731 : for (const auto & it : _indicator_field_to_error_vector)
416 2254 : _fe_problem.comm().sum((std::vector<float> &)*(it.second));
417 6477 : }
418 :
419 : bool
420 164935 : Adaptivity::isAdaptivityDue()
421 : {
422 164935 : return _mesh_refinement_on && (_start_time <= _t && _t < _stop_time) && _step % _interval == 0;
423 : }
424 :
425 : #endif // LIBMESH_ENABLE_AMR
|
