libMesh::AdjointResidualErrorEstimator Class Reference

#include <adjoint_residual_error_estimator.h>

Inheritance diagram for libMesh::AdjointResidualErrorEstimator:
Inheritance graph
[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 ()
const QoISet & qoi_set () const
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

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, (first proposed by Babuska and Miller in 1984) 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 55 of file adjoint_residual_error_estimator.h.

Member Typedef Documentation

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

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

Definition at line 107 of file error_estimator.h.

Constructor & Destructor Documentation

libMesh::AdjointResidualErrorEstimator::AdjointResidualErrorEstimator (  ) 

Constructor. Responsible for picking default subestimators.

Definition at line 40 of file adjoint_residual_error_estimator.C.

00040                                                               :  
00041   error_plot_suffix(),
00042   _primal_error_estimator(new PatchRecoveryErrorEstimator()),
00043   _dual_error_estimator(new PatchRecoveryErrorEstimator()),
00044   _qoi_set(QoISet())
00045 {
00046 }

libMesh::AdjointResidualErrorEstimator::~AdjointResidualErrorEstimator (  )  [inline]

Destructor.

Definition at line 67 of file adjoint_residual_error_estimator.h.

00067 {}

Member Function Documentation

AutoPtr<ErrorEstimator>& libMesh::AdjointResidualErrorEstimator::dual_error_estimator (  )  [inline]

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

Definition at line 79 of file adjoint_residual_error_estimator.h.

References _dual_error_estimator.

00079 { return _dual_error_estimator; }

void libMesh::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 libMesh::ErrorEstimator.

Definition at line 50 of file adjoint_residual_error_estimator.C.

References _dual_error_estimator, _primal_error_estimator, _qoi_set, libMesh::System::adjoint_solve(), libMesh::SystemNorm::calculate_norm(), libMesh::ErrorEstimator::error_norm, error_plot_suffix, libMesh::ErrorVectorReal, libMesh::AutoPtr< Tp >::get(), libMesh::System::get_adjoint_solution(), libMesh::System::get_equation_systems(), libMesh::System::get_mesh(), libMesh::QoISet::has_index(), libMesh::System::is_adjoint_already_solved(), libMesh::SystemNorm::is_identity(), libMesh::MeshBase::max_elem_id(), mesh, libMesh::System::n_vars(), n_vars, libMesh::ErrorVector::plot_error(), libMesh::System::qoi, libMesh::Real, libMesh::System::solution, and libMesh::QoISet::weight().

00054 {
00055   START_LOG("estimate_error()", "AdjointResidualErrorEstimator");
00056 
00057   // The current mesh
00058   const MeshBase& mesh = _system.get_mesh();
00059 
00060   // Resize the error_per_cell vector to be
00061   // the number of elements, initialize it to 0.
00062   error_per_cell.resize (mesh.max_elem_id());
00063   std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);
00064 
00065   // Get the number of variables in the system
00066   unsigned int n_vars = _system.n_vars();
00067 
00068   // We need to make a map of the pointer to the solution vector
00069   std::map<const System*, const NumericVector<Number>*>solutionvecs;
00070   solutionvecs[&_system] = _system.solution.get();
00071 
00072   // Solve the dual problem if we have to
00073   if (!_system.is_adjoint_already_solved())
00074     {
00075       // FIXME - we'll need to change a lot of APIs to make this trick
00076       // work with a const System...
00077       System&  system = const_cast<System&>(_system);
00078       system.adjoint_solve(_qoi_set);
00079     }
00080 
00081   // Flag to check whether we have not been asked to weight the variable error contributions in any specific manner
00082   bool error_norm_is_identity = error_norm.is_identity();
00083 
00084   // Create an ErrorMap/ErrorVector to store the primal, dual and total_dual variable errors
00085   ErrorMap primal_errors_per_cell;
00086   ErrorMap dual_errors_per_cell;
00087   ErrorMap total_dual_errors_per_cell;
00088   // Allocate ErrorVectors to this map if we're going to use it
00089   if (!error_norm_is_identity)
00090     for(unsigned int v = 0; v < n_vars; v++)
00091       {
00092         primal_errors_per_cell[std::make_pair(&_system, v)] = new ErrorVector;
00093         dual_errors_per_cell[std::make_pair(&_system, v)] = new ErrorVector;
00094         total_dual_errors_per_cell[std::make_pair(&_system, v)] = new ErrorVector;
00095       }
00096   ErrorVector primal_error_per_cell;
00097   ErrorVector dual_error_per_cell;
00098   ErrorVector total_dual_error_per_cell;
00099 
00100   // Have we been asked to weight the variable error contributions in any specific manner
00101   if(!error_norm_is_identity) // If we do
00102     {
00103       // Estimate the primal problem error for each variable
00104       _primal_error_estimator->estimate_errors
00105         (_system.get_equation_systems(), primal_errors_per_cell, &solutionvecs, estimate_parent_error);
00106     }
00107   else // If not
00108     {
00109       // Just get the combined error estimate
00110       _primal_error_estimator->estimate_error
00111         (_system, primal_error_per_cell, solution_vector, estimate_parent_error);
00112     }
00113 
00114   // Sum and weight the dual error estimate based on our QoISet
00115   for (unsigned int i = 0; i != _system.qoi.size(); ++i)
00116     {
00117       if (_qoi_set.has_index(i))
00118         {
00119           // Get the weight for the current QoI
00120           Real error_weight = _qoi_set.weight(i);
00121 
00122            // We need to make a map of the pointer to the adjoint solution vector
00123           std::map<const System*, const NumericVector<Number>*>adjointsolutionvecs;
00124           adjointsolutionvecs[&_system] = &_system.get_adjoint_solution(i);
00125 
00126           // Have we been asked to weight the variable error contributions in any specific manner
00127           if(!error_norm_is_identity) // If we have
00128             {
00129               _dual_error_estimator->estimate_errors
00130                 (_system.get_equation_systems(), dual_errors_per_cell, &adjointsolutionvecs,
00131                  estimate_parent_error);
00132             }
00133           else // If not
00134             {
00135             // Just get the combined error estimate
00136               _dual_error_estimator->estimate_error
00137                 (_system, dual_error_per_cell, &(_system.get_adjoint_solution(i)), estimate_parent_error);
00138             }
00139 
00140           std::size_t error_size;
00141 
00142           // Get the size of the first ErrorMap vector; this will give us the number of elements
00143           if(!error_norm_is_identity) // If in non default weights case
00144             {
00145               error_size = dual_errors_per_cell[std::make_pair(&_system, 0)]->size();
00146             }
00147           else // If in the standard default weights case
00148             {
00149               error_size = dual_error_per_cell.size();
00150             }
00151 
00152           // Resize the ErrorVector(s)
00153           if(!error_norm_is_identity)
00154             {
00155               // Loop over variables
00156               for(unsigned int v = 0; v < n_vars; v++)
00157                 {
00158                   libmesh_assert(!total_dual_errors_per_cell[std::make_pair(&_system, v)]->size() ||
00159                                  total_dual_errors_per_cell[std::make_pair(&_system, v)]->size() == error_size) ;
00160                   total_dual_errors_per_cell[std::make_pair(&_system, v)]->resize(error_size);
00161                 }
00162             }
00163           else
00164             {
00165               libmesh_assert(!total_dual_error_per_cell.size() ||
00166                              total_dual_error_per_cell.size() == error_size);
00167               total_dual_error_per_cell.resize(error_size);
00168             }
00169 
00170           for (std::size_t e = 0; e != error_size; ++e)
00171             {
00172               // Have we been asked to weight the variable error contributions in any specific manner
00173               if(!error_norm_is_identity) // If we have
00174                 {
00175                   // Loop over variables
00176                   for(unsigned int v = 0; v < n_vars; v++)
00177                     {
00178                       // Now fill in total_dual_error ErrorMap with the weight
00179                       (*total_dual_errors_per_cell[std::make_pair(&_system, v)])[e] +=
00180                         static_cast<ErrorVectorReal>
00181                         (error_weight * 
00182                          (*dual_errors_per_cell[std::make_pair(&_system, v)])[e]);
00183                     }
00184                 }
00185               else // If not
00186               {
00187                 total_dual_error_per_cell[e] +=
00188                   static_cast<ErrorVectorReal>(error_weight * dual_error_per_cell[e]);
00189               }
00190             }
00191         }
00192     }
00193 
00194   // Do some debugging plots if requested
00195   if (!error_plot_suffix.empty())
00196     {
00197       if(!error_norm_is_identity) // If we have
00198         {
00199           // Loop over variables
00200           for(unsigned int v = 0; v < n_vars; v++)
00201             {
00202               std::ostringstream primal_out;
00203               std::ostringstream dual_out;
00204               primal_out << "primal_" << error_plot_suffix << ".";
00205               dual_out << "dual_" << error_plot_suffix << ".";
00206 
00207               primal_out << std::setw(1)
00208                          << std::setprecision(0)
00209                          << std::setfill('0')
00210                          << std::right
00211                          << v;
00212 
00213               dual_out << std::setw(1)
00214                        << std::setprecision(0)
00215                        << std::setfill('0')
00216                        << std::right
00217                        << v;
00218 
00219               (*primal_errors_per_cell[std::make_pair(&_system, v)]).plot_error(primal_out.str(), _system.get_mesh());
00220               (*total_dual_errors_per_cell[std::make_pair(&_system, v)]).plot_error(dual_out.str(), _system.get_mesh());
00221 
00222               primal_out.clear();
00223               dual_out.clear();
00224             }
00225         }
00226       else // If not
00227         {
00228           std::ostringstream primal_out;
00229           std::ostringstream dual_out;
00230           primal_out << "primal_" << error_plot_suffix ;
00231           dual_out << "dual_" << error_plot_suffix ;
00232 
00233           primal_error_per_cell.plot_error(primal_out.str(), _system.get_mesh());
00234           total_dual_error_per_cell.plot_error(dual_out.str(), _system.get_mesh());
00235 
00236           primal_out.clear();
00237           dual_out.clear();
00238         }
00239     }
00240 
00241   // Weight the primal error by the dual error using the system norm object
00242   // FIXME: we ought to thread this
00243   for (unsigned int i=0; i != error_per_cell.size(); ++i)
00244     {
00245       // Have we been asked to weight the variable error contributions in any specific manner
00246       if(!error_norm_is_identity) // If we do
00247         {
00248           // Create Error Vectors to pass to calculate_norm
00249           std::vector<Real> cell_primal_error;
00250           std::vector<Real> cell_dual_error;
00251 
00252           for(unsigned int v = 0; v < n_vars; v++)
00253             {
00254               cell_primal_error.push_back((*primal_errors_per_cell[std::make_pair(&_system, v)])[i]);
00255               cell_dual_error.push_back((*total_dual_errors_per_cell[std::make_pair(&_system, v)])[i]);
00256             }
00257 
00258           error_per_cell[i] =
00259            static_cast<ErrorVectorReal>
00260              (error_norm.calculate_norm(cell_primal_error, cell_dual_error));
00261         }
00262       else // If not
00263         {
00264           error_per_cell[i] = primal_error_per_cell[i]*total_dual_error_per_cell[i];
00265         }
00266     }
00267 
00268   // Deallocate the ErrorMap contents if we allocated them earlier
00269   if (!error_norm_is_identity)
00270     for(unsigned int v = 0; v < n_vars; v++)
00271       {
00272         delete primal_errors_per_cell[std::make_pair(&_system, v)];
00273         delete dual_errors_per_cell[std::make_pair(&_system, v)];
00274         delete total_dual_errors_per_cell[std::make_pair(&_system, v)];
00275       }
00276 
00277   STOP_LOG("estimate_error()", "AdjointResidualErrorEstimator");
00278 }

void libMesh::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.

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(), libMesh::System::n_vars(), n_vars, and libMesh::SystemNorm::type().

00097 {
00098   SystemNorm old_error_norm = this->error_norm;
00099 
00100   // Find the requested error values from each system
00101   for (unsigned int s = 0; s != equation_systems.n_systems(); ++s)
00102     {
00103       const System &sys = equation_systems.get_system(s);
00104 
00105       unsigned int n_vars = sys.n_vars();
00106 
00107       for (unsigned int v = 0; v != n_vars; ++v)
00108         {
00109           // Only fill in ErrorVectors the user asks for
00110           if (errors_per_cell.find(std::make_pair(&sys, v)) ==
00111               errors_per_cell.end())
00112             continue;
00113 
00114           // Calculate error in only one variable
00115           std::vector<Real> weights(n_vars, 0.0);
00116           weights[v] = 1.0;
00117           this->error_norm =
00118             SystemNorm(std::vector<FEMNormType>(n_vars, old_error_norm.type(v)),
00119                        weights);
00120 
00121           const NumericVector<Number>* solution_vector = NULL;
00122           if (solution_vectors &&
00123               solution_vectors->find(&sys) != solution_vectors->end())
00124             solution_vector = solution_vectors->find(&sys)->second;
00125 
00126           this->estimate_error
00127             (sys, *errors_per_cell[std::make_pair(&sys, v)],
00128              solution_vector, estimate_parent_error);
00129         }
00130     }
00131 
00132   // Restore our old state before returning
00133   this->error_norm = old_error_norm;
00134 }

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 = 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 libMesh::UniformRefinementEstimator.

Definition at line 47 of file error_estimator.C.

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

00052 {
00053   SystemNorm old_error_norm = this->error_norm;
00054 
00055   // Sum the error values from each system
00056   for (unsigned int s = 0; s != equation_systems.n_systems(); ++s)
00057     {
00058       ErrorVector system_error_per_cell;
00059       const System &sys = equation_systems.get_system(s);
00060       if (error_norms.find(&sys) == error_norms.end())
00061         this->error_norm = old_error_norm;
00062       else
00063         this->error_norm = error_norms.find(&sys)->second;
00064 
00065       const NumericVector<Number>* solution_vector = NULL;
00066       if (solution_vectors &&
00067           solution_vectors->find(&sys) != solution_vectors->end())
00068         solution_vector = solution_vectors->find(&sys)->second;
00069 
00070       this->estimate_error(sys, system_error_per_cell,
00071                            solution_vector, estimate_parent_error);
00072 
00073       if (s)
00074         {
00075           libmesh_assert_equal_to (error_per_cell.size(), system_error_per_cell.size());
00076           for (unsigned int i=0; i != error_per_cell.size(); ++i)
00077             error_per_cell[i] += system_error_per_cell[i];
00078         }
00079       else
00080         error_per_cell = system_error_per_cell;
00081     }
00082 
00083   // Restore our old state before returning
00084   this->error_norm = old_error_norm;
00085 }

AutoPtr<ErrorEstimator>& libMesh::AdjointResidualErrorEstimator::primal_error_estimator (  )  [inline]

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

Definition at line 73 of file adjoint_residual_error_estimator.h.

References _primal_error_estimator.

00073 { return _primal_error_estimator; }

const QoISet& libMesh::AdjointResidualErrorEstimator::qoi_set (  )  const [inline]

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

Definition at line 91 of file adjoint_residual_error_estimator.h.

References _qoi_set.

00091 { return _qoi_set; }

QoISet& libMesh::AdjointResidualErrorEstimator::qoi_set (  )  [inline]

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

Definition at line 85 of file adjoint_residual_error_estimator.h.

References _qoi_set.

00085 { return _qoi_set; }

void libMesh::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 libMesh::WeightedPatchRecoveryErrorEstimator::estimate_error(), libMesh::PatchRecoveryErrorEstimator::estimate_error(), and libMesh::JumpErrorEstimator::estimate_error().

Member Data Documentation

AutoPtr<ErrorEstimator> libMesh::AdjointResidualErrorEstimator::_dual_error_estimator [protected]

An error estimator for the adjoint problem

Definition at line 125 of file adjoint_residual_error_estimator.h.

Referenced by dual_error_estimator(), and estimate_error().

AutoPtr<ErrorEstimator> libMesh::AdjointResidualErrorEstimator::_primal_error_estimator [protected]

An error estimator for the forward problem

Definition at line 120 of file adjoint_residual_error_estimator.h.

Referenced by estimate_error(), and primal_error_estimator().

QoISet libMesh::AdjointResidualErrorEstimator::_qoi_set [protected]

A QoISet to handle cases with multiple QoIs available

Definition at line 130 of file adjoint_residual_error_estimator.h.

Referenced by estimate_error(), and qoi_set().

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

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

std::string libMesh::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 100 of file adjoint_residual_error_estimator.h.

Referenced by estimate_error().

The documentation for this class was generated from the following files:
Site Created By: libMesh Developers
Last modified: February 05 2013 19:55:06 UTC

Hosted By:
SourceForge.net Logo