commit f3947228c46bcd49b18f7e05a85a9a3767086d24
parent 1c50b530ffe3e5dc4e75aa95e0b3e466408f3f17
Author: Fabian Wermelinger <fabianw@mavt.ethz.ch>
Date: Fri, 18 Dec 2020 20:38:44 +0100
Rename includes
Diffstat:
25 files changed, 1130 insertions(+), 1116 deletions(-)
diff --git a/GnuplotDump.h b/GnuplotDump.h
@@ -1,110 +0,0 @@
-/* File: GnuplotDump.h */
-/* Date: Tue Mar 1 21:07:49 2016 */
-/* Author: Fabian Wermelinger */
-/* Tag: Gnuplot Dumper */
-/* Copyright 2016 ETH Zurich. All Rights Reserved. */
-#ifndef GNUPLOTDUMP_H_2VNNTIUB
-#define GNUPLOTDUMP_H_2VNNTIUB
-
-#ifdef _FLOAT_PRECISION_
-using Real = float;
-#else
-using Real = double;
-#endif
-
-#ifdef _GNUPLOT_DOUBLE_
-using GPfloat = double;
-#else
-using GPfloat = float;
-#endif
-
-class GnuplotDump
-{
-private:
- const int m_ftype;
- std::string m_fname;
- std::ofstream m_file;
- size_t m_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())
- {
- m_file << step << '\t' << t << '\t' << dt;
- for (size_t i = 0; i < N; ++i)
- m_file << '\t' << src[i];
- m_file << std::endl;
- m_N = N+3;
- }
- }
-
- void _writeBIN(const size_t step, const GPfloat t, const GPfloat dt, const Real * const src, size_t N)
- {
- if (1024 < N)
- {
- std::cerr << "GnuplotDump: WARNING -- Max. Column width for BIN dump is 1024" << std::endl;
- N = 1024;
- }
- GPfloat buf[1024+3];
-
- buf[0] = static_cast<GPfloat>(step);
- buf[1] = t;
- buf[2] = dt;
- for (size_t i = 0; i < N; ++i)
- buf[i+3] = static_cast<GPfloat>(src[i]);
-
- if (m_file.good())
- {
- m_file.write((const char*)(&buf[0]), (N+3)*sizeof(GPfloat));
- m_N = N+3;
- }
- }
-
- void _writeGPscript() const
- {
- std::ofstream out(m_fname+".gp");
- if (m_ftype == ASCII)
- out << "plot '" << m_fname+".dat" << "' using 1:2 with lines" << std::endl;
- else
- {
- out << "plot '" << m_fname+".bin" << "' binary format='\%" << m_N;
- if (sizeof(GPfloat) == 4)
- out << "float";
- else
- out << "double";
- out << "' using 1:2 with lines" << std::endl;
- }
- out.close();
- }
-
-public:
- GnuplotDump(const std::string fname="gnuplot", const int type=BIN, const size_t pos=0) : m_ftype(type), m_fname(fname)
- {
- if (m_ftype == ASCII)
- {
- m_file.open(m_fname+".dat", std::ios_base::out|std::ios_base::app);
- m_file.setf(std::ios::scientific, std::ios::floatfield);
- }
- else
- m_file.open(m_fname+".bin", std::ios_base::out|std::ios_base::binary|std::ios_base::app);
-
- if (pos > 0) m_file.seekp(pos);
- }
- ~GnuplotDump()
- {
- m_file.close();
- _writeGPscript();
- }
-
- enum {ASCII, BIN};
-
- 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);
- else
- _writeBIN(step, t, dt, src, N);
- }
-};
-
-#endif /* GNUPLOTDUMP_H_2VNNTIUB */
diff --git a/TimeStepper/StepperBase.h b/TimeStepper/StepperBase.h
@@ -1,235 +0,0 @@
-// File : StepperBase.h
-// Date : Thu Sep 22 13:27:19 2016
-// Author : Fabian Wermelinger
-// Description: Time stepper base class
-// Copyright 2016 ETH Zurich. All Rights Reserved.
-#ifndef STEPPERBASE_H_FCCSM4HL
-#define STEPPERBASE_H_FCCSM4HL
-
-#include "KernelBase.h"
-
-// forward declarations
-template <typename Tinput, typename Trhs=Tinput>
-class StepperBase;
-template <typename Tinput, typename Trhs=Tinput>
-class Euler;
-template <typename Tinput, typename Trhs=Tinput>
-class RK4;
-template <typename Tinput, typename Trhs=Tinput>
-class LSRK3;
-template <typename Tinput, typename Trhs=Tinput>
-class RKF45;
-template <typename Tinput, typename Trhs=Tinput>
-class RKV56;
-#ifdef _USE_SUNDIALS_
-template <typename Tinput, typename Trhs=Tinput>
-class BDF;
-#endif /* _USE_SUNDIALS_ */
-
-#ifdef _FLOAT_PRECISION_
-using Real = float;
-#else
-using Real = double;
-#endif
-
-struct StepperSettings
-{
- StepperSettings(ArgumentParser& p)
- {
- aTol = p("-ts_aTol").asDouble(1.0e-8);
- rTol = p("-ts_rTol").asDouble(1.0e-8);
- minScale = p("-ts_minScale").asDouble(0.2);
- maxScale = p("-ts_maxScale").asDouble(10.0);
- alpha = p("-ts_alpha").asDouble(1.0);
- beta = p("-ts_beta").asDouble(0.0);
- safety = p("-ts_safety").asDouble(0.9);
- t = p("-ts_t0").asDouble(0.0);
- dt = p("-ts_dt").asDouble(1.0e-4);
- step = 0;
- nsteps = p("-ts_nsteps").asInt(0);
- p.set_strict_mode();
- tFinal = p("-ts_tend").asDouble();
- p.unset_strict_mode();
- dtDump = p("-ts_dtDump").asDouble(-1.0);
- tDump = t;
- writeGranularity = p("-ts_writeGranularity").asInt(1);
- writeCount = 0;
- reportGranularity= p("-ts_reportGranularity").asInt(100);
- bFixedStep = true;
-
- // restarts
- bRestart = p("-ts_restart").asBool(false);
- restartstep = p("-ts_restart_step").asInt(0);
-
- }
-
- Real aTol;
- Real rTol;
- Real minScale;
- Real maxScale;
- Real alpha;
- Real beta;
- Real safety;
- Real t;
- Real dt;
- size_t step;
- size_t nsteps;
- Real tFinal;
- Real dtDump;
- Real tDump;
- size_t writeGranularity;
- size_t writeCount;
- size_t reportGranularity;
- bool bFixedStep;
-
- // restarts
- bool bRestart;
- int restartstep;
-
-
- // helper
- void print() const
- {
- printf("Time Stepper :\n");
- printf("\tAbsolute tolerance = %e\n", aTol);
- printf("\tRelative tolerance = %e\n", rTol);
- printf("\tMinimum time step scale = %e\n", minScale);
- printf("\tMaximum time step scale = %e\n", maxScale);
- printf("\tError PI control alpha = %e\n", alpha);
- printf("\tError PI control beta = %e\n", beta);
- }
-
- template <typename Tinput, typename Trhs=Tinput>
- StepperBase<Tinput,Trhs>* stepperFactory(ArgumentParser& p, Tinput& U, KernelBase<Tinput,Trhs> * const kern=nullptr)
- {
- assert(kern != nullptr);
-
- if (p("-ts").asString("rkv56") == "euler")
- {
- cout << "Allocating Euler time stepper..." << endl;
- return new Euler<Tinput, Trhs>(U, *this, kern);
- }
- else if (p("-ts").asString("rkv56") == "rk4")
- {
- cout << "Allocating RK4 time stepper..." << endl;
- return new RK4<Tinput, Trhs>(U, *this, kern);
- }
- else if (p("-ts").asString("rkv56") == "lsrk3")
- {
- cout << "Allocating LSRK3 time stepper..." << endl;
- return new LSRK3<Tinput, Trhs>(U, *this, kern);
- }
- else if (p("-ts").asString("rkv56") == "rkf45")
- {
- cout << "Allocating RKF45 adaptive time stepper..." << endl;
- this->bFixedStep = false;
- this->print();
- return new RKF45<Tinput, Trhs>(U, *this, kern);
- }
- else if (p("-ts").asString("rkv56") == "rkv56")
- {
- cout << "Allocating RKV56 adaptive time stepper..." << endl;
- this->bFixedStep = false;
- this->print();
- return new RKV56<Tinput, Trhs>(U, *this, kern);
- }
-#ifdef _USE_SUNDIALS_
- else if (p("-ts").asString("rkv56") == "bdf")
- {
- cout << "Allocating Sundials BDF/Newton implicit time stepper..." << endl;
- return new BDF<Tinput, Trhs>(U, *this, kern);
- }
-#endif /* _USE_SUNDIALS_ */
- else
- return nullptr;
- }
-};
-
-
-template <typename Tinput, typename Trhs>
-class StepperBase
-{
- StepperBase(Tinput const& rhs) = delete;
- StepperBase& operator=(Tinput const& rhs) = delete;
-
-protected:
- StepperSettings& m_settings;
- Tinput& m_U;
- KernelBase<Tinput,Trhs> * const m_kernel;
-
-public:
- StepperBase(Tinput& U, StepperSettings& settings, KernelBase<Tinput,Trhs> * const kern=nullptr) : m_settings(settings), m_U(U), m_kernel(kern) { }
- virtual ~StepperBase() { }
-
- inline Tinput& U() { return m_U; }
- inline Tinput const& U() const { return m_U; }
-
- virtual void step(void const* const data=nullptr) = 0;
-
- inline void write(const size_t step, const Real t, const Real dt, const Tinput& U, const void * const data=nullptr)
- {
- m_kernel->write(step, t, dt, U, data);
- }
-
- inline void dispose() { delete m_kernel; }
-};
-
-// serializer
-struct _simulation_state_t
-{
- Real t;
- Real dt;
- Real tFinal;
- Real dtDump;
- Real tDump;
- size_t step;
- size_t nsteps;
- size_t writeCount;
-} __attribute__((packed));
-
-
-template <typename Tinput>
-void serialize(const Tinput& U, const StepperSettings& S, const std::string fname="restart.bin")
-{
- _simulation_state_t state;
- state.t = S.t;
- state.dt = S.dt;
- state.tFinal = S.tFinal;
- state.dtDump = S.dtDump;
- state.tDump = S.tDump;
- state.step = S.step;
- state.nsteps = S.nsteps;
- state.writeCount = S.writeCount;
-
- ofstream _stream(fname, ios::out | ios::binary);
- _stream.write((const char*)&state.t, sizeof(_simulation_state_t));
-
- const char* const _src = (const char*)U.data();
- _stream.write(_src, sizeof(typename Tinput::DataType)*U.size());
-
- _stream.close();
-}
-
-template <typename Tinput>
-void deserialize(Tinput& U, StepperSettings& S, const std::string fname="restart.bin")
-{
- ifstream _stream(fname, ios::in | ios::binary);
-
- _simulation_state_t state;
- _stream.read((char*)&state.t, sizeof(_simulation_state_t));
- S.t = state.t;
- S.dt = state.dt;
- S.tFinal = state.tFinal;
- S.dtDump = state.dtDump;
- S.tDump = state.tDump;
- S.step = state.step;
- S.nsteps = state.nsteps;
- S.writeCount = state.writeCount;
-
- char* const _dst = (char*)U.data();
- _stream.read(_dst, sizeof(typename Tinput::DataType)*U.size());
-
- _stream.close();
-}
-
-#endif /* STEPPERBASE_H_FCCSM4HL */
diff --git a/TimeStepper/TimeStepper.h b/TimeStepper/TimeStepper.h
@@ -1,23 +0,0 @@
-// File : TimeStepper.h
-// Date : Thu Sep 22 13:23:57 2016
-// Author : Fabian Wermelinger
-// Description: Structure to setup a time stepper
-// Copyright 2016 ETH Zurich. All Rights Reserved.
-#ifndef TIMESTEPPER_H_XFR2JEZJ
-#define TIMESTEPPER_H_XFR2JEZJ
-
-#include <cstdio>
-#include "ArgumentParser.h"
-#include "KernelBase.h"
-#include "StepperBase.h"
-
-#include "explicit/Euler.h"
-#include "explicit/RK4.h"
-#include "explicit/LSRK3.h"
-#include "explicit/variableStep/RKF45.h"
-#include "explicit/variableStep/RKV56.h"
-#ifdef _USE_SUNDIALS_
-#include "implicit/BDF.h"
-#endif /* _USE_SUNDIALS_ */
-
-#endif /* TIMESTEPPER_H_XFR2JEZJ */
diff --git a/TimeStepper/explicit/Euler.h b/TimeStepper/explicit/Euler.h
@@ -1,37 +0,0 @@
-/* File: Euler.h */
-/* Date: Fri Feb 5 18:57:57 2016 */
-/* Author: Fabian Wermelinger */
-/* Tag: 1st order Euler */
-/* Copyright 2016 ETH Zurich. All Rights Reserved. */
-#ifndef EULER_H_8R2CXYR0
-#define EULER_H_8R2CXYR0
-
-#include <cstdio>
-#include "StepperBase.h"
-
-template <typename Tinput, typename Trhs>
-class Euler : public StepperBase<Tinput,Trhs>
-{
-protected:
- using StepperBase<Tinput,Trhs>::m_settings;
- using StepperBase<Tinput,Trhs>::m_U;
- using StepperBase<Tinput,Trhs>::m_kernel;
- Trhs m_rhs;
-
-public:
- Euler(Tinput& U, StepperSettings& S, KernelBase<Tinput,Trhs> * const kern) : StepperBase<Tinput,Trhs>(U,S,kern), m_rhs(U.size()) { }
- virtual ~Euler() {}
-
- virtual void step(void const* const data=nullptr)
- {
- m_kernel->compute(m_U, m_rhs, m_settings.t, data);
- m_U += static_cast<Real>(m_settings.dt)*m_rhs;
- m_settings.t += m_settings.dt;
- ++(m_settings.step);
-
- if (m_settings.step % m_settings.reportGranularity == 0)
- std::printf("Time = %e;\tStep = %d\n", m_settings.t, m_settings.step);
- }
-};
-
-#endif /* EULER_H_8R2CXYR0 */
diff --git a/TimeStepper/explicit/LSRK3.h b/TimeStepper/explicit/LSRK3.h
@@ -1,56 +0,0 @@
-/* File: LSRK3.h */
-/* Date: Fri Feb 5 18:58:48 2016 */
-/* Author: Fabian Wermelinger */
-/* Tag: 3rd order low-storage Runge-Kutta. Requires a compatible kernel */
-/* Copyright 2016 ETH Zurich. All Rights Reserved. */
-#ifndef LSRK3_H_MJ6ACLH4
-#define LSRK3_H_MJ6ACLH4
-
-#include "explicit/Euler.h"
-
-template <typename Tinput, typename Trhs>
-class LSRK3 : public Euler<Tinput, Trhs>
-{
- const Real m_A[3];
- const Real m_B[3];
- using Euler<Tinput,Trhs>::m_settings;
- using Euler<Tinput,Trhs>::m_U;
- using Euler<Tinput,Trhs>::m_rhs;
- using Euler<Tinput,Trhs>::m_kernel;
-
-public:
- LSRK3(Tinput& U, StepperSettings& S, KernelBase<Tinput,Trhs> * const kern) :
- Euler<Tinput,Trhs>(U,S,kern),
- // m_A{0.0, -2.915492524638791, -0.000000093517376}, // Gottlieb & Shu (Optimal LSRK3 for CFL = 0.32)
- // m_B{0.924574000000000, 0.287713063186749, 0.626538109512740}
- m_A{ 0.0, -17./32, -32./27},
- m_B{1./4, 8./9, 3./4}
- { }
- virtual ~LSRK3() {}
-
- virtual void step(void const* const data=nullptr)
- {
- // stage 1
- Real lsrk_data[2] = {static_cast<Real>(m_settings.dt), m_A[0]};
- m_kernel->compute(m_U, m_rhs, m_settings.t, lsrk_data);
- m_U += m_B[0]*m_rhs;
-
- // stage 2
- lsrk_data[1] = m_A[1];
- m_kernel->compute(m_U, m_rhs, m_settings.t, lsrk_data);
- m_U += m_B[1]*m_rhs;
-
- // stage 3
- lsrk_data[1] = m_A[2];
- m_kernel->compute(m_U, m_rhs, m_settings.t, lsrk_data);
- m_U += m_B[2]*m_rhs;
-
- m_settings.t += m_settings.dt;
- ++(m_settings.step);
-
- if (m_settings.step % m_settings.reportGranularity == 0)
- std::printf("Time = %e;\tStep = %d\n", m_settings.t, m_settings.step);
- }
-};
-
-#endif /* LSRK3_H_MJ6ACLH4 */
diff --git a/TimeStepper/explicit/RK4.h b/TimeStepper/explicit/RK4.h
@@ -1,53 +0,0 @@
-/* File: RK4.h */
-/* Date: Fri Feb 5 19:14:44 2016 */
-/* Author: Fabian Wermelinger */
-/* Tag: Classical 4th order Runge-Kutta */
-/* Copyright 2016 ETH Zurich. All Rights Reserved. */
-#ifndef RK4_H_LROUVHDE
-#define RK4_H_LROUVHDE
-
-#include "explicit/Euler.h"
-
-template <typename Tinput, typename Trhs>
-class RK4 : public Euler<Tinput, Trhs>
-{
- using Euler<Tinput, Trhs>::m_settings;
- using Euler<Tinput, Trhs>::m_U;
- using Euler<Tinput, Trhs>::m_rhs;
- using Euler<Tinput, Trhs>::m_kernel;
-
-public:
- RK4(Tinput& U, StepperSettings& S, KernelBase<Tinput,Trhs> * const kern) : Euler<Tinput,Trhs>(U,S,kern) { }
- virtual ~RK4() {}
-
- virtual void step(void const* const data=nullptr)
- {
- const Real oneHalf = static_cast<Real>(0.5);
- const Real oneThird = static_cast<Real>(1./3.);
- Tinput tmpU = m_U;
-
- // stage 1
- m_kernel->compute(tmpU, m_rhs, m_settings.t, data);
- m_U += oneHalf*oneThird*m_settings.dt*m_rhs;
-
- // stage 2
- m_kernel->compute(tmpU + oneHalf*m_settings.dt*m_rhs, m_rhs, m_settings.t + oneHalf*m_settings.dt, data);
- m_U += oneThird*m_settings.dt*m_rhs;
-
- // stage 3
- m_kernel->compute(tmpU + oneHalf*m_settings.dt*m_rhs, m_rhs, m_settings.t + oneHalf*m_settings.dt, data);
- m_U += oneThird*m_settings.dt*m_rhs;
-
- // stage 4
- m_kernel->compute(tmpU + m_settings.dt*m_rhs, m_rhs, m_settings.t + m_settings.dt, data);
- m_U += oneHalf*oneThird*m_settings.dt*m_rhs;
-
- m_settings.t += m_settings.dt;
- ++(m_settings.step);
-
- if (m_settings.step % m_settings.reportGranularity == 0)
- std::printf("Time = %e;\tStep = %d\n", m_settings.t, m_settings.step);
- }
-};
-
-#endif /* RK4_H_LROUVHDE */
diff --git a/TimeStepper/explicit/variableStep/RKF45.h b/TimeStepper/explicit/variableStep/RKF45.h
@@ -1,187 +0,0 @@
-/* File: RKF45.h */
-/* Date: Wed Feb 10 17:33:42 2016 */
-/* Author: Fabian Wermelinger */
-/* Tag: Runge-Kutta-Fehlberg 4th-5th order variable step method */
-/* Copyright 2016 ETH Zurich. All Rights Reserved. */
-#ifndef RKF45_H_UWN6WM1V
-#define RKF45_H_UWN6WM1V
-
-#include <cmath>
-#include <cstdlib>
-#include <limits>
-
-#include "explicit/Euler.h"
-
-template <typename Tinput, typename Trhs>
-class RKF45 : public Euler<Tinput, Trhs>
-{
- using Euler<Tinput, Trhs>::m_settings;
- using Euler<Tinput, Trhs>::m_U;
- using Euler<Tinput, Trhs>::m_rhs;
- using Euler<Tinput, Trhs>::m_kernel;
-
- Real m_errold;
- const Real m_alpha;
- const Real m_beta;
- bool m_reject;
- bool m_firstPass;
-
- Trhs m_rhs2;
- Trhs m_rhs3;
- Trhs m_rhs4;
- Trhs m_rhs5;
- Trhs m_rhs6;
-
- Tinput m_output;
-
- static constexpr Real b1 = 16./135., b3 = 6656./12825., b4 = 28561./56430., b5 = -9./50., b6 = 2./55.;
- static constexpr Real c21 = 0.25, ct2 = 0.25;
- static constexpr Real c31 = 3./32., c32 = 9./32., ct3 = 3./8.;
- static constexpr Real c41 = 1932./2197., c42 = -7200./2197., c43 = 7296./2197., ct4 = 12./13.;
- static constexpr Real c51 = 439./216., c52 = -8., c53 = 3680./513., c54 = -845./4104., ct5 = 1.0;
- static constexpr Real c61 = -8./27., c62 = 2., c63 = -3544./2565., c64 = 1859./4104., c65 = -11./40., ct6 = 0.5;
- static constexpr Real e1 = 1./360., e3 = -128./4275., e4 = -2197./75240., e5 = 1./50., e6 = 2./55.;
-
- void _computeRHS(const Real t, const Real dt, void const* const data)
- {
- if (!m_reject) m_kernel->compute(m_U, m_rhs, t, data);
- m_kernel->compute(m_U + dt*(c21*m_rhs), m_rhs2, t + ct2*dt, data);
- m_kernel->compute(m_U + dt*(c31*m_rhs + c32*m_rhs2), m_rhs3, t + ct3*dt, data);
- m_kernel->compute(m_U + dt*(c41*m_rhs + c42*m_rhs2 + c43*m_rhs3), m_rhs4, t + ct4*dt, data);
- m_kernel->compute(m_U + dt*(c51*m_rhs + c52*m_rhs2 + c53*m_rhs3 + c54*m_rhs4), m_rhs5, t + ct5*dt, data);
- 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 Real _error() const
- {
- assert(m_U.size() == m_rhs.size());
- 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 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);
- }
- }
- assert(!isnan(maxerr));
- return maxerr;
- }
-
- bool _inBound(const Real t, Real& dt, void const* const data, Real& dt_next)
- {
- 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, static_cast<Real>(1.0e-4));
- m_reject = false;
- m_firstPass = true;
- return true;
- }
-
- const Real safety = m_settings.safety;
- const Real minScale = m_settings.minScale;
- const Real maxScale = m_settings.maxScale;
- Real scale;
- if (err <= 1.0)
- {
- if (err == 0.0)
- scale = maxScale;
- else
- {
- scale = safety*std::pow(err,-m_alpha)*std::pow(m_errold,m_beta);
- scale = (scale < minScale) ? minScale : ((scale > maxScale) ? maxScale : scale);
- }
- if (m_reject)
- dt_next = dt*std::min(scale, static_cast<Real>(1.0));
- else
- dt_next = dt*scale;
- m_errold = std::max(err, static_cast<Real>(1.0e-4));
- m_reject = false;
- m_firstPass = true;
- return true;
- }
- else
- {
- scale = safety*std::pow(err,-m_alpha);
- scale = std::max(scale, minScale);
- dt *= scale;
- m_reject = true;
- m_firstPass = false;
- _computeRHS(t, dt, data);
- return false;
- }
- }
-
- void _dumpOutput(const void * const data)
- {
- 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;
-
- m_kernel->compute(m_U, m_rhs, t+dt, data);
- m_firstPass = false;
-
- while(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);
- tDump += m_settings.dtDump;
- }
- }
-
-public:
- RKF45(Tinput& U, StepperSettings& S, KernelBase<Tinput,Trhs> * const kern) :
- Euler<Tinput,Trhs>(U,S,kern), m_errold(1.0e-4), m_alpha(0.25*(S.alpha-0.75*S.beta)), m_beta(0.25*S.beta), m_reject(false), m_firstPass(true),
- m_rhs2(U.size()), m_rhs3(U.size()), m_rhs4(U.size()), m_rhs5(U.size()), m_rhs6(U.size()), m_output(U.size()) {}
- virtual ~RKF45() {}
-
- virtual void step(void const* const data=nullptr)
- {
- size_t& step = m_settings.step;
- Real& t = m_settings.t;
- Real& dt = m_settings.dt;
- const Real& tFinal = m_settings.tFinal;
- Real& tDump = m_settings.tDump;
- if (!m_settings.bRestart && m_settings.dtDump > 0.0) tDump = t + m_settings.dtDump;
-
- while(t < tFinal)
- {
- Real dt_next;
- dt = (tFinal-t) < dt ? (tFinal-t) : dt;
- _computeRHS(t, dt, data);
- while (!_inBound(t, dt, data, dt_next)) {}
- m_output = m_U;
- m_U += dt*(b1*m_rhs + b3*m_rhs3 + b4*m_rhs4 + b5*m_rhs5 + b6*m_rhs6);
- if (m_settings.dtDump > 0.0)
- {
- if (tDump <= (t+dt)) _dumpOutput(data);
- }
- else
- if ((step+1) % m_settings.writeGranularity == 0)
- m_kernel->write((step+1), t+dt, dt, m_U, data);
- t += dt;
- dt = dt_next;
- ++step;
-
- if (m_settings.step % m_settings.reportGranularity == 0)
- std::printf("Time = %e;\tStep = %d\n", m_settings.t, m_settings.step);
-
- if (m_settings.restartstep > 0 && step % m_settings.restartstep == 0)
- serialize<Tinput>(m_U, m_settings);
- }
- }
-};
-
-#endif /* RKF45_H_UWN6WM1V */
diff --git a/TimeStepper/explicit/variableStep/RKV56.h b/TimeStepper/explicit/variableStep/RKV56.h
@@ -1,194 +0,0 @@
-/* File: RKV56.h */
-/* Date: Thu Feb 11 11:36:11 2016 */
-/* Author: Fabian Wermelinger */
-/* Tag: Runge-Kutta-Verner 5th-6th order variable step method */
-/* Copyright 2016 ETH Zurich. All Rights Reserved. */
-#ifndef RKV56_H_A6VN39JN
-#define RKV56_H_A6VN39JN
-
-#include <iostream>
-#include <cmath>
-#include <cstdlib>
-#include <limits>
-
-#include "explicit/Euler.h"
-
-template <typename Tinput, typename Trhs>
-class RKV56 : public Euler<Tinput, Trhs>
-{
- using Euler<Tinput, Trhs>::m_settings;
- using Euler<Tinput, Trhs>::m_U;
- using Euler<Tinput, Trhs>::m_rhs;
- using Euler<Tinput, Trhs>::m_kernel;
-
- Real m_errold;
- const Real m_alpha;
- const Real m_beta;
- bool m_reject;
- bool m_firstPass;
-
- Trhs m_rhs2;
- Trhs m_rhs3;
- Trhs m_rhs4;
- Trhs m_rhs5;
- Trhs m_rhs6;
- Trhs m_rhs7;
- Trhs m_rhs8;
-
- Tinput m_output;
-
- static constexpr Real b1 = 3./40., b3 = 875./2244., b4 = 23./72., b5 = 264./1955., b7 = 125./11592., b8 = 43./616.;
- static constexpr Real c21 = 1./6., ct2 = 1./6.;
- static constexpr Real c31 = 4./75., c32 = 16./75., ct3 = 4./15.;
- static constexpr Real c41 = 5./6., c42 = -8./3., c43 = 2.5, ct4 = 2./3.;
- static constexpr Real c51 = -165./64., c52 = 55./6., c53 = -425./64., c54 = 85./96., ct5 = 5./6.;
- static constexpr Real c61 = 2.4, c62 = -8., c63 = 4015./612., c64 = -11./36., c65 = 88./255., ct6 = 1.0;
- static constexpr Real c71 = -8263./15000., c72 = 124./75., c73 = -643./680., c74 = -81./250., c75 = 2484./10625., ct7 = 1./15.;
- static constexpr Real c81 = 3501./1720., c82 = -300./43., c83 = 297275./52632., c84 = -319./2322., c85 = 24068./84065., c87 = 3850./26703., ct8 = 1.0;
- static constexpr Real e1 = -1./160., e3 = -125./17952., e4 = 1./144., e5 = -12./1955., e6 = -3./44., e7 = 125./11592., e8 = 43./616.;
-
- void _computeRHS(const Real t, const Real dt, void const* const data)
- {
- if (!m_reject) m_kernel->compute(m_U, m_rhs, t, data);
- m_kernel->compute(m_U + dt*(c21*m_rhs), m_rhs2, t + ct2*dt, data);
- m_kernel->compute(m_U + dt*(c31*m_rhs + c32*m_rhs2), m_rhs3, t + ct3*dt, data);
- m_kernel->compute(m_U + dt*(c41*m_rhs + c42*m_rhs2 + c43*m_rhs3), m_rhs4, t + ct4*dt, data);
- m_kernel->compute(m_U + dt*(c51*m_rhs + c52*m_rhs2 + c53*m_rhs3 + c54*m_rhs4), m_rhs5, t + ct5*dt, data);
- 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);
- m_kernel->compute(m_U + dt*(c71*m_rhs + c72*m_rhs2 + c73*m_rhs3 + c74*m_rhs4 + c75*m_rhs5), m_rhs7, t + ct7*dt, data);
- 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 Real _error() const
- {
- assert(m_U.size() == m_rhs.size());
- 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 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);
- }
- }
- assert(!isnan(maxerr));
- return maxerr;
- }
-
- bool _inBound(const Real t, Real& dt, void const* const data, Real& dt_next)
- {
- 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, static_cast<Real>(1.0e-4));
- m_reject = false;
- m_firstPass = true;
- return true;
- }
-
- const Real safety = m_settings.safety;
- const Real minScale = m_settings.minScale;
- const Real maxScale = m_settings.maxScale;
- Real scale;
- if (err <= 1.0)
- {
- if (err == 0.0)
- scale = maxScale;
- else
- {
- scale = safety*std::pow(err,-m_alpha)*std::pow(m_errold,m_beta);
- scale = (scale < minScale) ? minScale : ((scale > maxScale) ? maxScale : scale);
- }
- if (m_reject)
- dt_next = dt*std::min(scale, static_cast<Real>(1.0));
- else
- dt_next = dt*scale;
- m_errold = std::max(err, static_cast<Real>(1.0e-4));
- m_reject = false;
- m_firstPass = true;
- return true;
- }
- else
- {
- scale = safety*std::pow(err,-m_alpha);
- scale = std::max(scale, minScale);
- dt *= scale;
- m_reject = true;
- m_firstPass = false;
- _computeRHS(t, dt, data);
- return false;
- }
- }
-
- void _dumpOutput(const void * const data)
- {
- 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;
-
- m_kernel->compute(m_U, m_rhs, t+dt, data);
- m_firstPass = false;
-
- while(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);
- tDump += m_settings.dtDump;
- }
- }
-
-public:
- RKV56(Tinput& U, StepperSettings& S, KernelBase<Tinput,Trhs> * const kern) :
- Euler<Tinput,Trhs>(U,S,kern), m_errold(1.0e-4), m_alpha(0.25*(S.alpha-0.75*S.beta)), m_beta(0.25*S.beta), m_reject(false), m_firstPass(true),
- m_rhs2(U.size()), m_rhs3(U.size()), m_rhs4(U.size()), m_rhs5(U.size()), m_rhs6(U.size()), m_rhs7(U.size()), m_rhs8(U.size()), m_output(U.size()) {}
- virtual ~RKV56() {}
-
- virtual void step(void const* const data=nullptr)
- {
- size_t& step = m_settings.step;
- Real& t = m_settings.t;
- Real& dt = m_settings.dt;
- const Real& tFinal = m_settings.tFinal;
- Real& tDump = m_settings.tDump;
- if (!m_settings.bRestart && m_settings.dtDump > 0.0) tDump = t + m_settings.dtDump;
-
- while(t < tFinal)
- {
- Real dt_next;
- dt = (tFinal-t) < dt ? (tFinal-t) : dt;
- _computeRHS(t, dt, data);
- while (!_inBound(t, dt, data, dt_next)) {}
- m_output = m_U;
- m_U += dt*(b1*m_rhs + b3*m_rhs3 + b4*m_rhs4 + b5*m_rhs5 + b7*m_rhs7 + b8*m_rhs8);
- if (m_settings.dtDump > 0.0)
- {
- if (tDump <= (t+dt)) _dumpOutput(data);
- }
- else
- if ((step+1) % m_settings.writeGranularity == 0)
- m_kernel->write((step+1), t+dt, dt, m_U, data);
- t += dt;
- dt = dt_next;
- ++step;
-
- if (m_settings.step % m_settings.reportGranularity == 0)
- std::printf("Time = %e;\tStep = %d\n", m_settings.t, m_settings.step);
-
- if (m_settings.restartstep > 0 && step % m_settings.restartstep == 0)
- serialize<Tinput>(m_U, m_settings);
- }
- }
-};
-
-#endif /* RKV56_H_A6VN39JN */
diff --git a/TimeStepper/implicit/BDF.h b/TimeStepper/implicit/BDF.h
@@ -1,76 +0,0 @@
-/* File: BDF.h */
-/* Date: Mon Feb 8 11:44:23 2016 */
-/* Author: Fabian Wermelinger */
-/* Tag: Backward differentiation formula using SUNDIALS library */
-/* Copyright 2016 ETH Zurich. All Rights Reserved. */
-#ifndef BDF_H_YF523CRX
-#define BDF_H_YF523CRX
-
-#ifdef _USE_SUNDIALS_
-#include <cassert>
-#include "StepperBase.h"
-
-#include "cvode/cvode.h"
-#include "cvode/cvode_dense.h"
-#include "nvector/nvector_serial.h"
-
-template <typename Tinput, typename Trhs>
-class BDF : public StepperBase<Tinput,Trhs>
-{
- Real m_a;
- void* m_cvode_mem;
- N_Vector m_y;
-
- void _init_sundials()
- {
- m_cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
- m_y = N_VMake_Serial(Tinput::SIZE, m_U.data());
- CVodeInit(m_cvode_mem, BDF::f, 0.0, m_y);
- CVodeSStolerances(m_cvode_mem, 1.0e-6, 1.0e-6);
- CVDense(m_cvode_mem, Trhs::SIZE);
- CVodeSetUserData(m_cvode_mem, static_cast<void*>(&m_a));
- }
-
-protected:
- using StepperBase<Tinput,Trhs>::m_settings;
- using StepperBase<Tinput,Trhs>::m_U;
- using StepperBase<Tinput,Trhs>::m_kernel;
- Trhs m_rhs;
-
-public:
- BDF(Tinput& U, StepperSettings& S, KernelBase<Tinput,Trhs> * const kern) :
- StepperBase<Tinput,Trhs>(U,S,kern), m_rhs(U.size())
- {
- _init_sundials();
- }
- ~BDF() { CVodeFree(&m_cvode_mem); }
-
- virtual void step(void const* const data=nullptr)
- {
- /* TODO: (fabianw; Fri Feb 12 13:26:02 2016) this does not work yet.
- * Do not compile with Sundials support */
- /* m_a = a; */
- Real tret; // interface dummy
- CVode(m_cvode_mem, t+dt, m_y, &tret, CV_ONE_STEP);
- t = tret;
-
- if (m_settings.step % m_settings.reportGranularity == 0)
- std::printf("Time = %e;\tStep = %d\n", m_settings.t, m_settings.step);
- }
-
- virtual inline void setIC(Real const ic, const Real t0=0)
- {
- m_U = ic;
- CVodeReInit(m_cvode_mem, t0, m_y);
- }
-
- static inline int f(Real t, N_Vector y, N_Vector ydot, void *a)
- {
- Tkernel kernel(Trhs::SIZE);
- kernel.compute(NV_DATA_S(y), NV_DATA_S(ydot), *static_cast<Real*>(a), t);
- return 0;
- }
-};
-#endif /* _USE_SUNDIALS_ */
-
-#endif /* BDF_H_YF523CRX */
diff --git a/TimeStepper/orderVerification/orderVerification.cpp b/TimeStepper/orderVerification/orderVerification.cpp
@@ -1,145 +0,0 @@
-/* File: orderVerification.cpp */
-/* Date: Sun Feb 7 22:48:30 2016 */
-/* Author: Fabian Wermelinger */
-/* Tag: Verify order of time stepper */
-/* Copyright 2016 ETH Zurich. All Rights Reserved. */
-#include <cassert>
-#include <iostream>
-#include <cmath>
-#include <vector>
-#include <string>
-#include <utility>
-#include <fstream>
-
-#include "common.h"
-#include "KernelBase.h"
-#include "TimeStepper.h"
-#include "GnuplotDump.h"
-
-using namespace std;
-
-constexpr Real DECAY = -1.0;
-
-template <typename Tinput, typename Trhs=Tinput>
-class Dahlquist : public KernelBase<Tinput,Trhs>
-{
- GnuplotDump m_gp;
-
-protected:
- const Real m_a;
-
-public:
- Dahlquist(const Real a=1.0) : m_gp("out", GnuplotDump::ASCII), m_a(a) {}
- ~Dahlquist() {}
-
- void compute(const Tinput& U, Trhs& rhs, const Real t, void const* const data=nullptr)
- {
- rhs = m_a*U;
- }
-
- void write(const size_t step, const Real t, const Real dt, const Tinput& U, const void * const data=nullptr)
- {
- m_gp.write(step, t, dt, U.data(), U.size()*Tinput::DataType::SIZE);
- }
-};
-
-template <typename Tinput, typename Trhs=Tinput>
-class DahlquistLSRK3 : public Dahlquist<Tinput,Trhs>
-{
- using Dahlquist<Tinput, Trhs>::m_a;
-
-public:
- DahlquistLSRK3(const Real a=1.0) : Dahlquist<Tinput,Trhs>(a) {}
- ~DahlquistLSRK3() {}
-
- void compute(const Tinput& U, Trhs& rhs, const Real t, void const* const data=nullptr)
- {
- assert(data != nullptr);
- const Real * const LSRKData = static_cast<const Real* const>(data);
- const Real dt = LSRKData[0];
- const Real A = LSRKData[1];
- rhs = A*rhs + dt*m_a*U;
- }
-};
-
-// types we use
-typedef State<Real,1> State1;
-using vec_t = LightVector<State1>;
-using baseKernel_t = KernelBase<vec_t>*;
-using baseStepper_t = StepperBase<vec_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) {}
- 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)
-{
- vector<Error> error(stages);
- S.dt = dt;
- for (int s = 0; s < stages; ++s)
- {
- size_t const N = static_cast<size_t>(S.tFinal/S.dt)+1;
- S.t = 0.0;
- stepper->U() = 1.0;
- Real& numerical = stepper->U()[0][0];
- for (size_t i = 0; i < N; ++i)
- {
- stepper->step();
- 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);
- }
- error[s].L1 /= N;
- error[s].L2 = sqrt(error[s].L2/N);
- error[s].dt = S.dt;
- S.dt *= scale;
- }
- return error;
-}
-
-int main(int argc, const char** argv)
-{
- // set up solution vector
- vec_t U;
-
- // kernel instances
- baseKernel_t dahlquist = new Dahlquist<vec_t>(DECAY);
- baseKernel_t dahlquistLSRK3 = new DahlquistLSRK3<vec_t>(DECAY);
-
- // setup time steppers
- ArgumentParser parser(argc, argv);
- StepperSettings S(parser);
-
- vector<pair<pair<string,Real>, baseStepper_t> > timeStepperList;
- timeStepperList.push_back(make_pair(make_pair("Euler",1.0), new Euler<vec_t>(U, S, dahlquist)));
- timeStepperList.push_back(make_pair(make_pair("RK4",4.0), new RK4<vec_t>(U, S, dahlquist)));
- timeStepperList.push_back(make_pair(make_pair("LSRK3",3.0), new LSRK3<vec_t>(U, S, dahlquistLSRK3)));
- // timeStepperList.push_back(make_pair(make_pair("RKF45",5.0), new RKF45<vec_t>(U, S, dahlquist)));
- // timeStepperList.push_back(make_pair(make_pair("RKV56",6.0), new RKV56<vec_t>(U, S, dahlquist)));
-#ifdef _USE_SUNDIALS_
- timeStepperList.push_back(make_pair(make_pair("BDF",5.0), new BDF<vec_t>(U, S, dahlquist)));
-#endif /* _USE_SUNDIALS_ */
-
- for (size_t i = 0; i < timeStepperList.size(); ++i)
- {
- 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 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();
- }
-
- for (size_t i = 0; i < timeStepperList.size(); ++i)
- delete timeStepperList[i].second;
- delete dahlquist;
- delete dahlquistLSRK3;
-
- return 0;
-}
diff --git a/ArgumentParser.h b/include/ArgumentParser.h
diff --git a/include/GnuplotDump.h b/include/GnuplotDump.h
@@ -0,0 +1,115 @@
+/* File: GnuplotDump.h */
+/* Date: Tue Mar 1 21:07:49 2016 */
+/* Author: Fabian Wermelinger */
+/* Tag: Gnuplot Dumper */
+/* Copyright 2016 ETH Zurich. All Rights Reserved. */
+#ifndef GNUPLOTDUMP_H_2VNNTIUB
+#define GNUPLOTDUMP_H_2VNNTIUB
+
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+#include <string>
+
+#ifdef _FLOAT_PRECISION_
+using Real = float;
+#else
+using Real = double;
+#endif
+
+#ifdef _GNUPLOT_DOUBLE_
+using GPfloat = double;
+#else
+using GPfloat = float;
+#endif
+
+class GnuplotDump
+{
+private:
+ const int m_ftype;
+ std::string m_fname;
+ std::ofstream m_file;
+ size_t m_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())
+ {
+ m_file << step << '\t' << t << '\t' << dt;
+ for (size_t i = 0; i < N; ++i)
+ m_file << '\t' << src[i];
+ m_file << std::endl;
+ m_N = N+3;
+ }
+ }
+
+ void _writeBIN(const size_t step, const GPfloat t, const GPfloat dt, const Real * const src, size_t N)
+ {
+ if (1024 < N)
+ {
+ std::cerr << "GnuplotDump: WARNING -- Max. Column width for BIN dump is 1024" << std::endl;
+ N = 1024;
+ }
+ GPfloat buf[1024+3];
+
+ buf[0] = static_cast<GPfloat>(step);
+ buf[1] = t;
+ buf[2] = dt;
+ for (size_t i = 0; i < N; ++i)
+ buf[i+3] = static_cast<GPfloat>(src[i]);
+
+ if (m_file.good())
+ {
+ m_file.write((const char*)(&buf[0]), (N+3)*sizeof(GPfloat));
+ m_N = N+3;
+ }
+ }
+
+ void _writeGPscript() const
+ {
+ std::ofstream out(m_fname+".gp");
+ if (m_ftype == ASCII)
+ out << "plot '" << m_fname+".dat" << "' using 1:2 with lines" << std::endl;
+ else
+ {
+ out << "plot '" << m_fname+".bin" << "' binary format='\%" << m_N;
+ if (sizeof(GPfloat) == 4)
+ out << "float";
+ else
+ out << "double";
+ out << "' using 1:2 with lines" << std::endl;
+ }
+ out.close();
+ }
+
+public:
+ GnuplotDump(const std::string fname="gnuplot", const int type=BIN, const size_t pos=0) : m_ftype(type), m_fname(fname)
+ {
+ if (m_ftype == ASCII)
+ {
+ m_file.open(m_fname+".dat", std::ios_base::out|std::ios_base::app);
+ m_file.setf(std::ios::scientific, std::ios::floatfield);
+ }
+ else
+ m_file.open(m_fname+".bin", std::ios_base::out|std::ios_base::binary|std::ios_base::app);
+
+ if (pos > 0) m_file.seekp(pos);
+ }
+ ~GnuplotDump()
+ {
+ m_file.close();
+ _writeGPscript();
+ }
+
+ enum {ASCII, BIN};
+
+ 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);
+ else
+ _writeBIN(step, t, dt, src, N);
+ }
+};
+
+#endif /* GNUPLOTDUMP_H_2VNNTIUB */
diff --git a/TimeStepper/KernelBase.h b/include/TimeStepper/KernelBase.h
diff --git a/include/TimeStepper/StepperBase.h b/include/TimeStepper/StepperBase.h
@@ -0,0 +1,235 @@
+// File : StepperBase.h
+// Date : Thu Sep 22 13:27:19 2016
+// Author : Fabian Wermelinger
+// Description: Time stepper base class
+// Copyright 2016 ETH Zurich. All Rights Reserved.
+#ifndef STEPPERBASE_H_FCCSM4HL
+#define STEPPERBASE_H_FCCSM4HL
+
+#include <TimeStepper/KernelBase.h>
+
+// forward declarations
+template <typename Tinput, typename Trhs=Tinput>
+class StepperBase;
+template <typename Tinput, typename Trhs=Tinput>
+class Euler;
+template <typename Tinput, typename Trhs=Tinput>
+class RK4;
+template <typename Tinput, typename Trhs=Tinput>
+class LSRK3;
+template <typename Tinput, typename Trhs=Tinput>
+class RKF45;
+template <typename Tinput, typename Trhs=Tinput>
+class RKV56;
+#ifdef _USE_SUNDIALS_
+template <typename Tinput, typename Trhs=Tinput>
+class BDF;
+#endif /* _USE_SUNDIALS_ */
+
+#ifdef _FLOAT_PRECISION_
+using Real = float;
+#else
+using Real = double;
+#endif
+
+struct StepperSettings
+{
+ StepperSettings(ArgumentParser& p)
+ {
+ aTol = p("-ts_aTol").asDouble(1.0e-8);
+ rTol = p("-ts_rTol").asDouble(1.0e-8);
+ minScale = p("-ts_minScale").asDouble(0.2);
+ maxScale = p("-ts_maxScale").asDouble(10.0);
+ alpha = p("-ts_alpha").asDouble(1.0);
+ beta = p("-ts_beta").asDouble(0.0);
+ safety = p("-ts_safety").asDouble(0.9);
+ t = p("-ts_t0").asDouble(0.0);
+ dt = p("-ts_dt").asDouble(1.0e-4);
+ step = 0;
+ nsteps = p("-ts_nsteps").asInt(0);
+ p.set_strict_mode();
+ tFinal = p("-ts_tend").asDouble();
+ p.unset_strict_mode();
+ dtDump = p("-ts_dtDump").asDouble(-1.0);
+ tDump = t;
+ writeGranularity = p("-ts_writeGranularity").asInt(1);
+ writeCount = 0;
+ reportGranularity= p("-ts_reportGranularity").asInt(100);
+ bFixedStep = true;
+
+ // restarts
+ bRestart = p("-ts_restart").asBool(false);
+ restartstep = p("-ts_restart_step").asInt(0);
+
+ }
+
+ Real aTol;
+ Real rTol;
+ Real minScale;
+ Real maxScale;
+ Real alpha;
+ Real beta;
+ Real safety;
+ Real t;
+ Real dt;
+ size_t step;
+ size_t nsteps;
+ Real tFinal;
+ Real dtDump;
+ Real tDump;
+ size_t writeGranularity;
+ size_t writeCount;
+ size_t reportGranularity;
+ bool bFixedStep;
+
+ // restarts
+ bool bRestart;
+ int restartstep;
+
+
+ // helper
+ void print() const
+ {
+ printf("Time Stepper :\n");
+ printf("\tAbsolute tolerance = %e\n", aTol);
+ printf("\tRelative tolerance = %e\n", rTol);
+ printf("\tMinimum time step scale = %e\n", minScale);
+ printf("\tMaximum time step scale = %e\n", maxScale);
+ printf("\tError PI control alpha = %e\n", alpha);
+ printf("\tError PI control beta = %e\n", beta);
+ }
+
+ template <typename Tinput, typename Trhs=Tinput>
+ StepperBase<Tinput,Trhs>* stepperFactory(ArgumentParser& p, Tinput& U, KernelBase<Tinput,Trhs> * const kern=nullptr)
+ {
+ assert(kern != nullptr);
+
+ if (p("-ts").asString("rkv56") == "euler")
+ {
+ cout << "Allocating Euler time stepper..." << endl;
+ return new Euler<Tinput, Trhs>(U, *this, kern);
+ }
+ else if (p("-ts").asString("rkv56") == "rk4")
+ {
+ cout << "Allocating RK4 time stepper..." << endl;
+ return new RK4<Tinput, Trhs>(U, *this, kern);
+ }
+ else if (p("-ts").asString("rkv56") == "lsrk3")
+ {
+ cout << "Allocating LSRK3 time stepper..." << endl;
+ return new LSRK3<Tinput, Trhs>(U, *this, kern);
+ }
+ else if (p("-ts").asString("rkv56") == "rkf45")
+ {
+ cout << "Allocating RKF45 adaptive time stepper..." << endl;
+ this->bFixedStep = false;
+ this->print();
+ return new RKF45<Tinput, Trhs>(U, *this, kern);
+ }
+ else if (p("-ts").asString("rkv56") == "rkv56")
+ {
+ cout << "Allocating RKV56 adaptive time stepper..." << endl;
+ this->bFixedStep = false;
+ this->print();
+ return new RKV56<Tinput, Trhs>(U, *this, kern);
+ }
+#ifdef _USE_SUNDIALS_
+ else if (p("-ts").asString("rkv56") == "bdf")
+ {
+ cout << "Allocating Sundials BDF/Newton implicit time stepper..." << endl;
+ return new BDF<Tinput, Trhs>(U, *this, kern);
+ }
+#endif /* _USE_SUNDIALS_ */
+ else
+ return nullptr;
+ }
+};
+
+
+template <typename Tinput, typename Trhs>
+class StepperBase
+{
+ StepperBase(Tinput const& rhs) = delete;
+ StepperBase& operator=(Tinput const& rhs) = delete;
+
+protected:
+ StepperSettings& m_settings;
+ Tinput& m_U;
+ KernelBase<Tinput,Trhs> * const m_kernel;
+
+public:
+ StepperBase(Tinput& U, StepperSettings& settings, KernelBase<Tinput,Trhs> * const kern=nullptr) : m_settings(settings), m_U(U), m_kernel(kern) { }
+ virtual ~StepperBase() { }
+
+ inline Tinput& U() { return m_U; }
+ inline Tinput const& U() const { return m_U; }
+
+ virtual void step(void const* const data=nullptr) = 0;
+
+ inline void write(const size_t step, const Real t, const Real dt, const Tinput& U, const void * const data=nullptr)
+ {
+ m_kernel->write(step, t, dt, U, data);
+ }
+
+ inline void dispose() { delete m_kernel; }
+};
+
+// serializer
+struct _simulation_state_t
+{
+ Real t;
+ Real dt;
+ Real tFinal;
+ Real dtDump;
+ Real tDump;
+ size_t step;
+ size_t nsteps;
+ size_t writeCount;
+} __attribute__((packed));
+
+
+template <typename Tinput>
+void serialize(const Tinput& U, const StepperSettings& S, const std::string fname="restart.bin")
+{
+ _simulation_state_t state;
+ state.t = S.t;
+ state.dt = S.dt;
+ state.tFinal = S.tFinal;
+ state.dtDump = S.dtDump;
+ state.tDump = S.tDump;
+ state.step = S.step;
+ state.nsteps = S.nsteps;
+ state.writeCount = S.writeCount;
+
+ ofstream _stream(fname, ios::out | ios::binary);
+ _stream.write((const char*)&state.t, sizeof(_simulation_state_t));
+
+ const char* const _src = (const char*)U.data();
+ _stream.write(_src, sizeof(typename Tinput::DataType)*U.size());
+
+ _stream.close();
+}
+
+template <typename Tinput>
+void deserialize(Tinput& U, StepperSettings& S, const std::string fname="restart.bin")
+{
+ ifstream _stream(fname, ios::in | ios::binary);
+
+ _simulation_state_t state;
+ _stream.read((char*)&state.t, sizeof(_simulation_state_t));
+ S.t = state.t;
+ S.dt = state.dt;
+ S.tFinal = state.tFinal;
+ S.dtDump = state.dtDump;
+ S.tDump = state.tDump;
+ S.step = state.step;
+ S.nsteps = state.nsteps;
+ S.writeCount = state.writeCount;
+
+ char* const _dst = (char*)U.data();
+ _stream.read(_dst, sizeof(typename Tinput::DataType)*U.size());
+
+ _stream.close();
+}
+
+#endif /* STEPPERBASE_H_FCCSM4HL */
diff --git a/include/TimeStepper/TimeStepper.h b/include/TimeStepper/TimeStepper.h
@@ -0,0 +1,24 @@
+// File : TimeStepper.h
+// Date : Thu Sep 22 13:23:57 2016
+// Author : Fabian Wermelinger
+// Description: Structure to setup a time stepper
+// Copyright 2016 ETH Zurich. All Rights Reserved.
+#ifndef TIMESTEPPER_H_XFR2JEZJ
+#define TIMESTEPPER_H_XFR2JEZJ
+
+#include <ArgumentParser.h>
+#include <TimeStepper/KernelBase.h>
+#include <TimeStepper/StepperBase.h>
+
+#include <TimeStepper/explicit/Euler.h>
+#include <TimeStepper/explicit/LSRK3.h>
+#include <TimeStepper/explicit/RK4.h>
+#include <TimeStepper/explicit/variableStep/RKF45.h>
+#include <TimeStepper/explicit/variableStep/RKV56.h>
+#ifdef _USE_SUNDIALS_
+#include <TimeStepper/implicit/BDF.h>
+#endif /* _USE_SUNDIALS_ */
+
+#include <cstdio>
+
+#endif /* TIMESTEPPER_H_XFR2JEZJ */
diff --git a/include/TimeStepper/explicit/Euler.h b/include/TimeStepper/explicit/Euler.h
@@ -0,0 +1,38 @@
+/* File: Euler.h */
+/* Date: Fri Feb 5 18:57:57 2016 */
+/* Author: Fabian Wermelinger */
+/* Tag: 1st order Euler */
+/* Copyright 2016 ETH Zurich. All Rights Reserved. */
+#ifndef EULER_H_8R2CXYR0
+#define EULER_H_8R2CXYR0
+
+#include <TimeStepper/StepperBase.h>
+
+#include <cstdio>
+
+template <typename Tinput, typename Trhs>
+class Euler : public StepperBase<Tinput,Trhs>
+{
+protected:
+ using StepperBase<Tinput,Trhs>::m_settings;
+ using StepperBase<Tinput,Trhs>::m_U;
+ using StepperBase<Tinput,Trhs>::m_kernel;
+ Trhs m_rhs;
+
+public:
+ Euler(Tinput& U, StepperSettings& S, KernelBase<Tinput,Trhs> * const kern) : StepperBase<Tinput,Trhs>(U,S,kern), m_rhs(U.size()) { }
+ virtual ~Euler() {}
+
+ virtual void step(void const* const data=nullptr)
+ {
+ m_kernel->compute(m_U, m_rhs, m_settings.t, data);
+ m_U += static_cast<Real>(m_settings.dt)*m_rhs;
+ m_settings.t += m_settings.dt;
+ ++(m_settings.step);
+
+ if (m_settings.step % m_settings.reportGranularity == 0)
+ std::printf("Time = %e;\tStep = %d\n", m_settings.t, m_settings.step);
+ }
+};
+
+#endif /* EULER_H_8R2CXYR0 */
diff --git a/include/TimeStepper/explicit/LSRK3.h b/include/TimeStepper/explicit/LSRK3.h
@@ -0,0 +1,56 @@
+/* File: LSRK3.h */
+/* Date: Fri Feb 5 18:58:48 2016 */
+/* Author: Fabian Wermelinger */
+/* Tag: 3rd order low-storage Runge-Kutta. Requires a compatible kernel */
+/* Copyright 2016 ETH Zurich. All Rights Reserved. */
+#ifndef LSRK3_H_MJ6ACLH4
+#define LSRK3_H_MJ6ACLH4
+
+#include <TimeStepper/explicit/Euler.h>
+
+template <typename Tinput, typename Trhs>
+class LSRK3 : public Euler<Tinput, Trhs>
+{
+ const Real m_A[3];
+ const Real m_B[3];
+ using Euler<Tinput,Trhs>::m_settings;
+ using Euler<Tinput,Trhs>::m_U;
+ using Euler<Tinput,Trhs>::m_rhs;
+ using Euler<Tinput,Trhs>::m_kernel;
+
+public:
+ LSRK3(Tinput& U, StepperSettings& S, KernelBase<Tinput,Trhs> * const kern) :
+ Euler<Tinput,Trhs>(U,S,kern),
+ // m_A{0.0, -2.915492524638791, -0.000000093517376}, // Gottlieb & Shu (Optimal LSRK3 for CFL = 0.32)
+ // m_B{0.924574000000000, 0.287713063186749, 0.626538109512740}
+ m_A{ 0.0, -17./32, -32./27},
+ m_B{1./4, 8./9, 3./4}
+ { }
+ virtual ~LSRK3() {}
+
+ virtual void step(void const* const data=nullptr)
+ {
+ // stage 1
+ Real lsrk_data[2] = {static_cast<Real>(m_settings.dt), m_A[0]};
+ m_kernel->compute(m_U, m_rhs, m_settings.t, lsrk_data);
+ m_U += m_B[0]*m_rhs;
+
+ // stage 2
+ lsrk_data[1] = m_A[1];
+ m_kernel->compute(m_U, m_rhs, m_settings.t, lsrk_data);
+ m_U += m_B[1]*m_rhs;
+
+ // stage 3
+ lsrk_data[1] = m_A[2];
+ m_kernel->compute(m_U, m_rhs, m_settings.t, lsrk_data);
+ m_U += m_B[2]*m_rhs;
+
+ m_settings.t += m_settings.dt;
+ ++(m_settings.step);
+
+ if (m_settings.step % m_settings.reportGranularity == 0)
+ std::printf("Time = %e;\tStep = %d\n", m_settings.t, m_settings.step);
+ }
+};
+
+#endif /* LSRK3_H_MJ6ACLH4 */
diff --git a/include/TimeStepper/explicit/RK4.h b/include/TimeStepper/explicit/RK4.h
@@ -0,0 +1,53 @@
+/* File: RK4.h */
+/* Date: Fri Feb 5 19:14:44 2016 */
+/* Author: Fabian Wermelinger */
+/* Tag: Classical 4th order Runge-Kutta */
+/* Copyright 2016 ETH Zurich. All Rights Reserved. */
+#ifndef RK4_H_LROUVHDE
+#define RK4_H_LROUVHDE
+
+#include <TimeStepper/explicit/Euler.h>
+
+template <typename Tinput, typename Trhs>
+class RK4 : public Euler<Tinput, Trhs>
+{
+ using Euler<Tinput, Trhs>::m_settings;
+ using Euler<Tinput, Trhs>::m_U;
+ using Euler<Tinput, Trhs>::m_rhs;
+ using Euler<Tinput, Trhs>::m_kernel;
+
+public:
+ RK4(Tinput& U, StepperSettings& S, KernelBase<Tinput,Trhs> * const kern) : Euler<Tinput,Trhs>(U,S,kern) { }
+ virtual ~RK4() {}
+
+ virtual void step(void const* const data=nullptr)
+ {
+ const Real oneHalf = static_cast<Real>(0.5);
+ const Real oneThird = static_cast<Real>(1./3.);
+ Tinput tmpU = m_U;
+
+ // stage 1
+ m_kernel->compute(tmpU, m_rhs, m_settings.t, data);
+ m_U += oneHalf*oneThird*m_settings.dt*m_rhs;
+
+ // stage 2
+ m_kernel->compute(tmpU + oneHalf*m_settings.dt*m_rhs, m_rhs, m_settings.t + oneHalf*m_settings.dt, data);
+ m_U += oneThird*m_settings.dt*m_rhs;
+
+ // stage 3
+ m_kernel->compute(tmpU + oneHalf*m_settings.dt*m_rhs, m_rhs, m_settings.t + oneHalf*m_settings.dt, data);
+ m_U += oneThird*m_settings.dt*m_rhs;
+
+ // stage 4
+ m_kernel->compute(tmpU + m_settings.dt*m_rhs, m_rhs, m_settings.t + m_settings.dt, data);
+ m_U += oneHalf*oneThird*m_settings.dt*m_rhs;
+
+ m_settings.t += m_settings.dt;
+ ++(m_settings.step);
+
+ if (m_settings.step % m_settings.reportGranularity == 0)
+ std::printf("Time = %e;\tStep = %d\n", m_settings.t, m_settings.step);
+ }
+};
+
+#endif /* RK4_H_LROUVHDE */
diff --git a/include/TimeStepper/explicit/variableStep/RKF45.h b/include/TimeStepper/explicit/variableStep/RKF45.h
@@ -0,0 +1,188 @@
+/* File: RKF45.h */
+/* Date: Wed Feb 10 17:33:42 2016 */
+/* Author: Fabian Wermelinger */
+/* Tag: Runge-Kutta-Fehlberg 4th-5th order variable step method */
+/* Copyright 2016 ETH Zurich. All Rights Reserved. */
+#ifndef RKF45_H_UWN6WM1V
+#define RKF45_H_UWN6WM1V
+
+#include <TimeStepper/explicit/Euler.h>
+
+#include <cassert>
+#include <cmath>
+#include <cstdlib>
+#include <limits>
+
+template <typename Tinput, typename Trhs>
+class RKF45 : public Euler<Tinput, Trhs>
+{
+ using Euler<Tinput, Trhs>::m_settings;
+ using Euler<Tinput, Trhs>::m_U;
+ using Euler<Tinput, Trhs>::m_rhs;
+ using Euler<Tinput, Trhs>::m_kernel;
+
+ Real m_errold;
+ const Real m_alpha;
+ const Real m_beta;
+ bool m_reject;
+ bool m_firstPass;
+
+ Trhs m_rhs2;
+ Trhs m_rhs3;
+ Trhs m_rhs4;
+ Trhs m_rhs5;
+ Trhs m_rhs6;
+
+ Tinput m_output;
+
+ static constexpr Real b1 = 16./135., b3 = 6656./12825., b4 = 28561./56430., b5 = -9./50., b6 = 2./55.;
+ static constexpr Real c21 = 0.25, ct2 = 0.25;
+ static constexpr Real c31 = 3./32., c32 = 9./32., ct3 = 3./8.;
+ static constexpr Real c41 = 1932./2197., c42 = -7200./2197., c43 = 7296./2197., ct4 = 12./13.;
+ static constexpr Real c51 = 439./216., c52 = -8., c53 = 3680./513., c54 = -845./4104., ct5 = 1.0;
+ static constexpr Real c61 = -8./27., c62 = 2., c63 = -3544./2565., c64 = 1859./4104., c65 = -11./40., ct6 = 0.5;
+ static constexpr Real e1 = 1./360., e3 = -128./4275., e4 = -2197./75240., e5 = 1./50., e6 = 2./55.;
+
+ void _computeRHS(const Real t, const Real dt, void const* const data)
+ {
+ if (!m_reject) m_kernel->compute(m_U, m_rhs, t, data);
+ m_kernel->compute(m_U + dt*(c21*m_rhs), m_rhs2, t + ct2*dt, data);
+ m_kernel->compute(m_U + dt*(c31*m_rhs + c32*m_rhs2), m_rhs3, t + ct3*dt, data);
+ m_kernel->compute(m_U + dt*(c41*m_rhs + c42*m_rhs2 + c43*m_rhs3), m_rhs4, t + ct4*dt, data);
+ m_kernel->compute(m_U + dt*(c51*m_rhs + c52*m_rhs2 + c53*m_rhs3 + c54*m_rhs4), m_rhs5, t + ct5*dt, data);
+ 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 Real _error() const
+ {
+ assert(m_U.size() == m_rhs.size());
+ 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 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);
+ }
+ }
+ assert(!isnan(maxerr));
+ return maxerr;
+ }
+
+ bool _inBound(const Real t, Real& dt, void const* const data, Real& dt_next)
+ {
+ 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, static_cast<Real>(1.0e-4));
+ m_reject = false;
+ m_firstPass = true;
+ return true;
+ }
+
+ const Real safety = m_settings.safety;
+ const Real minScale = m_settings.minScale;
+ const Real maxScale = m_settings.maxScale;
+ Real scale;
+ if (err <= 1.0)
+ {
+ if (err == 0.0)
+ scale = maxScale;
+ else
+ {
+ scale = safety*std::pow(err,-m_alpha)*std::pow(m_errold,m_beta);
+ scale = (scale < minScale) ? minScale : ((scale > maxScale) ? maxScale : scale);
+ }
+ if (m_reject)
+ dt_next = dt*std::min(scale, static_cast<Real>(1.0));
+ else
+ dt_next = dt*scale;
+ m_errold = std::max(err, static_cast<Real>(1.0e-4));
+ m_reject = false;
+ m_firstPass = true;
+ return true;
+ }
+ else
+ {
+ scale = safety*std::pow(err,-m_alpha);
+ scale = std::max(scale, minScale);
+ dt *= scale;
+ m_reject = true;
+ m_firstPass = false;
+ _computeRHS(t, dt, data);
+ return false;
+ }
+ }
+
+ void _dumpOutput(const void * const data)
+ {
+ 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;
+
+ m_kernel->compute(m_U, m_rhs, t+dt, data);
+ m_firstPass = false;
+
+ while(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);
+ tDump += m_settings.dtDump;
+ }
+ }
+
+public:
+ RKF45(Tinput& U, StepperSettings& S, KernelBase<Tinput,Trhs> * const kern) :
+ Euler<Tinput,Trhs>(U,S,kern), m_errold(1.0e-4), m_alpha(0.25*(S.alpha-0.75*S.beta)), m_beta(0.25*S.beta), m_reject(false), m_firstPass(true),
+ m_rhs2(U.size()), m_rhs3(U.size()), m_rhs4(U.size()), m_rhs5(U.size()), m_rhs6(U.size()), m_output(U.size()) {}
+ virtual ~RKF45() {}
+
+ virtual void step(void const* const data=nullptr)
+ {
+ size_t& step = m_settings.step;
+ Real& t = m_settings.t;
+ Real& dt = m_settings.dt;
+ const Real& tFinal = m_settings.tFinal;
+ Real& tDump = m_settings.tDump;
+ if (!m_settings.bRestart && m_settings.dtDump > 0.0) tDump = t + m_settings.dtDump;
+
+ while(t < tFinal)
+ {
+ Real dt_next;
+ dt = (tFinal-t) < dt ? (tFinal-t) : dt;
+ _computeRHS(t, dt, data);
+ while (!_inBound(t, dt, data, dt_next)) {}
+ m_output = m_U;
+ m_U += dt*(b1*m_rhs + b3*m_rhs3 + b4*m_rhs4 + b5*m_rhs5 + b6*m_rhs6);
+ if (m_settings.dtDump > 0.0)
+ {
+ if (tDump <= (t+dt)) _dumpOutput(data);
+ }
+ else
+ if ((step+1) % m_settings.writeGranularity == 0)
+ m_kernel->write((step+1), t+dt, dt, m_U, data);
+ t += dt;
+ dt = dt_next;
+ ++step;
+
+ if (m_settings.step % m_settings.reportGranularity == 0)
+ std::printf("Time = %e;\tStep = %d\n", m_settings.t, m_settings.step);
+
+ if (m_settings.restartstep > 0 && step % m_settings.restartstep == 0)
+ serialize<Tinput>(m_U, m_settings);
+ }
+ }
+};
+
+#endif /* RKF45_H_UWN6WM1V */
diff --git a/include/TimeStepper/explicit/variableStep/RKV56.h b/include/TimeStepper/explicit/variableStep/RKV56.h
@@ -0,0 +1,195 @@
+/* File: RKV56.h */
+/* Date: Thu Feb 11 11:36:11 2016 */
+/* Author: Fabian Wermelinger */
+/* Tag: Runge-Kutta-Verner 5th-6th order variable step method */
+/* Copyright 2016 ETH Zurich. All Rights Reserved. */
+#ifndef RKV56_H_A6VN39JN
+#define RKV56_H_A6VN39JN
+
+#include <TimeStepper/explicit/Euler.h>
+
+#include <cassert>
+#include <cmath>
+#include <cstdlib>
+#include <iostream>
+#include <limits>
+
+template <typename Tinput, typename Trhs>
+class RKV56 : public Euler<Tinput, Trhs>
+{
+ using Euler<Tinput, Trhs>::m_settings;
+ using Euler<Tinput, Trhs>::m_U;
+ using Euler<Tinput, Trhs>::m_rhs;
+ using Euler<Tinput, Trhs>::m_kernel;
+
+ Real m_errold;
+ const Real m_alpha;
+ const Real m_beta;
+ bool m_reject;
+ bool m_firstPass;
+
+ Trhs m_rhs2;
+ Trhs m_rhs3;
+ Trhs m_rhs4;
+ Trhs m_rhs5;
+ Trhs m_rhs6;
+ Trhs m_rhs7;
+ Trhs m_rhs8;
+
+ Tinput m_output;
+
+ static constexpr Real b1 = 3./40., b3 = 875./2244., b4 = 23./72., b5 = 264./1955., b7 = 125./11592., b8 = 43./616.;
+ static constexpr Real c21 = 1./6., ct2 = 1./6.;
+ static constexpr Real c31 = 4./75., c32 = 16./75., ct3 = 4./15.;
+ static constexpr Real c41 = 5./6., c42 = -8./3., c43 = 2.5, ct4 = 2./3.;
+ static constexpr Real c51 = -165./64., c52 = 55./6., c53 = -425./64., c54 = 85./96., ct5 = 5./6.;
+ static constexpr Real c61 = 2.4, c62 = -8., c63 = 4015./612., c64 = -11./36., c65 = 88./255., ct6 = 1.0;
+ static constexpr Real c71 = -8263./15000., c72 = 124./75., c73 = -643./680., c74 = -81./250., c75 = 2484./10625., ct7 = 1./15.;
+ static constexpr Real c81 = 3501./1720., c82 = -300./43., c83 = 297275./52632., c84 = -319./2322., c85 = 24068./84065., c87 = 3850./26703., ct8 = 1.0;
+ static constexpr Real e1 = -1./160., e3 = -125./17952., e4 = 1./144., e5 = -12./1955., e6 = -3./44., e7 = 125./11592., e8 = 43./616.;
+
+ void _computeRHS(const Real t, const Real dt, void const* const data)
+ {
+ if (!m_reject) m_kernel->compute(m_U, m_rhs, t, data);
+ m_kernel->compute(m_U + dt*(c21*m_rhs), m_rhs2, t + ct2*dt, data);
+ m_kernel->compute(m_U + dt*(c31*m_rhs + c32*m_rhs2), m_rhs3, t + ct3*dt, data);
+ m_kernel->compute(m_U + dt*(c41*m_rhs + c42*m_rhs2 + c43*m_rhs3), m_rhs4, t + ct4*dt, data);
+ m_kernel->compute(m_U + dt*(c51*m_rhs + c52*m_rhs2 + c53*m_rhs3 + c54*m_rhs4), m_rhs5, t + ct5*dt, data);
+ 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);
+ m_kernel->compute(m_U + dt*(c71*m_rhs + c72*m_rhs2 + c73*m_rhs3 + c74*m_rhs4 + c75*m_rhs5), m_rhs7, t + ct7*dt, data);
+ 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 Real _error() const
+ {
+ assert(m_U.size() == m_rhs.size());
+ 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 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);
+ }
+ }
+ assert(!isnan(maxerr));
+ return maxerr;
+ }
+
+ bool _inBound(const Real t, Real& dt, void const* const data, Real& dt_next)
+ {
+ 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, static_cast<Real>(1.0e-4));
+ m_reject = false;
+ m_firstPass = true;
+ return true;
+ }
+
+ const Real safety = m_settings.safety;
+ const Real minScale = m_settings.minScale;
+ const Real maxScale = m_settings.maxScale;
+ Real scale;
+ if (err <= 1.0)
+ {
+ if (err == 0.0)
+ scale = maxScale;
+ else
+ {
+ scale = safety*std::pow(err,-m_alpha)*std::pow(m_errold,m_beta);
+ scale = (scale < minScale) ? minScale : ((scale > maxScale) ? maxScale : scale);
+ }
+ if (m_reject)
+ dt_next = dt*std::min(scale, static_cast<Real>(1.0));
+ else
+ dt_next = dt*scale;
+ m_errold = std::max(err, static_cast<Real>(1.0e-4));
+ m_reject = false;
+ m_firstPass = true;
+ return true;
+ }
+ else
+ {
+ scale = safety*std::pow(err,-m_alpha);
+ scale = std::max(scale, minScale);
+ dt *= scale;
+ m_reject = true;
+ m_firstPass = false;
+ _computeRHS(t, dt, data);
+ return false;
+ }
+ }
+
+ void _dumpOutput(const void * const data)
+ {
+ 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;
+
+ m_kernel->compute(m_U, m_rhs, t+dt, data);
+ m_firstPass = false;
+
+ while(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);
+ tDump += m_settings.dtDump;
+ }
+ }
+
+public:
+ RKV56(Tinput& U, StepperSettings& S, KernelBase<Tinput,Trhs> * const kern) :
+ Euler<Tinput,Trhs>(U,S,kern), m_errold(1.0e-4), m_alpha(0.25*(S.alpha-0.75*S.beta)), m_beta(0.25*S.beta), m_reject(false), m_firstPass(true),
+ m_rhs2(U.size()), m_rhs3(U.size()), m_rhs4(U.size()), m_rhs5(U.size()), m_rhs6(U.size()), m_rhs7(U.size()), m_rhs8(U.size()), m_output(U.size()) {}
+ virtual ~RKV56() {}
+
+ virtual void step(void const* const data=nullptr)
+ {
+ size_t& step = m_settings.step;
+ Real& t = m_settings.t;
+ Real& dt = m_settings.dt;
+ const Real& tFinal = m_settings.tFinal;
+ Real& tDump = m_settings.tDump;
+ if (!m_settings.bRestart && m_settings.dtDump > 0.0) tDump = t + m_settings.dtDump;
+
+ while(t < tFinal)
+ {
+ Real dt_next;
+ dt = (tFinal-t) < dt ? (tFinal-t) : dt;
+ _computeRHS(t, dt, data);
+ while (!_inBound(t, dt, data, dt_next)) {}
+ m_output = m_U;
+ m_U += dt*(b1*m_rhs + b3*m_rhs3 + b4*m_rhs4 + b5*m_rhs5 + b7*m_rhs7 + b8*m_rhs8);
+ if (m_settings.dtDump > 0.0)
+ {
+ if (tDump <= (t+dt)) _dumpOutput(data);
+ }
+ else
+ if ((step+1) % m_settings.writeGranularity == 0)
+ m_kernel->write((step+1), t+dt, dt, m_U, data);
+ t += dt;
+ dt = dt_next;
+ ++step;
+
+ if (m_settings.step % m_settings.reportGranularity == 0)
+ std::printf("Time = %e;\tStep = %d\n", m_settings.t, m_settings.step);
+
+ if (m_settings.restartstep > 0 && step % m_settings.restartstep == 0)
+ serialize<Tinput>(m_U, m_settings);
+ }
+ }
+};
+
+#endif /* RKV56_H_A6VN39JN */
diff --git a/include/TimeStepper/implicit/BDF.h b/include/TimeStepper/implicit/BDF.h
@@ -0,0 +1,77 @@
+/* File: BDF.h */
+/* Date: Mon Feb 8 11:44:23 2016 */
+/* Author: Fabian Wermelinger */
+/* Tag: Backward differentiation formula using SUNDIALS library */
+/* Copyright 2016 ETH Zurich. All Rights Reserved. */
+#ifndef BDF_H_YF523CRX
+#define BDF_H_YF523CRX
+
+#ifdef _USE_SUNDIALS_
+#include <TimeStepper/StepperBase.h>
+
+#include <cvode/cvode.h>
+#include <cvode/cvode_dense.h>
+#include <nvector/nvector_serial.h>
+
+#include <cassert>
+
+template <typename Tinput, typename Trhs>
+class BDF : public StepperBase<Tinput,Trhs>
+{
+ Real m_a;
+ void* m_cvode_mem;
+ N_Vector m_y;
+
+ void _init_sundials()
+ {
+ m_cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
+ m_y = N_VMake_Serial(Tinput::SIZE, m_U.data());
+ CVodeInit(m_cvode_mem, BDF::f, 0.0, m_y);
+ CVodeSStolerances(m_cvode_mem, 1.0e-6, 1.0e-6);
+ CVDense(m_cvode_mem, Trhs::SIZE);
+ CVodeSetUserData(m_cvode_mem, static_cast<void*>(&m_a));
+ }
+
+protected:
+ using StepperBase<Tinput,Trhs>::m_settings;
+ using StepperBase<Tinput,Trhs>::m_U;
+ using StepperBase<Tinput,Trhs>::m_kernel;
+ Trhs m_rhs;
+
+public:
+ BDF(Tinput& U, StepperSettings& S, KernelBase<Tinput,Trhs> * const kern) :
+ StepperBase<Tinput,Trhs>(U,S,kern), m_rhs(U.size())
+ {
+ _init_sundials();
+ }
+ ~BDF() { CVodeFree(&m_cvode_mem); }
+
+ virtual void step(void const* const data=nullptr)
+ {
+ /* TODO: (fabianw; Fri Feb 12 13:26:02 2016) this does not work yet.
+ * Do not compile with Sundials support */
+ /* m_a = a; */
+ Real tret; // interface dummy
+ CVode(m_cvode_mem, t+dt, m_y, &tret, CV_ONE_STEP);
+ t = tret;
+
+ if (m_settings.step % m_settings.reportGranularity == 0)
+ std::printf("Time = %e;\tStep = %d\n", m_settings.t, m_settings.step);
+ }
+
+ virtual inline void setIC(Real const ic, const Real t0=0)
+ {
+ m_U = ic;
+ CVodeReInit(m_cvode_mem, t0, m_y);
+ }
+
+ static inline int f(Real t, N_Vector y, N_Vector ydot, void *a)
+ {
+ Tkernel kernel(Trhs::SIZE);
+ kernel.compute(NV_DATA_S(y), NV_DATA_S(ydot), *static_cast<Real*>(a), t);
+ return 0;
+ }
+};
+#endif /* _USE_SUNDIALS_ */
+
+#endif /* BDF_H_YF523CRX */
diff --git a/common.h b/include/common.h
diff --git a/TimeStepper/orderVerification/.gitignore b/test/orderVerification/.gitignore
diff --git a/TimeStepper/orderVerification/Makefile b/test/orderVerification/Makefile
diff --git a/test/orderVerification/orderVerification.cpp b/test/orderVerification/orderVerification.cpp
@@ -0,0 +1,149 @@
+/* File: orderVerification.cpp */
+/* Date: Sun Feb 7 22:48:30 2016 */
+/* Author: Fabian Wermelinger */
+/* Tag: Verify order of time stepper */
+/* Copyright 2016 ETH Zurich. All Rights Reserved. */
+#include <GnuplotDump.h>
+#include <TimeStepper/KernelBase.h>
+#include <TimeStepper/TimeStepper.h>
+#include <common.h>
+
+#include <cassert>
+#include <iostream>
+#include <cmath>
+#include <vector>
+#include <string>
+#include <utility>
+#include <fstream>
+
+using namespace std;
+
+constexpr Real DECAY = -1.0;
+
+template <typename Tinput, typename Trhs=Tinput>
+class Dahlquist : public KernelBase<Tinput,Trhs>
+{
+ GnuplotDump m_gp;
+
+protected:
+ const Real m_a;
+
+public:
+ Dahlquist(const Real a=1.0) : m_gp("out", GnuplotDump::ASCII), m_a(a) {}
+ ~Dahlquist() {}
+
+ void compute(const Tinput& U, Trhs& rhs, const Real t, void const* const data=nullptr)
+ {
+ rhs = m_a*U;
+ }
+
+ void write(const size_t step, const Real t, const Real dt, const Tinput& U, const void * const data=nullptr)
+ {
+ m_gp.write(step,
+ t,
+ dt,
+ (const Real *)U.data(),
+ U.size() * Tinput::DataType::SIZE);
+ }
+};
+
+template <typename Tinput, typename Trhs=Tinput>
+class DahlquistLSRK3 : public Dahlquist<Tinput,Trhs>
+{
+ using Dahlquist<Tinput, Trhs>::m_a;
+
+public:
+ DahlquistLSRK3(const Real a=1.0) : Dahlquist<Tinput,Trhs>(a) {}
+ ~DahlquistLSRK3() {}
+
+ void compute(const Tinput& U, Trhs& rhs, const Real t, void const* const data=nullptr)
+ {
+ assert(data != nullptr);
+ const Real * const LSRKData = static_cast<const Real* const>(data);
+ const Real dt = LSRKData[0];
+ const Real A = LSRKData[1];
+ rhs = A*rhs + dt*m_a*U;
+ }
+};
+
+// types we use
+typedef State<Real,1> State1;
+using vec_t = LightVector<State1>;
+using baseKernel_t = KernelBase<vec_t>*;
+using baseStepper_t = StepperBase<vec_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) {}
+ 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)
+{
+ vector<Error> error(stages);
+ S.dt = dt;
+ for (int s = 0; s < stages; ++s)
+ {
+ size_t const N = static_cast<size_t>(S.tFinal/S.dt)+1;
+ S.t = 0.0;
+ stepper->U() = 1.0;
+ Real& numerical = stepper->U()[0][0];
+ for (size_t i = 0; i < N; ++i)
+ {
+ stepper->step();
+ 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);
+ }
+ error[s].L1 /= N;
+ error[s].L2 = sqrt(error[s].L2/N);
+ error[s].dt = S.dt;
+ S.dt *= scale;
+ }
+ return error;
+}
+
+int main(int argc, const char** argv)
+{
+ // set up solution vector
+ vec_t U;
+
+ // kernel instances
+ baseKernel_t dahlquist = new Dahlquist<vec_t>(DECAY);
+ baseKernel_t dahlquistLSRK3 = new DahlquistLSRK3<vec_t>(DECAY);
+
+ // setup time steppers
+ ArgumentParser parser(argc, argv);
+ StepperSettings S(parser);
+
+ vector<pair<pair<string,Real>, baseStepper_t> > timeStepperList;
+ timeStepperList.push_back(make_pair(make_pair("Euler",1.0), new Euler<vec_t>(U, S, dahlquist)));
+ timeStepperList.push_back(make_pair(make_pair("RK4",4.0), new RK4<vec_t>(U, S, dahlquist)));
+ timeStepperList.push_back(make_pair(make_pair("LSRK3",3.0), new LSRK3<vec_t>(U, S, dahlquistLSRK3)));
+ // timeStepperList.push_back(make_pair(make_pair("RKF45",5.0), new RKF45<vec_t>(U, S, dahlquist)));
+ // timeStepperList.push_back(make_pair(make_pair("RKV56",6.0), new RKV56<vec_t>(U, S, dahlquist)));
+#ifdef _USE_SUNDIALS_
+ timeStepperList.push_back(make_pair(make_pair("BDF",5.0), new BDF<vec_t>(U, S, dahlquist)));
+#endif /* _USE_SUNDIALS_ */
+
+ for (size_t i = 0; i < timeStepperList.size(); ++i)
+ {
+ 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 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();
+ }
+
+ for (size_t i = 0; i < timeStepperList.size(); ++i)
+ delete timeStepperList[i].second;
+ delete dahlquist;
+ delete dahlquistLSRK3;
+
+ return 0;
+}
