Added possibility for varying time step size

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1289 e10b68d5-8a6e-419e-a041-bce267b0401d

src/ASM/ASMs2D.h

@@ -38,7 +38,7 @@ class ASMs2D : public ASMstruct, public ASM2D
     int J; //!< Index in second parameter direction
   };
 
-  typedef std::vector<IJ> IndexVec; //! Node index container
+  typedef std::vector<IJ> IndexVec; //!< Node index container
 
   //! \brief Struct for edge node definitions.
   struct Edge

src/ASM/ASMs3D.h

@@ -38,7 +38,7 @@ class ASMs3D : public ASMstruct
     int K; //!< Index in third parameter direction
   };
 
-  typedef std::vector<IJK> IndexVec; //! Node index container
+  typedef std::vector<IJK> IndexVec; //!< Node index container
 
   //! \brief Struct for edge node definitions.
   struct Edge

src/ASM/SAMpatchPara.h

@@ -29,7 +29,7 @@ class SAMpatchPara : public SAMpatch
 {
 public:
   //! \brief The constructor initializes the \a l2gn array.
-  //! \param[in] Global-to-local node numbers for this processor
+  //! \param[in] g2ln Global-to-local node numbers for this processor
   SAMpatchPara(const std::map<int,int>& g2ln);
   //! \brief The destructor destroys the index set arrays.
   virtual ~SAMpatchPara();

src/SIM/NonLinSIM.C

@@ -48,18 +48,31 @@ bool NonLinSIM::parse (char* keyWord, std::istream& is)
 {
   if (!strncasecmp(keyWord,"TIME_STEPPING",13))
   {
-    double dt;
+    int nstep = atoi(keyWord+13);
-    tInc.clear();
+    if (nstep < 1) nstep = 1;
-    tInc.reserve(5);
-    std::istringstream cline(utl::readLine(is));
 
-    cline >> startTime >> stopTime >> dt;
+    double dt;
-    if (cline.fail() || cline.bad()) return false;
+    steps.resize(nstep);
-    while (!cline.fail() && !cline.bad())
+    for (int i = 0; i < nstep; i++)
     {
-      tInc.push_back(dt);
+      std::istringstream cline(utl::readLine(is));
-      cline >> dt;
+      if (i == 0) cline >> startTime;
+      cline >> steps[i].second >> dt;
+      if (cline.fail() || cline.bad()) return false;
+      if (dt > 1.0 && ceil(dt) == dt)
+      {
+	// The number of time steps are specified
+	dt = (steps[i].second - (i == 0 ? startTime : steps[i-1].second))/dt;
+	steps[i].first.push_back(dt);
+      }
+      else while (!cline.fail() && !cline.bad())
+      {
+	// The time step size(s) is/are specified
+	steps[i].first.push_back(dt);
+	cline >> dt;
+      }
     }
+    stopTime = steps.back().second;
   }
   else if (!strncasecmp(keyWord,"NONLINEAR_SOLVER",16))
   {
@@ -100,11 +113,7 @@ void NonLinSIM::initSystem (SystemMatrix::Type mType, size_t nGauss)
 
 void NonLinSIM::init (SolvePrm& param, const RealArray& initVal)
 {
-  param.startTime = startTime;
+  param.initTime(startTime,stopTime,steps);
-  param.stopTime = stopTime;
-  param.tInc    = tInc;
-  param.time.t  = startTime;
-  param.time.dt = tInc.empty() ? stopTime-startTime : tInc.front();
   param.maxit   = maxit;
   param.nupdat  = nupdat;
   param.convTol = convTol;
@@ -140,8 +149,14 @@ bool NonLinSIM::solveStep (SolvePrm& param, SIM::SolutionMode mode,
   PROFILE1("NonLinSIM::solveStep");
 
   if (msgLevel >= 0 && myPid == 0)
+  {
+    std::streamsize oldPrec = 0;
+    double digits = log10(param.time.t)-log10(param.time.dt);
+    if (digits > 6.0) oldPrec = std::cout.precision(ceil(digits));
     std::cout <<"\n  step="<< param.step
 	      <<"  time="<< param.time.t << std::endl;
+    if (oldPrec > 0) std::cout.precision(oldPrec);
+  }
 
   param.iter = 0;
   param.alpha = 1.0;

src/SIM/NonLinSIM.h

@@ -181,7 +181,8 @@ private:
   // Nonlinear solution algorithm parameters
   double startTime; //!< Start time of the simulation
   double stopTime;  //!< Stop time of the simulation
-  RealArray tInc;   //!< Time increment size(s)
+  SIMparameters::TimeSteps steps; //!< Time increment specifications
+
   double convTol;   //!< Relative convergence tolerance
   double divgLim;   //!< Relative divergence limit
   double eta;       //!< Line search tolerance

src/SIM/SIMparameters.C

@@ -14,36 +14,63 @@
 #include "SIMparameters.h"
 
 
+void SIMparameters::initTime (double start, double stop, const TimeSteps& steps)
+{
+  starTime = start;
+  stopTime = stop;
+  mySteps  = steps;
+  time.t   = start;
+  time.dt  = steps.empty() ? stop-start : steps.front().first.front();
+  stepIt   = mySteps.begin();
+}
+
+
+static const double epsT = 1.0e-6; //!< Tolerance parameter
+
+
 bool SIMparameters::multiSteps () const
 {
-  if (tInc.empty()) return false;
+  if (mySteps.empty()) return false;
 
-  const double epsT = 1.0e-6;
+  return (starTime+(1.0+epsT)*mySteps.front().first.front() < stopTime);
-  return (startTime+(1.0+epsT)*tInc.front() < stopTime);
 }
 
 
 bool SIMparameters::increment ()
 {
-  const double epsT = 1.0e-6;
+  if (stepIt != mySteps.end())
-
+    if (++lstep <= stepIt->first.size())
-  if (++step <= (int)tInc.size() && step > 0)
+      time.dt = stepIt->first[lstep-1];
-    time.dt = tInc[step-1];
 
   if (time.t+time.dt*epsT >= stopTime)
-    return false;
+    return false; // We've reached the end of the simulation
-
-  if (tInc.size() <= (size_t)step)
-    tInc.push_back(time.dt);
 
+  ++step;
   time.t += time.dt;
 
+  if (stepIt != mySteps.end())
+  {
+    if (stepIt->first.size() <= lstep)
+      stepIt->first.push_back(time.dt);
+    if (time.t+time.dt*epsT >= stepIt->second)
+    {
+      if (time.t != stepIt->second)
+      {
+	// Adjust the size of the last time step
+	time.dt += stepIt->second - time.t;
+	time.t = stepIt->second;
+	stepIt->first.back() = time.dt;
+      }
+      lstep = 0;
+      ++stepIt;
+    }
+  }
+
   if (time.t > stopTime)
   {
     // Adjust the size of the last time step
     time.dt += stopTime - time.t;
     time.t = stopTime;
-    tInc.back() = time.dt;
   }
 
   return true;

src/SIM/SIMparameters.h

@@ -16,6 +16,7 @@
 
 #include "TimeDomain.h"
 #include <vector>
+#include <cstddef>
 
 
 /*!
@@ -26,11 +27,17 @@ class SIMparameters
 {
 public:
   //! \brief The constructor initializes the counters to zero.
-  SIMparameters() : step(0), iter(time.it) {}
+  SIMparameters() : step(0), iter(time.it), lstep(0) { stepIt = mySteps.end(); }
 
   //! \brief Empty destructor.
   virtual ~SIMparameters() {}
 
+  //! \brief Time stepping definition container.
+  typedef std::vector< std::pair<std::vector<double>,double> > TimeSteps;
+
+  //! \brief Initializes the time step definitions.
+  void initTime(double start, double stop, const TimeSteps& steps);
+
   //! \brief Returns \e true if the simulation consists of several load steps.
   bool multiSteps() const;
 
@@ -38,13 +45,18 @@ public:
   //! \return \e true, if we have reached the end of the simulation.
   bool increment();
 
-  int    step;      //!< Load/time step counter
+  int        step; //!< Load/time step counter
-  int&   iter;      //!< Iteration counter
+  int&       iter; //!< Iteration counter
-  double startTime; //!< Start (pseudo)time of simulation
+  TimeDomain time; //!< Time domain data
-  double stopTime;  //!< Stop (pseudo)time of simulation
 
-  std::vector<double> tInc; //!< Time (or pseudo time) increment size(s)
+  double starTime; //!< Start (pseudo)time of simulation
-  TimeDomain          time; //!< Time domain data
+  double stopTime; //!< Stop (pseudo)time of simulation
+
+private:
+  size_t lstep; //!< Local step counter, i.e., within current \a *stepIt
+
+  TimeSteps           mySteps; //!< Time step definitions
+  TimeSteps::iterator stepIt;  //!< Running iterator over the time steps
 };
 
 #endif