SIMCoupledSI.h

Go to the documentation of this file.

// $Id$
//==============================================================================
//==============================================================================
 
#ifndef SIM_COUPLED_SI_H_
#define SIM_COUPLED_SI_H_
 
#include "SIMCoupled.h"
#include "SIMenums.h"
#include "MatVec.h"
#include "TimeStep.h"
#include "Functions.h"
#include "Utilities.h"
#include "IFEM.h"
#include <tinyxml2.h>
 
 
template<class T1, class T2>
class SIMCoupledSI : public SIMCoupled<T1,T2>
{
public:
  SIMCoupledSI(T1& s1, T2& s2) : SIMCoupled<T1,T2>(s1,s2)
  {
    maxIter = 50;
    maxItFunc = nullptr;
    noStg = aitken = false;
    omega = omega0 = 0.0;
  }
 
  virtual ~SIMCoupledSI() { delete maxItFunc; }
 
  virtual void enableStaggering(bool enable = true) { noStg = !enable; }
 
  virtual const Vector& getAitkenResidual() const
  {
    static Vector empty;
    return empty;
  }
 
  virtual const Vector& getRelaxationVector() const
  {
    static Vector empty;
    return empty;
  }
 
  virtual void setRelaxedSolution(const Vector&) {}
 
  virtual bool solveStep(TimeStep& tp, bool firstS1 = true)
  {
    if (tp.multiSteps())
      this->S1.getProcessAdm().cout <<"\n  step="<< tp.step
                                    <<"  time="<< tp.time.t << std::endl;
 
    int maxit = this->getMaxit(tp.step);
    SIM::ConvStatus conv = SIM::OK;
    for (tp.iter = 0; tp.iter <= maxit && conv != SIM::CONVERGED; tp.iter++)
    {
      SIM::ConvStatus status1 = SIM::OK, status2 = SIM::OK;
      if (firstS1 && (status1 = this->S1.solveIteration(tp)) <= SIM::DIVERGED)
        return false;
 
      if ((status2 = this->S2.solveIteration(tp)) <= SIM::DIVERGED)
        return false;
 
      if (!firstS1 && (status1 = this->S1.solveIteration(tp)) <= SIM::DIVERGED)
        return false;
 
      if ((conv = this->checkConvergence(tp,status1,status2)) <= SIM::DIVERGED)
        return false;
 
      if (omega0 != 0.0) {
        if (tp.iter > 0) {
          // Calculation of Aitken acceleration factor
          if (aitken) {
            Vector r1 = this->getAitkenResidual();
            r1 -= prevRes;
            omega *= -prevRes.dot(r1) / r1.dot(r1);
            if (fabs(omega) < 1e-6) {
              std::cerr <<"\n  ** Relaxation weight "<< omega <<" too small,"
                        <<" resetting to default "<< omega0 << std::endl;
              omega = omega0;
            }
          }
 
          // Perform relaxation
          IFEM::cout <<", omega="<< omega;
          prevSol *= 1.0 - omega;
          prevSol.add(this->getRelaxationVector(), omega);
          this->setRelaxedSolution(prevSol);
        } else {
          omega = omega0;
          prevSol = this->getRelaxationVector();
        }
        if (aitken)
          prevRes = this->getAitkenResidual();
      }
      IFEM::cout << std::endl;
    }
 
    this->S1.postSolve(tp);
    this->S2.postSolve(tp);
    tp.time.first = false;
 
    return true;
  }
 
  virtual SIM::ConvStatus checkConvergence(const TimeStep&,
                                           SIM::ConvStatus status1,
                                           SIM::ConvStatus status2)
  {
    if (status1 == status2)
      return status1;
 
    if (status1 == SIM::FAILURE || status2 == SIM::FAILURE)
      return SIM::FAILURE;
 
    if (status1 == SIM::DIVERGED || status2 == SIM::DIVERGED)
      return SIM::DIVERGED;
 
    return SIM::OK;
  }
 
protected:
  void parseIterations(const tinyxml2::XMLElement* elem)
  {
    IFEM::cout <<"\tUsing sub-iterations\n";
 
    std::string func;
    if (utl::getAttribute(elem,"maxFunc",func) ||
        utl::getAttribute(elem,"max",func))
      if (func.find_first_of('t') != std::string::npos)
      {
        IFEM::cout <<"\t\tmax = ";
        maxItFunc = utl::parseTimeFunc(func.c_str(),"expression");
      }
 
    if (!maxItFunc && utl::getAttribute(elem,"max",maxIter))
      IFEM::cout <<"\t\tmax = "<< maxIter << std::endl;
 
    if ((utl::getAttribute(elem,"relax",omega) ||
         utl::getAttribute(elem,"omega",omega)) && omega != 0.0) {
      IFEM::cout <<"\t\trelaxation = "<< omega;
      if (utl::getAttribute(elem,"aitken",aitken) && aitken)
        IFEM::cout <<" (aitken)";
      IFEM::cout << std::endl;
      omega0 = omega;
    }
  }
 
  int getMaxit(int iStep = 0) const
  {
    if (noStg)
      return 0; // staggering is disabled
    else if (maxItFunc)
      return static_cast<int>((*maxItFunc)(iStep));
 
    return maxIter;
  }
 
private:
  int maxIter; 
  ScalarFunc* maxItFunc; 
  bool noStg;  
 
  double omega0; 
  double omega;  
  bool   aitken; 
 
  Vector prevSol; 
  Vector prevRes; 
};
 
#endif