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