Skip to content
Snippets Groups Projects
Commit 954fa396 authored by Oliver Sander's avatar Oliver Sander Committed by sander@PCPOOL.MI.FU-BERLIN.DE
Browse files

moved to new module dune-solvers

[[Imported from SVN: r2349]]
parent 7430b4d8
No related branches found
No related tags found
No related merge requests found
#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";
}
}
#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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment