https://mooseframework.inl.gov
Adaptivity.C
Go to the documentation of this file.
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"
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 
34  : ConsoleStreamInterface(fe_problem.getMooseApp()),
35  PerfGraphInterface(fe_problem.getMooseApp().perfGraph(), "Adaptivity"),
36  ParallelObject(fe_problem.getMooseApp()),
37  _fe_problem(fe_problem),
38  _mesh(_fe_problem.mesh()),
39  _mesh_refinement_on(false),
40  _initialized(false),
41  _initial_steps(0),
42  _steps(0),
43  _print_mesh_changed(false),
44  _t(_fe_problem.time()),
45  _step(_fe_problem.timeStep()),
46  _interval(1),
47  _start_time(-std::numeric_limits<Real>::max()),
48  _stop_time(std::numeric_limits<Real>::max()),
49  _cycles_per_step(1),
50  _use_new_system(false),
51  _max_h_level(0),
52  _recompute_markers_during_cycles(false)
53 {
54 }
55 
57 
58 void
59 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.
67 
68  _mesh_refinement = std::make_unique<MeshRefinement>(_mesh);
69  _error = std::make_unique<ErrorVector>();
70 
72  es.parameters.set<bool>("adaptivity") = true;
73 
74  _initial_steps = initial_steps;
75  _steps = steps;
76  _p_refinement_flag = p_refinement;
77  _mesh_refinement_on = true;
78 
80  _mesh.doingPRefinement(true);
81 
82  _mesh_refinement->set_periodic_boundaries_ptr(
84 
85  // displaced problem
86  if (_displaced_problem != nullptr)
87  {
88  EquationSystems & displaced_es = _displaced_problem->es();
89  displaced_es.parameters.set<bool>("adaptivity") = true;
90 
92  _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  _displaced_mesh_refinement->set_periodic_boundaries_ptr(
100 
101  // TODO: This is currently an empty function on the DisplacedProblem... could it be removed?
102  _displaced_problem->initAdaptivity();
103  }
104 
105  // indicate the Adaptivity system has been initialized
106  _initialized = true;
107 }
108 
109 void
110 Adaptivity::setErrorEstimator(const MooseEnum & error_estimator_name)
111 {
112  if (error_estimator_name == "KellyErrorEstimator")
113  _error_estimator = std::make_unique<KellyErrorEstimator>();
114  else if (error_estimator_name == "LaplacianErrorEstimator")
115  _error_estimator = std::make_unique<LaplacianErrorEstimator>();
116  else if (error_estimator_name == "PatchRecoveryErrorEstimator")
117  _error_estimator = std::make_unique<PatchRecoveryErrorEstimator>();
118  else
119  mooseError(std::string("Unknown error_estimator selection: ") +
120  std::string(error_estimator_name));
121 }
122 
123 void
125 {
126  mooseAssert(_error_estimator, "error_estimator not initialized. Did you call init_adaptivity()?");
127  _error_estimator->error_norm = sys_norm;
128 }
129 
130 bool
131 Adaptivity::adaptMesh(std::string marker_name /*=std::string()*/)
132 {
133  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  if (marker_name.empty())
137  marker_name = _marker_variable_name;
138 
139  bool mesh_changed = false;
140 
141  // If mesh adaptivity is carried out in a distributed (scalable) way
142  bool distributed_adaptivity = false;
143 
144  if (_use_new_system)
145  {
146  if (!marker_name.empty()) // Only flag if a marker variable name has been set
147  {
148  _mesh_refinement->clean_refinement_flags();
149 
150  std::vector<Number> serialized_solution;
151 
152  auto distributed_mesh = dynamic_cast<DistributedMesh *>(&_fe_problem.mesh().getMesh());
153 
154  // Element range
155  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  if (distributed_mesh && !distributed_mesh->is_serial_on_zero())
161  {
162  // We update here to make sure local solution is up-to-date
164  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  all_elems = std::make_unique<ConstElemRange>(
177  _fe_problem.mesh().getMesh().active_local_elements_begin(),
178  _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  {
183  _fe_problem.getAuxiliarySystem().solution().localize(serialized_solution);
184  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  std::make_unique<ConstElemRange>(_fe_problem.mesh().getMesh().active_elements_begin(),
194  _fe_problem.mesh().getMesh().active_elements_end());
195  }
196 
197  FlagElementsThread fet(
198  _fe_problem, serialized_solution, _max_h_level, marker_name, !distributed_adaptivity);
199  Threads::parallel_reduce(*all_elems, fet);
201  }
202  }
203  else
204  {
205  // Compute the error for each active element
206  _error_estimator->estimate_error(_fe_problem.getNonlinearSystemBase(/*nl_sys=*/0).system(),
207  *_error);
208 
209  // Flag elements to be refined and coarsened
210  _mesh_refinement->flag_elements_by_error_fraction(*_error);
211 
212  if (_displaced_problem)
213  // Reuse the error vector and refine the displaced mesh
214  _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  if (_displaced_problem)
225  _displaced_problem->undisplaceMesh();
226 
227  // If markers are added to only local elements,
228  // we sync them here.
229  if (distributed_adaptivity)
230  _mesh_refinement->make_flags_parallel_consistent();
231 
232  if (_p_refinement_flag)
233  _mesh_refinement->switch_h_to_p_refinement();
234 
235  // Perform refinement and coarsening
236  mesh_changed = _mesh_refinement->refine_and_coarsen_elements();
237 
238  if (_displaced_problem && mesh_changed)
239  {
240  // If markers are added to only local elements,
241  // we sync them here.
242  if (distributed_adaptivity)
243  _displaced_mesh_refinement->make_flags_parallel_consistent();
244 
245  if (_p_refinement_flag)
246  _displaced_mesh_refinement->switch_h_to_p_refinement();
247 
248 #ifndef NDEBUG
249  bool displaced_mesh_changed =
250 #endif
251  _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  if (mesh_changed && _print_mesh_changed)
258  {
259  _console << "\nMesh Changed:\n";
260  _mesh.printInfo();
261  _console << std::flush;
262  }
263 
264  return mesh_changed;
265 }
266 
267 bool
269 {
271 }
272 
273 void
274 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  MeshRefinement mesh_refinement(*mesh);
281  if (level == libMesh::invalid_uint)
282  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  if (mesh->skipDeletionRepartitionAfterRefine())
291  {
292  mesh->getMesh().skip_partitioning(true);
293  mesh->getMesh().allow_remote_element_removal(false);
294  mesh->needsRemoteElemDeletion(false);
295  }
296 
297  mesh_refinement.uniformly_refine(level);
298 }
299 
300 void
302 {
303  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  MeshRefinement mesh_refinement(_mesh);
308  unsigned int level = _mesh.uniformRefineLevel();
309  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  for (unsigned int i = 0; i < level; i++)
313  {
314  // See comment above about why refining the displaced mesh is potentially unsafe.
315  if (_displaced_problem)
316  _displaced_problem->undisplaceMesh();
317 
318  mesh_refinement.uniformly_refine(1);
319 
320  if (_displaced_problem)
321  displaced_mesh_refinement.uniformly_refine(1);
323  /*intermediate_step=*/false, /*contract_mesh=*/true, /*clean_refinement_flags=*/true);
324  }
325 }
326 
327 void
329 {
330  // check if Adaptivity has been initialized before turning on
331  if (state == true && !_initialized)
332  mooseError("Mesh adaptivity system not available");
333 
334  _mesh_refinement_on = state;
335 }
336 
337 void
338 Adaptivity::setTimeActive(Real start_time, Real stop_time)
339 {
340  _start_time = start_time;
341  _stop_time = stop_time;
342 }
343 
344 void
346 {
347  _use_new_system = true;
348 }
349 
350 void
351 Adaptivity::setMarkerVariableName(std::string marker_field)
352 {
353  _marker_variable_name = marker_field;
354 }
355 
356 void
358 {
359  _initial_marker_variable_name = marker_field;
360 }
361 
362 ErrorVector &
363 Adaptivity::getErrorVector(const std::string & indicator_field)
364 {
365  // Insert or retrieve error vector
366  auto insert_pair = moose_try_emplace(
367  _indicator_field_to_error_vector, indicator_field, std::make_unique<ErrorVector>());
368  return *insert_pair.first->second;
369 }
370 
371 void
373 {
374  TIME_SECTION("updateErrorVectors", 5, "Updating Error Vectors");
375 
376  // Resize all of the ErrorVectors in case the mesh has changed
377  for (const auto & it : _indicator_field_to_error_vector)
378  {
379  ErrorVector & vec = *(it.second);
380  vec.assign(_mesh.getMesh().max_elem_id(), 0);
381  }
382 
383  // Fill the vectors with the local contributions
385  Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), uevt);
386 
387  // Now sum across all processors
388  for (const auto & it : _indicator_field_to_error_vector)
389  _fe_problem.comm().sum((std::vector<float> &)*(it.second));
390 }
391 
392 bool
394 {
395  return _mesh_refinement_on && (_start_time <= _t && _t < _stop_time) && _step % _interval == 0;
396 }
397 
398 #endif // LIBMESH_ENABLE_AMR
bool adaptMesh(std::string marker_name=std::string())
Adapts the mesh based on the error estimator used.
Definition: Adaptivity.C:131
bool initialAdaptMesh()
Used during initial adaptivity.
Definition: Adaptivity.C:268
virtual void meshChanged(bool intermediate_change, bool contract_mesh, bool clean_refinement_flags)
Update data after a mesh change.
std::unique_ptr< libMesh::ErrorVector > _error
Error vector for use with the error estimator.
Definition: Adaptivity.h:269
libMesh::ConstElemRange * getActiveLocalElementRange()
Return pointers to range objects for various types of ranges (local nodes, boundary elems...
Definition: MooseMesh.C:1235
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
const unsigned int invalid_uint
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:195
static void uniformRefine(MooseMesh *mesh, unsigned int level=libMesh::invalid_uint)
Performs uniform refinement of the passed Mesh object.
Definition: Adaptivity.C:274
void skip_partitioning(bool skip)
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
virtual ~Adaptivity()
Definition: Adaptivity.C:56
std::unique_ptr< libMesh::ErrorEstimator > _error_estimator
Error estimator to be used by the apps.
Definition: Adaptivity.h:267
void updateErrorVectors()
Update the ErrorVectors that have been requested through calls to getErrorVector().
Definition: Adaptivity.C:372
std::unique_ptr< libMesh::MeshRefinement > _mesh_refinement
A mesh refinement object to be used either with initial refinement or with Adaptivity.
Definition: Adaptivity.h:265
void setTimeActive(Real start_time, Real stop_time)
Sets the time when the adaptivity is active.
Definition: Adaptivity.C:338
MeshBase & mesh
void init(const unsigned int steps, const unsigned int initial_steps, const bool p_refinement)
Initialize and turn on adaptivity for the simulation.
Definition: Adaptivity.C:59
const Parallel::Communicator & comm() const
void setUseNewSystem()
Tells this object we&#39;re using the "new" adaptivity system.
Definition: Adaptivity.C:345
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
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:93
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
auto max(const L &left, const R &right)
Real & _t
Time.
Definition: Adaptivity.h:285
void update()
Update the system (doing libMesh magic)
Definition: SystemBase.C:1245
bool _print_mesh_changed
True if we want to print out info when mesh has changed.
Definition: Adaptivity.h:282
virtual libMesh::DofMap & dofMap()
Gets writeable reference to the dof map.
Definition: SystemBase.C:1165
void allow_remote_element_removal(bool allow)
MeshBase & getMesh()
Accessor for the underlying libMesh Mesh object.
Definition: MooseMesh.C:3443
PeriodicBoundaries * get_periodic_boundaries()
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:33
virtual libMesh::EquationSystems & es() override
An inteface for the _console for outputting to the Console object.
virtual dof_id_type max_elem_id() const=0
MooseMesh wraps a libMesh::Mesh object and enhances its capabilities by caching additional data and s...
Definition: MooseMesh.h:88
void setAdaptivityOn(bool state)
Allow adaptivity to be toggled programatically.
Definition: Adaptivity.C:328
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
Definition: MooseEnum.h:33
std::unique_ptr< libMesh::MeshRefinement > _displaced_mesh_refinement
A mesh refinement object for displaced mesh.
Definition: Adaptivity.h:274
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:3211
void setErrorEstimator(const MooseEnum &error_estimator_name)
Set the error estimator.
Definition: Adaptivity.C:110
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:301
void printInfo(std::ostream &os=libMesh::out, const unsigned int verbosity=0) const
Calls print_info() on the underlying Mesh.
Definition: MooseMesh.C:3457
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
std::map< std::string, std::unique_ptr< libMesh::ErrorVector > > _indicator_field_to_error_vector
Stores pointers to ErrorVectors associated with indicator field names.
Definition: Adaptivity.h:313
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
libMesh::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:363
T & set(const std::string &)
virtual MooseMesh & mesh() override
Real _start_time
When adaptivity start.
Definition: Adaptivity.h:291
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:393
const ConsoleStream _console
An instance of helper class to write streams to the Console objects.
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:351
void doingPRefinement(bool doing_p_refinement)
Indicate whether the kind of adaptivity we&#39;re doing is p-refinement.
Definition: MooseMesh.h:1347
unsigned int _max_h_level
The maximum number of refinement levels.
Definition: Adaptivity.h:307
void setErrorNorm(libMesh::SystemNorm &sys_norm)
Set the error norm (FIXME: improve description)
Definition: Adaptivity.C:124
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:357
FEProblemBase & _fe_problem
Definition: Adaptivity.h:257
void uniformly_refine(unsigned int n=1)
virtual void localize(std::vector< T > &v_local) const=0
virtual libMesh::System & system() override
Get the reference to the libMesh system.