AdjointResidualErrorEstimator Class Reference

#include <adjoint_residual_error_estimator.h>

Inheritance diagram for AdjointResidualErrorEstimator:

[legend]

List of all members.

Public Types
typedef std::map< std::pair < const System *, unsigned int > , ErrorVector * >	ErrorMap
Public Member Functions
	AdjointResidualErrorEstimator ()
	~AdjointResidualErrorEstimator ()
AutoPtr< ErrorEstimator > &	primal_error_estimator ()
AutoPtr< ErrorEstimator > &	dual_error_estimator ()
QoISet &	qoi_set ()
void	estimate_error (const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=NULL, bool estimate_parent_error=false)
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=NULL, 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=NULL, bool estimate_parent_error=false)
Public Attributes
bool	adjoint_already_solved
std::string	error_plot_suffix
SystemNorm	error_norm
Protected Member Functions
void	reduce_error (std::vector< float > &error_per_cell) const
Protected Attributes
AutoPtr< ErrorEstimator >	_primal_error_estimator
AutoPtr< ErrorEstimator >	_dual_error_estimator
QoISet	_qoi_set

---

Detailed Description

This class implements a goal oriented error indicator, by weighting residual-based estimates from the primal problem against estimates from the adjoint problem.

This is based on a trick suggested by Brian Carnes, but if it doesn't actually work then the misunderstanding or misimplementation will be the fault of Roy Stogner. It's also Roy's fault there's no literature reference here yet.

Author:

Roy H. Stogner, 2009.

Definition at line 52 of file adjoint_residual_error_estimator.h.

---

Member Typedef Documentation

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

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

Definition at line 105 of file error_estimator.h.

---

Constructor & Destructor Documentation

AdjointResidualErrorEstimator::AdjointResidualErrorEstimator

(

)

Constructor. Responsible for picking default subestimators.

Definition at line 36 of file adjoint_residual_error_estimator.C.

                                                              :
  adjoint_already_solved(false),
  error_plot_suffix(),
  _primal_error_estimator(new PatchRecoveryErrorEstimator()),
  _dual_error_estimator(new PatchRecoveryErrorEstimator()),
  _qoi_set(QoISet())
{
}

AdjointResidualErrorEstimator::~AdjointResidualErrorEstimator

(

)

[inline]

Destructor.

Definition at line 64 of file adjoint_residual_error_estimator.h.

{}

---

Member Function Documentation

AutoPtr<ErrorEstimator>& AdjointResidualErrorEstimator::dual_error_estimator

(

)

[inline]

Access to the "subestimator" (default: PatchRecovery) to use on the dual/adjoint solution

Definition at line 76 of file adjoint_residual_error_estimator.h.

References _dual_error_estimator.

{ return _dual_error_estimator; }

void AdjointResidualErrorEstimator::estimate_error	(	const System &	system,
		ErrorVector &	error_per_cell,
		const NumericVector< Number > *	solution_vector = NULL,
		bool	estimate_parent_error = false
	)			[virtual]

Compute the adjoint-weighted error on each element and place it in the error_per_cell vector. Note that this->error_norm is ignored; the error estimate is in the seminorm given by the absolute value of the error in the quantity of interest functional. The primal and dual subestimator error_norm values are used, and should be chosen appropriately for your model.

Implements ErrorEstimator.

Definition at line 47 of file adjoint_residual_error_estimator.C.

References _dual_error_estimator, _primal_error_estimator, _qoi_set, adjoint_already_solved, System::adjoint_solve(), error_plot_suffix, System::get_adjoint_solution(), System::get_mesh(), QoISet::has_index(), ErrorVector::plot_error(), System::qoi, and QoISet::weight().

{
  START_LOG("estimate_error()", "AdjointResidualErrorEstimator");

  // Start by estimating the primal problem error
  _primal_error_estimator->estimate_error
    (_system, error_per_cell, solution_vector, estimate_parent_error);

  // Solve the dual problem if we have to
  if (!adjoint_already_solved)
    {
      // FIXME - we'll need to change a lot of APIs to make this trick
      // work with a const System...
      System&  system = const_cast<System&>(_system);
      system.adjoint_solve(_qoi_set);
    }

  // This bookkeeping should now be taken care of in subestimators
  // when they see a non-default solution_vector
  //
  // // Swap the system solution and dual solution, so our next error
  // // estimate will use the latter
  // // system.solution->swap(system.get_adjoint_solution());

  // // Don't forget to keep the current_local_solution consistent, as
  // // error estimators are likely to use that!
  // // system.update();

  // Get a separate estimate of the dual problem error
  ErrorVector total_dual_error_per_cell;

  // Sum and weight this estimate based on our QoISet
  for (unsigned int i = 0; i != _system.qoi.size(); ++i)
    if (_qoi_set.has_index(i))
      {
        Real error_weight = _qoi_set.weight(i);

        ErrorVector dual_error_per_cell;
        _dual_error_estimator->estimate_error
          (_system, dual_error_per_cell, &(_system.get_adjoint_solution(i)),
           estimate_parent_error);

        unsigned int error_size = dual_error_per_cell.size();
        libmesh_assert(!total_dual_error_per_cell.size() ||
                       total_dual_error_per_cell.size() == error_size);
        total_dual_error_per_cell.resize(error_size);
        for (unsigned int e = 0; e != error_size; ++e)
          total_dual_error_per_cell[e] += 
            error_weight * dual_error_per_cell[e];
      }

  // Do some debugging plots if requested
  if (!error_plot_suffix.empty())
    {
      std::string primal_out = "primal_";
      std::string dual_out = "dual_";
      primal_out += error_plot_suffix;
      dual_out += error_plot_suffix;
      error_per_cell.plot_error(primal_out, _system.get_mesh());
      total_dual_error_per_cell.plot_error(dual_out, _system.get_mesh());
    }

  // More bookkeeping for subestimators to do
  //
  // // Swap the system solution and dual solution back to their proper
  // // places
  // system.solution->swap(system.get_adjoint_solution());

  // // Don't forget to fix the current_local_solution
  // system.update();

  // Weight the primal error by the dual error
  // FIXME: we ought to thread this
  for (unsigned int i=0; i != error_per_cell.size(); ++i)
    error_per_cell[i] *= total_dual_error_per_cell[i];

  STOP_LOG("estimate_error()", "AdjointResidualErrorEstimator");
}

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

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.

Definition at line 92 of file error_estimator.C.

References ErrorEstimator::error_norm, ErrorEstimator::estimate_error(), EquationSystems::get_system(), EquationSystems::n_systems(), System::n_vars(), n_vars, and 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(0)),
                       weights);

          const NumericVector<Number>* solution_vector = NULL;
          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;
}

void 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 = NULL,
		bool	estimate_parent_error = false
	)			[virtual, inherited]

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

Definition at line 46 of file error_estimator.C.

References ErrorEstimator::error_norm, ErrorEstimator::estimate_error(), EquationSystems::get_system(), and 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 = NULL;
      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(error_per_cell.size() == system_error_per_cell.size());
          for (unsigned int i=0; i != error_per_cell.size(); ++i)
            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;
}

AutoPtr<ErrorEstimator>& AdjointResidualErrorEstimator::primal_error_estimator

(

)

[inline]

Access to the "subestimator" (default: PatchRecovery) to use on the primal/forward solution

Definition at line 70 of file adjoint_residual_error_estimator.h.

References _primal_error_estimator.

{ return _primal_error_estimator; }

QoISet& AdjointResidualErrorEstimator::qoi_set

(

)

[inline]

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

Definition at line 82 of file adjoint_residual_error_estimator.h.

References _qoi_set.

{ return _qoi_set; }

void ErrorEstimator::reduce_error

(

std::vector< float > &

error_per_cell

)

const [protected, inherited]

This method takes the local error contributions in error_per_cell from each processor and combines them to get the global error vector.

Referenced by PatchRecoveryErrorEstimator::estimate_error(), and JumpErrorEstimator::estimate_error().

---

Member Data Documentation

AutoPtr<ErrorEstimator> AdjointResidualErrorEstimator::_dual_error_estimator [protected]

An error estimator for the adjoint problem

Definition at line 123 of file adjoint_residual_error_estimator.h.

Referenced by dual_error_estimator(), and estimate_error().

AutoPtr<ErrorEstimator> AdjointResidualErrorEstimator::_primal_error_estimator [protected]

An error estimator for the forward problem

Definition at line 118 of file adjoint_residual_error_estimator.h.

Referenced by estimate_error(), and primal_error_estimator().

QoISet AdjointResidualErrorEstimator::_qoi_set [protected]

A QoISet to handle cases with multiple QoIs available

Definition at line 128 of file adjoint_residual_error_estimator.h.

Referenced by estimate_error(), and qoi_set().

bool AdjointResidualErrorEstimator::adjoint_already_solved

Has the adjoint problem already been solved? If the user sets adjoint_already_solved to true, we won't waste time solving it again.

Definition at line 89 of file adjoint_residual_error_estimator.h.

Referenced by estimate_error().

SystemNorm 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 137 of file error_estimator.h.

Referenced by KellyErrorEstimator::boundary_side_integration(), DiscontinuityMeasure::boundary_side_integration(), DiscontinuityMeasure::DiscontinuityMeasure(), JumpErrorEstimator::estimate_error(), ErrorEstimator::estimate_errors(), ExactErrorEstimator::ExactErrorEstimator(), ExactErrorEstimator::find_squared_element_error(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), DiscontinuityMeasure::internal_side_integration(), KellyErrorEstimator::KellyErrorEstimator(), LaplacianErrorEstimator::LaplacianErrorEstimator(), PatchRecoveryErrorEstimator::EstimateError::operator()(), PatchRecoveryErrorEstimator::PatchRecoveryErrorEstimator(), and UniformRefinementEstimator::UniformRefinementEstimator().

std::string AdjointResidualErrorEstimator::error_plot_suffix

To aid in investigating error estimator behavior, set this string to a suffix with which to plot (prefixed by "primal_" or "dual_") the subestimator results. The suffix should end with a file extension (e.g. ".gmv") that the ErrorVector::plot_error recognizes.

Definition at line 98 of file adjoint_residual_error_estimator.h.

Referenced by estimate_error().

---

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

• adjoint_residual_error_estimator.h
• adjoint_residual_error_estimator.C