From 24cf8b3f9ede524e48bafb22b9aa723f62a02385 Mon Sep 17 00:00:00 2001 From: kmo Date: Sun, 13 Nov 2011 16:38:17 +0000 Subject: [PATCH] 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 | 2 +- src/ASM/ASMs3D.h | 2 +- src/ASM/SAMpatchPara.h | 2 +- src/SIM/NonLinSIM.C | 43 +++++++++++++++++++++++----------- src/SIM/NonLinSIM.h | 3 ++- src/SIM/SIMparameters.C | 51 +++++++++++++++++++++++++++++++---------- src/SIM/SIMparameters.h | 26 +++++++++++++++------ 7 files changed, 92 insertions(+), 37 deletions(-) diff --git a/src/ASM/ASMs2D.h b/src/ASM/ASMs2D.h index e93a68da..091c92fa 100644 --- a/src/ASM/ASMs2D.h +++ b/src/ASM/ASMs2D.h @@ -38,7 +38,7 @@ class ASMs2D : public ASMstruct, public ASM2D int J; //!< Index in second parameter direction }; - typedef std::vector IndexVec; //! Node index container + typedef std::vector IndexVec; //!< Node index container //! \brief Struct for edge node definitions. struct Edge diff --git a/src/ASM/ASMs3D.h b/src/ASM/ASMs3D.h index a226cd5b..a0a5aebc 100644 --- a/src/ASM/ASMs3D.h +++ b/src/ASM/ASMs3D.h @@ -38,7 +38,7 @@ class ASMs3D : public ASMstruct int K; //!< Index in third parameter direction }; - typedef std::vector IndexVec; //! Node index container + typedef std::vector IndexVec; //!< Node index container //! \brief Struct for edge node definitions. struct Edge diff --git a/src/ASM/SAMpatchPara.h b/src/ASM/SAMpatchPara.h index 1bf6ff64..977a37e2 100644 --- a/src/ASM/SAMpatchPara.h +++ b/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& g2ln); //! \brief The destructor destroys the index set arrays. virtual ~SAMpatchPara(); diff --git a/src/SIM/NonLinSIM.C b/src/SIM/NonLinSIM.C index 69b0c64e..107c0195 100644 --- a/src/SIM/NonLinSIM.C +++ b/src/SIM/NonLinSIM.C @@ -48,18 +48,31 @@ bool NonLinSIM::parse (char* keyWord, std::istream& is) { if (!strncasecmp(keyWord,"TIME_STEPPING",13)) { - double dt; - tInc.clear(); - tInc.reserve(5); - std::istringstream cline(utl::readLine(is)); + int nstep = atoi(keyWord+13); + if (nstep < 1) nstep = 1; - cline >> startTime >> stopTime >> dt; - if (cline.fail() || cline.bad()) return false; - while (!cline.fail() && !cline.bad()) + double dt; + steps.resize(nstep); + for (int i = 0; i < nstep; i++) { - tInc.push_back(dt); - cline >> dt; + std::istringstream cline(utl::readLine(is)); + 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.stopTime = stopTime; - param.tInc = tInc; - param.time.t = startTime; - param.time.dt = tInc.empty() ? stopTime-startTime : tInc.front(); + param.initTime(startTime,stopTime,steps); 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; diff --git a/src/SIM/NonLinSIM.h b/src/SIM/NonLinSIM.h index fbb0f4bd..c19f7bdf 100644 --- a/src/SIM/NonLinSIM.h +++ b/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 diff --git a/src/SIM/SIMparameters.C b/src/SIM/SIMparameters.C index 7c5643b3..00e2055f 100644 --- a/src/SIM/SIMparameters.C +++ b/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 (startTime+(1.0+epsT)*tInc.front() < stopTime); + return (starTime+(1.0+epsT)*mySteps.front().first.front() < stopTime); } bool SIMparameters::increment () { - const double epsT = 1.0e-6; - - if (++step <= (int)tInc.size() && step > 0) - time.dt = tInc[step-1]; + if (stepIt != mySteps.end()) + if (++lstep <= stepIt->first.size()) + time.dt = stepIt->first[lstep-1]; if (time.t+time.dt*epsT >= stopTime) - return false; - - if (tInc.size() <= (size_t)step) - tInc.push_back(time.dt); + return false; // We've reached the end of the simulation + ++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; diff --git a/src/SIM/SIMparameters.h b/src/SIM/SIMparameters.h index 888bc549..d5ad56e6 100644 --- a/src/SIM/SIMparameters.h +++ b/src/SIM/SIMparameters.h @@ -16,6 +16,7 @@ #include "TimeDomain.h" #include +#include /*! @@ -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,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& iter; //!< Iteration counter - double startTime; //!< Start (pseudo)time of simulation - double stopTime; //!< Stop (pseudo)time of simulation + int step; //!< Load/time step counter + int& iter; //!< Iteration counter + TimeDomain time; //!< Time domain data - std::vector tInc; //!< Time (or pseudo time) increment size(s) - TimeDomain time; //!< Time domain data + double starTime; //!< Start (pseudo)time of simulation + 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