#include <adjoint_refinement_estimator.h>

Inheritance diagram for libMesh::AdjointRefinementEstimator:

Public Types
typedef std::map< std::pair< const System , unsigned int >, ErrorVector >	ErrorMap

Public Member Functions
	AdjointRefinementEstimator ()

	AdjointRefinementEstimator (const AdjointRefinementEstimator &)=default

	AdjointRefinementEstimator (AdjointRefinementEstimator &&)=default

AdjointRefinementEstimator &	operator= (const AdjointRefinementEstimator &)=default

AdjointRefinementEstimator &	operator= (AdjointRefinementEstimator &&)=default

virtual	~AdjointRefinementEstimator ()=default

QoISet &	qoi_set ()

const QoISet &	qoi_set () const

virtual void	estimate_error (const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=nullptr, bool estimate_parent_error=false)

Number &	get_global_QoI_error_estimate (unsigned int qoi_index)

virtual ErrorEstimatorType	type () const

DifferentiablePhysics *	get_residual_evaluation_physics ()

void	set_residual_evaluation_physics (DifferentiablePhysics *set_physics)

virtual void	estimate_errors (const EquationSystems &equation_systems, ErrorVector &error_per_cell, const std::map< const System , SystemNorm > &error_norms, const std::map< const System , const NumericVector< Number > > solution_vectors=nullptr, bool estimate_parent_error=false)

virtual void	estimate_errors (const EquationSystems &equation_systems, ErrorMap &errors_per_cell, const std::map< const System , const NumericVector< Number > > *solution_vectors=nullptr, bool estimate_parent_error=false)

Public Attributes
unsigned char	number_h_refinements

unsigned char	number_p_refinements

SystemNorm	error_norm

Protected Member Functions
void	reduce_error (std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm) const

Protected Attributes
DifferentiablePhysics *	_residual_evaluation_physics

std::vector< Number >	computed_global_QoI_errors

QoISet	_qoi_set

Detailed Description

This class implements a "brute force" goal-oriented error estimator which computes an estimate of error in a quantity of interest based on the residual of the current coarse grid primal solution as weighted against an adjoint solution on a uniformly refined (in h and/or p, for an arbitrary number of levels) grid.

Author: Roy H. Stogner

Date: 2009

Definition at line 50 of file adjoint_refinement_estimator.h.

Member Typedef Documentation

◆ ErrorMap

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

                                                        :
   ErrorEstimator(),
   number_h_refinements(1),
   number_p_refinements(0),
   _residual_evaluation_physics(nullptr),
   _qoi_set(QoISet())
 {
   // We're not actually going to use error_norm; our norms are
   // absolute values of QoI error.
   error_norm = INVALID_NORM;
 }

◆ 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`
	)

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

 {
   // We have to break the rules here, because we can't refine a const System
   System & system = const_cast<System &>(_system);
 
   // An EquationSystems reference will be convenient.
   EquationSystems & es = system.get_equation_systems();
 
   // The current mesh
   MeshBase & mesh = es.get_mesh();
 
   // Get coarse grid adjoint solutions.  This should be a relatively
   // quick (especially with preconditioner reuse) way to get a good
   // initial guess for the fine grid adjoint solutions.  More
   // importantly, subtracting off a coarse adjoint approximation gives
   // us better local error indication, and subtracting off *some* lift
   // function is necessary for correctness if we have heterogeneous
   // adjoint Dirichlet conditions.
 
   // Solve the adjoint problem(s) on the coarse FE space
   // Only if the user didn't already solve it for us
   if (!system.is_adjoint_already_solved())
     system.adjoint_solve(_qoi_set);
 
   // Loop over all the adjoint problems and, if any have heterogenous
   // Dirichlet conditions, get the corresponding coarse lift
   // function(s)
   for (unsigned int j=0,
        n_qois = system.n_qois(); j != n_qois; j++)
     {
       // Skip this QoI if it is not in the QoI Set or if there are no
       // heterogeneous Dirichlet boundaries for it
       if (_qoi_set.has_index(j) &&
           system.get_dof_map().has_adjoint_dirichlet_boundaries(j))
         {
           // Next, we are going to build up the residual for evaluating the flux QoI
           NumericVector<Number> * coarse_residual = nullptr;
 
           // The definition of a flux QoI is R(u^h, L) where R is the residual as defined
           // by a conservation law. Therefore, if we are using stabilization, the
           // R should be specified by the user via the residual_evaluation_physics
 
           // If the residual physics pointer is not null, use it when
           // evaluating here.
           {
             const bool swapping_physics = _residual_evaluation_physics;
             if (swapping_physics)
               dynamic_cast<DifferentiableSystem &>(system).swap_physics(_residual_evaluation_physics);
 
             // Assemble without applying constraints, to capture the solution values on the boundary
             (dynamic_cast<ImplicitSystem &>(system)).assembly(true, false, false, true);
 
             // Get the residual vector (no constraints applied on boundary, so we wont blow away the lift)
             coarse_residual = &(dynamic_cast<ExplicitSystem &>(system)).get_vector("RHS Vector");
             coarse_residual->close();
 
             // Now build the lift function and add it to the system vectors
             std::ostringstream liftfunc_name;
             liftfunc_name << "adjoint_lift_function" << j;
 
             system.add_vector(liftfunc_name.str());
 
             // Initialize lift with coarse adjoint solve associate with this flux QoI to begin with
             system.get_vector(liftfunc_name.str()).init(system.get_adjoint_solution(j), false);
 
             // Build the actual lift using adjoint dirichlet conditions
             system.get_dof_map().enforce_adjoint_constraints_exactly
               (system.get_vector(liftfunc_name.str()), static_cast<unsigned int>(j));
 
             // Compute the flux R(u^h, L)
             std::cout<<"The flux QoI "<<static_cast<unsigned int>(j)<<" is: "<<coarse_residual->dot(system.get_vector(liftfunc_name.str()))<<std::endl<<std::endl;
 
             // Swap back if needed
             if (swapping_physics)
               dynamic_cast<DifferentiableSystem &>(system).swap_physics(_residual_evaluation_physics);
           }
         } // End if QoI in set and flux/dirichlet boundary QoI
     } // End loop over QoIs
 
   // We'll want to back up all coarse grid vectors
   std::map<std::string, std::unique_ptr<NumericVector<Number>>> coarse_vectors;
   for (const auto & pr : as_range(system.vectors_begin(), system.vectors_end()))
     {
       // The (string) name of this vector
       const std::string & var_name = pr.first;
 
       coarse_vectors[var_name] = pr.second->clone();
     }
 
   // Back up the coarse solution and coarse local solution
   std::unique_ptr<NumericVector<Number>> coarse_solution = system.solution->clone();
   std::unique_ptr<NumericVector<Number>> coarse_local_solution = system.current_local_solution->clone();
 
   // And we'll need to temporarily change solution projection settings
   bool old_projection_setting;
   old_projection_setting = system.project_solution_on_reinit();
 
   // Make sure the solution is projected when we refine the mesh
   system.project_solution_on_reinit() = true;
 
   // And it'll be best to avoid any repartitioning
   std::unique_ptr<Partitioner> old_partitioner(mesh.partitioner().release());
 
   // And we can't allow any renumbering
   const bool old_renumbering_setting = mesh.allow_renumbering();
   mesh.allow_renumbering(false);
 
   // Use a non-standard solution vector if necessary
   if (solution_vector && solution_vector != system.solution.get())
     {
       NumericVector<Number> * newsol =
         const_cast<NumericVector<Number> *> (solution_vector);
       newsol->swap(*system.solution);
       system.update();
     }
 
   // Resize the error_per_cell vector to be
   // the number of elements, initialized to 0.
   error_per_cell.clear();
   error_per_cell.resize (mesh.max_elem_id(), 0.);
 
 #ifndef NDEBUG
   // These variables are only used in assertions later so
   // avoid declaring them unless asserts are active.
   const dof_id_type n_coarse_elem = mesh.n_active_elem();
 
   dof_id_type local_dof_bearing_nodes = 0;
   const unsigned int sysnum = system.number();
   for (const auto * node : mesh.local_node_ptr_range())
     for (unsigned int v=0, nvars = node->n_vars(sysnum);
          v != nvars; ++v)
       if (node->n_comp(sysnum, v))
         {
           local_dof_bearing_nodes++;
           break;
         }
 #endif // NDEBUG
 
   // Uniformly refine the mesh
   MeshRefinement mesh_refinement(mesh);
 
   libmesh_assert (number_h_refinements > 0 || number_p_refinements > 0);
 
   // FIXME: this may break if there is more than one System
   // on this mesh but estimate_error was still called instead of
   // estimate_errors
   for (unsigned int i = 0; i != number_h_refinements; ++i)
     {
       mesh_refinement.uniformly_refine(1);
       es.reinit();
     }
 
   for (unsigned int i = 0; i != number_p_refinements; ++i)
     {
       mesh_refinement.uniformly_p_refine(1);
       es.reinit();
     }
 
   // Copy the projected coarse grid solutions, which will be
   // overwritten by solve()
   std::vector<std::unique_ptr<NumericVector<Number>>> coarse_adjoints;
   for (unsigned int j=0; j != system.n_qois(); j++)
     {
       if (_qoi_set.has_index(j))
         {
           auto coarse_adjoint = NumericVector<Number>::build(mesh.comm());
 
           // Can do "fast" init since we're overwriting this in a sec
           coarse_adjoint->init(system.get_adjoint_solution(j),
                                /* fast = */ true);
 
           *coarse_adjoint = system.get_adjoint_solution(j);
 
           coarse_adjoints.emplace_back(std::move(coarse_adjoint));
         }
       else
         coarse_adjoints.emplace_back(nullptr);
     }
 
   // Next, we are going to build up the residual for evaluating the
   // error estimate.
 
   // If the residual physics pointer is not null, use it when
   // evaluating here.
   {
     const bool swapping_physics = _residual_evaluation_physics;
     if (swapping_physics)
       dynamic_cast<DifferentiableSystem &>(system).swap_physics(_residual_evaluation_physics);
 
     // Rebuild the rhs with the projected primal solution, constraints
     // have to be applied to get the correct error estimate since error
     // on the Dirichlet boundary is zero
     (dynamic_cast<ImplicitSystem &>(system)).assembly(true, false);
 
     // Swap back if needed
     if (swapping_physics)
       dynamic_cast<DifferentiableSystem &>(system).swap_physics(_residual_evaluation_physics);
   }
 
   NumericVector<Number> & projected_residual = (dynamic_cast<ExplicitSystem &>(system)).get_vector("RHS Vector");
   projected_residual.close();
 
   // Solve the adjoint problem(s) on the refined FE space
   system.adjoint_solve(_qoi_set);
 
   // Now that we have the refined adjoint solution and the projected primal solution,
   // we first compute the global QoI error estimate
 
   // Resize the computed_global_QoI_errors vector to hold the error estimates for each QoI
   computed_global_QoI_errors.resize(system.n_qois());
 
   // Loop over all the adjoint solutions and get the QoI error
   // contributions from all of them.  While we're looping anyway we'll
   // pull off the coarse adjoints
   for (unsigned int j=0; j != system.n_qois(); j++)
     {
       // Skip this QoI if not in the QoI Set
       if (_qoi_set.has_index(j))
         {
 
           // If the adjoint solution has heterogeneous dirichlet
           // values, then to get a proper error estimate here we need
           // to subtract off a coarse grid lift function.
           // |Q(u) - Q(u^h)| = |R([u^h]+, z^h+ - [L]+)| + HOT
           if(system.get_dof_map().has_adjoint_dirichlet_boundaries(j))
             {
               // Need to create a string with current loop index to retrieve
               // the correct vector from the liftvectors map
               std::ostringstream liftfunc_name;
               liftfunc_name << "adjoint_lift_function" << j;
 
               // Subtract off the corresponding lift vector from the adjoint solution
               system.get_adjoint_solution(j) -= system.get_vector(liftfunc_name.str());
 
               // Now evaluate R(u^h, z^h+ - lift)
               computed_global_QoI_errors[j] = projected_residual.dot(system.get_adjoint_solution(j));
 
               // Add the lift back to get back the adjoint
               system.get_adjoint_solution(j) += system.get_vector(liftfunc_name.str());
 
             }
           else
             {
               // Usual dual weighted residual error estimate
               // |Q(u) - Q(u^h)| = |R([u^h]+, z^h+)| + HOT
               computed_global_QoI_errors[j] = projected_residual.dot(system.get_adjoint_solution(j));
             }
 
         }
     }
 
   // Done with the global error estimates, now construct the element wise error indicators
 
   // To get a better element wise breakdown of the error estimate,
   // we subtract off a coarse representation of the adjoint solution.
   // |Q(u) - Q(u^h)| = |R([u^h]+, z^h+ - [z^h]+)|
   // This remains valid for all combinations of heterogenous adjoint bcs and
   // stabilized/non-stabilized formulations, except for the case where we not using a
   // heterogenous adjoint bc and have a stabilized formulation.
   // Then, R(u^h_s, z^h_s)  != 0 (no Galerkin orthogonality w.r.t the non-stabilized residual)
   for (unsigned int j=0; j != system.n_qois(); j++)
     {
       // Skip this QoI if not in the QoI Set
       if (_qoi_set.has_index(j))
         {
           // If we have a nullptr residual evaluation physics pointer, we
           // assume the user's formulation is consistent from mesh to
           // mesh, so we have Galerkin orthogonality and we can get
           // better indicator results by subtracting a coarse adjoint.
 
           // If we have a residual evaluation physics pointer, but we
           // also have heterogeneous adjoint dirichlet boundaries,
           // then we have to subtract off *some* lift function for
           // consistency, and we choose the coarse adjoint in lieu of
           // having any better ideas.
 
           // If we have a residual evaluation physics pointer and we
           // have homogeneous adjoint dirichlet boundaries, then we
           // don't have to subtract off anything, and with stabilized
           // formulations we get the best results if we don't.
           if(system.get_dof_map().has_adjoint_dirichlet_boundaries(j)
              || !_residual_evaluation_physics)
             {
               // z^h+ -> z^h+ - [z^h]+
               system.get_adjoint_solution(j) -= *coarse_adjoints[j];
             }
         }
     }
 
   // We ought to account for 'spill-over' effects while computing the
   // element error indicators This happens because the same dof is
   // shared by multiple elements, one way of mitigating this is to
   // scale the contribution from each dof by the number of elements it
   // belongs to We first obtain the number of elements each node
   // belongs to
 
   // A map that relates a node id to an int that will tell us how many elements it is a node of
   std::unordered_map<dof_id_type, unsigned int> shared_element_count;
 
   // To fill this map, we will loop over elements, and then in each element, we will loop
   // over the nodes each element contains, and then query it for the number of coarse
   // grid elements it was a node of
 
   // Keep track of which nodes we have already dealt with
   std::unordered_set<dof_id_type> processed_node_ids;
 
   // We will be iterating over all the active elements in the fine mesh that live on
   // this processor.
   for (const auto & elem : mesh.active_local_element_ptr_range())
     for (unsigned int n=0; n != elem->n_nodes(); ++n)
       {
         // Get a reference to the current node
         const Node & node = elem->node_ref(n);
 
         // Get the id of this node
         dof_id_type node_id = node.id();
 
         // If we havent already processed this node, do so now
         if (processed_node_ids.find(node_id) == processed_node_ids.end())
           {
             // Declare a neighbor_set to be filled by the find_point_neighbors
             std::set<const Elem *> fine_grid_neighbor_set;
 
             // Call find_point_neighbors to fill the neighbor_set
             elem->find_point_neighbors(node, fine_grid_neighbor_set);
 
             // A vector to hold the coarse grid parents neighbors
             std::vector<dof_id_type> coarse_grid_neighbors;
 
             // Loop over all the fine neighbors of this node
             for (const auto & fine_elem : fine_grid_neighbor_set)
               {
                 // Find the element id for the corresponding coarse grid element
                 const Elem * coarse_elem = fine_elem;
                 for (unsigned int j = 0; j != number_h_refinements; ++j)
                   {
                     libmesh_assert (coarse_elem->parent());
 
                     coarse_elem = coarse_elem->parent();
                   }
 
                 // Loop over the existing coarse neighbors and check if this one is
                 // already in there
                 const dof_id_type coarse_id = coarse_elem->id();
                 std::size_t j = 0;
 
                 // If the set already contains this element break out of the loop
                 for (; j != coarse_grid_neighbors.size(); j++)
                   if (coarse_grid_neighbors[j] == coarse_id)
                     break;
 
                 // If we didn't leave the loop even at the last element,
                 // this is a new neighbour, put in the coarse_grid_neighbor_set
                 if (j == coarse_grid_neighbors.size())
                   coarse_grid_neighbors.push_back(coarse_id);
               } // End loop over fine neighbors
 
             // Set the shared_neighbour index for this node to the
             // size of the coarse grid neighbor set
             shared_element_count[node_id] =
               cast_int<unsigned int>(coarse_grid_neighbors.size());
 
             // Add this node to processed_node_ids vector
             processed_node_ids.insert(node_id);
           } // End if not processed node
       } // End loop over nodes
 
   // Get a DoF map, will be used to get the nodal dof_indices for each element
   DofMap & dof_map = system.get_dof_map();
 
   // The global DOF indices, we will use these later on when we compute the element wise indicators
   std::vector<dof_id_type> dof_indices;
 
   // Localize the global rhs and adjoint solution vectors (which might be shared on multiple processors) onto a
   // local ghosted vector, this ensures each processor has all the dof_indices to compute an error indicator for
   // an element it owns
   std::unique_ptr<NumericVector<Number>> localized_projected_residual = NumericVector<Number>::build(system.comm());
   localized_projected_residual->init(system.n_dofs(), system.n_local_dofs(), system.get_dof_map().get_send_list(), false, GHOSTED);
   projected_residual.localize(*localized_projected_residual, system.get_dof_map().get_send_list());
 
   // Each adjoint solution will also require ghosting; for efficiency we'll reuse the same memory
   std::unique_ptr<NumericVector<Number>> localized_adjoint_solution = NumericVector<Number>::build(system.comm());
   localized_adjoint_solution->init(system.n_dofs(), system.n_local_dofs(), system.get_dof_map().get_send_list(), false, GHOSTED);
 
   // We will loop over each adjoint solution, localize that adjoint
   // solution and then loop over local elements
   for (unsigned int i=0; i != system.n_qois(); i++)
     {
       // Skip this QoI if not in the QoI Set
       if (_qoi_set.has_index(i))
         {
           // Get the weight for the current QoI
           Real error_weight = _qoi_set.weight(i);
 
           (system.get_adjoint_solution(i)).localize(*localized_adjoint_solution, system.get_dof_map().get_send_list());
 
           // Loop over elements
           for (const auto & elem : mesh.active_local_element_ptr_range())
             {
               // Go up number_h_refinements levels up to find the coarse parent
               const Elem * coarse = elem;
 
               for (unsigned int j = 0; j != number_h_refinements; ++j)
                 {
                   libmesh_assert (coarse->parent());
 
                   coarse = coarse->parent();
                 }
 
               const dof_id_type e_id = coarse->id();
 
               // Get the local to global degree of freedom maps for this element
               dof_map.dof_indices (elem, dof_indices);
 
               // We will have to manually do the dot products.
               Number local_contribution = 0.;
 
               // Sum the contribution to the error indicator for each element from the current QoI
               for (const auto & dof : dof_indices)
                 local_contribution += (*localized_projected_residual)(dof) * (*localized_adjoint_solution)(dof);
 
               // Multiply by the error weight for this QoI
               local_contribution *= error_weight;
 
               // FIXME: we're throwing away information in the
               // --enable-complex case
               error_per_cell[e_id] += static_cast<ErrorVectorReal>
                 (std::abs(local_contribution));
 
             } // End loop over elements
         } // End if belong to QoI set
     } // End loop over QoIs
 
   // Don't bother projecting the solution; we'll restore from backup
   // after coarsening
   system.project_solution_on_reinit() = false;
 
   // Uniformly coarsen the mesh, without projecting the solution
   libmesh_assert (number_h_refinements > 0 || number_p_refinements > 0);
 
   for (unsigned int i = 0; i != number_h_refinements; ++i)
     {
       mesh_refinement.uniformly_coarsen(1);
       // FIXME - should the reinits here be necessary? - RHS
       es.reinit();
     }
 
   for (unsigned int i = 0; i != number_p_refinements; ++i)
     {
       mesh_refinement.uniformly_p_coarsen(1);
       es.reinit();
     }
 
   // We should have the same number of active elements as when we started,
   // but not necessarily the same number of elements since reinit() doesn't
   // always call contract()
   libmesh_assert_equal_to (n_coarse_elem, mesh.n_active_elem());
 
   // We should have the same number of dof-bearing nodes as when we
   // started
 #ifndef NDEBUG
   dof_id_type final_local_dof_bearing_nodes = 0;
   for (const auto * node : mesh.local_node_ptr_range())
     for (unsigned int v=0, nvars = node->n_vars(sysnum);
          v != nvars; ++v)
       if (node->n_comp(sysnum, v))
         {
           final_local_dof_bearing_nodes++;
           break;
         }
   libmesh_assert_equal_to (local_dof_bearing_nodes,
                            final_local_dof_bearing_nodes);
 #endif // NDEBUG
 
   // Restore old solutions and clean up the heap
   system.project_solution_on_reinit() = old_projection_setting;
 
   // Restore the coarse solution vectors and delete their copies
   *system.solution = *coarse_solution;
   *system.current_local_solution = *coarse_local_solution;
 
   for (const auto & pr : as_range(system.vectors_begin(), system.vectors_end()))
     {
       // The (string) name of this vector
       const std::string & var_name = pr.first;
 
       // If it's a vector we already had (and not a newly created
       // vector like an adjoint rhs), we need to restore it.
       auto it = coarse_vectors.find(var_name);
       if (it != coarse_vectors.end())
         {
           NumericVector<Number> * coarsevec = it->second.get();
           system.get_vector(var_name) = *coarsevec;
 
           coarsevec->clear();
         }
     }
 
   // Restore old partitioner and renumbering settings
   mesh.partitioner().reset(old_partitioner.release());
   mesh.allow_renumbering(old_renumbering_setting);
 
   // Finally sum the vector of estimated error values.
   this->reduce_error(error_per_cell, system.comm());
 
   // We don't take a square root here; this is a goal-oriented
   // estimate not a Hilbert norm estimate.
 } // end estimate_error function

◆ 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(), and libMesh::EquationSystems::n_systems().

 {
   SystemNorm old_error_norm = this->error_norm;
 
   // Sum the error values from each system
   for (unsigned int s = 0; s != equation_systems.n_systems(); ++s)
     {
       ErrorVector system_error_per_cell;
       const System & sys = equation_systems.get_system(s);
       if (error_norms.find(&sys) == error_norms.end())
         this->error_norm = old_error_norm;
       else
         this->error_norm = error_norms.find(&sys)->second;
 
       const NumericVector<Number> * solution_vector = nullptr;
       if (solution_vectors &&
           solution_vectors->find(&sys) != solution_vectors->end())
         solution_vector = solution_vectors->find(&sys)->second;
 
       this->estimate_error(sys, system_error_per_cell,
                            solution_vector, estimate_parent_error);
 
       if (s)
         {
           libmesh_assert_equal_to (error_per_cell.size(), system_error_per_cell.size());
           for (auto i : index_range(error_per_cell))
             error_per_cell[i] += system_error_per_cell[i];
         }
       else
         error_per_cell = system_error_per_cell;
     }
 
   // Restore our old state before returning
   this->error_norm = old_error_norm;
 }

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

Currently this function ignores the error_norm.weight() values because it calculates each variable's error individually, unscaled.

The user selects which errors get computed by filling a map with error vectors: If errors_per_cell[&system][v] exists, it will be filled with the error values in variable v of system

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

Reimplemented in libMesh::UniformRefinementEstimator.

Definition at line 93 of file error_estimator.C.

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

 {
   SystemNorm old_error_norm = this->error_norm;
 
   // Find the requested error values from each system
   for (unsigned int s = 0; s != equation_systems.n_systems(); ++s)
     {
       const System & sys = equation_systems.get_system(s);
 
       unsigned int n_vars = sys.n_vars();
 
       for (unsigned int v = 0; v != n_vars; ++v)
         {
           // Only fill in ErrorVectors the user asks for
           if (errors_per_cell.find(std::make_pair(&sys, v)) ==
               errors_per_cell.end())
             continue;
 
           // Calculate error in only one variable
           std::vector<Real> weights(n_vars, 0.0);
           weights[v] = 1.0;
           this->error_norm =
             SystemNorm(std::vector<FEMNormType>(n_vars, old_error_norm.type(v)),
                        weights);
 
           const NumericVector<Number> * solution_vector = nullptr;
           if (solution_vectors &&
               solution_vectors->find(&sys) != solution_vectors->end())
             solution_vector = solution_vectors->find(&sys)->second;
 
           this->estimate_error
             (sys, *errors_per_cell[std::make_pair(&sys, v)],
              solution_vector, estimate_parent_error);
         }
     }
 
   // Restore our old state before returning
   this->error_norm = old_error_norm;
 }

◆ get_global_QoI_error_estimate()

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

inline

This is an accessor function to access the computed global QoI error estimates

Definition at line 106 of file adjoint_refinement_estimator.h.

References computed_global_QoI_errors.

   {
     return computed_global_QoI_errors[qoi_index];
   }

◆ 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; }

libMesh::AdjointRefinementEstimator::_residual_evaluation_physics

DifferentiablePhysics * _residual_evaluation_physics

Definition: adjoint_refinement_estimator.h:142

◆ 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

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

Definition at line 73 of file adjoint_refinement_estimator.h.

References _qoi_set.

73 { return _qoi_set; }

libMesh::AdjointRefinementEstimator::_qoi_set

QoISet _qoi_set

Definition: adjoint_refinement_estimator.h:150

◆ 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; }

libMesh::AdjointRefinementEstimator::_qoi_set

QoISet _qoi_set

Definition: adjoint_refinement_estimator.h:150

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

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

◆ 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; }

libMesh::AdjointRefinementEstimator::_residual_evaluation_physics

DifferentiablePhysics * _residual_evaluation_physics

Definition: adjoint_refinement_estimator.h:142

◆ type()

ErrorEstimatorType libMesh::AdjointRefinementEstimator::type ( ) const

virtual

Returns: The type for the ErrorEstimator subclass.

Implements libMesh::ErrorEstimator.

Definition at line 85 of file adjoint_refinement_estimator.C.

References libMesh::ADJOINT_REFINEMENT.

 {
   return ADJOINT_REFINEMENT;
 }

Member Data Documentation

◆ _qoi_set

QoISet libMesh::AdjointRefinementEstimator::_qoi_set

protected

A QoISet to handle cases with multiple QoIs available

Definition at line 150 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 142 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

Definition at line 145 of file adjoint_refinement_estimator.h.

Referenced by estimate_error(), and get_global_QoI_error_estimate().

◆ 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 161 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:

Public Types

Public Member Functions

Public Attributes

Protected Member Functions

Protected Attributes

Detailed Description

Member Typedef Documentation

◆ ErrorMap

Constructor & Destructor Documentation

◆ AdjointRefinementEstimator() [1/3]

◆ AdjointRefinementEstimator() [2/3]

◆ AdjointRefinementEstimator() [3/3]

◆ ~AdjointRefinementEstimator()

Member Function Documentation

◆ estimate_error()

◆ estimate_errors() [1/2]

◆ estimate_errors() [2/2]

◆ get_global_QoI_error_estimate()

◆ get_residual_evaluation_physics()

◆ operator=() [1/2]

◆ operator=() [2/2]

◆ qoi_set() [1/2]

◆ qoi_set() [2/2]

◆ reduce_error()

◆ set_residual_evaluation_physics()

◆ type()

Member Data Documentation

◆ _qoi_set

◆ _residual_evaluation_physics

◆ computed_global_QoI_errors

◆ error_norm

◆ number_h_refinements

◆ number_p_refinements