diff --git a/ugbase/bridge/disc_bridges/algebra_bridge.cpp b/ugbase/bridge/disc_bridges/algebra_bridge.cpp index c13faa218..87c9a32a7 100644 --- a/ugbase/bridge/disc_bridges/algebra_bridge.cpp +++ b/ugbase/bridge/disc_bridges/algebra_bridge.cpp @@ -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); } diff --git a/ugbase/lib_disc/operator/linear_operator/nested_iteration/nested_iteration.h b/ugbase/lib_disc/operator/linear_operator/nested_iteration/nested_iteration.h index 275df9eec..244fa6e41 100644 --- a/ugbase/lib_disc/operator/linear_operator/nested_iteration/nested_iteration.h +++ b/ugbase/lib_disc/operator/linear_operator/nested_iteration/nested_iteration.h @@ -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" @@ -138,7 +139,13 @@ class NestedIterationSolver //! set marking strategy void set_refinement_marking(SmartPtr > m) {m_spRefinementMarking =m; } void set_coarsening_marking(SmartPtr > 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 < > spSpace) + { m_spAssociatedSpace = spSpace;} + + void set_absolute_tolerance(number atol) + { + + m_absTOL = atol; + UG_LOG("NI::set_absolute_tolerance" << m_absTOL < > spMarking, bool bClearMarks = true); number refine_domain(const grid_function_type& u); number coarsen_domain(const grid_function_type& u); @@ -191,10 +209,18 @@ class NestedIterationSolver SmartPtr > m_spCoarseningMarking; bool m_bUseAdaptiveRefinement; - number m_TOL; int m_maxSteps; number m_lastError; + /// tolerance + number m_TOL; + + /// associated norm (for relative error) + SmartPtr > m_spAssociatedSpace; + + /// absolute tolerance (for relative error) + number m_absTOL; + /// (optional) debug info for adaptive refinement SmartPtr m_spElemError; diff --git a/ugbase/lib_disc/operator/linear_operator/nested_iteration/nested_iteration_impl.h b/ugbase/lib_disc/operator/linear_operator/nested_iteration/nested_iteration_impl.h index 1ae4ce08d..875fea665 100644 --- a/ugbase/lib_disc/operator/linear_operator/nested_iteration/nested_iteration_impl.h +++ b/ugbase/lib_disc/operator/linear_operator/nested_iteration/nested_iteration_impl.h @@ -62,7 +62,7 @@ NestedIterationSolver(SmartPtr > 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 @@ -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 @@ -86,7 +87,7 @@ NestedIterationSolver(SmartPtr > 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); }; @@ -102,7 +103,7 @@ NestedIterationSolver(SmartPtr > 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 >(new AssembledOperator(m_spAss)); }; @@ -150,12 +151,12 @@ bool NestedIterationSolver::prepare(vector_type& u) //! Refines domain and provides error estimate. /*! Values depend on m_spDomErr */ template -number NestedIterationSolver::refine_domain(const grid_function_type& u) +void NestedIterationSolver::estimate_and_mark_domain(const grid_function_type& u, SmartPtr > spMarking, bool bClearMarks) { UG_LOG("Computing error... "<< std::endl); typedef DomainDiscretization domain_disc_type; - // typedef IDomainErrorIndicator domain_indicator_type; + // typedef IDomainErrorIndicator domain_indicator_type; SmartPtr spDomainEstimator = m_spDomErr.template cast_dynamic(); @@ -165,19 +166,21 @@ number NestedIterationSolver::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; } @@ -286,14 +289,14 @@ bool NestedIterationSolver::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. @@ -315,25 +318,34 @@ bool NestedIterationSolver::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 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 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 < "<< m_TOL << std::endl); + // Refine and continue. + UG_LOG("FAILED: Error (still) greater than tolerance: " << err << " > "<< desiredTOL << std::endl); + m_spRefiner->refine(); } }