www.mooseframework.org
Adaptivity.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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"
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 #ifdef LIBMESH_ENABLE_AMR
30 
32  : ConsoleStreamInterface(fe_problem.getMooseApp()),
33  PerfGraphInterface(fe_problem.getMooseApp().perfGraph(), "Adaptivity"),
34  ParallelObject(fe_problem.getMooseApp()),
35  _fe_problem(fe_problem),
36  _mesh(_fe_problem.mesh()),
37  _mesh_refinement_on(false),
38  _initialized(false),
39  _initial_steps(0),
40  _steps(0),
41  _print_mesh_changed(false),
42  _t(_fe_problem.time()),
43  _step(_fe_problem.timeStep()),
44  _interval(1),
45  _start_time(-std::numeric_limits<Real>::max()),
46  _stop_time(std::numeric_limits<Real>::max()),
47  _cycles_per_step(1),
48  _use_new_system(false),
49  _max_h_level(0),
50  _recompute_markers_during_cycles(false)
51 {
52 }
53 
55 
56 void
57 Adaptivity::init(unsigned int steps, unsigned int initial_steps)
58 {
59  // Get the pointer to the DisplacedProblem, this cannot be done at construction because
60  // DisplacedProblem
61  // does not exist at that point.
63 
64  _mesh_refinement = std::make_unique<MeshRefinement>(_mesh);
65  _error = std::make_unique<ErrorVector>();
66 
67  EquationSystems & es = _fe_problem.es();
68  es.parameters.set<bool>("adaptivity") = true;
69 
70  _initial_steps = initial_steps;
71  _steps = steps;
72  _mesh_refinement_on = true;
73 
74  _mesh_refinement->set_periodic_boundaries_ptr(
75  _fe_problem.getNonlinearSystemBase(/*nl_sys=*/0).dofMap().get_periodic_boundaries());
76 
77  // displaced problem
78  if (_displaced_problem != nullptr)
79  {
80  EquationSystems & displaced_es = _displaced_problem->es();
81  displaced_es.parameters.set<bool>("adaptivity") = true;
82 
84  _displaced_mesh_refinement = std::make_unique<MeshRefinement>(_displaced_problem->mesh());
85 
86  // The periodic boundaries pointer allows the MeshRefinement
87  // object to determine elements which are "topological" neighbors,
88  // i.e. neighbors across periodic boundaries, for the purposes of
89  // refinement.
90  _displaced_mesh_refinement->set_periodic_boundaries_ptr(
91  _fe_problem.getNonlinearSystemBase(/*nl_sys=*/0).dofMap().get_periodic_boundaries());
92 
93  // TODO: This is currently an empty function on the DisplacedProblem... could it be removed?
94  _displaced_problem->initAdaptivity();
95  }
96 
97  // indicate the Adaptivity system has been initialized
98  _initialized = true;
99 }
100 
101 void
102 Adaptivity::setErrorEstimator(const MooseEnum & error_estimator_name)
103 {
104  if (error_estimator_name == "KellyErrorEstimator")
105  _error_estimator = std::make_unique<KellyErrorEstimator>();
106  else if (error_estimator_name == "LaplacianErrorEstimator")
107  _error_estimator = std::make_unique<LaplacianErrorEstimator>();
108  else if (error_estimator_name == "PatchRecoveryErrorEstimator")
109  _error_estimator = std::make_unique<PatchRecoveryErrorEstimator>();
110  else
111  mooseError(std::string("Unknown error_estimator selection: ") +
112  std::string(error_estimator_name));
113 }
114 
115 void
116 Adaptivity::setErrorNorm(SystemNorm & sys_norm)
117 {
118  mooseAssert(_error_estimator, "error_estimator not initialized. Did you call init_adaptivity()?");
119  _error_estimator->error_norm = sys_norm;
120 }
121 
122 bool
123 Adaptivity::adaptMesh(std::string marker_name /*=std::string()*/)
124 {
125  TIME_SECTION("adaptMesh", 3, "Adapting Mesh");
126 
127  // If the marker name is supplied, use it. Otherwise, use the one in _marker_variable_name
128  if (marker_name.empty())
129  marker_name = _marker_variable_name;
130 
131  bool mesh_changed = false;
132 
133  // If mesh adaptivity is carried out in a distributed (scalable) way
134  bool distributed_adaptivity = false;
135 
136  if (_use_new_system)
137  {
138  if (!marker_name.empty()) // Only flag if a marker variable name has been set
139  {
140  _mesh_refinement->clean_refinement_flags();
141 
142  std::vector<Number> serialized_solution;
143 
144  auto distributed_mesh = dynamic_cast<DistributedMesh *>(&_fe_problem.mesh().getMesh());
145 
146  // Element range
147  std::unique_ptr<ConstElemRange> all_elems;
148  // If the mesh is distributed and we do not do "gather to zero" or "allgather".
149  // Then it is safe to not serialize solution.
150  // Some output idiom (Exodus) will do "gather to zero". That being said,
151  // if you have exodus output on, mesh adaptivty is not scalable.
152  if (distributed_mesh && !distributed_mesh->is_serial_on_zero())
153  {
154  // We update here to make sure local solution is up-to-date
156  distributed_adaptivity = true;
157 
158  // We can not assume that geometric and algebraic ghosting functors cover
159  // the same set of elements/nodes. That being said, in general,
160  // we would expect G(e) > A(e). Here G(e) is the set of elements reserved
161  // by the geometric ghosting functors, and A(e) corresponds to
162  // the one covered by the algebraic ghosting functors.
163  // Therefore, we have to work only on local elements instead of
164  // ghosted + local elements. The ghosted solution might not be enough
165  // for ghosted+local elements. But it is always sufficient for local elements.
166  // After we set markers for all local elements, we will do a global
167  // communication to sync markers for ghosted elements from their owners.
168  all_elems = std::make_unique<ConstElemRange>(
169  _fe_problem.mesh().getMesh().active_local_elements_begin(),
170  _fe_problem.mesh().getMesh().active_local_elements_end());
171  }
172  else // This is not scalable but it might be useful for small-size problems
173  {
175  _fe_problem.getAuxiliarySystem().solution().localize(serialized_solution);
176  distributed_adaptivity = false;
177 
178  // For a replicated mesh or a serialized distributed mesh, the solution
179  // is serialized to everyone. Then we update markers for all active elements.
180  // In this case, we can avoid a global communication to update mesh.
181  // I do not know if it is a good idea, but it the old code behavior.
182  // We might not care about much since a replicated mesh
183  // or a serialized distributed mesh is not scalable anyway.
184  all_elems =
185  std::make_unique<ConstElemRange>(_fe_problem.mesh().getMesh().active_elements_begin(),
186  _fe_problem.mesh().getMesh().active_elements_end());
187  }
188 
189  FlagElementsThread fet(
190  _fe_problem, serialized_solution, _max_h_level, marker_name, !distributed_adaptivity);
191  Threads::parallel_reduce(*all_elems, fet);
193  }
194  }
195  else
196  {
197  // Compute the error for each active element
198  _error_estimator->estimate_error(_fe_problem.getNonlinearSystemBase(/*nl_sys=*/0).system(),
199  *_error);
200 
201  // Flag elements to be refined and coarsened
202  _mesh_refinement->flag_elements_by_error_fraction(*_error);
203 
204  if (_displaced_problem)
205  // Reuse the error vector and refine the displaced mesh
206  _displaced_mesh_refinement->flag_elements_by_error_fraction(*_error);
207  }
208 
209  // If the DisplacedProblem is active, undisplace the DisplacedMesh
210  // in preparation for refinement. We can't safely refine the
211  // DisplacedMesh directly, since the Hilbert keys computed on the
212  // inconsistenly-displaced Mesh are different on different
213  // processors, leading to inconsistent Hilbert keys. We must do
214  // this before the undisplaced Mesh is refined, so that the
215  // element and node numbering is still consistent.
216  if (_displaced_problem)
217  _displaced_problem->undisplaceMesh();
218 
219  // If markers are added to only local elements,
220  // we sync them here.
221  if (distributed_adaptivity)
222  _mesh_refinement->make_flags_parallel_consistent();
223 
224  if (_p_refinement_flag)
225  _mesh_refinement->switch_h_to_p_refinement();
226 
227  // Perform refinement and coarsening
228  mesh_changed = _mesh_refinement->refine_and_coarsen_elements();
229 
230  if (_displaced_problem && mesh_changed)
231  {
232  // If markers are added to only local elements,
233  // we sync them here.
234  if (distributed_adaptivity)
235  _displaced_mesh_refinement->make_flags_parallel_consistent();
236 
237  if (_p_refinement_flag)
238  _displaced_mesh_refinement->switch_h_to_p_refinement();
239 
240 #ifndef NDEBUG
241  bool displaced_mesh_changed =
242 #endif
243  _displaced_mesh_refinement->refine_and_coarsen_elements();
244 
245  // Since the undisplaced mesh changed, the displaced mesh better have changed!
246  mooseAssert(displaced_mesh_changed, "Undisplaced mesh changed, but displaced mesh did not!");
247  }
248 
249  if (mesh_changed && _print_mesh_changed)
250  {
251  _console << "\nMesh Changed:\n";
252  _mesh.printInfo();
253  _console << std::flush;
254  }
255 
256  return mesh_changed;
257 }
258 
259 bool
261 {
263 }
264 
265 void
266 Adaptivity::uniformRefine(MooseMesh * mesh, unsigned int level /*=libMesh::invalid_uint*/)
267 {
268  mooseAssert(mesh, "Mesh pointer must not be NULL");
269 
270  // NOTE: we are using a separate object here, since adaptivity may not be on, but we need to be
271  // able to do refinements
272  MeshRefinement mesh_refinement(*mesh);
273  if (level == libMesh::invalid_uint)
274  level = mesh->uniformRefineLevel();
275 
276  // Skip deletion and repartition will make uniform refinements run more
277  // efficiently, but at the same time, there might be extra ghosting elements.
278  // The number of layers of additional ghosting elements depends on the number
279  // of uniform refinement levels. This should happen only when you have a "fine enough"
280  // coarse mesh and want to refine the mesh by a few levels. Otherwise, it might
281  // introduce an unbalanced workload and too large ghosting domain.
282  if (mesh->skipDeletionRepartitionAfterRefine())
283  {
284  mesh->getMesh().skip_partitioning(true);
285  mesh->getMesh().allow_remote_element_removal(false);
286  mesh->needsRemoteElemDeletion(false);
287  }
288 
289  mesh_refinement.uniformly_refine(level);
290 }
291 
292 void
294 {
295  TIME_SECTION("uniformRefineWithProjection", 2, "Uniformly Refining and Reprojecting");
296 
297  // NOTE: we are using a separate object here, since adaptivity may not be on, but we need to be
298  // able to do refinements
299  MeshRefinement mesh_refinement(_mesh);
300  unsigned int level = _mesh.uniformRefineLevel();
301  MeshRefinement displaced_mesh_refinement(_displaced_problem ? _displaced_problem->mesh() : _mesh);
302 
303  // we have to go step by step so EquationSystems::reinit() won't freak out
304  for (unsigned int i = 0; i < level; i++)
305  {
306  // See comment above about why refining the displaced mesh is potentially unsafe.
307  if (_displaced_problem)
308  _displaced_problem->undisplaceMesh();
309 
310  mesh_refinement.uniformly_refine(1);
311 
312  if (_displaced_problem)
313  displaced_mesh_refinement.uniformly_refine(1);
315  }
316 }
317 
318 void
320 {
321  // check if Adaptivity has been initialized before turning on
322  if (state == true && !_initialized)
323  mooseError("Mesh adaptivity system not available");
324 
325  _mesh_refinement_on = state;
326 }
327 
328 void
329 Adaptivity::setTimeActive(Real start_time, Real stop_time)
330 {
331  _start_time = start_time;
332  _stop_time = stop_time;
333 }
334 
335 void
337 {
338  _use_new_system = true;
339 }
340 
341 void
342 Adaptivity::setMarkerVariableName(std::string marker_field)
343 {
344  _marker_variable_name = marker_field;
345 }
346 
347 void
349 {
350  _initial_marker_variable_name = marker_field;
351 }
352 
353 ErrorVector &
354 Adaptivity::getErrorVector(const std::string & indicator_field)
355 {
356  // Insert or retrieve error vector
357  auto insert_pair = moose_try_emplace(
358  _indicator_field_to_error_vector, indicator_field, std::make_unique<ErrorVector>());
359  return *insert_pair.first->second;
360 }
361 
362 void
364 {
365  TIME_SECTION("updateErrorVectors", 5, "Updating Error Vectors");
366 
367  // Resize all of the ErrorVectors in case the mesh has changed
368  for (const auto & it : _indicator_field_to_error_vector)
369  {
370  ErrorVector & vec = *(it.second);
371  vec.assign(_mesh.getMesh().max_elem_id(), 0);
372  }
373 
374  // Fill the vectors with the local contributions
376  Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), uevt);
377 
378  // Now sum across all processors
379  for (const auto & it : _indicator_field_to_error_vector)
380  _fe_problem.comm().sum((std::vector<float> &)*(it.second));
381 }
382 
383 bool
385 {
386  return _mesh_refinement_on && (_start_time <= _t && _t < _stop_time) && _step % _interval == 0;
387 }
388 
389 void
390 Adaptivity::doingPRefinement(const bool doing_p_refinement,
391  const MultiMooseEnum & disable_p_refinement_for_families)
392 {
393  _p_refinement_flag = doing_p_refinement;
394  _fe_problem.doingPRefinement(doing_p_refinement, disable_p_refinement_for_families);
395 }
396 
397 #endif // LIBMESH_ENABLE_AMR
virtual void update(bool update_libmesh_system=true)
Update the system (doing libMesh magic)
Definition: SystemBase.C:1211
bool adaptMesh(std::string marker_name=std::string())
Adapts the mesh based on the error estimator used.
Definition: Adaptivity.C:123
bool initialAdaptMesh()
Used during initial adaptivity.
Definition: Adaptivity.C:260
void doingPRefinement(bool doing_p_refinement, const MultiMooseEnum &disable_p_refinement_for_families)
Indicate whether the kind of adaptivity we&#39;re doing is p-refinement.
Definition: Adaptivity.C:390
ConstElemRange * getActiveLocalElementRange()
Return pointers to range objects for various types of ranges (local nodes, boundary elems...
Definition: MooseMesh.C:1040
std::shared_ptr< DisplacedProblem > _displaced_problem
Definition: Adaptivity.h:271
Real _stop_time
When adaptivity stops.
Definition: Adaptivity.h:293
std::string _marker_variable_name
Name of the marker variable if using the new adaptivity system.
Definition: Adaptivity.h:301
std::unique_ptr< MeshRefinement > _displaced_mesh_refinement
A mesh refinement object for displaced mesh.
Definition: Adaptivity.h:274
const unsigned int invalid_uint
std::unique_ptr< ErrorEstimator > _error_estimator
Error estimator to be used by the apps.
Definition: Adaptivity.h:267
bool _p_refinement_flag
Definition: Adaptivity.h:315
bool _initialized
on/off flag reporting if the adaptivity system has been initialized
Definition: Adaptivity.h:263
NumericVector< Number > & solution()
Definition: SystemBase.h:176
static void uniformRefine(MooseMesh *mesh, unsigned int level=libMesh::invalid_uint)
Performs uniform refinement of the passed Mesh object.
Definition: Adaptivity.C:266
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:284
virtual ~Adaptivity()
Definition: Adaptivity.C:54
void updateErrorVectors()
Update the ErrorVectors that have been requested through calls to getErrorVector().
Definition: Adaptivity.C:363
void init(unsigned int steps, unsigned int initial_steps)
Initialize and turn on adaptivity for the simulation.
Definition: Adaptivity.C:57
void setTimeActive(Real start_time, Real stop_time)
Sets the time when the adaptivity is active.
Definition: Adaptivity.C:329
MeshBase & mesh
virtual void doingPRefinement(bool doing_p_refinement, const MultiMooseEnum &disable_p_refinement_for_families) override
Indicate whether the kind of adaptivity we&#39;re doing is p-refinement.
const Parallel::Communicator & comm() const
std::unique_ptr< ErrorVector > _error
Error vector for use with the error estimator.
Definition: Adaptivity.h:269
void setUseNewSystem()
Tells this object we&#39;re using the "new" adaptivity system.
Definition: Adaptivity.C:336
std::pair< typename M::iterator, bool > moose_try_emplace(M &m, const typename M::key_type &k, Args &&... args)
Function to mirror the behavior of the C++17 std::map::try_emplace() method (no hint).
Definition: Moose.h:94
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
std::string _initial_marker_variable_name
Name of the initial marker variable if using the new adaptivity system.
Definition: Adaptivity.h:304
MooseMesh & _mesh
Definition: Adaptivity.h:258
std::map< std::string, std::unique_ptr< ErrorVector > > _indicator_field_to_error_vector
Stores pointers to ErrorVectors associated with indicator field names.
Definition: Adaptivity.h:313
auto max(const L &left, const R &right)
Real & _t
Time.
Definition: Adaptivity.h:285
virtual EquationSystems & es() override
bool _print_mesh_changed
True if we want to print out info when mesh has changed.
Definition: Adaptivity.h:282
virtual DofMap & dofMap()
Gets writeable reference to the dof map.
Definition: SystemBase.C:1131
std::unique_ptr< MeshRefinement > _mesh_refinement
A mesh refinement object to be used either with initial refinement or with Adaptivity.
Definition: Adaptivity.h:265
MeshBase & getMesh()
Accessor for the underlying libMesh Mesh object.
Definition: MooseMesh.C:3198
bool _mesh_refinement_on
on/off flag reporting if the adaptivity is being used
Definition: Adaptivity.h:261
Adaptivity(FEProblemBase &fe_problem)
Definition: Adaptivity.C:31
An inteface for the _console for outputting to the Console object.
MooseMesh wraps a libMesh::Mesh object and enhances its capabilities by caching additional data and s...
Definition: MooseMesh.h:88
virtual void meshChanged() override
Update data after a mesh change.
void setAdaptivityOn(bool state)
Allow adaptivity to be toggled programatically.
Definition: Adaptivity.C:319
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
Definition: MooseEnum.h:31
NonlinearSystemBase & getNonlinearSystemBase(const unsigned int sys_num)
unsigned int uniformRefineLevel() const
Returns the level of uniform refinement requested (zero if AMR is disabled).
Definition: MooseMesh.C:2966
void setErrorEstimator(const MooseEnum &error_estimator_name)
Set the error estimator.
Definition: Adaptivity.C:102
AuxiliarySystem & getAuxiliarySystem()
unsigned int _steps
steps of adaptivity to perform
Definition: Adaptivity.h:279
Interface for objects interacting with the PerfGraph.
virtual void close()=0
void uniformRefineWithProjection()
Performs uniform refinement on the meshes in the current object.
Definition: Adaptivity.C:293
void printInfo(std::ostream &os=libMesh::out, const unsigned int verbosity=0) const
Calls print_info() on the underlying Mesh.
Definition: MooseMesh.C:3212
unsigned int _initial_steps
the number of adaptivity steps to do at the beginning of simulation
Definition: Adaptivity.h:277
bool _use_new_system
Whether or not to use the "new" adaptivity system.
Definition: Adaptivity.h:298
unsigned int _interval
intreval between adaptivity runs
Definition: Adaptivity.h:289
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual std::shared_ptr< const DisplacedProblem > getDisplacedProblem() const
ErrorVector & getErrorVector(const std::string &indicator_field)
Get an ErrorVector that will be filled up with values corresponding to the indicator field name passe...
Definition: Adaptivity.C:354
virtual System & system() override
Get the reference to the libMesh system.
virtual MooseMesh & mesh() override
Real _start_time
When adaptivity start.
Definition: Adaptivity.h:291
void setErrorNorm(SystemNorm &sys_norm)
Set the error norm (FIXME: improve description)
Definition: Adaptivity.C:116
int & _step
Time Step.
Definition: Adaptivity.h:287
bool isAdaptivityDue()
Query if an adaptivity step should be performed at the current time / time step.
Definition: Adaptivity.C:384
const ConsoleStream _console
An instance of helper class to write streams to the Console objects.
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
void setMarkerVariableName(std::string marker_field)
Sets the name of the field variable to actually use to flag elements for refinement / coarsening...
Definition: Adaptivity.C:342
unsigned int _max_h_level
The maximum number of refinement levels.
Definition: Adaptivity.h:307
void setInitialMarkerVariableName(std::string marker_field)
Sets the name of the field variable to actually use to flag elements for initial refinement / coarsen...
Definition: Adaptivity.C:348
FEProblemBase & _fe_problem
Definition: Adaptivity.h:257
virtual void localize(std::vector< Number > &v_local) const =0