libMesh
adjoint_refinement_estimator.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 // C++ includes
19 #include <algorithm> // for std::fill
20 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
21 #include <cmath> // for sqrt
22 #include <set>
23 #include <sstream> // for ostringstream
24 #include <unordered_map>
25 #include <unordered_set>
26 
27 // Local Includes
28 #include "libmesh/dof_map.h"
29 #include "libmesh/elem.h"
30 #include "libmesh/equation_systems.h"
31 #include "libmesh/error_vector.h"
32 #include "libmesh/fe.h"
33 #include "libmesh/fe_interface.h"
34 #include "libmesh/libmesh_common.h"
35 #include "libmesh/libmesh_logging.h"
36 #include "libmesh/mesh_base.h"
37 #include "libmesh/mesh_refinement.h"
38 #include "libmesh/numeric_vector.h"
39 #include "libmesh/quadrature.h"
40 #include "libmesh/system.h"
41 #include "libmesh/diff_physics.h"
42 #include "libmesh/fem_system.h"
43 #include "libmesh/implicit_system.h"
44 #include "libmesh/partitioner.h"
45 #include "libmesh/adjoint_refinement_estimator.h"
46 #include "libmesh/enum_error_estimator_type.h"
47 #include "libmesh/enum_norm_type.h"
48 #include "libmesh/utility.h"
49 
50 #ifdef LIBMESH_ENABLE_AMR
51 
52 namespace libMesh
53 {
54 
55 //-----------------------------------------------------------------
56 // AdjointRefinementErrorEstimator
57 
58 // As of 10/2/2012, this function implements a 'brute-force' adjoint based QoI
59 // error estimator, using the following relationship:
60 // Q(u) - Q(u_h) \approx - R( (u_h)_(h/2), z_(h/2) ) .
61 // where: Q(u) is the true QoI
62 // u_h is the approximate primal solution on the current FE space
63 // Q(u_h) is the approximate QoI
64 // z_(h/2) is the adjoint corresponding to Q, on a richer FE space
65 // (u_h)_(h/2) is a projection of the primal solution on the richer FE space
66 // By richer FE space, we mean a grid that has been refined once and a polynomial order
67 // that has been increased once, i.e. one h and one p refinement
68 
69 // Both a global QoI error estimate and element wise error indicators are included
70 // Note that the element wise error indicators slightly over estimate the error in
71 // each element
72 
75  number_h_refinements(1),
76  number_p_refinements(0),
77  _residual_evaluation_physics(nullptr),
78  _adjoint_evaluation_physics(nullptr),
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 }
85 
87 {
88  return ADJOINT_REFINEMENT;
89 }
90 
92  ErrorVector & error_per_cell,
93  const NumericVector<Number> * solution_vector,
94  bool /*estimate_parent_error*/)
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) &&
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
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
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
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.
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  if (auto it = coarse_vectors.find(var_name);
668  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
687 
688 }// namespace libMesh
689 
690 #endif // #ifdef LIBMESH_ENABLE_AMR
virtual ErrorEstimatorType type() const override
vectors_iterator vectors_end()
End of vectors container.
Definition: system.h:2576
void uniformly_p_refine(unsigned int n=1)
Uniformly p refines the mesh n times.
This is the EquationSystems class.
const Elem * parent() const
Definition: elem.h:3030
A Node is like a Point, but with more information.
Definition: node.h:52
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
virtual void get(const std::vector< numeric_index_type > &index, T *values) const
Access multiple components at once.
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:2164
Data structure for specifying which Quantities of Interest should be calculated in an adjoint or a pa...
Definition: qoi_set.h:45
DifferentiablePhysics * _adjoint_evaluation_physics
Pointer to object to use for adjoint assembly.
unsigned int n_qois() const
Number of currently active quantities of interest.
Definition: system.h:2621
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
Definition: error_vector.h:50
ErrorEstimatorType
Defines an enum for the different types of error estimators which are available.
const EquationSystems & get_equation_systems() const
Definition: system.h:721
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
void uniformly_p_coarsen(unsigned int n=1)
Attempts to uniformly p coarsen the mesh n times.
const Parallel::Communicator & comm() const
NumericVector< Number > & add_vector(std::string_view vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Adds the additional vector vec_name to this system.
Definition: system.C:768
The libMesh namespace provides an interface to certain functionality in the library.
vectors_iterator vectors_begin()
Beginning of vectors container.
Definition: system.h:2564
QoISet _qoi_set
A QoISet to handle cases with multiple QoIs available.
virtual T dot(const NumericVector< T > &v) const =0
dof_id_type n_local_dofs() const
Definition: system.C:158
virtual std::pair< unsigned int, Real > adjoint_solve(const QoISet &qoi_indices=QoISet()) override
Assembles & solves the linear system (dR/du)^T*z = dq/du, for those quantities of interest q specifie...
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
dof_id_type n_dofs() const
Definition: system.C:121
This is the MeshBase class.
Definition: mesh_base.h:75
void enforce_adjoint_constraints_exactly(NumericVector< Number > &v, unsigned int q) const
Heterogeneously constrains the numeric vector v, which represents an adjoint solution defined on the ...
Definition: dof_map.h:2354
This class provides a specific system class.
Definition: diff_system.h:54
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
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
unsigned int number() const
Definition: system.h:2350
virtual void assembly(bool, bool, bool=false, bool=false)
Assembles a residual in rhs and/or a jacobian in matrix, as requested.
void remove_vector(std::string_view vec_name)
Removes the additional vector vec_name from this system.
Definition: system.C:873
Implements (adaptive) mesh refinement algorithms for a MeshBase.
bool has_adjoint_dirichlet_boundaries(unsigned int q) const
dof_id_type id() const
Definition: dof_object.h:828
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:96
virtual void reinit()
Handle any mesh changes and reinitialize all the systems on the updated mesh.
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 uniformly_coarsen(unsigned int n=1)
Attempts to uniformly coarsen the mesh n times.
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1593
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
This class holds functions that will estimate the error in a finite element solution on a given mesh...
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...
void set_project_with_constraints(bool _project_with_constraints)
Definition: system.h:1800
DifferentiablePhysics * _residual_evaluation_physics
Pointer to object to use for physics assembly evaluations.
bool is_adjoint_already_solved() const
Accessor for the adjoint_already_solved boolean.
Definition: system.h:409
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
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...
Real weight(std::size_t) const
Get the weight for this index (default 1.0)
Definition: qoi_set.h:243
const std::string & vector_name(const unsigned int vec_num) const
Definition: system.C:983
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:510
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void set_vector_preservation(const std::string &vec_name, bool preserve)
Allows one to set the boolean controlling whether the vector identified by vec_name should be "preser...
Definition: system.C:1138
virtual void swap(NumericVector< T > &v)
Swaps the contents of this with v.
virtual std::unique_ptr< DifferentiableQoI > clone() override
We don&#39;t allow systems to be attached to each other.
Definition: diff_system.h:168
bool & project_solution_on_reinit(void)
Tells the System whether or not to project the solution vector onto new grids when the system is rein...
Definition: system.h:838
const MeshBase & get_mesh() const
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:140
static std::unique_ptr< NumericVector< T > > build(const Parallel::Communicator &comm, SolverPackage solver_package=libMesh::default_solver_package(), ParallelType parallel_type=AUTOMATIC)
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
NumericVector< Number > & get_adjoint_solution(unsigned int i=0)
Definition: system.C:1245
std::unique_ptr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1605
virtual void clear()
Restores the NumericVector<T> to a pristine state.
bool get_project_with_constraints()
Setter and getter functions for project_with_constraints boolean.
Definition: system.h:1795
bool vector_preservation(std::string_view vec_name) const
Definition: system.C:1148
const DofMap & get_dof_map() const
Definition: system.h:2374
unsigned char number_h_refinements
How many h refinements to perform to get the fine grid.
const std::vector< dof_id_type > & get_send_list() const
Definition: dof_map.h:526
Manages consistently variables, degrees of freedom, and coefficient vectors for explicit systems...
const NumericVector< Number > & get_vector(std::string_view vec_name) const
Definition: system.C:943
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.
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.
Manages consistently variables, degrees of freedom, coefficient vectors, and matrices for implicit sy...
virtual void localize(std::vector< T > &v_local) const =0
Creates a copy of the global vector in the local vector v_local.