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 >, std::unique_ptr< 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 (const AdjointRefinementEstimator &)=default
 Copy/move ctor, copy/move assignment operator, and destructor are all explicitly defaulted for this class. More...
 
 AdjointRefinementEstimator (AdjointRefinementEstimator &&)=default
 
AdjointRefinementEstimatoroperator= (const AdjointRefinementEstimator &)=default
 
AdjointRefinementEstimatoroperator= (AdjointRefinementEstimator &&)=default
 
virtual ~AdjointRefinementEstimator ()=default
 
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=nullptr, bool estimate_parent_error=false) override
 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 override
 
DifferentiablePhysicsget_residual_evaluation_physics ()
 
void set_residual_evaluation_physics (DifferentiablePhysics *set_physics)
 Set the _residual_evaluation_physics member to argument. More...
 
DifferentiablePhysicsget_adjoint_evaluation_physics ()
 
void set_adjoint_evaluation_physics (DifferentiablePhysics *set_physics)
 Set the _adjoint_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=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=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) 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...
 
DifferentiablePhysics_adjoint_evaluation_physics
 Pointer to object to use for adjoint assembly. 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

◆ ErrorMap

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

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

Definition at line 121 of file error_estimator.h.

Constructor & Destructor Documentation

◆ AdjointRefinementEstimator() [1/3]

libMesh::AdjointRefinementEstimator::AdjointRefinementEstimator ( )

Constructor.

Sets the most common default parameter values.

Definition at line 73 of file adjoint_refinement_estimator.C.

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

73  :
79  _qoi_set(QoISet())
80 {
81  // We're not actually going to use error_norm; our norms are
82  // absolute values of QoI error.
84 }
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
DifferentiablePhysics * _adjoint_evaluation_physics
Pointer to object to use for adjoint assembly.
QoISet _qoi_set
A QoISet to handle cases with multiple QoIs available.
ErrorEstimator()=default
Constructor.
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.

◆ AdjointRefinementEstimator() [2/3]

libMesh::AdjointRefinementEstimator::AdjointRefinementEstimator ( const AdjointRefinementEstimator )
default

Copy/move ctor, copy/move assignment operator, and destructor are all explicitly defaulted for this class.

◆ AdjointRefinementEstimator() [3/3]

libMesh::AdjointRefinementEstimator::AdjointRefinementEstimator ( AdjointRefinementEstimator &&  )
default

◆ ~AdjointRefinementEstimator()

virtual libMesh::AdjointRefinementEstimator::~AdjointRefinementEstimator ( )
virtualdefault

Member Function Documentation

◆ estimate_error()

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

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

References _adjoint_evaluation_physics, _qoi_set, _residual_evaluation_physics, std::abs(), libMesh::System::add_vector(), libMesh::ImplicitSystem::adjoint_solve(), libMesh::as_range(), libMesh::ImplicitSystem::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::System::get_project_with_constraints(), 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::libmesh_assert(), libMesh::NumericVector< T >::localize(), libMesh::make_range(), mesh, libMesh::System::n_dofs(), libMesh::System::n_local_dofs(), libMesh::System::n_qois(), libMesh::System::number(), number_h_refinements, number_p_refinements, libMesh::Elem::parent(), libMesh::System::project_solution_on_reinit(), libMesh::Real, libMesh::ErrorEstimator::reduce_error(), libMesh::EquationSystems::reinit(), libMesh::System::remove_vector(), libMesh::System::set_project_with_constraints(), libMesh::System::set_vector_preservation(), 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::vector_name(), libMesh::System::vector_preservation(), libMesh::System::vectors_begin(), libMesh::System::vectors_end(), and libMesh::QoISet::weight().

Referenced by libMesh::Euler2Solver::integrate_adjoint_refinement_error_estimate(), libMesh::EulerSolver::integrate_adjoint_refinement_error_estimate(), and libMesh::SteadySolver::integrate_adjoint_refinement_error_estimate().

95 {
96  // We have to break the rules here, because we can't refine a const System
97  System & mutable_system = const_cast<System &>(_system);
98 
99  // We really have to break the rules, because we can't do an
100  // adjoint_solve without a matrix.
101  ImplicitSystem & system = dynamic_cast<ImplicitSystem &>(mutable_system);
102 
103  // An EquationSystems reference will be convenient.
104  EquationSystems & es = system.get_equation_systems();
105 
106  // The current mesh
107  MeshBase & mesh = es.get_mesh();
108 
109  // Get coarse grid adjoint solutions. This should be a relatively
110  // quick (especially with preconditioner reuse) way to get a good
111  // initial guess for the fine grid adjoint solutions. More
112  // importantly, subtracting off a coarse adjoint approximation gives
113  // us better local error indication, and subtracting off *some* lift
114  // function is necessary for correctness if we have heterogeneous
115  // adjoint Dirichlet conditions.
116 
117  // Solve the adjoint problem(s) on the coarse FE space
118  // Only if the user didn't already solve it for us
119  // If _adjoint_evaluation_physics pointer is not null, swap
120  // the current physics with the one held by _adjoint_evaluation_physics
121  // before assembling the adjoint problem
122  if (!system.is_adjoint_already_solved())
123  {
124  // Swap in different physics if needed
126  dynamic_cast<DifferentiableSystem &>(system).push_physics(*_adjoint_evaluation_physics);
127 
128  // Solve the adjoint problem, remember physics swap also resets the cache, so
129  // we will assemble again, otherwise we just take the transpose
130  system.adjoint_solve(_qoi_set);
131 
133  dynamic_cast<DifferentiableSystem &>(system).pop_physics();
134  }
135 
136  // Loop over all the adjoint problems and, if any have heterogenous
137  // Dirichlet conditions, get the corresponding coarse lift
138  // function(s)
139  for (unsigned int j=0,
140  n_qois = system.n_qois(); j != n_qois; j++)
141  {
142  // Skip this QoI if it is not in the QoI Set or if there are no
143  // heterogeneous Dirichlet boundaries for it
144  if (_qoi_set.has_index(j) &&
145  system.get_dof_map().has_adjoint_dirichlet_boundaries(j))
146  {
147  // Next, we are going to build up the residual for evaluating the flux QoI
148  NumericVector<Number> * coarse_residual = nullptr;
149 
150  // The definition of a flux QoI is R(u^h, L) where R is the residual as defined
151  // by a conservation law. Therefore, if we are using stabilization, the
152  // R should be specified by the user via the residual_evaluation_physics
153 
154  // If the residual physics pointer is not null, use it when
155  // evaluating here.
156  {
158  dynamic_cast<DifferentiableSystem &>(system).push_physics(*_residual_evaluation_physics);
159 
160  // Assemble without applying constraints, to capture the solution values on the boundary
161  system.assembly(true, false, false, true);
162 
163  // Get the residual vector (no constraints applied on boundary, so we wont blow away the lift)
164  coarse_residual = &system.get_vector("RHS Vector");
165  coarse_residual->close();
166 
167  // Now build the lift function and add it to the system vectors
168  std::ostringstream liftfunc_name;
169  liftfunc_name << "adjoint_lift_function" << j;
170 
171  system.add_vector(liftfunc_name.str());
172 
173  // Initialize lift with coarse adjoint solve associate with this flux QoI to begin with
174  system.get_vector(liftfunc_name.str()).init(system.get_adjoint_solution(j), false);
175 
176  // Build the actual lift using adjoint dirichlet conditions
177  system.get_dof_map().enforce_adjoint_constraints_exactly
178  (system.get_vector(liftfunc_name.str()), static_cast<unsigned int>(j));
179 
180  // Compute the flux R(u^h, L)
181  std::cout<<"The flux QoI "<<static_cast<unsigned int>(j)<<" is: "<<coarse_residual->dot(system.get_vector(liftfunc_name.str()))<<std::endl<<std::endl;
182 
183  // Restore the original physics if needed
185  dynamic_cast<DifferentiableSystem &>(system).pop_physics();
186  }
187  } // End if QoI in set and flux/dirichlet boundary QoI
188  } // End loop over QoIs
189 
190  // We'll want to back up all coarse grid vectors
191  std::map<std::string, std::unique_ptr<NumericVector<Number>>> coarse_vectors;
192  for (const auto & [var_name, vec] : as_range(system.vectors_begin(), system.vectors_end()))
193  coarse_vectors[var_name] = vec->clone();
194 
195  // Back up the coarse solution and coarse local solution
196  std::unique_ptr<NumericVector<Number>> coarse_solution = system.solution->clone();
197  std::unique_ptr<NumericVector<Number>> coarse_local_solution = system.current_local_solution->clone();
198 
199  // We need to make sure that the coarse adjoint vectors used in the
200  // calculations below are preserved during reinit, regardless of how
201  // the user is treating them in their code
202  // The adjoint lift function we have defined above is set to be preserved
203  // by default
204  std::vector<bool> old_adjoints_projection_settings(system.n_qois());
205  for (auto j : make_range(system.n_qois()))
206  {
207  if (_qoi_set.has_index(j))
208  {
209  // Get the vector preservation setting for this adjoint vector
210  auto adjoint_vector_name = system.vector_name(system.get_adjoint_solution(j));
211  auto old_adjoint_vector_projection_setting = system.vector_preservation(adjoint_vector_name);
212 
213  // Save for restoration later on
214  old_adjoints_projection_settings[j] = old_adjoint_vector_projection_setting;
215 
216  // Set the preservation to true for the upcoming reinits
217  system.set_vector_preservation(adjoint_vector_name, true);
218  }
219  }
220 
221  // And we'll need to temporarily change solution projection settings
222  bool old_projection_setting;
223  old_projection_setting = system.project_solution_on_reinit();
224 
225  // Make sure the solution is projected when we refine the mesh
226  system.project_solution_on_reinit() = true;
227 
228  // And we need to make sure we dont reapply constraints after refining the mesh
229  bool old_project_with_constraints_setting;
230  old_project_with_constraints_setting = system.get_project_with_constraints();
231 
232  system.set_project_with_constraints(false);
233 
234  // And it'll be best to avoid any repartitioning
235  std::unique_ptr<Partitioner> old_partitioner(mesh.partitioner().release());
236 
237  // And we can't allow any renumbering
238  const bool old_renumbering_setting = mesh.allow_renumbering();
239  mesh.allow_renumbering(false);
240 
241  // Use a non-standard solution vector if necessary
242  if (solution_vector && solution_vector != system.solution.get())
243  {
244  NumericVector<Number> * newsol =
245  const_cast<NumericVector<Number> *> (solution_vector);
246  newsol->swap(*system.solution);
247  system.update();
248  }
249 
250  // Resize the error_per_cell vector to be
251  // the number of elements, initialized to 0.
252  error_per_cell.clear();
253  error_per_cell.resize (mesh.max_elem_id(), 0.);
254 
255 #ifndef NDEBUG
256  // These variables are only used in assertions later so
257  // avoid declaring them unless asserts are active.
258  const dof_id_type n_coarse_elem = mesh.n_active_elem();
259 
260  dof_id_type local_dof_bearing_nodes = 0;
261  const unsigned int sysnum = system.number();
262  for (const auto * node : mesh.local_node_ptr_range())
263  for (unsigned int v=0, nvars = node->n_vars(sysnum);
264  v != nvars; ++v)
265  if (node->n_comp(sysnum, v))
266  {
267  local_dof_bearing_nodes++;
268  break;
269  }
270 #endif // NDEBUG
271 
272  // Uniformly refine the mesh
273  MeshRefinement mesh_refinement(mesh);
274 
275  // We only need to worry about Galerkin orthogonality if we
276  // are estimating discretization error in a single model setting
277  {
278  const bool swapping_adjoint_physics = _adjoint_evaluation_physics;
279  if(!swapping_adjoint_physics)
281  }
282 
283  // FIXME: this may break if there is more than one System
284  // on this mesh but estimate_error was still called instead of
285  // estimate_errors
286  for (unsigned int i = 0; i != number_h_refinements; ++i)
287  {
288  mesh_refinement.uniformly_refine(1);
289  es.reinit();
290  }
291 
292  for (unsigned int i = 0; i != number_p_refinements; ++i)
293  {
294  mesh_refinement.uniformly_p_refine(1);
295  es.reinit();
296  }
297 
298  // Copy the projected coarse grid solutions, which will be
299  // overwritten by solve()
300  std::vector<std::unique_ptr<NumericVector<Number>>> coarse_adjoints;
301  for (auto j : make_range(system.n_qois()))
302  {
303  if (_qoi_set.has_index(j))
304  {
305  auto coarse_adjoint = NumericVector<Number>::build(mesh.comm());
306 
307  // Can do "fast" init since we're overwriting this in a sec
308  coarse_adjoint->init(system.get_adjoint_solution(j),
309  /* fast = */ true);
310 
311  *coarse_adjoint = system.get_adjoint_solution(j);
312 
313  coarse_adjoints.emplace_back(std::move(coarse_adjoint));
314  }
315  else
316  coarse_adjoints.emplace_back(nullptr);
317  }
318 
319  // Next, we are going to build up the residual for evaluating the
320  // error estimate.
321 
322  // If the residual physics pointer is not null, use it when
323  // evaluating here.
324  {
326  dynamic_cast<DifferentiableSystem &>(system).push_physics(*_residual_evaluation_physics);
327 
328  // Rebuild the rhs with the projected primal solution, do not apply constraints
329  system.assembly(true, false, false, true);
330 
331  // Restore the original physics if needed
333  dynamic_cast<DifferentiableSystem &>(system).pop_physics();
334  }
335 
336  NumericVector<Number> & projected_residual = (dynamic_cast<ExplicitSystem &>(system)).get_vector("RHS Vector");
337  projected_residual.close();
338 
340  dynamic_cast<DifferentiableSystem &>(system).push_physics(*_adjoint_evaluation_physics);
341 
342  // Solve the adjoint problem(s) on the refined FE space
343  // The matrix will be reassembled again because we have refined the mesh
344  // If we have no h or p refinements, no need to solve for a fine adjoint
346  system.adjoint_solve(_qoi_set);
347 
348  // Swap back if needed, recall that _adjoint_evaluation_physics now holds the pointer
349  // to the pre-swap physics
351  dynamic_cast<DifferentiableSystem &>(system).pop_physics();
352 
353  // Now that we have the refined adjoint solution and the projected primal solution,
354  // we first compute the global QoI error estimate
355 
356  // Resize the computed_global_QoI_errors vector to hold the error estimates for each QoI
357  computed_global_QoI_errors.resize(system.n_qois());
358 
359  // Loop over all the adjoint solutions and get the QoI error
360  // contributions from all of them. While we're looping anyway we'll
361  // pull off the coarse adjoints
362  for (auto j : make_range(system.n_qois()))
363  {
364  // Skip this QoI if not in the QoI Set
365  if (_qoi_set.has_index(j))
366  {
367 
368  // If the adjoint solution has heterogeneous dirichlet
369  // values, then to get a proper error estimate here we need
370  // to subtract off a coarse grid lift function.
371  // |Q(u) - Q(u^h)| = |R([u^h]+, z^h+ - [L]+)| + HOT
372  if(system.get_dof_map().has_adjoint_dirichlet_boundaries(j))
373  {
374  // Need to create a string with current loop index to retrieve
375  // the correct vector from the liftvectors map
376  std::ostringstream liftfunc_name;
377  liftfunc_name << "adjoint_lift_function" << j;
378 
379  // Subtract off the corresponding lift vector from the adjoint solution
380  system.get_adjoint_solution(j) -= system.get_vector(liftfunc_name.str());
381 
382  // Now evaluate R(u^h, z^h+ - lift)
383  computed_global_QoI_errors[j] = projected_residual.dot(system.get_adjoint_solution(j));
384 
385  // Add the lift back to get back the adjoint
386  system.get_adjoint_solution(j) += system.get_vector(liftfunc_name.str());
387 
388  }
389  else
390  {
391  // Usual dual weighted residual error estimate
392  // |Q(u) - Q(u^h)| = |R([u^h]+, z^h+)| + HOT
393  computed_global_QoI_errors[j] = projected_residual.dot(system.get_adjoint_solution(j));
394  }
395 
396  }
397  }
398 
399  // We are all done with Dirichlet lift vectors and they should be removed, lest we run into I/O issues later
400  for (auto j : make_range(system.n_qois()))
401  {
402  // Skip this QoI if not in the QoI Set
403  if (_qoi_set.has_index(j))
404  {
405  // Lifts are created only for adjoint dirichlet QoIs
406  if(system.get_dof_map().has_adjoint_dirichlet_boundaries(j))
407  {
408  // Need to create a string with current loop index to retrieve
409  // the correct vector from the liftvectors map
410  std::ostringstream liftfunc_name;
411  liftfunc_name << "adjoint_lift_function" << j;
412 
413  // Remove the lift vector from system since we did not write it to file and it cannot be retrieved
414  system.remove_vector(liftfunc_name.str());
415  }
416  }
417  }
418 
419  // Done with the global error estimates, now construct the element wise error indicators
420 
421  // To get a better element wise breakdown of the error estimate,
422  // we subtract off a coarse representation of the adjoint solution.
423  // |Q(u) - Q(u^h)| = |R([u^h]+, z^h+ - [z^h]+)|
424  // This remains valid for all combinations of heterogenous adjoint bcs and
425  // stabilized/non-stabilized formulations, except for the case where we not using a
426  // heterogenous adjoint bc and have a stabilized formulation.
427  // Then, R(u^h_s, z^h_s) != 0 (no Galerkin orthogonality w.r.t the non-stabilized residual)
428  for (auto j : make_range(system.n_qois()))
429  {
430  // Skip this QoI if not in the QoI Set
431  if (_qoi_set.has_index(j))
432  {
433  // If we have a nullptr residual evaluation physics pointer, we
434  // assume the user's formulation is consistent from mesh to
435  // mesh, so we have Galerkin orthogonality and we can get
436  // better indicator results by subtracting a coarse adjoint.
437 
438  // If we have a residual evaluation physics pointer, but we
439  // also have heterogeneous adjoint dirichlet boundaries,
440  // then we have to subtract off *some* lift function for
441  // consistency, and we choose the coarse adjoint in lieu of
442  // having any better ideas.
443 
444  // If we have a residual evaluation physics pointer and we
445  // have homogeneous adjoint dirichlet boundaries, then we
446  // don't have to subtract off anything, and with stabilized
447  // formulations we get the best results if we don't.
448  if(system.get_dof_map().has_adjoint_dirichlet_boundaries(j)
450  {
451  // z^h+ -> z^h+ - [z^h]+
452  system.get_adjoint_solution(j) -= *coarse_adjoints[j];
453  }
454  }
455  }
456 
457  // We ought to account for 'spill-over' effects while computing the
458  // element error indicators This happens because the same dof is
459  // shared by multiple elements, one way of mitigating this is to
460  // scale the contribution from each dof by the number of elements it
461  // belongs to We first obtain the number of elements each node
462  // belongs to
463 
464  // A map that relates a node id to an int that will tell us how many elements it is a node of
465  std::unordered_map<dof_id_type, unsigned int> shared_element_count;
466 
467  // To fill this map, we will loop over elements, and then in each element, we will loop
468  // over the nodes each element contains, and then query it for the number of coarse
469  // grid elements it was a node of
470 
471  // Keep track of which nodes we have already dealt with
472  std::unordered_set<dof_id_type> processed_node_ids;
473 
474  // We will be iterating over all the active elements in the fine mesh that live on
475  // this processor.
476  for (const auto & elem : mesh.active_local_element_ptr_range())
477  for (const Node & node : elem->node_ref_range())
478  {
479  // Get the id of this node
480  dof_id_type node_id = node.id();
481 
482  // If we haven't already processed this node, do so now
483  if (processed_node_ids.find(node_id) == processed_node_ids.end())
484  {
485  // Declare a neighbor_set to be filled by the find_point_neighbors
486  std::set<const Elem *> fine_grid_neighbor_set;
487 
488  // Call find_point_neighbors to fill the neighbor_set
489  elem->find_point_neighbors(node, fine_grid_neighbor_set);
490 
491  // A vector to hold the coarse grid parents neighbors
492  std::vector<dof_id_type> coarse_grid_neighbors;
493 
494  // Loop over all the fine neighbors of this node
495  for (const auto & fine_elem : fine_grid_neighbor_set)
496  {
497  // Find the element id for the corresponding coarse grid element
498  const Elem * coarse_elem = fine_elem;
499  for (unsigned int j = 0; j != number_h_refinements; ++j)
500  {
501  libmesh_assert (coarse_elem->parent());
502 
503  coarse_elem = coarse_elem->parent();
504  }
505 
506  // Loop over the existing coarse neighbors and check if this one is
507  // already in there
508  const dof_id_type coarse_id = coarse_elem->id();
509  std::size_t j = 0;
510 
511  // If the set already contains this element break out of the loop
512  for (std::size_t cgns = coarse_grid_neighbors.size(); j != cgns; j++)
513  if (coarse_grid_neighbors[j] == coarse_id)
514  break;
515 
516  // If we didn't leave the loop even at the last element,
517  // this is a new neighbour, put in the coarse_grid_neighbor_set
518  if (j == coarse_grid_neighbors.size())
519  coarse_grid_neighbors.push_back(coarse_id);
520  } // End loop over fine neighbors
521 
522  // Set the shared_neighbour index for this node to the
523  // size of the coarse grid neighbor set
524  shared_element_count[node_id] =
525  cast_int<unsigned int>(coarse_grid_neighbors.size());
526 
527  // Add this node to processed_node_ids vector
528  processed_node_ids.insert(node_id);
529  } // End if not processed node
530  } // End loop over nodes
531 
532  // Get a DoF map, will be used to get the nodal dof_indices for each element
533  DofMap & dof_map = system.get_dof_map();
534 
535  // The global DOF indices, we will use these later on when we compute the element wise indicators
536  std::vector<dof_id_type> dof_indices;
537 
538  // Localize the global rhs and adjoint solution vectors (which might be shared on multiple processors) onto a
539  // local ghosted vector, this ensures each processor has all the dof_indices to compute an error indicator for
540  // an element it owns
541  std::unique_ptr<NumericVector<Number>> localized_projected_residual = NumericVector<Number>::build(system.comm());
542  localized_projected_residual->init(system.n_dofs(), system.n_local_dofs(), system.get_dof_map().get_send_list(), false, GHOSTED);
543  projected_residual.localize(*localized_projected_residual, system.get_dof_map().get_send_list());
544 
545  // Each adjoint solution will also require ghosting; for efficiency we'll reuse the same memory
546  std::unique_ptr<NumericVector<Number>> localized_adjoint_solution = NumericVector<Number>::build(system.comm());
547  localized_adjoint_solution->init(system.n_dofs(), system.n_local_dofs(), system.get_dof_map().get_send_list(), false, GHOSTED);
548 
549  // We will loop over each adjoint solution, localize that adjoint
550  // solution and then loop over local elements
551  for (auto i : make_range(system.n_qois()))
552  {
553  // Skip this QoI if not in the QoI Set
554  if (_qoi_set.has_index(i))
555  {
556  // Get the weight for the current QoI
557  Real error_weight = _qoi_set.weight(i);
558 
559  (system.get_adjoint_solution(i)).localize(*localized_adjoint_solution, system.get_dof_map().get_send_list());
560 
561  // Loop over elements
562  for (const auto & elem : mesh.active_local_element_ptr_range())
563  {
564  // Go up number_h_refinements levels up to find the coarse parent
565  const Elem * coarse = elem;
566 
567  for (unsigned int j = 0; j != number_h_refinements; ++j)
568  {
569  libmesh_assert (coarse->parent());
570 
571  coarse = coarse->parent();
572  }
573 
574  const dof_id_type e_id = coarse->id();
575 
576  // Get the local to global degree of freedom maps for this element
577  dof_map.dof_indices (elem, dof_indices);
578 
579  // We will have to manually do the dot products.
580  Number local_contribution = 0.;
581 
582  // Sum the contribution to the error indicator for each element from the current QoI
583  for (const auto & dof : dof_indices)
584  local_contribution += (*localized_projected_residual)(dof) * (*localized_adjoint_solution)(dof);
585 
586  // Multiply by the error weight for this QoI
587  local_contribution *= error_weight;
588 
589  // FIXME: we're throwing away information in the
590  // --enable-complex case
591  error_per_cell[e_id] += static_cast<ErrorVectorReal>
592  (std::abs(local_contribution));
593 
594  } // End loop over elements
595  } // End if belong to QoI set
596  } // End loop over QoIs
597 
598  // Don't bother projecting the solution; we'll restore from backup
599  // after coarsening
600  system.project_solution_on_reinit() = false;
601 
602  // Uniformly coarsen the mesh, without projecting the solution
603  // Only need to do this if we are estimating discretization error
604  // with a single physics residual
607 
608  for (unsigned int i = 0; i != number_h_refinements; ++i)
609  {
610  mesh_refinement.uniformly_coarsen(1);
611  // FIXME - should the reinits here be necessary? - RHS
612  es.reinit();
613  }
614 
615  for (unsigned int i = 0; i != number_p_refinements; ++i)
616  {
617  mesh_refinement.uniformly_p_coarsen(1);
618  es.reinit();
619  }
620 
621  // We should have the same number of active elements as when we started,
622  // but not necessarily the same number of elements since reinit() doesn't
623  // always call contract()
624  libmesh_assert_equal_to (n_coarse_elem, mesh.n_active_elem());
625 
626  // We should have the same number of dof-bearing nodes as when we
627  // started
628 #ifndef NDEBUG
629  dof_id_type final_local_dof_bearing_nodes = 0;
630  for (const auto * node : mesh.local_node_ptr_range())
631  for (unsigned int v=0, nvars = node->n_vars(sysnum);
632  v != nvars; ++v)
633  if (node->n_comp(sysnum, v))
634  {
635  final_local_dof_bearing_nodes++;
636  break;
637  }
638  libmesh_assert_equal_to (local_dof_bearing_nodes,
639  final_local_dof_bearing_nodes);
640 #endif // NDEBUG
641 
642  // Restore old solutions and clean up the heap
643  system.project_solution_on_reinit() = old_projection_setting;
644  system.set_project_with_constraints(old_project_with_constraints_setting);
645 
646  // Restore the adjoint vector preservation settings
647  for (auto j : make_range(system.n_qois()))
648  {
649  if (_qoi_set.has_index(j))
650  {
651  auto adjoint_vector_name = system.vector_name(system.get_adjoint_solution(j));
652  system.set_vector_preservation(adjoint_vector_name, old_adjoints_projection_settings[j]);
653  }
654  }
655 
656  // Restore the coarse solution vectors and delete their copies
657  *system.solution = *coarse_solution;
658  *system.current_local_solution = *coarse_local_solution;
659 
660  for (const auto & pr : as_range(system.vectors_begin(), system.vectors_end()))
661  {
662  // The (string) name of this vector
663  const std::string & var_name = pr.first;
664 
665  // If it's a vector we already had (and not a newly created
666  // vector like an adjoint rhs), we need to restore it.
667  auto it = coarse_vectors.find(var_name);
668  if (it != coarse_vectors.end())
669  {
670  NumericVector<Number> * coarsevec = it->second.get();
671  system.get_vector(var_name) = *coarsevec;
672 
673  coarsevec->clear();
674  }
675  }
676 
677  // Restore old partitioner and renumbering settings
678  mesh.partitioner().reset(old_partitioner.release());
679  mesh.allow_renumbering(old_renumbering_setting);
680 
681  // Finally sum the vector of estimated error values.
682  this->reduce_error(error_per_cell, system.comm());
683 
684  // We don't take a square root here; this is a goal-oriented
685  // estimate not a Hilbert norm estimate.
686 } // end estimate_error function
DifferentiablePhysics * _adjoint_evaluation_physics
Pointer to object to use for adjoint assembly.
MeshBase & mesh
QoISet _qoi_set
A QoISet to handle cases with multiple QoIs available.
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
Definition: type_vector.h:57
bool has_index(std::size_t) const
Return whether or not this index is in the set to be calculated.
Definition: qoi_set.h:224
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
Helper function that allows us to treat a homogenous pair as a range.
Definition: simple_range.h:57
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
libmesh_assert(ctx)
void reduce_error(std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm) const
This method takes the local error contributions in error_per_cell from each processor and combines th...
DifferentiablePhysics * _residual_evaluation_physics
Pointer to object to use for physics assembly evaluations.
static std::unique_ptr< 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(std::size_t) const
Get the weight for this index (default 1.0)
Definition: qoi_set.h:243
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:134
template class LIBMESH_EXPORT NumericVector< Number >
unsigned char number_h_refinements
How many h refinements to perform to get the fine grid.
uint8_t dof_id_type
Definition: id_types.h:67
unsigned char number_p_refinements
How many p refinements to perform to get the fine grid.

◆ estimate_errors() [1/2]

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 = 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 47 of file error_estimator.C.

References libMesh::ErrorEstimator::error_norm, libMesh::ErrorEstimator::estimate_error(), libMesh::EquationSystems::get_system(), libMesh::index_range(), libMesh::make_range(), and libMesh::EquationSystems::n_systems().

52 {
53  SystemNorm old_error_norm = this->error_norm;
54 
55  // Sum the error values from each system
56  for (auto s : make_range(equation_systems.n_systems()))
57  {
58  ErrorVector system_error_per_cell;
59  const System & sys = equation_systems.get_system(s);
60  if (error_norms.find(&sys) == error_norms.end())
61  this->error_norm = old_error_norm;
62  else
63  this->error_norm = error_norms.find(&sys)->second;
64 
65  const NumericVector<Number> * solution_vector = nullptr;
66  if (solution_vectors &&
67  solution_vectors->find(&sys) != solution_vectors->end())
68  solution_vector = solution_vectors->find(&sys)->second;
69 
70  this->estimate_error(sys, system_error_per_cell,
71  solution_vector, estimate_parent_error);
72 
73  if (s)
74  {
75  libmesh_assert_equal_to (error_per_cell.size(), system_error_per_cell.size());
76  for (auto i : index_range(error_per_cell))
77  error_per_cell[i] += system_error_per_cell[i];
78  }
79  else
80  error_per_cell = system_error_per_cell;
81  }
82 
83  // Restore our old state before returning
84  this->error_norm = old_error_norm;
85 }
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=nullptr, bool estimate_parent_error=false)=0
This pure virtual function must be redefined in derived classes to compute the error for each active ...
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:134
template class LIBMESH_EXPORT NumericVector< Number >
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:111

◆ estimate_errors() [2/2]

void libMesh::ErrorEstimator::estimate_errors ( const EquationSystems equation_systems,
ErrorMap errors_per_cell,
const std::map< const System *, const NumericVector< Number > *> *  solution_vectors = 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 93 of file error_estimator.C.

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

97 {
98  SystemNorm old_error_norm = this->error_norm;
99 
100  // Find the requested error values from each system
101  for (auto s : make_range(equation_systems.n_systems()))
102  {
103  const System & sys = equation_systems.get_system(s);
104 
105  unsigned int n_vars = sys.n_vars();
106 
107  for (unsigned int v = 0; v != n_vars; ++v)
108  {
109  // Only fill in ErrorVectors the user asks for
110  if (!errors_per_cell.count(std::make_pair(&sys, v)))
111  continue;
112 
113  // Calculate error in only one variable
114  std::vector<Real> weights(n_vars, 0.0);
115  weights[v] = 1.0;
116  this->error_norm =
117  SystemNorm(std::vector<FEMNormType>(n_vars, old_error_norm.type(v)),
118  weights);
119 
120  const NumericVector<Number> * solution_vector = nullptr;
121  if (solution_vectors &&
122  solution_vectors->find(&sys) != solution_vectors->end())
123  solution_vector = solution_vectors->find(&sys)->second;
124 
125  this->estimate_error
126  (sys, *errors_per_cell[std::make_pair(&sys, v)],
127  solution_vector, estimate_parent_error);
128  }
129  }
130 
131  // Restore our old state before returning
132  this->error_norm = old_error_norm;
133 }
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=nullptr, bool estimate_parent_error=false)=0
This pure virtual function must be redefined in derived classes to compute the error for each active ...
unsigned int n_vars
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:134
template class LIBMESH_EXPORT NumericVector< Number >

◆ get_adjoint_evaluation_physics()

DifferentiablePhysics* libMesh::AdjointRefinementEstimator::get_adjoint_evaluation_physics ( )
inline
Returns
A pointer to the DifferentiablePhysics object or nullptr if no external Physics object is attached.

Definition at line 140 of file adjoint_refinement_estimator.h.

References _adjoint_evaluation_physics.

141  { return this->_adjoint_evaluation_physics; }
DifferentiablePhysics * _adjoint_evaluation_physics
Pointer to object to use for adjoint assembly.

◆ get_global_QoI_error_estimate()

Number& libMesh::AdjointRefinementEstimator::get_global_QoI_error_estimate ( unsigned int  qoi_index)
inline

◆ get_residual_evaluation_physics()

DifferentiablePhysics* libMesh::AdjointRefinementEstimator::get_residual_evaluation_physics ( )
inline
Returns
A pointer to the DifferentiablePhysics object or nullptr if no external Physics object is attached.

Definition at line 127 of file adjoint_refinement_estimator.h.

References _residual_evaluation_physics.

128  { return this->_residual_evaluation_physics; }
DifferentiablePhysics * _residual_evaluation_physics
Pointer to object to use for physics assembly evaluations.

◆ operator=() [1/2]

AdjointRefinementEstimator& libMesh::AdjointRefinementEstimator::operator= ( const AdjointRefinementEstimator )
default

◆ operator=() [2/2]

AdjointRefinementEstimator& libMesh::AdjointRefinementEstimator::operator= ( AdjointRefinementEstimator &&  )
default

◆ qoi_set() [1/2]

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

◆ qoi_set() [2/2]

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 79 of file adjoint_refinement_estimator.h.

References _qoi_set.

79 { return _qoi_set; }
QoISet _qoi_set
A QoISet to handle cases with multiple QoIs available.

◆ reduce_error()

void libMesh::ErrorEstimator::reduce_error ( std::vector< ErrorVectorReal > &  error_per_cell,
const Parallel::Communicator comm 
) 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 32 of file error_estimator.C.

References TIMPI::Communicator::sum().

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

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

◆ set_adjoint_evaluation_physics()

void libMesh::AdjointRefinementEstimator::set_adjoint_evaluation_physics ( DifferentiablePhysics set_physics)
inline

Set the _adjoint_evaluation_physics member to argument.

Definition at line 146 of file adjoint_refinement_estimator.h.

References _adjoint_evaluation_physics.

147  { this->_adjoint_evaluation_physics = set_physics; }
DifferentiablePhysics * _adjoint_evaluation_physics
Pointer to object to use for adjoint assembly.

◆ set_residual_evaluation_physics()

void libMesh::AdjointRefinementEstimator::set_residual_evaluation_physics ( DifferentiablePhysics set_physics)
inline

Set the _residual_evaluation_physics member to argument.

Definition at line 133 of file adjoint_refinement_estimator.h.

References _residual_evaluation_physics.

134  { this->_residual_evaluation_physics = set_physics; }
DifferentiablePhysics * _residual_evaluation_physics
Pointer to object to use for physics assembly evaluations.

◆ type()

ErrorEstimatorType libMesh::AdjointRefinementEstimator::type ( ) const
overridevirtual
Returns
The type for the ErrorEstimator subclass.

Implements libMesh::ErrorEstimator.

Definition at line 86 of file adjoint_refinement_estimator.C.

References libMesh::ADJOINT_REFINEMENT.

Member Data Documentation

◆ _adjoint_evaluation_physics

DifferentiablePhysics* libMesh::AdjointRefinementEstimator::_adjoint_evaluation_physics
protected

Pointer to object to use for adjoint assembly.

Defaults to nullptr for backwards compatibility.

Definition at line 161 of file adjoint_refinement_estimator.h.

Referenced by estimate_error(), get_adjoint_evaluation_physics(), and set_adjoint_evaluation_physics().

◆ _qoi_set

QoISet libMesh::AdjointRefinementEstimator::_qoi_set
protected

A QoISet to handle cases with multiple QoIs available.

Definition at line 169 of file adjoint_refinement_estimator.h.

Referenced by estimate_error(), and qoi_set().

◆ _residual_evaluation_physics

DifferentiablePhysics* libMesh::AdjointRefinementEstimator::_residual_evaluation_physics
protected

Pointer to object to use for physics assembly evaluations.

Defaults to nullptr for backwards compatibility.

Definition at line 155 of file adjoint_refinement_estimator.h.

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

◆ computed_global_QoI_errors

std::vector<Number> libMesh::AdjointRefinementEstimator::computed_global_QoI_errors
protected

◆ error_norm

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 158 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::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::PatchRecoveryErrorEstimator::PatchRecoveryErrorEstimator(), libMesh::JumpErrorEstimator::reinit_sides(), and libMesh::UniformRefinementEstimator::UniformRefinementEstimator().

◆ number_h_refinements

unsigned char libMesh::AdjointRefinementEstimator::number_h_refinements

How many h refinements to perform to get the fine grid.

Definition at line 116 of file adjoint_refinement_estimator.h.

Referenced by estimate_error().

◆ number_p_refinements

unsigned char libMesh::AdjointRefinementEstimator::number_p_refinements

How many p refinements to perform to get the fine grid.

Definition at line 121 of file adjoint_refinement_estimator.h.

Referenced by estimate_error().


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