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
 
DifferentiablePhysicsget_residual_evaluation_physics ()
 
void set_residual_evaluation_physics (DifferentiablePhysics *set_physics)
 
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

DifferentiablePhysics_residual_evaluation_physics
 
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 50 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 113 of file error_estimator.h.

Constructor & Destructor Documentation

libMesh::AdjointRefinementEstimator::AdjointRefinementEstimator ( )
inline

Constructor. Sets the most common default parameter values.

Definition at line 57 of file adjoint_refinement_estimator.h.

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

57  :
62  _qoi_set(QoISet())
63  {
64  // We're not actually going to use error_norm; our norms are
65  // absolute values of QoI error.
67  }
const class libmesh_nullptr_t libmesh_nullptr
libMesh::AdjointRefinementEstimator::~AdjointRefinementEstimator ( )
inline

Destructor.

Definition at line 72 of file adjoint_refinement_estimator.h.

72 {}

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 71 of file adjoint_refinement_estimator.C.

References _qoi_set, _residual_evaluation_physics, 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::DifferentiableSystem::assembly(), libMesh::NumericVector< T >::build(), libMesh::NumericVector< T >::clear(), libMesh::DifferentiableSystem::clone(), 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_adjoint_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::TriangleWrapper::init(), 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::MeshBase::n_active_elem(), libMesh::System::n_dofs(), 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().

75 {
76  // We have to break the rules here, because we can't refine a const System
77  System & system = const_cast<System &>(_system);
78 
79  // An EquationSystems reference will be convenient.
80  EquationSystems & es = system.get_equation_systems();
81 
82  // The current mesh
83  MeshBase & mesh = es.get_mesh();
84 
85  // Get coarse grid adjoint solutions. This should be a relatively
86  // quick (especially with preconditioner reuse) way to get a good
87  // initial guess for the fine grid adjoint solutions. More
88  // importantly, subtracting off a coarse adjoint approximation gives
89  // us better local error indication, and subtracting off *some* lift
90  // function is necessary for correctness if we have heterogeneous
91  // adjoint Dirichlet conditions.
92 
93  // Solve the adjoint problem(s) on the coarse FE space
94  // Only if the user didn't already solve it for us
95  if (!system.is_adjoint_already_solved())
96  system.adjoint_solve(_qoi_set);
97 
98  // Loop over all the adjoint problems and, if any have heterogenous
99  // Dirichlet conditions, get the corresponding coarse lift
100  // function(s)
101  for (std::size_t j=0; j != system.qoi.size(); j++)
102  {
103  // Skip this QoI if it is not in the QoI Set or if there are no
104  // heterogeneous Dirichlet boundaries for it
105  if (_qoi_set.has_index(j) &&
106  system.get_dof_map().has_adjoint_dirichlet_boundaries(j))
107  {
108  // Next, we are going to build up the residual for evaluating the flux QoI
109  NumericVector<Number> * coarse_residual = libmesh_nullptr;
110 
111  // The definition of a flux QoI is R(u^h, L) where R is the residual as defined
112  // by a conservation law. Therefore, if we are using stabilization, the
113  // R should be specified by the user via the residual_evaluation_physics
114 
115  // If the residual physics pointer is not null, use it when
116  // evaluating here.
117  {
118  const bool swapping_physics = _residual_evaluation_physics;
119  if (swapping_physics)
120  dynamic_cast<DifferentiableSystem &>(system).swap_physics(_residual_evaluation_physics);
121 
122  // Assemble without applying constraints, to capture the solution values on the boundary
123  (dynamic_cast<ImplicitSystem &>(system)).assembly(true, false, false, true);
124 
125  // Get the residual vector (no constraints applied on boundary, so we wont blow away the lift)
126  coarse_residual = &(dynamic_cast<ExplicitSystem &>(system)).get_vector("RHS Vector");
127  coarse_residual->close();
128 
129  // Now build the lift function and add it to the system vectors
130  std::ostringstream liftfunc_name;
131  liftfunc_name << "adjoint_lift_function" << j;
132 
133  system.add_vector(liftfunc_name.str());
134 
135  // Initialize lift with coarse adjoint solve associate with this flux QoI to begin with
136  system.get_vector(liftfunc_name.str()).init(system.get_adjoint_solution(j), false);
137 
138  // Build the actual lift using adjoint dirichlet conditions
139  system.get_dof_map().enforce_adjoint_constraints_exactly
140  (system.get_vector(liftfunc_name.str()), static_cast<unsigned int>(j));
141 
142  // Compute the flux R(u^h, L)
143  std::cout<<"The flux QoI "<<static_cast<unsigned int>(j)<<" is: "<<coarse_residual->dot(system.get_vector(liftfunc_name.str()))<<std::endl<<std::endl;
144 
145  // Swap back if needed
146  if (swapping_physics)
147  dynamic_cast<DifferentiableSystem &>(system).swap_physics(_residual_evaluation_physics);
148  }
149  } // End if QoI in set and flux/dirichlet boundary QoI
150  } // End loop over QoIs
151 
152  // We'll want to back up all coarse grid vectors
153  std::map<std::string, NumericVector<Number> *> coarse_vectors;
154  for (System::vectors_iterator vec = system.vectors_begin(); vec !=
155  system.vectors_end(); ++vec)
156  {
157  // The (string) name of this vector
158  const std::string & var_name = vec->first;
159 
160  coarse_vectors[var_name] = vec->second->clone().release();
161  }
162 
163  // Back up the coarse solution and coarse local solution
164  NumericVector<Number> * coarse_solution =
165  system.solution->clone().release();
166  NumericVector<Number> * coarse_local_solution =
167  system.current_local_solution->clone().release();
168 
169  // And we'll need to temporarily change solution projection settings
170  bool old_projection_setting;
171  old_projection_setting = system.project_solution_on_reinit();
172 
173  // Make sure the solution is projected when we refine the mesh
174  system.project_solution_on_reinit() = true;
175 
176  // And it'll be best to avoid any repartitioning
177  UniquePtr<Partitioner> old_partitioner(mesh.partitioner().release());
178 
179  // And we can't allow any renumbering
180  const bool old_renumbering_setting = mesh.allow_renumbering();
181  mesh.allow_renumbering(false);
182 
183  // Use a non-standard solution vector if necessary
184  if (solution_vector && solution_vector != system.solution.get())
185  {
186  NumericVector<Number> * newsol =
187  const_cast<NumericVector<Number> *> (solution_vector);
188  newsol->swap(*system.solution);
189  system.update();
190  }
191 
192  // Resize the error_per_cell vector to be
193  // the number of elements, initialized to 0.
194  error_per_cell.clear();
195  error_per_cell.resize (mesh.max_elem_id(), 0.);
196 
197 #ifndef NDEBUG
198  // n_coarse_elem is only used in an assertion later so
199  // avoid declaring it unless asserts are active.
200  const dof_id_type n_coarse_elem = mesh.n_active_elem();
201 #endif
202 
203  // Uniformly refine the mesh
204  MeshRefinement mesh_refinement(mesh);
205 
207 
208  // FIXME: this may break if there is more than one System
209  // on this mesh but estimate_error was still called instead of
210  // estimate_errors
211  for (unsigned int i = 0; i != number_h_refinements; ++i)
212  {
213  mesh_refinement.uniformly_refine(1);
214  es.reinit();
215  }
216 
217  for (unsigned int i = 0; i != number_p_refinements; ++i)
218  {
219  mesh_refinement.uniformly_p_refine(1);
220  es.reinit();
221  }
222 
223  // Copy the projected coarse grid solutions, which will be
224  // overwritten by solve()
225  std::vector<NumericVector<Number> *> coarse_adjoints;
226  for (std::size_t j=0; j != system.qoi.size(); j++)
227  {
228  if (_qoi_set.has_index(j))
229  {
230  NumericVector<Number> * coarse_adjoint =
231  NumericVector<Number>::build(mesh.comm()).release();
232 
233  // Can do "fast" init since we're overwriting this in a sec
234  coarse_adjoint->init(system.get_adjoint_solution(j),
235  /* fast = */ true);
236 
237  *coarse_adjoint = system.get_adjoint_solution(j);
238 
239  coarse_adjoints.push_back(coarse_adjoint);
240  }
241  else
242  coarse_adjoints.push_back(static_cast<NumericVector<Number> *>(libmesh_nullptr));
243  }
244 
245  // Next, we are going to build up the residual for evaluating the
246  // error estimate.
247 
248  // If the residual physics pointer is not null, use it when
249  // evaluating here.
250  {
251  const bool swapping_physics = _residual_evaluation_physics;
252  if (swapping_physics)
253  dynamic_cast<DifferentiableSystem &>(system).swap_physics(_residual_evaluation_physics);
254 
255  // Rebuild the rhs with the projected primal solution, constraints
256  // have to be applied to get the correct error estimate since error
257  // on the Dirichlet boundary is zero
258  (dynamic_cast<ImplicitSystem &>(system)).assembly(true, false);
259 
260  // Swap back if needed
261  if (swapping_physics)
262  dynamic_cast<DifferentiableSystem &>(system).swap_physics(_residual_evaluation_physics);
263  }
264 
265  NumericVector<Number> & projected_residual = (dynamic_cast<ExplicitSystem &>(system)).get_vector("RHS Vector");
266  projected_residual.close();
267 
268  // Solve the adjoint problem(s) on the refined FE space
269  system.adjoint_solve(_qoi_set);
270 
271  // Now that we have the refined adjoint solution and the projected primal solution,
272  // we first compute the global QoI error estimate
273 
274  // Resize the computed_global_QoI_errors vector to hold the error estimates for each QoI
275  computed_global_QoI_errors.resize(system.qoi.size());
276 
277  // Loop over all the adjoint solutions and get the QoI error
278  // contributions from all of them. While we're looping anyway we'll
279  // pull off the coarse adjoints
280  for (std::size_t j=0; j != system.qoi.size(); j++)
281  {
282  // Skip this QoI if not in the QoI Set
283  if (_qoi_set.has_index(j))
284  {
285 
286  // If the adjoint solution has heterogeneous dirichlet
287  // values, then to get a proper error estimate here we need
288  // to subtract off a coarse grid lift function.
289  // |Q(u) - Q(u^h)| = |R([u^h]+, z^h+ - [L]+)| + HOT
290  if(system.get_dof_map().has_adjoint_dirichlet_boundaries(j))
291  {
292  // Need to create a string with current loop index to retrieve
293  // the correct vector from the liftvectors map
294  std::ostringstream liftfunc_name;
295  liftfunc_name << "adjoint_lift_function" << j;
296 
297  // Subtract off the corresponding lift vector from the adjoint solution
298  system.get_adjoint_solution(j) -= system.get_vector(liftfunc_name.str());
299 
300  // Now evaluate R(u^h, z^h+ - lift)
301  computed_global_QoI_errors[j] = projected_residual.dot(system.get_adjoint_solution(j));
302 
303  // Add the lift back to get back the adjoint
304  system.get_adjoint_solution(j) += system.get_vector(liftfunc_name.str());
305 
306  }
307  else
308  {
309  // Usual dual weighted residual error estimate
310  // |Q(u) - Q(u^h)| = |R([u^h]+, z^h+)| + HOT
311  computed_global_QoI_errors[j] = projected_residual.dot(system.get_adjoint_solution(j));
312  }
313 
314  }
315  }
316 
317  // Done with the global error estimates, now construct the element wise error indicators
318 
319  // To get a better element wise breakdown of the error estimate,
320  // we subtract off a coarse representation of the adjoint solution.
321  // |Q(u) - Q(u^h)| = |R([u^h]+, z^h+ - [z^h]+)|
322  // This remains valid for all combinations of heterogenous adjoint bcs and
323  // stabilized/non-stabilized formulations, except for the case where we not using a
324  // heterogenous adjoint bc and have a stabilized formulation.
325  // Then, R(u^h_s, z^h_s) != 0 (no Galerkin orthogonality w.r.t the non-stabilized residual)
326  for (std::size_t j=0; j != system.qoi.size(); j++)
327  {
328  // Skip this QoI if not in the QoI Set
329  if (_qoi_set.has_index(j))
330  {
331  // If we have a NULL residual evaluation physics pointer, we
332  // assume the user's formulation is consistent from mesh to
333  // mesh, so we have Galerkin orthogonality and we can get
334  // better indicator results by subtracting a coarse adjoint.
335 
336  // If we have a residual evaluation physics pointer, but we
337  // also have heterogeneous adjoint dirichlet boundaries,
338  // then we have to subtract off *some* lift function for
339  // consistency, and we choose the coarse adjoint in lieu of
340  // having any better ideas.
341 
342  // If we have a residual evaluation physics pointer and we
343  // have homogeneous adjoint dirichlet boundaries, then we
344  // don't have to subtract off anything, and with stabilized
345  // formulations we get the best results if we don't.
346  if(system.get_dof_map().has_adjoint_dirichlet_boundaries(j)
348  {
349  // z^h+ -> z^h+ - [z^h]+
350  system.get_adjoint_solution(j) -= *coarse_adjoints[j];
351  }
352  }
353  }
354 
355  // We ought to account for 'spill-over' effects while computing the
356  // element error indicators This happens because the same dof is
357  // shared by multiple elements, one way of mitigating this is to
358  // scale the contribution from each dof by the number of elements it
359  // belongs to We first obtain the number of elements each node
360  // belongs to
361 
362  // A map that relates a node id to an int that will tell us how many elements it is a node of
363  LIBMESH_BEST_UNORDERED_MAP<dof_id_type, unsigned int>shared_element_count;
364 
365  // To fill this map, we will loop over elements, and then in each element, we will loop
366  // over the nodes each element contains, and then query it for the number of coarse
367  // grid elements it was a node of
368 
369  // Keep track of which nodes we have already dealt with
370  LIBMESH_BEST_UNORDERED_SET<dof_id_type> processed_node_ids;
371 
372  // We will be iterating over all the active elements in the fine mesh that live on
373  // this processor
374  {
375  MeshBase::const_element_iterator elem_it = mesh.active_local_elements_begin();
376  const MeshBase::const_element_iterator elem_end = mesh.active_local_elements_end();
377 
378  // Start loop over elems
379  for (; elem_it != elem_end; ++elem_it)
380  {
381  // Pointer to this element
382  const Elem * elem = *elem_it;
383 
384  // Loop over the nodes in the element
385  for (unsigned int n=0; n != elem->n_nodes(); ++n)
386  {
387  // Get a reference to the current node
388  const Node & node = elem->node_ref(n);
389 
390  // Get the id of this node
391  dof_id_type node_id = node.id();
392 
393  // If we havent already processed this node, do so now
394  if (processed_node_ids.find(node_id) == processed_node_ids.end())
395  {
396  // Declare a neighbor_set to be filled by the find_point_neighbors
397  std::set<const Elem *> fine_grid_neighbor_set;
398 
399  // Call find_point_neighbors to fill the neighbor_set
400  elem->find_point_neighbors(node, fine_grid_neighbor_set);
401 
402  // A vector to hold the coarse grid parents neighbors
403  std::vector<dof_id_type> coarse_grid_neighbors;
404 
405  // Iterators over the fine grid neighbors set
406  std::set<const Elem *>::iterator fine_neighbor_it = fine_grid_neighbor_set.begin();
407  const std::set<const Elem *>::iterator fine_neighbor_end = fine_grid_neighbor_set.end();
408 
409  // Loop over all the fine neighbors of this node
410  for (; fine_neighbor_it != fine_neighbor_end ; ++fine_neighbor_it)
411  {
412  // Pointer to the current fine neighbor element
413  const Elem * fine_elem = *fine_neighbor_it;
414 
415  // Find the element id for the corresponding coarse grid element
416  const Elem * coarse_elem = fine_elem;
417  for (unsigned int j = 0; j != number_h_refinements; ++j)
418  {
419  libmesh_assert (coarse_elem->parent());
420 
421  coarse_elem = coarse_elem->parent();
422  }
423 
424  // Loop over the existing coarse neighbors and check if this one is
425  // already in there
426  const dof_id_type coarse_id = coarse_elem->id();
427  std::size_t j = 0;
428  for (; j != coarse_grid_neighbors.size(); j++)
429  {
430  // If the set already contains this element break out of the loop
431  if (coarse_grid_neighbors[j] == coarse_id)
432  {
433  break;
434  }
435  }
436 
437  // If we didn't leave the loop even at the last element,
438  // this is a new neighbour, put in the coarse_grid_neighbor_set
439  if (j == coarse_grid_neighbors.size())
440  {
441  coarse_grid_neighbors.push_back(coarse_id);
442  }
443 
444  } // End loop over fine neighbors
445 
446  // Set the shared_neighbour index for this node to the
447  // size of the coarse grid neighbor set
448  shared_element_count[node_id] =
449  cast_int<unsigned int>(coarse_grid_neighbors.size());
450 
451  // Add this node to processed_node_ids vector
452  processed_node_ids.insert(node_id);
453 
454  } // End if not processed node
455 
456  } // End loop over nodes
457 
458  } // End loop over elems
459  }
460 
461  // Get a DoF map, will be used to get the nodal dof_indices for each element
462  DofMap & dof_map = system.get_dof_map();
463 
464  // The global DOF indices, we will use these later on when we compute the element wise indicators
465  std::vector<dof_id_type> dof_indices;
466 
467  // Localize the global rhs and adjoint solution vectors (which might be shared on multiple processors) onto a
468  // local ghosted vector, this ensures each processor has all the dof_indices to compute an error indicator for
469  // an element it owns
470  UniquePtr<NumericVector<Number> > localized_projected_residual = NumericVector<Number>::build(system.comm());
471  localized_projected_residual->init(system.n_dofs(), system.n_local_dofs(), system.get_dof_map().get_send_list(), false, GHOSTED);
472  projected_residual.localize(*localized_projected_residual, system.get_dof_map().get_send_list());
473 
474  // Each adjoint solution will also require ghosting; for efficiency we'll reuse the same memory
475  UniquePtr<NumericVector<Number> > localized_adjoint_solution = NumericVector<Number>::build(system.comm());
476  localized_adjoint_solution->init(system.n_dofs(), system.n_local_dofs(), system.get_dof_map().get_send_list(), false, GHOSTED);
477 
478  // We will loop over each adjoint solution, localize that adjoint
479  // solution and then loop over local elements
480  for (std::size_t i=0; i != system.qoi.size(); i++)
481  {
482  // Skip this QoI if not in the QoI Set
483  if (_qoi_set.has_index(i))
484  {
485  // Get the weight for the current QoI
486  Real error_weight = _qoi_set.weight(i);
487 
488  (system.get_adjoint_solution(i)).localize(*localized_adjoint_solution, system.get_dof_map().get_send_list());
489 
490  // Loop over elements
491  MeshBase::const_element_iterator elem_it = mesh.active_local_elements_begin();
492  const MeshBase::const_element_iterator elem_end = mesh.active_local_elements_end();
493 
494  for (; elem_it != elem_end; ++elem_it)
495  {
496  // Pointer to the element
497  const Elem * elem = *elem_it;
498 
499  // Go up number_h_refinements levels up to find the coarse parent
500  const Elem * coarse = elem;
501 
502  for (unsigned int j = 0; j != number_h_refinements; ++j)
503  {
504  libmesh_assert (coarse->parent());
505 
506  coarse = coarse->parent();
507  }
508 
509  const dof_id_type e_id = coarse->id();
510 
511  // Get the local to global degree of freedom maps for this element
512  dof_map.dof_indices (elem, dof_indices);
513 
514  // We will have to manually do the dot products.
515  Number local_contribution = 0.;
516 
517  for (std::size_t j=0; j != dof_indices.size(); j++)
518  {
519  // The contribution to the error indicator for this element from the current QoI
520  local_contribution += (*localized_projected_residual)(dof_indices[j]) * (*localized_adjoint_solution)(dof_indices[j]);
521  }
522 
523  // Multiply by the error weight for this QoI
524  local_contribution *= error_weight;
525 
526  // FIXME: we're throwing away information in the
527  // --enable-complex case
528  error_per_cell[e_id] += static_cast<ErrorVectorReal>
529  (std::abs(local_contribution));
530 
531  } // End loop over elements
532 
533  } // End if belong to QoI set
534 
535  } // End loop over QoIs
536 
537  for (std::size_t j=0; j != system.qoi.size(); j++)
538  {
539  if (_qoi_set.has_index(j))
540  {
541  delete coarse_adjoints[j];
542  }
543  }
544 
545  // Don't bother projecting the solution; we'll restore from backup
546  // after coarsening
547  system.project_solution_on_reinit() = false;
548 
549  // Uniformly coarsen the mesh, without projecting the solution
551 
552  for (unsigned int i = 0; i != number_h_refinements; ++i)
553  {
554  mesh_refinement.uniformly_coarsen(1);
555  // FIXME - should the reinits here be necessary? - RHS
556  es.reinit();
557  }
558 
559  for (unsigned int i = 0; i != number_p_refinements; ++i)
560  {
561  mesh_refinement.uniformly_p_coarsen(1);
562  es.reinit();
563  }
564 
565  // We should have the same number of active elements as when we started,
566  // but not necessarily the same number of elements since reinit() doesn't
567  // always call contract()
568  libmesh_assert_equal_to (n_coarse_elem, mesh.n_active_elem());
569 
570  // Restore old solutions and clean up the heap
571  system.project_solution_on_reinit() = old_projection_setting;
572 
573  // Restore the coarse solution vectors and delete their copies
574  *system.solution = *coarse_solution;
575  delete coarse_solution;
576  *system.current_local_solution = *coarse_local_solution;
577  delete coarse_local_solution;
578 
579  for (System::vectors_iterator vec = system.vectors_begin(); vec !=
580  system.vectors_end(); ++vec)
581  {
582  // The (string) name of this vector
583  const std::string & var_name = vec->first;
584 
585  // If it's a vector we already had (and not a newly created
586  // vector like an adjoint rhs), we need to restore it.
587  std::map<std::string, NumericVector<Number> *>::iterator it =
588  coarse_vectors.find(var_name);
589  if (it != coarse_vectors.end())
590  {
591  NumericVector<Number> * coarsevec = it->second;
592  system.get_vector(var_name) = *coarsevec;
593 
594  coarsevec->clear();
595  delete coarsevec;
596  }
597  }
598 
599  // Restore old partitioner and renumbering settings
600  mesh.partitioner().reset(old_partitioner.release());
601  mesh.allow_renumbering(old_renumbering_setting);
602 
603  // Finally sum the vector of estimated error values.
604  this->reduce_error(error_per_cell, system.comm());
605 
606  // We don't take a square root here; this is a goal-oriented
607  // estimate not a Hilbert norm estimate.
608 } // 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 init(triangulateio &t)
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:748
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 111 of file adjoint_refinement_estimator.h.

References computed_global_QoI_errors.

112  {
113  return computed_global_QoI_errors[qoi_index];
114  }
DifferentiablePhysics* libMesh::AdjointRefinementEstimator::get_residual_evaluation_physics ( )
inline
Returns
A pointer to the DifferentiablePhysics object or NULL if no external Physics object is attached.

Definition at line 133 of file adjoint_refinement_estimator.h.

References _residual_evaluation_physics.

134  { return this->_residual_evaluation_physics; }
QoISet& libMesh::AdjointRefinementEstimator::qoi_set ( )
inline

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

Definition at line 78 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 84 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 contributions
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 }
void libMesh::AdjointRefinementEstimator::set_residual_evaluation_physics ( DifferentiablePhysics set_physics)
inline

Set the _residual_evaluation_physics member to argument

Definition at line 139 of file adjoint_refinement_estimator.h.

References _residual_evaluation_physics.

140  { this->_residual_evaluation_physics = set_physics; }
virtual ErrorEstimatorType libMesh::AdjointRefinementEstimator::type ( ) const
inlinevirtual
Returns
The type for the ErrorEstimator subclass.

Implements libMesh::ErrorEstimator.

Definition at line 116 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 156 of file adjoint_refinement_estimator.h.

Referenced by estimate_error(), and qoi_set().

DifferentiablePhysics* libMesh::AdjointRefinementEstimator::_residual_evaluation_physics
protected

Pointer to object to use for physics assembly evaluations. Defaults to libmesh_nullptr for backwards compatibility.

Definition at line 148 of file adjoint_refinement_estimator.h.

Referenced by estimate_error(), get_residual_evaluation_physics(), and set_residual_evaluation_physics().

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 150 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 122 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 127 of file adjoint_refinement_estimator.h.

Referenced by estimate_error().


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