Skip to content

Commit

Permalink
Update for nested iteration (space-time-adaptive prototype)
Browse files Browse the repository at this point in the history
  • Loading branch information
@anaegel
anaegel committed Mar 29, 2018
1 parent 794b4d0 commit 36f3860
Show file tree
Hide file tree
Showing 3 changed files with 76 additions and 35 deletions.
3 changes: 3 additions & 0 deletions ugbase/bridge/disc_bridges/algebra_bridge.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -549,6 +549,9 @@ static void DomainAlgebra(Registry& reg, string parentGroup)
.add_method("enable_adaptive_refinement", &T::enable_adaptive_refinement)
.add_method("disable_adaptive_refinement", &T::disable_adaptive_refinement)
.add_method("config_string", &T::config_string)

.add_method("set_associated_space", &T::set_associated_space)
.add_method("set_absolute_tolerance", &T::set_absolute_tolerance)
.set_construct_as_smart_pointer(true);
reg.add_class_to_group(name, "NestedIterationSolver", tag);
}
Expand Down
32 changes: 29 additions & 3 deletions ugbase/lib_disc/operator/linear_operator/nested_iteration/nested_iteration.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@
#include "lib_disc/assemble_interface.h"

#include "lib_disc/spatial_disc/domain_disc.h"
#include "lib_disc/function_spaces/metric_spaces.h"
#include "lib_disc/function_spaces/error_elem_marking_strategy.h"

#include "lib_algebra/operator/debug_writer.h"
Expand Down Expand Up @@ -138,7 +139,13 @@ class NestedIterationSolver
//! set marking strategy
void set_refinement_marking(SmartPtr<IElementMarkingStrategy<TDomain> > m) {m_spRefinementMarking =m; }
void set_coarsening_marking(SmartPtr<IElementMarkingStrategy<TDomain> > m) {m_spCoarseningMarking =m; }
void set_tolerance(number tol) {m_TOL = tol;}
void set_tolerance(number tol)
{
m_TOL = tol;
UG_LOG("NI::set_tolerance" << m_TOL <<std::endl);
}



//! indicates, if grids should be refined further (if desired accuracy has not been achieved)
bool use_adaptive_refinement() const {return m_bUseAdaptiveRefinement;}
Expand All @@ -150,6 +157,17 @@ class NestedIterationSolver

void set_max_steps(int steps) {m_maxSteps = steps;}
number last_error() const {return m_lastError;}

void set_associated_space(SmartPtr<IGridFunctionSpace<grid_function_type> > spSpace)
{ m_spAssociatedSpace = spSpace;}

void set_absolute_tolerance(number atol)
{

m_absTOL = atol;
UG_LOG("NI::set_absolute_tolerance" << m_absTOL <<std::endl);
}

///@}

/// for debug output
Expand All @@ -159,7 +177,7 @@ class NestedIterationSolver
{m_spElemError = spErrEta;}

protected:
//number estimate_error(const grid_function_type& u);
void estimate_and_mark_domain(const grid_function_type& u, SmartPtr<IElementMarkingStrategy<TDomain> > spMarking, bool bClearMarks = true);
number refine_domain(const grid_function_type& u);
number coarsen_domain(const grid_function_type& u);

Expand Down Expand Up @@ -191,10 +209,18 @@ class NestedIterationSolver
SmartPtr<IElementMarkingStrategy<TDomain> > m_spCoarseningMarking;
bool m_bUseAdaptiveRefinement;

number m_TOL;
int m_maxSteps;
number m_lastError;

/// tolerance
number m_TOL;

/// associated norm (for relative error)
SmartPtr<IGridFunctionSpace<grid_function_type> > m_spAssociatedSpace;

/// absolute tolerance (for relative error)
number m_absTOL;

/// (optional) debug info for adaptive refinement
SmartPtr<grid_function_type> m_spElemError;

Expand Down
76 changes: 44 additions & 32 deletions ugbase/lib_disc/operator/linear_operator/nested_iteration/nested_iteration_impl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ NestedIterationSolver(SmartPtr<ILinearOperatorInverse<vector_type> > LinearSolve
m_spAss(NULL),
m_dgbCall(0),
m_lastNumSteps(0),
m_bUseAdaptiveRefinement(true)
m_bUseAdaptiveRefinement(true), m_TOL(1e-3), m_absTOL(1e-12), m_maxSteps(100)
{};

template <typename TDomain, typename TAlgebra>
Expand All @@ -74,7 +74,8 @@ NestedIterationSolver() :
m_spAss(NULL),
m_dgbCall(0),
m_lastNumSteps(0),
m_bUseAdaptiveRefinement(true)
m_bUseAdaptiveRefinement(true),
m_TOL(1e-3), m_absTOL(1e-12), m_maxSteps(100)
{};

template <typename TDomain, typename TAlgebra>
Expand All @@ -86,7 +87,7 @@ NestedIterationSolver(SmartPtr<IOperator<vector_type> > N) :
m_spAss(NULL),
m_dgbCall(0),
m_lastNumSteps(0),
m_bUseAdaptiveRefinement(true)
m_bUseAdaptiveRefinement(true), m_TOL(1e-3), m_absTOL(1e-12), m_maxSteps(100)
{
init(N);
};
Expand All @@ -102,7 +103,7 @@ NestedIterationSolver(SmartPtr<IAssemble<TAlgebra> > spAss) :
m_dgbCall(0),
m_lastNumSteps(0),
m_baseLevel(0),
m_bUseAdaptiveRefinement(true)
m_bUseAdaptiveRefinement(true), m_TOL(1e-3), m_absTOL(1e-12), m_maxSteps(100)
{
m_N = SmartPtr<AssembledOperator<TAlgebra> >(new AssembledOperator<TAlgebra>(m_spAss));
};
Expand Down Expand Up @@ -150,12 +151,12 @@ bool NestedIterationSolver<TDomain,TAlgebra>::prepare(vector_type& u)
//! Refines domain and provides error estimate.
/*! Values depend on m_spDomErr */
template <typename TDomain, typename TAlgebra>
number NestedIterationSolver<TDomain,TAlgebra>::refine_domain(const grid_function_type& u)
void NestedIterationSolver<TDomain,TAlgebra>::estimate_and_mark_domain(const grid_function_type& u, SmartPtr<IElementMarkingStrategy<TDomain> > spMarking, bool bClearMarks)
{
UG_LOG("Computing error... "<< std::endl);

typedef DomainDiscretization<TDomain, TAlgebra> domain_disc_type;
// typedef IDomainErrorIndicator<TAlgebra> domain_indicator_type;
// typedef IDomainErrorIndicator<TAlgebra> domain_indicator_type;

SmartPtr<domain_disc_type> spDomainEstimator = m_spDomErr.template cast_dynamic<domain_disc_type>();

Expand All @@ -165,19 +166,21 @@ number NestedIterationSolver<TDomain,TAlgebra>::refine_domain(const grid_functio

// compute error
if (m_spElemError.valid()) {
// debug version
spDomainEstimator->calc_error(u, u.dd(), &(*m_spElemError));
} else {
// standard version
spDomainEstimator->calc_error(u, u.dd());
}

// set (new) marks
m_spRefiner->clear_marks();
spDomainEstimator->mark_with_strategy(*m_spRefiner, m_spRefinementMarking);
UG_LOG("Estimated error: " << m_spRefinementMarking->global_estimated_error());
if (bClearMarks) m_spRefiner->clear_marks();
spDomainEstimator->mark_with_strategy(*m_spRefiner, spMarking);
UG_LOG("Estimated error: " << spMarking->global_estimated_error());



const number err = m_spRefinementMarking->global_estimated_error();
if(err >= m_TOL) {m_spRefiner->refine();}

return err;
}


Expand Down Expand Up @@ -286,14 +289,14 @@ bool NestedIterationSolver<TDomain,TAlgebra>::apply(vector_type& u)

// Solve linearized system.
try{
NESTED_ITER_PROFILE_BEGIN(NewtonApplyLinSolver);
if(!m_spLinearSolver->apply(u, *spD))
{
UG_LOG("ERROR in 'NestedIterationSolver::apply': Cannot apply Inverse Linear "
"Operator for Jacobi-Operator.\n");
return false;
}
NESTED_ITER_PROFILE_END();
NESTED_ITER_PROFILE_BEGIN(NewtonApplyLinSolver);
if(!m_spLinearSolver->apply(u, *spD))
{
UG_LOG("ERROR in 'NestedIterationSolver::apply': Cannot apply Inverse Linear "
"Operator for Jacobi-Operator.\n");
return false;
}
NESTED_ITER_PROFILE_END();
}UG_CATCH_THROW("NestedIterationSolver::apply: Application of Linear Solver failed.");

// Adjust solution.
Expand All @@ -315,25 +318,34 @@ bool NestedIterationSolver<TDomain,TAlgebra>::apply(vector_type& u)
break;
}

// Refine domain and check error.
number err = this->refine_domain(usol);
// Estimate and mark domain.
this->estimate_and_mark_domain(usol, m_spRefinementMarking);

if (m_spElemError.valid())
{
// write defect for debug
std::string name("NESTED_ITER_Error");
name.append(ext);
VTKOutput<TDomain::dim> outError;
outError.print(name.c_str(), *m_spElemError);
// OPTIONAL: write defect for debug
if (m_spElemError.valid())
{
std::string name("NESTED_ITER_Error");
name.append(ext);
VTKOutput<TDomain::dim> outError;
outError.print(name.c_str(), *m_spElemError);
}

// Check error.
const number err = m_spRefinementMarking->global_estimated_error();
double desiredTOL = (m_spAssociatedSpace.valid()) ? m_TOL*m_spAssociatedSpace->norm(usol) + m_absTOL : m_TOL;

UG_LOG("NI *** Desired tolerance: " << m_TOL << " * "<< m_spAssociatedSpace->norm(usol) << " + "<< m_absTOL <<std::endl);

// Quit or continue?
if(err < m_TOL)
if(err < desiredTOL)
{
UG_LOG("SUCCESS: Error smaller than tolerance: " << err << " < "<< m_TOL << std::endl);
// Quit.
UG_LOG("SUCCESS: Error smaller than tolerance: " << err << " < "<< desiredTOL << std::endl);
break;
} else {
UG_LOG("FAILED: Error (still) greater than tolerance: " << err << " > "<< m_TOL << std::endl);
// Refine and continue.
UG_LOG("FAILED: Error (still) greater than tolerance: " << err << " > "<< desiredTOL << std::endl);
m_spRefiner->refine();
}
}

Expand Down

0 comments on commit 36f3860

Please sign in to comment.