commit 74baa0d17070e41260cd2c26ed16ed2e47ba1457
parent 7e780512a4d3f7ab709394ffcaa96d61434bba24
Author: Fabian Wermelinger <fabianw@mavt.ethz.ch>
Date: Sun, 25 Sep 2016 13:45:19 +0200
removed mixed precision computations in variable step integrators
Diffstat:
6 files changed, 81 insertions(+), 67 deletions(-)
diff --git a/GnuplotDump.h b/GnuplotDump.h
@@ -26,7 +26,7 @@ private:
std::ofstream m_file;
size_t m_N;
- void _writeASCII(const size_t step, const double t, const double dt, const Real * const src, const size_t N)
+ void _writeASCII(const size_t step, const Real t, const Real dt, const Real * const src, const size_t N)
{
if (m_file.good())
{
@@ -96,7 +96,7 @@ public:
enum {ASCII, BIN};
- void write(const size_t step, const double t, const double dt, const Real * const src, const size_t N)
+ void write(const size_t step, const Real t, const Real dt, const Real * const src, const size_t N)
{
if (m_ftype == ASCII)
_writeASCII(step, t, dt, src, N);
diff --git a/TimeStepper/StepperBase.h b/TimeStepper/StepperBase.h
@@ -54,20 +54,20 @@ struct StepperSettings
bFixedStep = true;
}
- double aTol;
- double rTol;
- double minScale;
- double maxScale;
- double alpha;
- double beta;
- double safety;
- double t;
- double dt;
+ Real aTol;
+ Real rTol;
+ Real minScale;
+ Real maxScale;
+ Real alpha;
+ Real beta;
+ Real safety;
+ Real t;
+ Real dt;
size_t step;
size_t nsteps;
- double tFinal;
- double dtDump;
- double tDump;
+ Real tFinal;
+ Real dtDump;
+ Real tDump;
size_t writeGranularity;
size_t writeCount;
size_t reportGranularity;
diff --git a/TimeStepper/explicit/variableStep/RKF45.h b/TimeStepper/explicit/variableStep/RKF45.h
@@ -20,9 +20,9 @@ class RKF45 : public Euler<Tinput, Trhs>
using Euler<Tinput, Trhs>::m_rhs;
using Euler<Tinput, Trhs>::m_kernel;
- double m_errold;
- const double m_alpha;
- const double m_beta;
+ Real m_errold;
+ const Real m_alpha;
+ const Real m_beta;
bool m_reject;
bool m_firstPass;
@@ -52,18 +52,18 @@ class RKF45 : public Euler<Tinput, Trhs>
m_kernel->compute(m_U + dt*(c61*m_rhs + c62*m_rhs2 + c63*m_rhs3 + c64*m_rhs4 + c65*m_rhs5), m_rhs6, t + ct6*dt, data);
}
- inline double _error() const
+ inline Real _error() const
{
assert(m_U.size() == m_rhs.size());
- double maxerr = 0.0;
+ Real maxerr = 0.0;
// using inf-norm (could also do L2-norm)
for (size_t i = 0; i < m_U.size(); ++i)
{
const typename Trhs::DataType E = e1*m_rhs[i] + e3*m_rhs3[i] + e4*m_rhs4[i] + e5*m_rhs5[i] + e6*m_rhs6[i];
for (int j = 0; j < Trhs::DataType::SIZE; ++j)
{
- const double IepsI = std::abs(E[j]);
- const double scale = m_settings.aTol + std::abs(m_U[i][j])*m_settings.rTol;
+ const Real IepsI = std::abs(E[j]);
+ const Real scale = m_settings.aTol + std::abs(m_U[i][j])*m_settings.rTol;
maxerr = std::max(maxerr, IepsI/scale);
}
}
@@ -73,19 +73,19 @@ class RKF45 : public Euler<Tinput, Trhs>
bool _inBound(const Real t, Real& dt, void const* const data, Real& dt_next)
{
- const double err = _error();
+ const Real err = _error();
if (dt < std::numeric_limits<Real>::epsilon())
{
// if you request the impossible
dt_next = 1.5*std::numeric_limits<Real>::epsilon();
- m_errold = std::max(err, 1.0e-4);
+ m_errold = std::max(err, static_cast<Real>(1.0e-4));
m_reject = false;
m_firstPass = true;
return true;
}
- const double safety = m_settings.safety;
+ const Real safety = m_settings.safety;
const Real minScale = m_settings.minScale;
const Real maxScale = m_settings.maxScale;
Real scale;
@@ -102,7 +102,7 @@ class RKF45 : public Euler<Tinput, Trhs>
dt_next = dt*std::min(scale, static_cast<Real>(1.0));
else
dt_next = dt*scale;
- m_errold = std::max(err, 1.0e-4);
+ m_errold = std::max(err, static_cast<Real>(1.0e-4));
m_reject = false;
m_firstPass = true;
return true;
@@ -121,9 +121,9 @@ class RKF45 : public Euler<Tinput, Trhs>
void _dumpOutput(const void * const data)
{
- const double t = m_settings.t;
- const double dt = m_settings.dt;
- double& tDump = m_settings.tDump;
+ const Real t = m_settings.t;
+ const Real dt = m_settings.dt;
+ Real& tDump = m_settings.tDump;
const Tinput Uold = m_output;
const Trhs rhsOld = m_rhs;
@@ -133,7 +133,7 @@ class RKF45 : public Euler<Tinput, Trhs>
while(tDump <= (t+dt))
{
- const double phi = (tDump - t)/dt;
+ const Real phi = (tDump - t)/dt;
// Hermite 3rd order
m_output = (1.0-phi)*Uold + phi*m_U + phi*(phi-1.0)*((1.0-2.0*phi)*(m_U - Uold) + (phi-1.0)*dt*rhsOld + phi*dt*m_rhs);
m_kernel->write(m_settings.step, t+phi*dt, phi*dt, m_output, data);
@@ -150,10 +150,10 @@ public:
virtual void step(void const* const data=nullptr)
{
size_t& step = m_settings.step;
- double& t = m_settings.t;
- double& dt = m_settings.dt;
- const double& tFinal = m_settings.tFinal;
- double& tDump = m_settings.tDump;
+ Real& t = m_settings.t;
+ Real& dt = m_settings.dt;
+ const Real& tFinal = m_settings.tFinal;
+ Real& tDump = m_settings.tDump;
if (m_settings.dtDump > 0.0) tDump = t + m_settings.dtDump;
while(t < tFinal)
diff --git a/TimeStepper/explicit/variableStep/RKV56.h b/TimeStepper/explicit/variableStep/RKV56.h
@@ -21,9 +21,9 @@ class RKV56 : public Euler<Tinput, Trhs>
using Euler<Tinput, Trhs>::m_rhs;
using Euler<Tinput, Trhs>::m_kernel;
- double m_errold;
- const double m_alpha;
- const double m_beta;
+ Real m_errold;
+ const Real m_alpha;
+ const Real m_beta;
bool m_reject;
bool m_firstPass;
@@ -59,18 +59,18 @@ class RKV56 : public Euler<Tinput, Trhs>
m_kernel->compute(m_U + dt*(c81*m_rhs + c82*m_rhs2 + c83*m_rhs3 + c84*m_rhs4 + c85*m_rhs5 + c87*m_rhs7), m_rhs8, t + ct8*dt, data);
}
- inline double _error() const
+ inline Real _error() const
{
assert(m_U.size() == m_rhs.size());
- double maxerr = 0.0;
+ Real maxerr = 0.0;
// using inf-norm (could also do L2-norm)
for (size_t i = 0; i < m_U.size(); ++i)
{
const typename Trhs::DataType E = e1*m_rhs[i] + e3*m_rhs3[i] + e4*m_rhs4[i] + e5*m_rhs5[i] + e6*m_rhs6[i] + e7*m_rhs7[i] + e8*m_rhs8[i];
for (int j = 0; j < Trhs::DataType::SIZE; ++j)
{
- const double IepsI = std::abs(E[j]);
- const double scale = m_settings.aTol + std::abs(m_U[i][j])*m_settings.rTol;
+ const Real IepsI = std::abs(E[j]);
+ const Real scale = m_settings.aTol + std::abs(m_U[i][j])*m_settings.rTol;
maxerr = std::max(maxerr, IepsI/scale);
}
}
@@ -80,19 +80,19 @@ class RKV56 : public Euler<Tinput, Trhs>
bool _inBound(const Real t, Real& dt, void const* const data, Real& dt_next)
{
- const double err = _error();
+ const Real err = _error();
if (dt < std::numeric_limits<Real>::epsilon())
{
// if you request the impossible
dt_next = 1.5*std::numeric_limits<Real>::epsilon();
- m_errold = std::max(err, 1.0e-4);
+ m_errold = std::max(err, static_cast<Real>(1.0e-4));
m_reject = false;
m_firstPass = true;
return true;
}
- const double safety = m_settings.safety;
+ const Real safety = m_settings.safety;
const Real minScale = m_settings.minScale;
const Real maxScale = m_settings.maxScale;
Real scale;
@@ -109,7 +109,7 @@ class RKV56 : public Euler<Tinput, Trhs>
dt_next = dt*std::min(scale, static_cast<Real>(1.0));
else
dt_next = dt*scale;
- m_errold = std::max(err, 1.0e-4);
+ m_errold = std::max(err, static_cast<Real>(1.0e-4));
m_reject = false;
m_firstPass = true;
return true;
@@ -128,9 +128,9 @@ class RKV56 : public Euler<Tinput, Trhs>
void _dumpOutput(const void * const data)
{
- const double t = m_settings.t;
- const double dt = m_settings.dt;
- double& tDump = m_settings.tDump;
+ const Real t = m_settings.t;
+ const Real dt = m_settings.dt;
+ Real& tDump = m_settings.tDump;
const Tinput Uold = m_output;
const Trhs rhsOld = m_rhs;
@@ -140,7 +140,7 @@ class RKV56 : public Euler<Tinput, Trhs>
while(tDump <= (t+dt))
{
- const double phi = (tDump - t)/dt;
+ const Real phi = (tDump - t)/dt;
// Hermite 3rd order
m_output = (1.0-phi)*Uold + phi*m_U + phi*(phi-1.0)*((1.0-2.0*phi)*(m_U - Uold) + (phi-1.0)*dt*rhsOld + phi*dt*m_rhs);
m_kernel->write(m_settings.step, t+phi*dt, phi*dt, m_output, data);
@@ -157,10 +157,10 @@ public:
virtual void step(void const* const data=nullptr)
{
size_t& step = m_settings.step;
- double& t = m_settings.t;
- double& dt = m_settings.dt;
- const double& tFinal = m_settings.tFinal;
- double& tDump = m_settings.tDump;
+ Real& t = m_settings.t;
+ Real& dt = m_settings.dt;
+ const Real& tFinal = m_settings.tFinal;
+ Real& tDump = m_settings.tDump;
if (m_settings.dtDump > 0.0) tDump = t + m_settings.dtDump;
while(t < tFinal)
diff --git a/TimeStepper/orderVerification/orderVerification.cpp b/TimeStepper/orderVerification/orderVerification.cpp
@@ -18,7 +18,7 @@
using namespace std;
-constexpr double DECAY = -1.0;
+constexpr Real DECAY = -1.0;
template <typename Tinput, typename Trhs>
class Dahlquist : public KernelBase<Tinput,Trhs>
@@ -67,12 +67,12 @@ using vec_t = StateVector<State1>;
using baseKernel_t = KernelBase<vec_t,vec_t>*;
using baseStepper_t = StepperBase<vec_t,vec_t>*;
-inline double exact(const double t, const double c) { return exp(c*t); }
+inline Real exact(const Real t, const Real c) { return exp(c*t); }
struct Error
{
Error() : dt(0), L1(0), L2(0), Linf(0) {}
- double dt, L1, L2, Linf;
+ Real dt, L1, L2, Linf;
};
vector<Error> verify(baseStepper_t stepper, int const stages, StepperSettings& S, const Real dt=1.0, const Real scale=0.5)
@@ -88,7 +88,7 @@ vector<Error> verify(baseStepper_t stepper, int const stages, StepperSettings& S
for (size_t i = 0; i < N; ++i)
{
stepper->step();
- const double absErr = abs(numerical - exact(S.t, DECAY));
+ const Real absErr = abs(numerical - exact(S.t, DECAY));
error[s].L1 += absErr;
error[s].L2 += absErr*absErr;
error[s].Linf = (error[s].Linf < absErr ? absErr : error[s].Linf);
@@ -114,7 +114,7 @@ int main(int argc, const char** argv)
ArgumentParser parser(argc, argv);
StepperSettings S(parser);
- vector<pair<pair<string,double>, baseStepper_t> > timeStepperList;
+ vector<pair<pair<string,Real>, baseStepper_t> > timeStepperList;
timeStepperList.push_back(make_pair(make_pair("Euler",1.0), new Euler<vec_t, vec_t>(U, S, dahlquist)));
timeStepperList.push_back(make_pair(make_pair("RK4",4.0), new RK4<vec_t, vec_t>(U, S, dahlquist)));
timeStepperList.push_back(make_pair(make_pair("LSRK3",3.0), new LSRK3<vec_t, vec_t>(U, S, dahlquistLSRK3)));
@@ -129,7 +129,7 @@ int main(int argc, const char** argv)
vector<Error> error = verify(timeStepperList[i].second, 7, S);
ofstream save(timeStepperList[i].first.first + ".dat");
save.setf(std::ios::scientific, std::ios::floatfield);
- auto expected = [&](const double x){ return error[0].L1*pow(x, timeStepperList[i].first.second); };
+ auto expected = [&](const Real x){ return error[0].L1*pow(x, timeStepperList[i].first.second); };
for (size_t j = 0; j < error.size(); ++j)
save << error[j].dt << '\t' << 1.0/error[j].dt << '\t' << error[j].L1 << '\t' << error[j].L2 << '\t' << error[j].Linf << '\t' << expected(error[j].dt) << endl;
save.close();
diff --git a/common.h b/common.h
@@ -53,7 +53,7 @@ struct State
s[i] = rhs.s[i];
return *this;
}
- inline State& operator=(const double c)
+ inline State& operator=(const Real c)
{
for (int i = 0; i < _SIZE; ++i)
s[i] = c;
@@ -133,17 +133,29 @@ typedef State<Real,2> State2;
template <typename T, int _SS=0, int _SE=0>
class LightVector
{
- const size_t _N;
+ bool _bAllocated;
+ size_t _N;
T* _data;
inline T* _alloc(const size_t n)
{
- void* pmem;
- posix_memalign(&pmem, _ALIGN_, n*sizeof(T));
- return (T*)pmem;
+ if (!_data)
+ {
+ void* pmem;
+ posix_memalign(&pmem, _ALIGN_, n*sizeof(T));
+ _bAllocated = true;
+ return (T*)pmem;
+ }
+ else
+ return _data;
}
- inline void _dealloc() { free(_data); _data=0; }
+ inline void _dealloc()
+ {
+ if (_data) free(_data);
+ _bAllocated = false;
+ _data = nullptr;
+ }
inline void _copy(const T* const src)
{
@@ -160,8 +172,9 @@ class LightVector
}
public:
- LightVector(const int n) : _N(n), _data(_alloc(n+(_SE-_SS))) { _clear(); }
- LightVector(const LightVector& rhs) : _N(rhs._N), _data(_alloc(rhs._N+(_SE-_SS))) { _copy(rhs._data); }
+ LightVector() : _bAllocated(false), _N(0), _data(nullptr) {}
+ LightVector(const int n) : _bAllocated(false), _N(n), _data(_alloc(n+(_SE-_SS))) { _clear(); }
+ LightVector(const LightVector& rhs) : _bAllocated(false), _N(rhs._N), _data(_alloc(rhs._N+(_SE-_SS))) { _copy(rhs._data); }
virtual ~LightVector() { _dealloc(); }
static const int SS = _SS;
@@ -202,6 +215,7 @@ public:
return *this;
}
+ inline void allocate(const size_t n) { _N=n; _alloc(n+(_SE-_SS)); _clear(); }
inline size_t size() const { return _N; }
inline size_t size_all() const { return _N+_SE-_SS; }
inline void clear() { _clear(); }
@@ -255,7 +269,7 @@ public:
if (this != &rhs) _copy(rhs);
return *this;
}
- StateVector& operator=(double const rhs)
+ StateVector& operator=(Real const rhs)
{
for (size_t i = 0; i < m_state.size(); ++i)
m_state[i] = rhs;
