stationaryadaptivealgorithm.hh

#ifndef __DUNE_ACFEM_ALGORITHMS_STATIONARY_ADAPTIVE_ALGORITHM_HH__
#define __DUNE_ACFEM_ALGORITHMS_STATIONARY_ADAPTIVE_ALGORITHM_HH__
 
#include <dune/fem/solver/timeprovider.hh>
#include "../algorithms/femschemeinterface.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    static inline
    std::pair<double, double>
    adaptiveAlgorithm(AdaptiveFemScheme& scheme, const std::string& prefix_ = "")
    {
      std::string prefix = prefix_ == std::string("") ? prefix_ : prefix_ + ".";
 
      const bool verbose = Dune::Fem::Parameter::getValue<int>("fem.verboserank", -1) != -1;
 
      // initialize discrete solution
      scheme.initialize();
 
      // solve once
      scheme.solve();
 
      // write initial solve
      scheme.output();
 
      // calculate error
      double error = -1., oldError;
      double estimate = 1., oldEstimate;
      size_t size = 0, oldSize;
 
      // get tolerance for adaptive algorithm
      const double tolerance = Dune::Fem::Parameter::getValue<double>(prefix+"adaptation.tolerance", 0.1);
 
      // get error estimate
      estimate = scheme.estimate();
 
      // TODO: maybe trigger computation of errors by a dedicated
      // flag parameter (also second error computation below)
      if (verbose) {
        error = scheme.error();
        size_t size = scheme.size(); // attention: communication, NOT inside verbose()!!!
        if (Fem::Parameter::verbose()) {
          std::cout << "Error: " << error << std::endl;
          std::cout << "#DoFs: " << size << std::endl;
        }
      }
 
      // until estimator is below tolerance for *THE LOOP*
      while (estimate > tolerance) {
        // mark element for adaptation
        scheme.mark(tolerance);
 
        // adapt grid
        scheme.adapt();
 
        // solve again
        scheme.solve();
 
        // data I/O
        scheme.output();
 
        // save old estimate
        oldEstimate = estimate;
 
        // calculate new error
        estimate = scheme.estimate();
 
        // TODO: maybe trigger computation of errors by a dedicated
        // flag parameter (also second error computation below)
        if (verbose) {
          oldError = error;
          error = scheme.error();
          oldSize = size;
          size = scheme.size(); // attention: communication, NOT inside Parameter::verbose()!!!
          if (Fem::Parameter::verbose()) {
            std::cout << "Estimate/Error: " << estimate << ", " << error << std::endl;
            std::cout << "#DoFs:          " << size << std::endl;
            if (oldEstimate > 0) {
              double estimatorRatio = oldEstimate/estimate;
              double errorRatio = oldError/error;
              double sizeRatio = (double)size/(double)oldSize;
 
              std::cout << "Efficiency:     "
                        << std::log(estimatorRatio)/std::log(sizeRatio)
                        << ", "
                        << std::log(errorRatio)/std::log(sizeRatio)
                        << std::endl;
            }
          }
        }
      }
 
      if (!verbose) {
        error = scheme.error();
      }
 
      return std::pair<double,double>(estimate, error);
    }
 
 
  } // ACFem
 
} // Dune
 
#endif //  __DUNE_ACFEM_ALGORITHMS_STATIONARY_ADAPTIVE_ALGORITHM_HH__