diff --git a/dune-solvers/solvers/loopsolver.cc b/dune-solvers/solvers/loopsolver.cc new file mode 100644 index 0000000000000000000000000000000000000000..2e027f4c7188427f772a232a20f0165f466f3495 --- /dev/null +++ b/dune-solvers/solvers/loopsolver.cc @@ -0,0 +1,156 @@ +#include <cmath> +#include <limits> +#include <iostream> +#include <iomanip> + +template <class VectorType, class BitVectorType> +void LoopSolver<VectorType, BitVectorType>::check() const +{ + if (!iterationStep_) + DUNE_THROW(SolverError, "You need to supply an iteration step to an iterative solver!"); + + iterationStep_->check(); + + // check base class + IterativeSolver<VectorType,BitVectorType>::check(); +} + +template <class VectorType, class BitVectorType> +void LoopSolver<VectorType, BitVectorType>::solve() +{ + + int i; + + // Check whether the solver is set up properly + this->check(); + + this->iterationStep_->preprocess(); + + if (this->verbosity_ != NumProc::QUIET) + std::cout << "--- LoopSolver ---\n"; + + if (this->verbosity_ == NumProc::FULL) + { + std::cout << " iter"; + if (referenceSolution_) + { + if (this->useRelativeError_) + std::cout << " error"; + std::cout << " abs error"; + std::cout << " abs correction"; + } + else + { + if (this->useRelativeError_) + std::cout << " correction"; + std::cout << " abs correction"; + } + std::cout << " rate"; + std::string header = this->iterationStep_->getOutput(); + std::cout << header; + std::cout << std::endl; + std::cout << "-----"; + if (this->useRelativeError_) + std::cout << "---------------"; + if (referenceSolution_) + std::cout << "---------------"; + std::cout << "---------------"; + std::cout << "---------"; + for(int i=0; i<header.size(); ++i) + std::cout << "-"; + std::cout << std::endl; + } + + double error = std::numeric_limits<double>::max(); + + double normOfOldCorrection = 0; + double normOfOldError = 0; + double totalConvRate = 1; + this->maxTotalConvRate_ = 0; + int convRateCounter = 0; + + // Loop until desired tolerance or maximum number of iterations is reached + for (i=0; i<this->maxIterations_ && (error>this->tolerance_ || std::isnan(error)); i++) { + + // Backup of the current solution for the error computation later on + VectorType oldSolution = iterationStep_->getSol(); + + // Perform one iteration step + iterationStep_->iterate(); + + // write iteration to file, if requested + if (this->historyBuffer_!="") + this->writeIterate(iterationStep_->getSol(), i); + + // Compute error + double oldNorm = this->errorNorm_->operator()(oldSolution); + oldSolution -= iterationStep_->getSol(); + + double normOfCorrection; + double normOfError; + double convRate; + + normOfCorrection = this->errorNorm_->operator()(oldSolution); + convRate = normOfCorrection / normOfOldCorrection; + error = normOfCorrection; + normOfOldCorrection = normOfCorrection; + + if (referenceSolution_) + { + normOfError = this->errorNorm_->diff(iterationStep_->getSol(), *referenceSolution_); + convRate = normOfError / normOfOldError; + error = normOfError; + normOfOldError = normOfError; + } + + if (this->useRelativeError_) + error = error / oldNorm; + else + error = error; + + if (!isinf(convRate) && !isnan(convRate)) { + totalConvRate *= convRate; + this->maxTotalConvRate_ = std::max(this->maxTotalConvRate_, std::pow(totalConvRate, 1/((double)convRateCounter+1))); + convRateCounter++; + } + + // Output + if (this->verbosity_ == NumProc::FULL) { + std::cout << std::setw(5) << i; + + if (this->useRelativeError_) + { + std::cout << std::setiosflags(std::ios::scientific); + std::cout << std::setw(15) << std::setprecision(7) << error; + std::cout << std::resetiosflags(std::ios::scientific); + } + + if (referenceSolution_) + { + std::cout << std::setiosflags(std::ios::scientific); + std::cout << std::setw(15) << std::setprecision(7) << normOfError; + std::cout << std::resetiosflags(std::ios::scientific); + } + + std::cout << std::setiosflags(std::ios::scientific); + std::cout << std::setw(15) << std::setprecision(7) << normOfCorrection; + std::cout << std::resetiosflags(std::ios::scientific); + + std::cout << std::setiosflags(std::ios::fixed); + std::cout << std::setw(9) << std::setprecision(5) << convRate; + std::cout << std::resetiosflags(std::ios::fixed); + + std::cout << this->iterationStep_->getOutput(); + std::cout << std::endl; + } + + } + + + if (this->verbosity_ != NumProc::QUIET) { + std::cout << "maxTotalConvRate: " << this->maxTotalConvRate_ << ", " + << i << " iterations performed\n"; + std::cout << "--------------------\n"; + } + +} diff --git a/dune-solvers/solvers/loopsolver.hh b/dune-solvers/solvers/loopsolver.hh new file mode 100644 index 0000000000000000000000000000000000000000..fa6fa84d8fa36836dde78eb2873c8535b1a9443e --- /dev/null +++ b/dune-solvers/solvers/loopsolver.hh @@ -0,0 +1,53 @@ +#ifndef LOOP_SOLVER_HH +#define LOOP_SOLVER_HH + +#include <dune/ag-common/iterativesolver.hh> +#include <dune/ag-common/iterationstep.hh> +#include <dune/ag-common/norms/norm.hh> + +/** \brief A solver which consists of a single loop + * + This class basically implements a loop that calls an iteration procedure + (which is to be supplied by the user). It also monitors convergence. +*/ +template <class VectorType, class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension> > +class LoopSolver : public IterativeSolver<VectorType, BitVectorType> +{ +public: + + /** \brief Constructor taking all relevant data */ + LoopSolver(IterationStep<VectorType, BitVectorType>* iterationStep, + int maxIterations, + double tolerance, + Norm<VectorType>* errorNorm, + Solver::VerbosityMode verbosity, + bool useRelativeError=true, + const VectorType* referenceSolution=0) + : IterativeSolver<VectorType, BitVectorType>(maxIterations, + tolerance, + errorNorm, + verbosity, + useRelativeError), + iterationStep_(iterationStep), + referenceSolution_(referenceSolution) + {} + + /** \brief Checks whether all relevant member variables are set + * \exception SolverError if the iteration step is not set up properly + */ + virtual void check() const; + + /** \brief Loop, call the iteration procedure + * and monitor convergence */ + virtual void solve(); + + //! The iteration step used by the algorithm + IterationStep<VectorType, BitVectorType>* iterationStep_; + + const VectorType* referenceSolution_; + +}; + +#include "loopsolver.cc" + +#endif