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

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 }
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(), 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 }
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: