transientadaptivealgorithm.hh

#ifndef __DUNE_ACFEM_ALGORITHMS_TRANSIENT_ADAPTIVE_ALGORITHM_HH__
#define __DUNE_ACFEM_ALGORITHMS_TRANSIENT_ADAPTIVE_ALGORITHM_HH__
 
#include <dune/fem/solver/timeprovider.hh>
#include "../algorithms/femschemeinterface.hh"
 
// for log-file
#include <fstream>
 
namespace Dune {
 
  namespace ACFem {
 
    template<class TimeProvider>
    static inline
    double adaptiveAlgorithm(TransientAdaptiveFemScheme& scheme,
                             TimeProvider& timeProvider,
                             double endTime,
                             const std::string& prefix_ = "")
    {
      std::string prefix = prefix_ == std::string("") ? prefix_ : prefix_ + ".";
 
      const bool verbose = Dune::Fem::Parameter::getValue<int>("fem.verboserank", -1) != -1;
 
      bool restart = Dune::Fem::Parameter::getValue<bool>(prefix + "restart", false);
 
      // initialtimeStep from initilized TimeProvider.
      double timeStep = timeProvider.deltaT();
      // The time-step and mu are intentionally non-const. timeStep is
      // only the initial time-step size for the "implicit"
      // time-adaptation method
 
      const double tolerance = Dune::Fem::Parameter::getValue<double>(prefix + "adaptation.tolerance", 1e-2);
      const double initialTolerance =
        tolerance * Dune::Fem::Parameter::getValue<double>(prefix + "adaptation.initialToleranceFraction", 0.2);
      const double spaceTolerance =
        tolerance * Dune::Fem::Parameter::getValue<double>(prefix + "adaptation.spaceToleranceFraction", 0.4);
      const double timeTolerance =
        tolerance * Dune::Fem::Parameter::getValue<double>(prefix + "adaptation.timeToleranceFraction", 0.4);
      const double minDeltaT =
        Dune::Fem::Parameter::getValue<double>(prefix + "adaptation.time.minimalStep", -1);
      const double maxDeltaT =
        Dune::Fem::Parameter::getValue<double>(prefix + "adaptation.time.maximalStep", -1);
 
      const int maximumIterations =
        Dune::Fem::Parameter::getValue<int>(prefix + "adaptation.maximumIterations", 10);
      const int spaceMaximumIterations =
        Dune::Fem::Parameter::getValue<int>(prefix + "adaptation.space.maximumIterations", 10);
      const int initialMaximumIterations =
        Dune::Fem::Parameter::getValue<int>(prefix + "adaptation.initial.maximumIterations", 10);
 
      const std::string timeStrategyNames[] = { "explicit", "implicit" };
      enum { explicitTimeStrategy, implicitTimeStrategy }; // by purpose explicit == false
      const int timeStrategy = Dune::Fem::Parameter::getEnum(prefix + "adaptation.time.strategy", timeStrategyNames, explicitTimeStrategy);
      const double timeNarrowFactor = Dune::Fem::Parameter::getValue<double>(prefix + "adaptation.time.narrowFactor", M_SQRT1_2);
      const double timeRelaxFactor = Dune::Fem::Parameter::getValue<double>(prefix + "adaptation.time.relaxFactor", M_SQRT2);
      const double timeRelaxTolerance =
        timeTolerance * Dune::Fem::Parameter::getValue<double>(prefix + "adaptation.time.relaxToleranceFraction", 0.3);
 
      // construct a file-name template in order to dump statistical data to disk
     int writeCount = 0;
      std::string dataPath = Dune::Fem::Parameter::commonOutputPath() + "/";
      dataPath += Dune::Fem::Parameter::getValue<std::string>("fem.io.path", "./") + "/";
 
      /*
       *
       ****************************************************************************
       *
       * Algorithm starts from here. First try to create a good initial
       * approximation, then run the adaptive space-time method which is
       * also implemented in the ALBERTA toolbox
       *
       */
 
      /******************** Compute an initial approximation *********************/
 
      double maxError;
      double maxEstimate;
      double estimate;
 
      Dune::Fem::persistenceManager << maxError;
      Dune::Fem::persistenceManager << maxEstimate;
 
      if (restart) {
        scheme.restart(/* checkpoint name taken from parameter file */);
      } else {
        scheme.initialize();
 
        // Compute an initial estimate
        estimate = scheme.initialEstimate();
 
        int cnt = 0;
        while (estimate > initialTolerance && ++cnt < initialMaximumIterations) {
 
          // mark element for adaptation
          scheme.initialMarking(initialTolerance);
 
          // actual refine or coarsen
          scheme.adapt();
 
          // Compute initial approximation
          scheme.initialize();
 
          estimate = scheme.initialEstimate();
        }
        assert(cnt < initialMaximumIterations);
 
        maxError = estimate; // actually the initial error
        maxEstimate = 0.0;
 
        scheme.output();
      }
 
      if (timeStrategy == explicitTimeStrategy) {
 
        /****************** variant with fixed time step size ********************/
 
        // ??? really ???
        //estimate = scheme.estimate();
        //maxEstimate = std::max(maxEstimate, estimate);
 
        // time loop
        while (timeProvider.time() < endTime) {
 
          if (Dune::Fem::Parameter::verbose()) {
            std::cout << "**** Starting time-step " << timeProvider.timeStep() << "@" << timeProvider.time() << std::endl;
          }
 
          if ( timeProvider.timeStep() > 0 ) {
            // Mark for refinement and coarsening
            scheme.mark(spaceTolerance);
 
            // adapt the mesh
            scheme.adapt();
          }
 
          // solve on the new mesh
          scheme.solve();
 
          // estimate the error
          estimate = scheme.estimate();
 
          int cnt = 0;
 
          while (estimate > spaceTolerance && cnt < spaceMaximumIterations ) {
            // Mark for refinement and coarsening
            scheme.mark(spaceTolerance);
 
            // adapt the mesh
            scheme.adapt();
 
            // solve on the new mesh
            scheme.solve();
 
            // estimate the error
            estimate = scheme.estimate();
 
            ++cnt;
          }
 
          // Switch to the next time step in order to be able to
          // compute the exact "error".
          scheme.next();
          timeProvider.next(timeStep);
 
          // calculate standard error
          double error = scheme.error();
 
          if (Dune::Fem::Parameter::verbose()) {
            std::cout << "Error: " << error << std::endl;
            std::cout << "Estimate: " << estimate << std::endl;
          }
 
          maxError    = std::max(error, maxError);
          maxEstimate = std::max(estimate + scheme.timeEstimate(), maxEstimate);
 
          // write fem data
          if ((writeCount = scheme.output()) >= 0) {
            // FIXME: should in principle only generated on one machine
            std::string name = Dune::Fem::generateFilename(dataPath + prefix_ + "-errorlog", writeCount);
            std::ofstream f;
            f.open(name, std::ios::out|std::ios::trunc);
 
            f << std::scientific << std::setprecision( 18 );
            f <<
              "#TimeStep "
              "Time                          "
              "DeltaT                        "
              "Estimate                      "
              "Error                         "
              "Dofs" << std::endl;
            f << std::setw(10) << std::left << timeProvider.timeStep() - 1 // already at next step
              << std::setw(30) << timeProvider.time() // already at next step
              << std::setw(30) << timeStep // old == new step size
              << std::setw(30) << estimate
              << std::setw(30) << error
              << std::setw(30) << scheme.size()
              << std::endl;
 
            f.close();
          }
        }
 
      } else { // implicitTimeStrategy
 
        /***************** variant with adaptive time step size ******************/
 
        // time loop
        while (timeProvider.time() < endTime) {
 
          double spaceEstimate;
          double timeEstimate;
 
          if (Dune::Fem::Parameter::verbose()) {
            std::cout << "**** Starting time-step " << timeProvider.timeStep() << "@" << timeProvider.time() << std::endl;
          }
 
          int cnt = 0;
          do {
 
            if (Dune::Fem::Parameter::verbose()) {
              std::cout << "Solving with step-size " << timeStep << std::endl;
            }
 
            // solve on the new mesh
            scheme.solve();
 
            // estimate the error
            spaceEstimate = scheme.estimate();
            timeEstimate  = scheme.timeEstimate();
 
            if (Dune::Fem::Parameter::verbose()) {
              std::cout << "spaceEstimate " << spaceEstimate << " (tolerance: " << spaceTolerance << ")" << std::endl
                        << "timeEstimate  " << timeEstimate << " (tolerance: " << timeTolerance << ")" << std::endl;
            }
 
            if (++cnt > maximumIterations) {
              // or terminate with error?
              break;
            }
 
            // First reduce the time-error
            if (timeEstimate > timeTolerance && timeStep > minDeltaT) {
              timeStep *= timeNarrowFactor;
              timeStep = std::max(minDeltaT, timeStep);
              timeProvider.init(timeStep);
              continue;
            }
 
            // Space adaptivity loop
            int spaceCnt = 0;
            do {
              // Mark for refinement and coarsening
              if (scheme.mark(spaceTolerance)) {
                scheme.adapt();
                scheme.solve();
                spaceEstimate = scheme.estimate();
                timeEstimate  = scheme.timeEstimate();
                if (timeEstimate > timeTolerance && timeStep > minDeltaT) {
                  if (Dune::Fem::Parameter::verbose()) {
                    std::cout << "Estimate: " << timeEstimate
                              << " (tolerance: " << timeTolerance << ")"
                              << std::endl;
                    std::cout << "timeEstimate > timeTolerance, restart"
                              << std::endl;
                  }
                  timeStep *= timeNarrowFactor;
                  timeStep = std::max(minDeltaT, timeStep);
                  timeProvider.init(timeStep);
                  break;
                }
              }
              if (++spaceCnt > spaceMaximumIterations) {
                break;
              }
            } while (spaceEstimate > spaceTolerance);
          } while (timeEstimate > timeTolerance);
 
          if (timeEstimate <= timeRelaxTolerance) {
            timeStep *= timeRelaxFactor;
            if (maxDeltaT > 0) {
              timeStep = std::min(maxDeltaT, timeStep);
            }
          }
 
          // Switch to the next time step in order to be able to
          // compute the exact "error".
 
          scheme.next();
 
          // remember for statistics
          double oldDeltaT = timeProvider.deltaT();
          timeProvider.next(timeStep);
 
          // calculate error
          double error = scheme.error();
 
          if (verbose) {
            size_t size = scheme.size();
            if (Dune::Fem::Parameter::verbose()) {
              std::cout << "Error         : " << error << std::endl;
              std::cout << "Space Estimate: " << spaceEstimate << std::endl;
              std::cout << "Time Estimate : " << timeEstimate << std::endl;
              std::cout << "#DoFs         : " << size << std::endl;
            }
          }
 
          maxError    = std::max(error, maxError);
          maxEstimate = std::max(spaceEstimate + timeEstimate, maxEstimate);
 
          // write fem data
          if ((writeCount = scheme.output()) >= 0) {
            // FIXME: should in principle only generated on one machine
            std::string name = Dune::Fem::generateFilename(dataPath + prefix_ + "-errorlog", writeCount);
            std::ofstream f;
            f.open(name, std::ios::out|std::ios::trunc);
 
            f << std::scientific << std::setprecision( 18 );
            f <<
              "#TimeStep "
              "Time                          "
              "DeltaT                        "
              "spaceEstimate                 "
              "timeEstimate                  "
              "Estimate                      "
              "Error                         "
              "Dofs" << std::endl;
            f << std::setw(10) << std::left << timeProvider.timeStep() - 1  // already at next step
              << std::setw(30) << timeProvider.time() // already at next step
              << std::setw(30) << oldDeltaT  // already at next step
              << std::setw(30) << spaceEstimate
              << std::setw(30) << timeEstimate
              << std::setw(30) << spaceEstimate + timeEstimate
              << std::setw(30) << error
              << std::setw(30) << scheme.size()
              << std::endl;
 
            f.close();
          }
 
        }
 
      }
 
      return maxEstimate;
    }
 
 
  } // ACFem
 
} // Dune
 
#endif //  __DUNE_ACFEM_ALGORITHMS_TRANSIENT_ADAPTIVE_ALGORITHM_HH__