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 57974 : Adaptivity::Adaptivity(FEProblemBase & fe_problem)
34 : : ConsoleStreamInterface(fe_problem.getMooseApp()),
35 57974 : PerfGraphInterface(fe_problem.getMooseApp().perfGraph(), "Adaptivity"),
36 57974 : ParallelObject(fe_problem.getMooseApp()),
37 57974 : _fe_problem(fe_problem),
38 115948 : _mesh(_fe_problem.mesh()),
39 57974 : _mesh_refinement_on(false),
40 57974 : _initialized(false),
41 57974 : _initial_steps(0),
42 57974 : _steps(0),
43 57974 : _print_mesh_changed(false),
44 57974 : _t(_fe_problem.time()),
45 57974 : _step(_fe_problem.timeStep()),
46 57974 : _interval(1),
47 57974 : _start_time(-std::numeric_limits<Real>::max()),
48 57974 : _stop_time(std::numeric_limits<Real>::max()),
49 57974 : _cycles_per_step(1),
50 57974 : _use_new_system(false),
51 57974 : _max_h_level(0),
52 173922 : _recompute_markers_during_cycles(false)
53 : {
54 57974 : }
55 :
56 53701 : Adaptivity::~Adaptivity() {}
57 :
58 : void
59 2202 : 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 2202 : _displaced_problem = _fe_problem.getDisplacedProblem();
67 :
68 2202 : _mesh_refinement = std::make_unique<MeshRefinement>(_mesh);
69 2202 : _error = std::make_unique<ErrorVector>();
70 :
71 2202 : EquationSystems & es = _fe_problem.es();
72 2202 : es.parameters.set<bool>("adaptivity") = true;
73 :
74 2202 : _initial_steps = initial_steps;
75 2202 : _steps = steps;
76 2202 : _p_refinement_flag = p_refinement;
77 2202 : _mesh_refinement_on = true;
78 :
79 2202 : if (_p_refinement_flag)
80 204 : _mesh.doingPRefinement(true);
81 :
82 4404 : _mesh_refinement->set_periodic_boundaries_ptr(
83 2202 : _fe_problem.getNonlinearSystemBase(/*nl_sys=*/0).dofMap().get_periodic_boundaries());
84 :
85 : // displaced problem
86 2202 : if (_displaced_problem != nullptr)
87 : {
88 130 : EquationSystems & displaced_es = _displaced_problem->es();
89 130 : displaced_es.parameters.set<bool>("adaptivity") = true;
90 :
91 130 : if (!_displaced_mesh_refinement)
92 130 : _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 260 : _displaced_mesh_refinement->set_periodic_boundaries_ptr(
99 130 : _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 130 : _displaced_problem->initAdaptivity();
103 : }
104 :
105 : // indicate the Adaptivity system has been initialized
106 2202 : _initialized = true;
107 2202 : }
108 :
109 : void
110 426 : Adaptivity::setErrorEstimator(const MooseEnum & error_estimator_name)
111 : {
112 426 : if (error_estimator_name == "KellyErrorEstimator")
113 417 : _error_estimator = std::make_unique<KellyErrorEstimator>();
114 9 : else if (error_estimator_name == "LaplacianErrorEstimator")
115 0 : _error_estimator = std::make_unique<LaplacianErrorEstimator>();
116 9 : else if (error_estimator_name == "PatchRecoveryErrorEstimator")
117 9 : _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 426 : }
122 :
123 : void
124 12 : Adaptivity::setErrorNorm(SystemNorm & sys_norm)
125 : {
126 : mooseAssert(_error_estimator, "error_estimator not initialized. Did you call init_adaptivity()?");
127 12 : _error_estimator->error_norm = sys_norm;
128 12 : }
129 :
130 : bool
131 5725 : Adaptivity::adaptMesh(std::string marker_name /*=std::string()*/)
132 : {
133 5725 : 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 5725 : if (marker_name.empty())
137 4870 : marker_name = _marker_variable_name;
138 :
139 5725 : bool mesh_changed = false;
140 :
141 : // If mesh adaptivity is carried out in a distributed (scalable) way
142 5725 : bool distributed_adaptivity = false;
143 :
144 5725 : if (_use_new_system)
145 : {
146 4866 : if (!marker_name.empty()) // Only flag if a marker variable name has been set
147 : {
148 4833 : _mesh_refinement->clean_refinement_flags();
149 :
150 4833 : std::vector<Number> serialized_solution;
151 :
152 4833 : auto distributed_mesh = dynamic_cast<DistributedMesh *>(&_fe_problem.mesh().getMesh());
153 :
154 : // Element range
155 4833 : 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 4833 : if (distributed_mesh && !distributed_mesh->is_serial_on_zero())
161 : {
162 : // We update here to make sure local solution is up-to-date
163 894 : _fe_problem.getAuxiliarySystem().update();
164 894 : 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 1788 : all_elems = std::make_unique<ConstElemRange>(
177 1788 : _fe_problem.mesh().getMesh().active_local_elements_begin(),
178 2682 : _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 3939 : _fe_problem.getAuxiliarySystem().solution().close();
183 3939 : _fe_problem.getAuxiliarySystem().solution().localize(serialized_solution);
184 3939 : 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 7878 : std::make_unique<ConstElemRange>(_fe_problem.mesh().getMesh().active_elements_begin(),
194 11817 : _fe_problem.mesh().getMesh().active_elements_end());
195 : }
196 :
197 : FlagElementsThread fet(
198 4833 : _fe_problem, serialized_solution, _max_h_level, marker_name, !distributed_adaptivity);
199 4833 : Threads::parallel_reduce(*all_elems, fet);
200 4833 : _fe_problem.getAuxiliarySystem().solution().close();
201 4833 : }
202 : }
203 : else
204 : {
205 : // Compute the error for each active element
206 1718 : _error_estimator->estimate_error(_fe_problem.getNonlinearSystemBase(/*nl_sys=*/0).system(),
207 859 : *_error);
208 :
209 : // Flag elements to be refined and coarsened
210 859 : _mesh_refinement->flag_elements_by_error_fraction(*_error);
211 :
212 859 : if (_displaced_problem)
213 : // Reuse the error vector and refine the displaced mesh
214 31 : _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 5725 : if (_displaced_problem)
225 524 : _displaced_problem->undisplaceMesh();
226 :
227 : // If markers are added to only local elements,
228 : // we sync them here.
229 5725 : if (distributed_adaptivity)
230 894 : _mesh_refinement->make_flags_parallel_consistent();
231 :
232 5725 : if (_p_refinement_flag)
233 253 : _mesh_refinement->switch_h_to_p_refinement();
234 :
235 : // Perform refinement and coarsening
236 5725 : mesh_changed = _mesh_refinement->refine_and_coarsen_elements();
237 :
238 5725 : if (_displaced_problem && mesh_changed)
239 : {
240 : // If markers are added to only local elements,
241 : // we sync them here.
242 399 : if (distributed_adaptivity)
243 129 : _displaced_mesh_refinement->make_flags_parallel_consistent();
244 :
245 399 : 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 399 : _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 5725 : 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 5725 : return mesh_changed;
265 5725 : }
266 :
267 : bool
268 1036 : Adaptivity::initialAdaptMesh()
269 : {
270 1036 : return adaptMesh(_initial_marker_variable_name);
271 : }
272 :
273 : void
274 4803 : 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 4803 : MeshRefinement mesh_refinement(*mesh);
281 4803 : if (level == libMesh::invalid_uint)
282 4719 : 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 4803 : if (mesh->skipDeletionRepartitionAfterRefine())
291 : {
292 24 : mesh->getMesh().skip_partitioning(true);
293 24 : mesh->getMesh().allow_remote_element_removal(false);
294 24 : mesh->needsRemoteElemDeletion(false);
295 : }
296 :
297 4803 : mesh_refinement.uniformly_refine(level);
298 4803 : }
299 :
300 : void
301 10 : Adaptivity::uniformRefineWithProjection()
302 : {
303 10 : 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 10 : MeshRefinement mesh_refinement(_mesh);
308 10 : unsigned int level = _mesh.uniformRefineLevel();
309 10 : 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 20 : for (unsigned int i = 0; i < level; i++)
313 : {
314 : // See comment above about why refining the displaced mesh is potentially unsafe.
315 10 : if (_displaced_problem)
316 0 : _displaced_problem->undisplaceMesh();
317 :
318 10 : mesh_refinement.uniformly_refine(1);
319 :
320 10 : if (_displaced_problem)
321 0 : displaced_mesh_refinement.uniformly_refine(1);
322 10 : _fe_problem.meshChanged(
323 : /*intermediate_step=*/false, /*contract_mesh=*/true, /*clean_refinement_flags=*/true);
324 : }
325 10 : }
326 :
327 : void
328 203 : Adaptivity::setAdaptivityOn(bool state)
329 : {
330 : // check if Adaptivity has been initialized before turning on
331 203 : if (state == true && !_initialized)
332 0 : mooseError("Mesh adaptivity system not available");
333 :
334 203 : _mesh_refinement_on = state;
335 203 : }
336 :
337 : void
338 2202 : Adaptivity::setTimeActive(Real start_time, Real stop_time)
339 : {
340 2202 : _start_time = start_time;
341 2202 : _stop_time = stop_time;
342 2202 : }
343 :
344 : void
345 1776 : Adaptivity::setUseNewSystem()
346 : {
347 1776 : _use_new_system = true;
348 1776 : }
349 :
350 : void
351 1140 : Adaptivity::setMarkerVariableName(std::string marker_field)
352 : {
353 1140 : _marker_variable_name = marker_field;
354 1140 : }
355 :
356 : void
357 654 : Adaptivity::setInitialMarkerVariableName(std::string marker_field)
358 : {
359 654 : _initial_marker_variable_name = marker_field;
360 654 : }
361 :
362 : ErrorVector &
363 503 : Adaptivity::getErrorVector(const std::string & indicator_field)
364 : {
365 : // Insert or retrieve error vector
366 503 : auto insert_pair = moose_try_emplace(
367 503 : _indicator_field_to_error_vector, indicator_field, std::make_unique<ErrorVector>());
368 1006 : return *insert_pair.first->second;
369 : }
370 :
371 : void
372 6428 : Adaptivity::updateErrorVectors()
373 : {
374 6428 : TIME_SECTION("updateErrorVectors", 5, "Updating Error Vectors");
375 :
376 : // Resize all of the ErrorVectors in case the mesh has changed
377 8624 : for (const auto & it : _indicator_field_to_error_vector)
378 : {
379 2196 : ErrorVector & vec = *(it.second);
380 2196 : vec.assign(_mesh.getMesh().max_elem_id(), 0);
381 : }
382 :
383 : // Fill the vectors with the local contributions
384 6428 : UpdateErrorVectorsThread uevt(_fe_problem, _indicator_field_to_error_vector);
385 6428 : Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), uevt);
386 :
387 : // Now sum across all processors
388 8624 : for (const auto & it : _indicator_field_to_error_vector)
389 2196 : _fe_problem.comm().sum((std::vector<float> &)*(it.second));
390 6428 : }
391 :
392 : bool
393 163150 : Adaptivity::isAdaptivityDue()
394 : {
395 163150 : return _mesh_refinement_on && (_start_time <= _t && _t < _stop_time) && _step % _interval == 0;
396 : }
397 :
398 : #endif // LIBMESH_ENABLE_AMR
|
