#include <adjoint_residual_error_estimator.h>
typedef std::map< std::pair< const System *, unsigned int >, ErrorVector * > ErrorMap
AdjointResidualErrorEstimator ()
~AdjointResidualErrorEstimator ()
AutoPtr< ErrorEstimator > & primal_error_estimator ()
AutoPtr< ErrorEstimator > & dual_error_estimator ()
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)
bool adjoint_already_solved
std::string error_plot_suffix
SystemNorm error_norm
void reduce_error (std::vector< float > &error_per_cell) const
AutoPtr< ErrorEstimator > _primal_error_estimator
AutoPtr< ErrorEstimator > _dual_error_estimator
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:
Definition at line 51 of file adjoint_residual_error_estimator.h.
Definition at line 105 of file error_estimator.h.
Definition at line 35 of file adjoint_residual_error_estimator.C.
:
adjoint_already_solved(false),
error_plot_suffix(),
_primal_error_estimator(new PatchRecoveryErrorEstimator()),
_dual_error_estimator(new PatchRecoveryErrorEstimator())
{
}
Definition at line 63 of file adjoint_residual_error_estimator.h.
{}
Definition at line 73 of file adjoint_residual_error_estimator.h.
References _dual_error_estimator.
{ return _dual_error_estimator; }
Implements ErrorEstimator.
Definition at line 45 of file adjoint_residual_error_estimator.C.
References _dual_error_estimator, _primal_error_estimator, adjoint_already_solved, System::adjoint_solve(), error_plot_suffix, System::get_adjoint_solution(), System::get_mesh(), and ErrorVector::plot_error().
{
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();
}
// 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 dual_error_per_cell;
_dual_error_estimator->estimate_error
(_system, dual_error_per_cell, &(_system.get_adjoint_solution()),
estimate_parent_error);
// 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());
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] *= dual_error_per_cell[i];
STOP_LOG('estimate_error()', 'AdjointResidualErrorEstimator');
}
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 UniformRefinementEstimator.
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(), 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;
}
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;
}
Definition at line 68 of file adjoint_residual_error_estimator.h.
References _primal_error_estimator.
{ return _primal_error_estimator; }
Definition at line 32 of file error_estimator.C.
Referenced by UniformRefinementEstimator::_estimate_error(), PatchRecoveryErrorEstimator::estimate_error(), and JumpErrorEstimator::estimate_error().
{
// This function must be run on all processors at once
parallel_only();
// Each processor has now computed the error contribuions
// for its local elements. We may need to sum the vector to
// recover the error for each element.
Parallel::sum(error_per_cell);
}
Definition at line 114 of file adjoint_residual_error_estimator.h.
Referenced by dual_error_estimator(), and estimate_error().
Definition at line 109 of file adjoint_residual_error_estimator.h.
Referenced by estimate_error(), and primal_error_estimator().
Definition at line 80 of file adjoint_residual_error_estimator.h.
Referenced by estimate_error().
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 UniformRefinementEstimator::_estimate_error(), 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().
Definition at line 89 of file adjoint_residual_error_estimator.h.
Referenced by estimate_error().
Generated automatically by Doxygen for libMesh from the source code.