libMesh
Public Types | Public Member Functions | Public Attributes | Protected Member Functions | Protected Attributes | List of all members
libMesh::AdjointRefinementEstimator Class Reference

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. More...

#include <adjoint_refinement_estimator.h>

Inheritance diagram for libMesh::AdjointRefinementEstimator:
[legend]

Public Types

typedef std::map< std::pair< const System *, unsigned int >, ErrorVector * > ErrorMap
 When calculating many error vectors at once, we need a data structure to hold them all. More...
 

Public Member Functions

 AdjointRefinementEstimator ()
 Constructor. More...
 
 ~AdjointRefinementEstimator ()
 Destructor. More...
 
QoISetqoi_set ()
 Access to the QoISet (default: weight all QoIs equally) to use when computing errors. More...
 
const QoISetqoi_set () const
 Access to the QoISet (default: weight all QoIs equally) to use when computing errors. More...
 
virtual void estimate_error (const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=libmesh_nullptr, bool estimate_parent_error=false)
 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. More...
 
Numberget_global_QoI_error_estimate (unsigned int qoi_index)
 This is an accessor function to access the computed global QoI error estimates. More...
 
virtual ErrorEstimatorType type () const
 
DifferentiablePhysicsget_residual_evaluation_physics ()
 
void set_residual_evaluation_physics (DifferentiablePhysics *set_physics)
 Set the _residual_evaluation_physics member to argument. More...
 
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)
 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. More...
 
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)
 This virtual function can be redefined in derived classes, but by default it calls estimate_error repeatedly to calculate the requested error vectors. More...
 

Public Attributes

unsigned char number_h_refinements
 How many h refinements to perform to get the fine grid. More...
 
unsigned char number_p_refinements
 How many p refinements to perform to get the fine grid. More...
 
SystemNorm error_norm
 When estimating the error in a single system, the error_norm is used to control the scaling and norm choice for each variable. More...
 

Protected Member Functions

void reduce_error (std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm LIBMESH_CAN_DEFAULT_TO_COMMWORLD) const
 This method takes the local error contributions in error_per_cell from each processor and combines them to get the global error vector. More...
 

Protected Attributes

DifferentiablePhysics_residual_evaluation_physics
 Pointer to object to use for physics assembly evaluations. More...
 
std::vector< Numbercomputed_global_QoI_errors
 
QoISet _qoi_set
 A QoISet to handle cases with multiple QoIs available. More...
 

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 ( )

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  }
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
ErrorEstimator()
Constructor.
const class libmesh_nullptr_t libmesh_nullptr
QoISet _qoi_set
A QoISet to handle cases with multiple QoIs available.
DifferentiablePhysics * _residual_evaluation_physics
Pointer to object to use for physics assembly evaluations.
unsigned char number_h_refinements
How many h refinements to perform to get the fine grid.
unsigned char number_p_refinements
How many p refinements to perform to get the fine grid.
libMesh::AdjointRefinementEstimator::~AdjointRefinementEstimator ( )

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::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(), 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  std::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  std::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  for (const auto & elem : mesh.active_local_element_ptr_range())
375  for (unsigned int n=0; n != elem->n_nodes(); ++n)
376  {
377  // Get a reference to the current node
378  const Node & node = elem->node_ref(n);
379 
380  // Get the id of this node
381  dof_id_type node_id = node.id();
382 
383  // If we havent already processed this node, do so now
384  if (processed_node_ids.find(node_id) == processed_node_ids.end())
385  {
386  // Declare a neighbor_set to be filled by the find_point_neighbors
387  std::set<const Elem *> fine_grid_neighbor_set;
388 
389  // Call find_point_neighbors to fill the neighbor_set
390  elem->find_point_neighbors(node, fine_grid_neighbor_set);
391 
392  // A vector to hold the coarse grid parents neighbors
393  std::vector<dof_id_type> coarse_grid_neighbors;
394 
395  // Iterators over the fine grid neighbors set
396  std::set<const Elem *>::iterator fine_neighbor_it = fine_grid_neighbor_set.begin();
397  const std::set<const Elem *>::iterator fine_neighbor_end = fine_grid_neighbor_set.end();
398 
399  // Loop over all the fine neighbors of this node
400  for (; fine_neighbor_it != fine_neighbor_end ; ++fine_neighbor_it)
401  {
402  // Pointer to the current fine neighbor element
403  const Elem * fine_elem = *fine_neighbor_it;
404 
405  // Find the element id for the corresponding coarse grid element
406  const Elem * coarse_elem = fine_elem;
407  for (unsigned int j = 0; j != number_h_refinements; ++j)
408  {
409  libmesh_assert (coarse_elem->parent());
410 
411  coarse_elem = coarse_elem->parent();
412  }
413 
414  // Loop over the existing coarse neighbors and check if this one is
415  // already in there
416  const dof_id_type coarse_id = coarse_elem->id();
417  std::size_t j = 0;
418 
419  // If the set already contains this element break out of the loop
420  for (; j != coarse_grid_neighbors.size(); j++)
421  if (coarse_grid_neighbors[j] == coarse_id)
422  break;
423 
424  // If we didn't leave the loop even at the last element,
425  // this is a new neighbour, put in the coarse_grid_neighbor_set
426  if (j == coarse_grid_neighbors.size())
427  coarse_grid_neighbors.push_back(coarse_id);
428  } // End loop over fine neighbors
429 
430  // Set the shared_neighbour index for this node to the
431  // size of the coarse grid neighbor set
432  shared_element_count[node_id] =
433  cast_int<unsigned int>(coarse_grid_neighbors.size());
434 
435  // Add this node to processed_node_ids vector
436  processed_node_ids.insert(node_id);
437  } // End if not processed node
438  } // End loop over nodes
439 
440  // Get a DoF map, will be used to get the nodal dof_indices for each element
441  DofMap & dof_map = system.get_dof_map();
442 
443  // The global DOF indices, we will use these later on when we compute the element wise indicators
444  std::vector<dof_id_type> dof_indices;
445 
446  // Localize the global rhs and adjoint solution vectors (which might be shared on multiple processors) onto a
447  // local ghosted vector, this ensures each processor has all the dof_indices to compute an error indicator for
448  // an element it owns
449  UniquePtr<NumericVector<Number>> localized_projected_residual = NumericVector<Number>::build(system.comm());
450  localized_projected_residual->init(system.n_dofs(), system.n_local_dofs(), system.get_dof_map().get_send_list(), false, GHOSTED);
451  projected_residual.localize(*localized_projected_residual, system.get_dof_map().get_send_list());
452 
453  // Each adjoint solution will also require ghosting; for efficiency we'll reuse the same memory
454  UniquePtr<NumericVector<Number>> localized_adjoint_solution = NumericVector<Number>::build(system.comm());
455  localized_adjoint_solution->init(system.n_dofs(), system.n_local_dofs(), system.get_dof_map().get_send_list(), false, GHOSTED);
456 
457  // We will loop over each adjoint solution, localize that adjoint
458  // solution and then loop over local elements
459  for (std::size_t i=0; i != system.qoi.size(); i++)
460  {
461  // Skip this QoI if not in the QoI Set
462  if (_qoi_set.has_index(i))
463  {
464  // Get the weight for the current QoI
465  Real error_weight = _qoi_set.weight(i);
466 
467  (system.get_adjoint_solution(i)).localize(*localized_adjoint_solution, system.get_dof_map().get_send_list());
468 
469  // Loop over elements
470  for (const auto & elem : mesh.active_local_element_ptr_range())
471  {
472  // Go up number_h_refinements levels up to find the coarse parent
473  const Elem * coarse = elem;
474 
475  for (unsigned int j = 0; j != number_h_refinements; ++j)
476  {
477  libmesh_assert (coarse->parent());
478 
479  coarse = coarse->parent();
480  }
481 
482  const dof_id_type e_id = coarse->id();
483 
484  // Get the local to global degree of freedom maps for this element
485  dof_map.dof_indices (elem, dof_indices);
486 
487  // We will have to manually do the dot products.
488  Number local_contribution = 0.;
489 
490  for (std::size_t j=0; j != dof_indices.size(); j++)
491  {
492  // The contribution to the error indicator for this element from the current QoI
493  local_contribution += (*localized_projected_residual)(dof_indices[j]) * (*localized_adjoint_solution)(dof_indices[j]);
494  }
495 
496  // Multiply by the error weight for this QoI
497  local_contribution *= error_weight;
498 
499  // FIXME: we're throwing away information in the
500  // --enable-complex case
501  error_per_cell[e_id] += static_cast<ErrorVectorReal>
502  (std::abs(local_contribution));
503 
504  } // End loop over elements
505  } // End if belong to QoI set
506  } // End loop over QoIs
507 
508  for (std::size_t j=0; j != system.qoi.size(); j++)
509  if (_qoi_set.has_index(j))
510  delete coarse_adjoints[j];
511 
512  // Don't bother projecting the solution; we'll restore from backup
513  // after coarsening
514  system.project_solution_on_reinit() = false;
515 
516  // Uniformly coarsen the mesh, without projecting the solution
518 
519  for (unsigned int i = 0; i != number_h_refinements; ++i)
520  {
521  mesh_refinement.uniformly_coarsen(1);
522  // FIXME - should the reinits here be necessary? - RHS
523  es.reinit();
524  }
525 
526  for (unsigned int i = 0; i != number_p_refinements; ++i)
527  {
528  mesh_refinement.uniformly_p_coarsen(1);
529  es.reinit();
530  }
531 
532  // We should have the same number of active elements as when we started,
533  // but not necessarily the same number of elements since reinit() doesn't
534  // always call contract()
535  libmesh_assert_equal_to (n_coarse_elem, mesh.n_active_elem());
536 
537  // Restore old solutions and clean up the heap
538  system.project_solution_on_reinit() = old_projection_setting;
539 
540  // Restore the coarse solution vectors and delete their copies
541  *system.solution = *coarse_solution;
542  delete coarse_solution;
543  *system.current_local_solution = *coarse_local_solution;
544  delete coarse_local_solution;
545 
546  for (System::vectors_iterator vec = system.vectors_begin(); vec !=
547  system.vectors_end(); ++vec)
548  {
549  // The (string) name of this vector
550  const std::string & var_name = vec->first;
551 
552  // If it's a vector we already had (and not a newly created
553  // vector like an adjoint rhs), we need to restore it.
554  std::map<std::string, NumericVector<Number> *>::iterator it =
555  coarse_vectors.find(var_name);
556  if (it != coarse_vectors.end())
557  {
558  NumericVector<Number> * coarsevec = it->second;
559  system.get_vector(var_name) = *coarsevec;
560 
561  coarsevec->clear();
562  delete coarsevec;
563  }
564  }
565 
566  // Restore old partitioner and renumbering settings
567  mesh.partitioner().reset(old_partitioner.release());
568  mesh.allow_renumbering(old_renumbering_setting);
569 
570  // Finally sum the vector of estimated error values.
571  this->reduce_error(error_per_cell, system.comm());
572 
573  // We don't take a square root here; this is a goal-oriented
574  // estimate not a Hilbert norm estimate.
575 } // end estimate_error function
double abs(double a)
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
QoISet _qoi_set
A QoISet to handle cases with multiple QoIs available.
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())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
Real weight(unsigned int) const
Get the weight for this index (default 1.0)
Definition: qoi_set.h:240
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
DifferentiablePhysics * _residual_evaluation_physics
Pointer to object to use for physics assembly evaluations.
void reduce_error(std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm LIBMESH_CAN_DEFAULT_TO_COMMWORLD) const
This method takes the local error contributions in error_per_cell from each processor and combines th...
std::map< std::string, NumericVector< Number > * >::iterator vectors_iterator
Vector iterator typedefs.
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
Access multiple components at once.
unsigned char number_h_refinements
How many h refinements to perform to get the fine grid.
bool has_index(unsigned int) const
Return whether or not this index is in the set to be calculated.
Definition: qoi_set.h:221
uint8_t dof_id_type
Definition: id_types.h:64
unsigned char number_p_refinements
How many p refinements to perform to get the fine grid.
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 }
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
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
This pure virtual function must be redefined in derived classes to compute the error for each active ...
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.

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

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

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 }
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
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
This pure virtual function must be redefined in derived classes to compute the error for each active ...
Number& libMesh::AdjointRefinementEstimator::get_global_QoI_error_estimate ( unsigned int  qoi_index)

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 ( )
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; }
DifferentiablePhysics * _residual_evaluation_physics
Pointer to object to use for physics assembly evaluations.
QoISet& libMesh::AdjointRefinementEstimator::qoi_set ( )

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.

Referenced by build_adjoint_refinement_error_estimator().

78 { return _qoi_set; }
QoISet _qoi_set
A QoISet to handle cases with multiple QoIs available.
const QoISet& libMesh::AdjointRefinementEstimator::qoi_set ( ) const

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.

84 { return _qoi_set; }
QoISet _qoi_set
A QoISet to handle cases with multiple QoIs available.
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)

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; }
DifferentiablePhysics * _residual_evaluation_physics
Pointer to object to use for physics assembly evaluations.
virtual ErrorEstimatorType libMesh::AdjointRefinementEstimator::type ( ) const
virtual
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(), main(), 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 build_adjoint_refinement_error_estimator(), and 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: