libMesh::AdjointRefinementEstimator Class Reference

#include <adjoint_refinement_estimator.h>

Inheritance diagram for libMesh::AdjointRefinementEstimator:

Public Types

typedef std::map< std::pair< const System *, unsigned int >, ErrorVector * > ErrorMap
 

Public Member Functions

 AdjointRefinementEstimator ()
 
 ~AdjointRefinementEstimator ()
 
QoISetqoi_set ()
 
const QoISetqoi_set () const
 
virtual void estimate_error (const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=libmesh_nullptr, bool estimate_parent_error=false)
 
Numberget_global_QoI_error_estimate (unsigned int qoi_index)
 
virtual ErrorEstimatorType type () const
 
virtual void estimate_errors (const EquationSystems &equation_systems, ErrorVector &error_per_cell, const std::map< const System *, SystemNorm > &error_norms, const std::map< const System *, const NumericVector< Number > * > *solution_vectors=libmesh_nullptr, bool estimate_parent_error=false)
 
virtual void estimate_errors (const EquationSystems &equation_systems, ErrorMap &errors_per_cell, const std::map< const System *, const NumericVector< Number > * > *solution_vectors=libmesh_nullptr, bool estimate_parent_error=false)
 

Public Attributes

unsigned char number_h_refinements
 
unsigned char number_p_refinements
 
SystemNorm error_norm
 

Protected Member Functions

void reduce_error (std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm LIBMESH_CAN_DEFAULT_TO_COMMWORLD) const
 

Protected Attributes

std::vector< Numbercomputed_global_QoI_errors
 
QoISet _qoi_set
 

Detailed Description

This class implements a ``brute force'' goal-oriented error estimator which computes an estimate of error in a quantity of interest based on the residual of the current coarse grid primal solution as weighted against an adjoint solution on a uniformly refined (in h and/or p, for an arbitrary number of levels) grid.

Author
Roy H. Stogner
Date
2009

Definition at line 47 of file adjoint_refinement_estimator.h.

Member Typedef Documentation

typedef std::map<std::pair<const System *, unsigned int>, ErrorVector *> libMesh::ErrorEstimator::ErrorMap
inherited

When calculating many error vectors at once, we need a data structure to hold them all

Definition at line 112 of file error_estimator.h.

Constructor & Destructor Documentation

libMesh::AdjointRefinementEstimator::AdjointRefinementEstimator ( )
inline

Constructor. Sets the most common default parameter values.

Definition at line 54 of file adjoint_refinement_estimator.h.

References libMesh::ErrorEstimator::error_norm, and libMesh::INVALID_NORM.

54  :
58  _qoi_set(QoISet())
59  {
60  // We're not actually going to use error_norm; our norms are
61  // absolute values of QoI error.
63  }
libMesh::AdjointRefinementEstimator::~AdjointRefinementEstimator ( )
inline

Destructor.

Definition at line 68 of file adjoint_refinement_estimator.h.

68 {}

Member Function Documentation

void libMesh::AdjointRefinementEstimator::estimate_error ( const System system,
ErrorVector error_per_cell,
const NumericVector< Number > *  solution_vector = libmesh_nullptr,
bool  estimate_parent_error = false 
)
virtual

This function does uniform refinements and an adjoint solve to get an adjoint solution on each cell, then estimates the error by finding the weighted residual of the coarse solution with the fine adjoint solution.

system.solve() and system.assembly() must be called, and so should have no side effects.

Only the provided system is solved on the refined mesh; we don't support adjoint solves on loosely coupled collections of Systems.

The estimated error is output in the vector error_per_cell

Implements libMesh::ErrorEstimator.

Definition at line 69 of file adjoint_refinement_estimator.C.

References _qoi_set, std::abs(), libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), libMesh::System::add_vector(), libMesh::System::adjoint_solve(), libMesh::MeshBase::allow_renumbering(), libMesh::NumericVector< T >::build(), libMesh::NumericVector< T >::clear(), libMesh::NumericVector< T >::close(), libMesh::ParallelObject::comm(), computed_global_QoI_errors, libMesh::System::current_local_solution, libMesh::DofMap::dof_indices(), libMesh::NumericVector< T >::dot(), libMesh::DofMap::enforce_constraints_exactly(), libMesh::ErrorVectorReal, libMesh::Elem::find_point_neighbors(), libMesh::System::get_adjoint_solution(), libMesh::System::get_dof_map(), libMesh::System::get_equation_systems(), libMesh::EquationSystems::get_mesh(), libMesh::DofMap::get_send_list(), libMesh::System::get_vector(), libMesh::GHOSTED, libMesh::DofMap::has_adjoint_dirichlet_boundaries(), libMesh::QoISet::has_index(), libMesh::DofObject::id(), libMesh::NumericVector< T >::init(), libMesh::System::is_adjoint_already_solved(), libMesh::libmesh_assert(), libmesh_nullptr, libMesh::NumericVector< T >::localize(), libMesh::MeshBase::max_elem_id(), mesh, libMesh::System::n_dofs(), libMesh::MeshBase::n_elem(), libMesh::System::n_local_dofs(), libMesh::Elem::n_nodes(), libMesh::Elem::node_ref(), number_h_refinements, number_p_refinements, libMesh::Elem::parent(), libMesh::MeshBase::partitioner(), libMesh::System::project_solution_on_reinit(), libMesh::System::qoi, libMesh::Real, libMesh::ErrorEstimator::reduce_error(), libMesh::EquationSystems::reinit(), libMesh::System::solution, libMesh::NumericVector< T >::swap(), libMesh::MeshRefinement::uniformly_coarsen(), libMesh::MeshRefinement::uniformly_p_coarsen(), libMesh::MeshRefinement::uniformly_p_refine(), libMesh::MeshRefinement::uniformly_refine(), libMesh::System::update(), libMesh::System::vectors_begin(), libMesh::System::vectors_end(), and libMesh::QoISet::weight().

Referenced by qoi_set().

73 {
74  // We have to break the rules here, because we can't refine a const System
75  System & system = const_cast<System &>(_system);
76 
77  // An EquationSystems reference will be convenient.
78  EquationSystems & es = system.get_equation_systems();
79 
80  // The current mesh
81  MeshBase & mesh = es.get_mesh();
82 
83  // Get coarse grid adjoint solutions. This should be a relatively
84  // quick (especially with preconditioner reuse) way to get a good
85  // initial guess for the fine grid adjoint solutions. More
86  // importantly, subtracting off a coarse adjoint approximation gives
87  // us better local error indication, and subtracting off *some* lift
88  // function is necessary for correctness if we have heterogeneous
89  // adjoint Dirichlet conditions.
90 
91  // Solve the adjoint problem(s) on the coarse FE space
92  // Only if the user didn't already solve it for us
93  if (!system.is_adjoint_already_solved())
94  system.adjoint_solve(_qoi_set);
95 
96  // Loop over all the adjoint problems and, if any have heterogenous
97  // Dirichlet conditions, get the corresponding coarse lift
98  // function(s)
99  for (std::size_t j=0; j != system.qoi.size(); j++)
100  {
101  // Skip this QoI if it is not in the QoI Set or if there are no
102  // heterogeneous Dirichlet boundaries for it
103  if (_qoi_set.has_index(j) &&
104  system.get_dof_map().has_adjoint_dirichlet_boundaries(j))
105  {
106  std::ostringstream liftfunc_name;
107  liftfunc_name << "adjoint_lift_function" << j;
108  NumericVector<Number> & liftvec =
109  system.add_vector(liftfunc_name.str());
110 
111  system.get_dof_map().enforce_constraints_exactly
112  (system, &liftvec, true);
113  }
114  }
115 
116  // We'll want to back up all coarse grid vectors
117  std::map<std::string, NumericVector<Number> *> coarse_vectors;
118  for (System::vectors_iterator vec = system.vectors_begin(); vec !=
119  system.vectors_end(); ++vec)
120  {
121  // The (string) name of this vector
122  const std::string & var_name = vec->first;
123 
124  coarse_vectors[var_name] = vec->second->clone().release();
125  }
126  // Back up the coarse solution and coarse local solution
127  NumericVector<Number> * coarse_solution =
128  system.solution->clone().release();
129  NumericVector<Number> * coarse_local_solution =
130  system.current_local_solution->clone().release();
131 
132  // And we'll need to temporarily change solution projection settings
133  bool old_projection_setting;
134  old_projection_setting = system.project_solution_on_reinit();
135 
136  // Make sure the solution is projected when we refine the mesh
137  system.project_solution_on_reinit() = true;
138 
139  // And it'll be best to avoid any repartitioning
140  UniquePtr<Partitioner> old_partitioner(mesh.partitioner().release());
141 
142  // And we can't allow any renumbering
143  const bool old_renumbering_setting = mesh.allow_renumbering();
144  mesh.allow_renumbering(false);
145 
146  // Use a non-standard solution vector if necessary
147  if (solution_vector && solution_vector != system.solution.get())
148  {
149  NumericVector<Number> * newsol =
150  const_cast<NumericVector<Number> *> (solution_vector);
151  newsol->swap(*system.solution);
152  system.update();
153  }
154 
155  // Resize the error_per_cell vector to be
156  // the number of elements, initialized to 0.
157  error_per_cell.clear();
158  error_per_cell.resize (mesh.max_elem_id(), 0.);
159 
160 #ifndef NDEBUG
161  // n_coarse_elem is only used in an assertion later so
162  // avoid declaring it unless asserts are active.
163  const dof_id_type n_coarse_elem = mesh.n_elem();
164 #endif
165 
166  // Uniformly refine the mesh
167  MeshRefinement mesh_refinement(mesh);
168 
170 
171  // FIXME: this may break if there is more than one System
172  // on this mesh but estimate_error was still called instead of
173  // estimate_errors
174  for (unsigned int i = 0; i != number_h_refinements; ++i)
175  {
176  mesh_refinement.uniformly_refine(1);
177  es.reinit();
178  }
179 
180  for (unsigned int i = 0; i != number_p_refinements; ++i)
181  {
182  mesh_refinement.uniformly_p_refine(1);
183  es.reinit();
184  }
185 
186  // Copy the projected coarse grid solutions, which will be
187  // overwritten by solve()
188  std::vector<NumericVector<Number> *> coarse_adjoints;
189  for (std::size_t j=0; j != system.qoi.size(); j++)
190  {
191  if (_qoi_set.has_index(j))
192  {
193  NumericVector<Number> * coarse_adjoint =
194  NumericVector<Number>::build(mesh.comm()).release();
195 
196  // Can do "fast" init since we're overwriting this in a sec
197  coarse_adjoint->init(system.get_adjoint_solution(j),
198  /* fast = */ true);
199 
200  *coarse_adjoint = system.get_adjoint_solution(j);
201 
202  coarse_adjoints.push_back(coarse_adjoint);
203  }
204  else
205  coarse_adjoints.push_back(static_cast<NumericVector<Number> *>(libmesh_nullptr));
206  }
207 
208  // Rebuild the rhs with the projected primal solution
209  (dynamic_cast<ImplicitSystem &>(system)).assembly(true, false);
210  NumericVector<Number> & projected_residual = (dynamic_cast<ExplicitSystem &>(system)).get_vector("RHS Vector");
211  projected_residual.close();
212 
213  // Solve the adjoint problem(s) on the refined FE space
214  system.adjoint_solve(_qoi_set);
215 
216  // Now that we have the refined adjoint solution and the projected primal solution,
217  // we first compute the global QoI error estimate
218 
219  // Resize the computed_global_QoI_errors vector to hold the error estimates for each QoI
220  computed_global_QoI_errors.resize(system.qoi.size());
221 
222  // Loop over all the adjoint solutions and get the QoI error
223  // contributions from all of them. While we're looping anyway we'll
224  // pull off the coarse adjoints
225  for (std::size_t j=0; j != system.qoi.size(); j++)
226  {
227  // Skip this QoI if not in the QoI Set
228  if (_qoi_set.has_index(j))
229  {
230  // If the adjoint solution has heterogeneous dirichlet
231  // values, then to get a proper error estimate here we need
232  // to subtract off a coarse grid lift function. In any case
233  // we can get a better error estimate by separating off a
234  // coarse representation of the adjoint solution, so we'll
235  // use that for our lift function.
236  system.get_adjoint_solution(j) -= *coarse_adjoints[j];
237 
238  computed_global_QoI_errors[j] = projected_residual.dot(system.get_adjoint_solution(j));
239  }
240  }
241 
242  // Done with the global error estimates, now construct the element wise error indicators
243 
244  // We ought to account for 'spill-over' effects while computing the
245  // element error indicators This happens because the same dof is
246  // shared by multiple elements, one way of mitigating this is to
247  // scale the contribution from each dof by the number of elements it
248  // belongs to We first obtain the number of elements each node
249  // belongs to
250 
251  // A map that relates a node id to an int that will tell us how many elements it is a node of
252  LIBMESH_BEST_UNORDERED_MAP<dof_id_type, unsigned int>shared_element_count;
253 
254  // To fill this map, we will loop over elements, and then in each element, we will loop
255  // over the nodes each element contains, and then query it for the number of coarse
256  // grid elements it was a node of
257 
258  // Keep track of which nodes we have already dealt with
259  LIBMESH_BEST_UNORDERED_SET<dof_id_type> processed_node_ids;
260 
261  // We will be iterating over all the active elements in the fine mesh that live on
262  // this processor
263  {
264  MeshBase::const_element_iterator elem_it = mesh.active_local_elements_begin();
265  const MeshBase::const_element_iterator elem_end = mesh.active_local_elements_end();
266 
267  // Start loop over elems
268  for(; elem_it != elem_end; ++elem_it)
269  {
270  // Pointer to this element
271  const Elem * elem = *elem_it;
272 
273  // Loop over the nodes in the element
274  for(unsigned int n=0; n != elem->n_nodes(); ++n)
275  {
276  // Get a reference to the current node
277  const Node & node = elem->node_ref(n);
278 
279  // Get the id of this node
280  dof_id_type node_id = node.id();
281 
282  // If we havent already processed this node, do so now
283  if(processed_node_ids.find(node_id) == processed_node_ids.end())
284  {
285  // Declare a neighbor_set to be filled by the find_point_neighbors
286  std::set<const Elem *> fine_grid_neighbor_set;
287 
288  // Call find_point_neighbors to fill the neighbor_set
289  elem->find_point_neighbors(node, fine_grid_neighbor_set);
290 
291  // A vector to hold the coarse grid parents neighbors
292  std::vector<dof_id_type> coarse_grid_neighbors;
293 
294  // Iterators over the fine grid neighbors set
295  std::set<const Elem *>::iterator fine_neighbor_it = fine_grid_neighbor_set.begin();
296  const std::set<const Elem *>::iterator fine_neighbor_end = fine_grid_neighbor_set.end();
297 
298  // Loop over all the fine neighbors of this node
299  for(; fine_neighbor_it != fine_neighbor_end ; ++fine_neighbor_it)
300  {
301  // Pointer to the current fine neighbor element
302  const Elem * fine_elem = *fine_neighbor_it;
303 
304  // Find the element id for the corresponding coarse grid element
305  const Elem * coarse_elem = fine_elem;
306  for (unsigned int j = 0; j != number_h_refinements; ++j)
307  {
308  libmesh_assert (coarse_elem->parent());
309 
310  coarse_elem = coarse_elem->parent();
311  }
312 
313  // Loop over the existing coarse neighbors and check if this one is
314  // already in there
315  const dof_id_type coarse_id = coarse_elem->id();
316  std::size_t j = 0;
317  for (; j != coarse_grid_neighbors.size(); j++)
318  {
319  // If the set already contains this element break out of the loop
320  if(coarse_grid_neighbors[j] == coarse_id)
321  {
322  break;
323  }
324  }
325 
326  // If we didn't leave the loop even at the last element,
327  // this is a new neighbour, put in the coarse_grid_neighbor_set
328  if(j == coarse_grid_neighbors.size())
329  {
330  coarse_grid_neighbors.push_back(coarse_id);
331  }
332 
333  } // End loop over fine neighbors
334 
335  // Set the shared_neighbour index for this node to the
336  // size of the coarse grid neighbor set
337  shared_element_count[node_id] =
338  cast_int<unsigned int>(coarse_grid_neighbors.size());
339 
340  // Add this node to processed_node_ids vector
341  processed_node_ids.insert(node_id);
342 
343  } // End if not processed node
344 
345  } // End loop over nodes
346 
347  } // End loop over elems
348  }
349 
350  // Get a DoF map, will be used to get the nodal dof_indices for each element
351  DofMap & dof_map = system.get_dof_map();
352 
353  // The global DOF indices, we will use these later on when we compute the element wise indicators
354  std::vector<dof_id_type> dof_indices;
355 
356  // Localize the global rhs and adjoint solution vectors (which might be shared on multiple processsors) onto a
357  // local ghosted vector, this ensures each processor has all the dof_indices to compute an error indicator for
358  // an element it owns
359  UniquePtr<NumericVector<Number> > localized_projected_residual = NumericVector<Number>::build(system.comm());
360  localized_projected_residual->init(system.n_dofs(), system.n_local_dofs(), system.get_dof_map().get_send_list(), false, GHOSTED);
361  projected_residual.localize(*localized_projected_residual, system.get_dof_map().get_send_list());
362 
363  // Each adjoint solution will also require ghosting; for efficiency we'll reuse the same memory
364  UniquePtr<NumericVector<Number> > localized_adjoint_solution = NumericVector<Number>::build(system.comm());
365  localized_adjoint_solution->init(system.n_dofs(), system.n_local_dofs(), system.get_dof_map().get_send_list(), false, GHOSTED);
366 
367  // We will loop over each adjoint solution, localize that adjoint
368  // solution and then loop over local elements
369  for (std::size_t i=0; i != system.qoi.size(); i++)
370  {
371  // Skip this QoI if not in the QoI Set
372  if (_qoi_set.has_index(i))
373  {
374  // Get the weight for the current QoI
375  Real error_weight = _qoi_set.weight(i);
376 
377  (system.get_adjoint_solution(i)).localize(*localized_adjoint_solution, system.get_dof_map().get_send_list());
378 
379  // Loop over elements
380  MeshBase::const_element_iterator elem_it = mesh.active_local_elements_begin();
381  const MeshBase::const_element_iterator elem_end = mesh.active_local_elements_end();
382 
383  for(; elem_it != elem_end; ++elem_it)
384  {
385  // Pointer to the element
386  const Elem * elem = *elem_it;
387 
388  // Go up number_h_refinements levels up to find the coarse parent
389  const Elem * coarse = elem;
390 
391  for (unsigned int j = 0; j != number_h_refinements; ++j)
392  {
393  libmesh_assert (coarse->parent());
394 
395  coarse = coarse->parent();
396  }
397 
398  const dof_id_type e_id = coarse->id();
399 
400  // Get the local to global degree of freedom maps for this element
401  dof_map.dof_indices (elem, dof_indices);
402 
403  // We will have to manually do the dot products.
404  Number local_contribution = 0.;
405 
406  for (std::size_t j=0; j != dof_indices.size(); j++)
407  {
408  // The contribution to the error indicator for this element from the current QoI
409  local_contribution += (*localized_projected_residual)(dof_indices[j]) * (*localized_adjoint_solution)(dof_indices[j]);
410  }
411 
412  // Multiply by the error weight for this QoI
413  local_contribution *= error_weight;
414 
415  // FIXME: we're throwing away information in the
416  // --enable-complex case
417  error_per_cell[e_id] += static_cast<ErrorVectorReal>
418  (std::abs(local_contribution));
419 
420  } // End loop over elements
421 
422  } // End if belong to QoI set
423 
424  } // End loop over QoIs
425 
426  for (std::size_t j=0; j != system.qoi.size(); j++)
427  {
428  if (_qoi_set.has_index(j))
429  {
430  delete coarse_adjoints[j];
431  }
432  }
433 
434  // Don't bother projecting the solution; we'll restore from backup
435  // after coarsening
436  system.project_solution_on_reinit() = false;
437 
438  // Uniformly coarsen the mesh, without projecting the solution
440 
441  for (unsigned int i = 0; i != number_h_refinements; ++i)
442  {
443  mesh_refinement.uniformly_coarsen(1);
444  // FIXME - should the reinits here be necessary? - RHS
445  es.reinit();
446  }
447 
448  for (unsigned int i = 0; i != number_p_refinements; ++i)
449  {
450  mesh_refinement.uniformly_p_coarsen(1);
451  es.reinit();
452  }
453 
454  // We should be back where we started
455  libmesh_assert_equal_to (n_coarse_elem, mesh.n_elem());
456 
457  // Restore old solutions and clean up the heap
458  system.project_solution_on_reinit() = old_projection_setting;
459 
460  // Restore the coarse solution vectors and delete their copies
461  *system.solution = *coarse_solution;
462  delete coarse_solution;
463  *system.current_local_solution = *coarse_local_solution;
464  delete coarse_local_solution;
465 
466  for (System::vectors_iterator vec = system.vectors_begin(); vec !=
467  system.vectors_end(); ++vec)
468  {
469  // The (string) name of this vector
470  const std::string & var_name = vec->first;
471 
472  // If it's a vector we already had (and not a newly created
473  // vector like an adjoint rhs), we need to restore it.
474  std::map<std::string, NumericVector<Number> *>::iterator it =
475  coarse_vectors.find(var_name);
476  if (it != coarse_vectors.end())
477  {
478  NumericVector<Number> * coarsevec = it->second;
479  system.get_vector(var_name) = *coarsevec;
480 
481  coarsevec->clear();
482  delete coarsevec;
483  }
484  }
485 
486  // Restore old partitioner and renumbering settings
487  mesh.partitioner().reset(old_partitioner.release());
488  mesh.allow_renumbering(old_renumbering_setting);
489 
490  // Fiinally sum the vector of estimated error values.
491  this->reduce_error(error_per_cell, system.comm());
492 
493  // We don't take a square root here; this is a goal-oriented
494  // estimate not a Hilbert norm estimate.
495 } // end estimate_error function
double abs(double a)
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
libmesh_assert(j)
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Real weight(unsigned int) const
Definition: qoi_set.h:240
void reduce_error(std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm LIBMESH_CAN_DEFAULT_TO_COMMWORLD) const
std::map< std::string, NumericVector< Number > * >::iterator vectors_iterator
Definition: system.h:747
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void get(const std::vector< numeric_index_type > &index, T *values) const
bool has_index(unsigned int) const
Definition: qoi_set.h:221
uint8_t dof_id_type
Definition: id_types.h:64
void libMesh::ErrorEstimator::estimate_errors ( const EquationSystems equation_systems,
ErrorVector error_per_cell,
const std::map< const System *, SystemNorm > &  error_norms,
const std::map< const System *, const NumericVector< Number > * > *  solution_vectors = libmesh_nullptr,
bool  estimate_parent_error = false 
)
virtualinherited

This virtual function can be redefined in derived classes, but by default computes the sum of the error_per_cell for each system in the equation_systems.

Currently this function ignores the error_norm member variable, and uses the function argument error_norms instead.

This function is named estimate_errors instead of estimate_error because otherwise C++ can get confused.

Reimplemented in libMesh::UniformRefinementEstimator.

Definition at line 48 of file error_estimator.C.

References libMesh::ErrorEstimator::error_norm, libMesh::ErrorEstimator::estimate_error(), libMesh::EquationSystems::get_system(), libmesh_nullptr, libMesh::EquationSystems::n_systems(), and libMesh::sys.

Referenced by libMesh::ErrorEstimator::~ErrorEstimator().

53 {
54  SystemNorm old_error_norm = this->error_norm;
55 
56  // Sum the error values from each system
57  for (unsigned int s = 0; s != equation_systems.n_systems(); ++s)
58  {
59  ErrorVector system_error_per_cell;
60  const System & sys = equation_systems.get_system(s);
61  if (error_norms.find(&sys) == error_norms.end())
62  this->error_norm = old_error_norm;
63  else
64  this->error_norm = error_norms.find(&sys)->second;
65 
66  const NumericVector<Number> * solution_vector = libmesh_nullptr;
67  if (solution_vectors &&
68  solution_vectors->find(&sys) != solution_vectors->end())
69  solution_vector = solution_vectors->find(&sys)->second;
70 
71  this->estimate_error(sys, system_error_per_cell,
72  solution_vector, estimate_parent_error);
73 
74  if (s)
75  {
76  libmesh_assert_equal_to (error_per_cell.size(), system_error_per_cell.size());
77  for (std::size_t i=0; i != error_per_cell.size(); ++i)
78  error_per_cell[i] += system_error_per_cell[i];
79  }
80  else
81  error_per_cell = system_error_per_cell;
82  }
83 
84  // Restore our old state before returning
85  this->error_norm = old_error_norm;
86 }
ImplicitSystem & sys
const class libmesh_nullptr_t libmesh_nullptr
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=libmesh_nullptr, bool estimate_parent_error=false)=0
void libMesh::ErrorEstimator::estimate_errors ( const EquationSystems equation_systems,
ErrorMap errors_per_cell,
const std::map< const System *, const NumericVector< Number > * > *  solution_vectors = libmesh_nullptr,
bool  estimate_parent_error = false 
)
virtualinherited

This virtual function can be redefined in derived classes, but by default it calls estimate_error repeatedly to calculate the requested error vectors.

Currently this function ignores the error_norm.weight() values because it calculates each variable's error individually, unscaled.

The user selects which errors get computed by filling a map with error vectors: If errors_per_cell[&system][v] exists, it will be filled with the error values in variable v of system

FIXME: This is a default implementation - derived classes should reimplement it for efficiency.

Reimplemented in libMesh::UniformRefinementEstimator.

Definition at line 94 of file error_estimator.C.

References libMesh::ErrorEstimator::error_norm, libMesh::ErrorEstimator::estimate_error(), libMesh::EquationSystems::get_system(), libmesh_nullptr, libMesh::EquationSystems::n_systems(), n_vars, libMesh::System::n_vars(), libMesh::sys, and libMesh::SystemNorm::type().

98 {
99  SystemNorm old_error_norm = this->error_norm;
100 
101  // Find the requested error values from each system
102  for (unsigned int s = 0; s != equation_systems.n_systems(); ++s)
103  {
104  const System & sys = equation_systems.get_system(s);
105 
106  unsigned int n_vars = sys.n_vars();
107 
108  for (unsigned int v = 0; v != n_vars; ++v)
109  {
110  // Only fill in ErrorVectors the user asks for
111  if (errors_per_cell.find(std::make_pair(&sys, v)) ==
112  errors_per_cell.end())
113  continue;
114 
115  // Calculate error in only one variable
116  std::vector<Real> weights(n_vars, 0.0);
117  weights[v] = 1.0;
118  this->error_norm =
119  SystemNorm(std::vector<FEMNormType>(n_vars, old_error_norm.type(v)),
120  weights);
121 
122  const NumericVector<Number> * solution_vector = libmesh_nullptr;
123  if (solution_vectors &&
124  solution_vectors->find(&sys) != solution_vectors->end())
125  solution_vector = solution_vectors->find(&sys)->second;
126 
127  this->estimate_error
128  (sys, *errors_per_cell[std::make_pair(&sys, v)],
129  solution_vector, estimate_parent_error);
130  }
131  }
132 
133  // Restore our old state before returning
134  this->error_norm = old_error_norm;
135 }
ImplicitSystem & sys
const class libmesh_nullptr_t libmesh_nullptr
const unsigned int n_vars
Definition: tecplot_io.C:68
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=libmesh_nullptr, bool estimate_parent_error=false)=0
Number& libMesh::AdjointRefinementEstimator::get_global_QoI_error_estimate ( unsigned int  qoi_index)
inline

This is an accessor function to access the computed global QoI error estimates

Definition at line 107 of file adjoint_refinement_estimator.h.

References computed_global_QoI_errors.

108  {
109  return computed_global_QoI_errors[qoi_index];
110  }
QoISet& libMesh::AdjointRefinementEstimator::qoi_set ( )
inline

Access to the QoISet (default: weight all QoIs equally) to use when computing errors

Definition at line 74 of file adjoint_refinement_estimator.h.

References _qoi_set.

const QoISet& libMesh::AdjointRefinementEstimator::qoi_set ( ) const
inline

Access to the QoISet (default: weight all QoIs equally) to use when computing errors

Definition at line 80 of file adjoint_refinement_estimator.h.

References _qoi_set, estimate_error(), and libmesh_nullptr.

void libMesh::ErrorEstimator::reduce_error ( std::vector< ErrorVectorReal > &  error_per_cell,
const Parallel::Communicator &comm  LIBMESH_CAN_DEFAULT_TO_COMMWORLD 
) const
protectedinherited

This method takes the local error contributions in error_per_cell from each processor and combines them to get the global error vector.

Definition at line 33 of file error_estimator.C.

References libMesh::Parallel::Communicator::sum().

Referenced by libMesh::UniformRefinementEstimator::_estimate_error(), libMesh::WeightedPatchRecoveryErrorEstimator::estimate_error(), libMesh::PatchRecoveryErrorEstimator::estimate_error(), libMesh::JumpErrorEstimator::estimate_error(), estimate_error(), and libMesh::ExactErrorEstimator::estimate_error().

35 {
36  // This function must be run on all processors at once
37  // parallel_object_only();
38 
39  // Each processor has now computed the error contribuions
40  // for its local elements. We may need to sum the vector to
41  // recover the error for each element.
42 
43  comm.sum(error_per_cell);
44 }
virtual ErrorEstimatorType libMesh::AdjointRefinementEstimator::type ( ) const
inlinevirtual

Returns the type for the ErrorEstimator subclass.

Implements libMesh::ErrorEstimator.

Definition at line 112 of file adjoint_refinement_estimator.h.

References libMesh::ADJOINT_REFINEMENT.

Member Data Documentation

QoISet libMesh::AdjointRefinementEstimator::_qoi_set
protected

A QoISet to handle cases with multiple QoIs available

Definition at line 133 of file adjoint_refinement_estimator.h.

Referenced by estimate_error(), and qoi_set().

std::vector<Number> libMesh::AdjointRefinementEstimator::computed_global_QoI_errors
protected
SystemNorm libMesh::ErrorEstimator::error_norm
inherited

When estimating the error in a single system, the error_norm is used to control the scaling and norm choice for each variable. Not all estimators will support all norm choices. The default scaling is for all variables to be weighted equally. The default norm choice depends on the error estimator.

Part of this functionality was supported via component_scale and sobolev_order in older libMesh versions, and a small part was supported via component_mask in even older versions. Hopefully the encapsulation here will allow us to avoid changing this API again.

Definition at line 149 of file error_estimator.h.

Referenced by libMesh::UniformRefinementEstimator::_estimate_error(), AdjointRefinementEstimator(), libMesh::DiscontinuityMeasure::boundary_side_integration(), libMesh::KellyErrorEstimator::boundary_side_integration(), libMesh::DiscontinuityMeasure::DiscontinuityMeasure(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::AdjointResidualErrorEstimator::estimate_error(), libMesh::ExactErrorEstimator::estimate_error(), libMesh::ErrorEstimator::estimate_errors(), libMesh::ExactErrorEstimator::ExactErrorEstimator(), libMesh::ExactErrorEstimator::find_squared_element_error(), libMesh::LaplacianErrorEstimator::init_context(), libMesh::DiscontinuityMeasure::init_context(), libMesh::KellyErrorEstimator::init_context(), libMesh::LaplacianErrorEstimator::internal_side_integration(), libMesh::DiscontinuityMeasure::internal_side_integration(), libMesh::KellyErrorEstimator::internal_side_integration(), libMesh::KellyErrorEstimator::KellyErrorEstimator(), libMesh::LaplacianErrorEstimator::LaplacianErrorEstimator(), libMesh::PatchRecoveryErrorEstimator::PatchRecoveryErrorEstimator(), libMesh::JumpErrorEstimator::reinit_sides(), and libMesh::UniformRefinementEstimator::UniformRefinementEstimator().

unsigned char libMesh::AdjointRefinementEstimator::number_h_refinements

How many h refinements to perform to get the fine grid

Definition at line 118 of file adjoint_refinement_estimator.h.

Referenced by estimate_error().

unsigned char libMesh::AdjointRefinementEstimator::number_p_refinements

How many p refinements to perform to get the fine grid

Definition at line 123 of file adjoint_refinement_estimator.h.

Referenced by estimate_error().


The documentation for this class was generated from the following files: