commit 1bc8a2d5cbaf162ce4f10aa90342b7bbc3008f5e
Author: Fabian Wermelinger <fabianw@mavt.ethz.ch>
Date: Thu, 22 Sep 2016 18:27:17 +0200
initial
Diffstat:
13 files changed, 1662 insertions(+), 0 deletions(-)
diff --git a/.gitignore b/.gitignore
@@ -0,0 +1,4 @@
+*~
+*.bak
+*.swp
+*.swo
diff --git a/ArgumentParser.h b/ArgumentParser.h
@@ -0,0 +1,375 @@
+/*
+ * ArgumentParser.h
+ * Cubism
+ *
+ * This argument parser assumes that all arguments are optional ie, each of the argument names is preceded by a '-'
+ * all arguments are however NOT optional to avoid a mess with default values and returned values when not found!
+ *
+ * More converter could be required:
+ * add as needed
+ * TypeName as{TypeName}() in Value
+ *
+ * Created by Christian Conti on 6/7/10.
+ * Copyright 2010 ETH Zurich. All rights reserved.
+ *
+ */
+
+#pragma once
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <map>
+#include <vector>
+#include <string>
+#include <iostream>
+#include <sstream>
+#include <iomanip>
+#include <fstream>
+#include <limits>
+
+
+using namespace std;
+
+class Value
+{
+private:
+ string content;
+
+public:
+
+ Value() : content("") {}
+
+ Value(string content_) : content(content_) { /*printf("%s\n",content.c_str());*/ }
+
+ double asDouble(double def=0)
+ {
+ if (content == "")
+ {
+ ostringstream sbuf;
+ sbuf << def;
+ content = sbuf.str();
+ }
+ return (double) atof(content.c_str());
+ }
+
+ int asInt(int def=0)
+ {
+ if (content == "")
+ {
+ ostringstream sbuf;
+ sbuf << def;
+ content = sbuf.str();
+ }
+ return atoi(content.c_str());
+ }
+
+ bool asBool(bool def=false)
+ {
+ if (content == "")
+ {
+ if (def) content = "true";
+ else content = "false";
+ }
+ if (content == "0") return false;
+ if (content == "false") return false;
+
+ return true;
+ }
+
+ string asString(string def="")
+ {
+ if (content == "") content = def;
+
+ return content;
+ }
+};
+
+class CommandlineParser
+{
+private:
+ const int iArgC;
+ const char** vArgV;
+ bool bStrictMode, bVerbose;
+
+protected:
+ map<string,Value> mapArguments;
+
+public:
+
+ Value& operator()(const string arg)
+ {
+ if (bStrictMode)
+ {
+ map<string,Value>::const_iterator it = mapArguments.find(arg);
+
+ if (it == mapArguments.end())
+ {
+ printf("Runtime option NOT SPECIFIED! ABORTING! name: %s\n",arg.data());
+ abort();
+ }
+ }
+
+ if (bVerbose) printf("%s is %s\n", arg.data(), mapArguments[arg].asString().data());
+ return mapArguments[arg];
+ }
+
+ bool check(const string arg) const
+ {
+ return mapArguments.find(arg) != mapArguments.end();
+ }
+
+ CommandlineParser(const int argc, const char ** argv) : mapArguments(), iArgC(argc), vArgV(argv), bStrictMode(false), bVerbose(true)
+ {
+ for (int i=1; i<argc; i++)
+ if (argv[i][0] == '-')
+ {
+ string values = "";
+ int itemCount = 0;
+
+ for (int j=i+1; j<argc; j++)
+ if (argv[j][0] == '-')
+ break;
+ else
+ {
+ if (strcmp(values.c_str(), ""))
+ values += ' ';
+
+ values += argv[j];
+ itemCount++;
+ }
+
+ if (itemCount == 0)
+ values = "true";
+
+ mapArguments[argv[i]] = Value(values);
+ i += itemCount;
+ }
+
+ mute();
+ //printf("found %ld arguments of %d\n",mapArguments.size(),argc);
+ }
+
+ int getargc() const { return iArgC; }
+
+ const char** getargv() const { return vArgV; }
+
+ void set_strict_mode()
+ {
+ bStrictMode = true;
+ }
+
+ void unset_strict_mode()
+ {
+ bStrictMode = false;
+ }
+
+ void mute()
+ {
+ bVerbose = false;
+ }
+
+ void loud()
+ {
+ bVerbose = true;
+ }
+
+ void save_options(string path=".")
+ {
+ string options;
+ for(map<string,Value>::iterator it=mapArguments.begin(); it!=mapArguments.end(); it++)
+ {
+ options+= it->first + " " + it->second.asString() + " ";
+ }
+ string filepath = (path + "/" + string("argumentparser.log"));
+ FILE * f = fopen(filepath.data(), "a");
+ if (f == NULL)
+ {
+ printf("impossible to write %s.\n", filepath.data());
+ return;
+ }
+ fprintf(f, "%s\n", options.data());
+ fclose(f);
+ }
+
+ void print_args()
+ {
+ for(map<string,Value>::iterator it=mapArguments.begin(); it!=mapArguments.end(); it++)
+ {
+ std::cout.width(50);
+ std::cout.fill('.');
+ std::cout << std::left << it->first;
+ std::cout << ": " << it->second.asString() << std::endl;
+ }
+ }
+};
+
+
+class ArgumentParser: public CommandlineParser
+{
+ typedef std::map<std::string, Value> ArgMap;
+ typedef std::map<std::string, Value*> pArgMap;
+ typedef std::map<std::string, ArgMap* > FileMap;
+
+ // keep a reference form option origin
+ ArgMap from_commandline;
+ FileMap from_files;
+ pArgMap from_code;
+
+ // helper
+ void _ignoreComments(std::ifstream& stream, const char commentChar = '#')
+ {
+ stream >> std::ws;
+ int nextchar = stream.peek();
+ while (nextchar == commentChar)
+ {
+ stream.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
+ stream >> std::ws;
+ nextchar = stream.peek();
+ }
+ }
+
+ bool _existKey(std::string& key) const
+ {
+ const std::string og_key = key;
+ bool bExist = true;
+
+ // look for both possible keys (i.e. with leading "-" and without)
+ ArgMap::const_iterator it = mapArguments.find(key);
+ if (it == mapArguments.end())
+ {
+ if (key[0] == '-') key = key.erase(0, 1);
+ else key = "-" + key;
+ it = mapArguments.find(key);
+ if (it == mapArguments.end())
+ {
+ key = og_key;
+ bExist = false;
+ }
+ }
+ return bExist;
+ }
+
+ inline std::string _stripKey(std::string key) const
+ {
+ if (key[0] == '-') key = key.erase(0, 1);
+ return key;
+ }
+
+
+public:
+ ArgumentParser(const int _argc, const char ** _argv):
+ CommandlineParser(_argc, _argv) { from_commandline = mapArguments; }
+
+ virtual ~ArgumentParser()
+ {
+ for (FileMap::iterator it = from_files.begin(); it != from_files.end(); it++)
+ delete it->second;
+ }
+
+ void readFile(const std::string filepath)
+ {
+ from_files[filepath] = new ArgMap;
+ ArgMap& myFMap = *(from_files[filepath]);
+
+ std::ifstream confFile(filepath.c_str());
+ if (confFile.is_open())
+ {
+ // read (key value) pairs from input file, ignore comments
+ // beginning with "#"
+ std::string key, val;
+ while (!confFile.eof())
+ {
+ _ignoreComments(confFile);
+ confFile >> key >> val;
+ if (_existKey(key)) continue;
+ std::pair<string, Value> item(key, Value(val));
+ mapArguments.insert(item); // add to parent container
+ myFMap.insert(item); // add to private container
+ }
+ }
+ }
+
+ Value& operator()(std::string key)
+ {
+ const bool bDefaultInCode = !_existKey(key);
+ Value& retval = CommandlineParser::operator()(key);
+ if (bDefaultInCode) from_code[key] = &retval;
+ return retval;
+ }
+
+ inline bool check(std::string key) const { return _existKey(key); }
+ inline bool exist(std::string key) const { return _existKey(key); }
+
+ void print_args()
+ {
+ std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
+ std::cout << "* Summary:" << std::endl;
+ std::cout << "* Parameter read from command line: " << from_commandline.size() << std::endl;
+ size_t nFiles = 0;
+ size_t nFileParameter = 0;
+ for (FileMap::const_iterator it=from_files.begin(); it!=from_files.end(); ++it)
+ {
+ if (it->second->size() > 0)
+ {
+ ++nFiles;
+ nFileParameter += it->second->size();
+ }
+ }
+ std::cout << "* Parameter read from " << std::setw(3) << std::right << nFiles << " file(s): " << nFileParameter << std::endl;
+ std::cout << "* Parameter read from defaults in code: " << from_code.size() << std::endl;
+ std::cout << "* Total number of parameter read from all sources: " << mapArguments.size() << std::endl;
+ std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
+
+ // command line given arguments
+ if (!from_commandline.empty())
+ {
+ std::cout << "* Command Line:" << std::endl;
+ std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
+ for(ArgMap::iterator it=from_commandline.begin(); it!=from_commandline.end(); it++)
+ {
+ std::cout.width(50);
+ std::cout.fill('.');
+ std::cout << std::left << _stripKey(it->first);
+ std::cout << ": " << it->second.asString() << std::endl;
+ }
+ std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
+ }
+
+ // options read from input files
+ if (!from_files.empty())
+ {
+ for (FileMap::iterator itFile=from_files.begin(); itFile!=from_files.end(); itFile++)
+ {
+ if (!itFile->second->empty())
+ {
+ std::cout << "* File: " << itFile->first << std::endl;
+ std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
+ ArgMap& fileArgs = *(itFile->second);
+ for(ArgMap::iterator it=fileArgs.begin(); it!=fileArgs.end(); it++)
+ {
+ std::cout.width(50);
+ std::cout.fill('.');
+ std::cout << std::left << _stripKey(it->first);
+ std::cout << ": " << it->second.asString() << std::endl;
+ }
+ std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
+ }
+ }
+ }
+
+ // defaults defined in code
+ if (!from_code.empty())
+ {
+ std::cout << "* Defined in Code:" << std::endl;
+ std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
+ for(pArgMap::iterator it=from_code.begin(); it!=from_code.end(); it++)
+ {
+ std::cout.width(50);
+ std::cout.fill('.');
+ std::cout << std::left << _stripKey(it->first);
+ std::cout << ": " << it->second->asString() << std::endl;
+ }
+ std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
+ }
+ }
+};
diff --git a/GnuplotDump.h b/GnuplotDump.h
@@ -0,0 +1,108 @@
+/* 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 double t, const double 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) : m_ftype(type), m_fname(fname)
+ {
+ if (m_ftype == ASCII)
+ {
+ m_file.open(m_fname+".dat", std::ios_base::out);
+ m_file.setf(std::ios::scientific, std::ios::floatfield);
+ }
+ else
+ m_file.open(m_fname+".bin", std::ios_base::out|std::ios_base::binary);
+ }
+ ~GnuplotDump()
+ {
+ m_file.close();
+ _writeGPscript();
+ }
+
+ enum {ASCII, BIN};
+
+ void write(const size_t step, const double t, const double 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/TimeStepper/KernelBase.h
@@ -0,0 +1,31 @@
+/* File: KernelBase.h */
+/* Date: Fri Feb 12 13:08:27 2016 */
+/* Author: Fabian Wermelinger */
+/* Tag: Kernel base class. Required for the computation of the
+ * right-hand-side in the timestepper and defines what happens when
+ * writing data. */
+/* Copyright 2016 ETH Zurich. All Rights Reserved. */
+#ifndef KERNELBASE_H_OTQSHL94
+#define KERNELBASE_H_OTQSHL94
+
+#ifdef _FLOAT_PRECISION_
+using Real = float;
+#else
+using Real = double;
+#endif
+
+template <typename Tinput, typename Trhs>
+class KernelBase
+{
+ KernelBase(KernelBase const& c) = delete;
+ KernelBase& operator=(KernelBase const& c) = delete;
+
+public:
+ KernelBase() = default;
+ virtual ~KernelBase() = default;
+
+ virtual void compute(const Tinput& U, Trhs& rhs, const Real t=0.0, void const* const data=nullptr) = 0;
+ virtual void write(const size_t step, const Real t, const Real dt, const Tinput& U, const void * const data=nullptr) = 0;
+};
+
+#endif /* KERNELBASE_H_OTQSHL94 */
diff --git a/TimeStepper/StepperBase.h b/TimeStepper/StepperBase.h
@@ -0,0 +1,163 @@
+// 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>
+class StepperBase;
+template <typename Tinput, typename Trhs>
+class Euler;
+template <typename Tinput, typename Trhs>
+class RK4;
+template <typename Tinput, typename Trhs>
+class LSRK3;
+template <typename Tinput, typename Trhs>
+class RKF45;
+template <typename Tinput, typename Trhs>
+class RKV56;
+
+#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;
+ }
+
+ double aTol;
+ double rTol;
+ double minScale;
+ double maxScale;
+ double alpha;
+ double beta;
+ double safety;
+ double t;
+ double dt;
+ size_t step;
+ size_t nsteps;
+ double tFinal;
+ double dtDump;
+ double tDump;
+ size_t writeGranularity;
+ size_t writeCount;
+ size_t reportGranularity;
+ bool bFixedStep;
+
+ // 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>
+ 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; }
+};
+
+#endif /* STEPPERBASE_H_FCCSM4HL */
diff --git a/TimeStepper/TimeStepper.h b/TimeStepper/TimeStepper.h
@@ -0,0 +1,23 @@
+// 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
@@ -0,0 +1,37 @@
+/* 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
@@ -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 "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
@@ -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 "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
@@ -0,0 +1,184 @@
+/* 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;
+
+ double m_errold;
+ const double m_alpha;
+ const double 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 double _error() const
+ {
+ assert(m_U.size() == m_rhs.size());
+ double 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;
+ 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 double 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_reject = false;
+ m_firstPass = true;
+ return true;
+ }
+
+ const double 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, 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 double t = m_settings.t;
+ const double dt = m_settings.dt;
+ double& 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 double 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;
+ double& t = m_settings.t;
+ double& dt = m_settings.dt;
+ const double& tFinal = m_settings.tFinal;
+ double& tDump = m_settings.tDump;
+ if (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);
+ }
+ }
+};
+
+#endif /* RKF45_H_UWN6WM1V */
diff --git a/TimeStepper/explicit/variableStep/RKV56.h b/TimeStepper/explicit/variableStep/RKV56.h
@@ -0,0 +1,191 @@
+/* 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;
+
+ double m_errold;
+ const double m_alpha;
+ const double 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 double _error() const
+ {
+ assert(m_U.size() == m_rhs.size());
+ double 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;
+ 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 double 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_reject = false;
+ m_firstPass = true;
+ return true;
+ }
+
+ const double 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, 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 double t = m_settings.t;
+ const double dt = m_settings.dt;
+ double& 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 double 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;
+ double& t = m_settings.t;
+ double& dt = m_settings.dt;
+ const double& tFinal = m_settings.tFinal;
+ double& tDump = m_settings.tDump;
+ if (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);
+ }
+ }
+};
+
+#endif /* RKV56_H_A6VN39JN */
diff --git a/TimeStepper/implicit/BDF.h b/TimeStepper/implicit/BDF.h
@@ -0,0 +1,76 @@
+/* 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/common.h b/common.h
@@ -0,0 +1,361 @@
+/* File: common.h */
+/* Date: Mon Dec 21 18:00:32 2015 */
+/* Author: Fabian Wermelinger */
+/* Tag: common stuff */
+/* Copyright 2015 ETH Zurich. All Rights Reserved. */
+#ifndef COMMON_H_4JKNXGVM
+#define COMMON_H_4JKNXGVM
+
+#include <cassert>
+#include <ostream>
+#include <algorithm>
+#include <cstdlib>
+#include <cstring>
+#include <cmath>
+#include <omp.h>
+
+#define LOOPSCHED
+
+#ifdef _FLOAT_PRECISION_
+using Real = float;
+#else
+using Real = double;
+#endif
+
+template <typename T>
+inline T mysqrt(T a) { abort(); return a; }
+
+template <>
+inline float mysqrt(float a) { return sqrtf(a); }
+template <>
+inline double mysqrt(double a) { return sqrt(a); }
+
+
+template <typename T, int _SIZE>
+struct State
+{
+ static constexpr int SIZE = _SIZE;
+ T s[_SIZE];
+
+ State() : s{0} {}
+ State(const State& rhs)
+ {
+ for (int i = 0; i < _SIZE; ++i)
+ s[i] = rhs.s[i];
+ }
+
+ inline T& operator[](const size_t i) { assert(i<_SIZE); return s[i]; }
+ inline T operator[](const size_t i) const { assert(i<_SIZE); return s[i]; }
+
+ inline State& operator=(const State& rhs)
+ {
+ for (int i = 0; i < _SIZE; ++i)
+ s[i] = rhs.s[i];
+ return *this;
+ }
+ inline State& operator=(const double c)
+ {
+ for (int i = 0; i < _SIZE; ++i)
+ s[i] = c;
+ return *this;
+ }
+ inline State& operator+=(const State& rhs)
+ {
+ for (int i = 0; i < _SIZE; ++i)
+ s[i] += rhs.s[i];
+ return *this;
+ }
+ inline State& operator-=(const State& rhs)
+ {
+ for (int i = 0; i < _SIZE; ++i)
+ s[i] -= rhs.s[i];
+ return *this;
+ }
+ inline State& operator*=(const State& rhs)
+ {
+ for (int i = 0; i < _SIZE; ++i)
+ s[i] *= rhs.s[i];
+ return *this;
+ }
+ inline State& operator/=(const State& rhs)
+ {
+ for (int i = 0; i < _SIZE; ++i)
+ {
+ assert(rhs.s[i] != 0);
+ s[i] /= rhs.s[i];
+ }
+ return *this;
+ }
+ inline State operator+(const State& rhs) const
+ {
+ State thisCopy(*this);
+ return thisCopy += rhs;
+ }
+ inline State operator-(const State& rhs) const
+ {
+ State thisCopy(*this);
+ return thisCopy -= rhs;
+ }
+ inline State operator*(const State& rhs) const
+ {
+ State thisCopy(*this);
+ return thisCopy *= rhs;
+ }
+ inline State operator/(const State& rhs) const
+ {
+ State thisCopy(*this);
+ return thisCopy /= rhs;
+ }
+ inline State& operator*=(const Real c)
+ {
+ for (int i = 0; i < _SIZE; ++i)
+ s[i] *= c;
+ return *this;
+ }
+ inline State operator*(const Real c) const
+ {
+ State thisCopy(*this);
+ return thisCopy *= c;
+ }
+};
+
+template <typename T, int _SIZE>
+inline State<T,_SIZE> operator*(Real const c, State<T,_SIZE> const& SV)
+{
+ State<T,_SIZE> SVcopy(SV);
+ return SVcopy *= c;
+}
+
+typedef State<Real,1> State1;
+typedef State<Real,2> State2;
+
+
+template <typename T, int _SS=0, int _SE=0>
+class LightVector
+{
+ const size_t _N;
+ T* _data;
+
+ inline T* _alloc(const size_t n)
+ {
+ void* pmem;
+ posix_memalign(&pmem, _ALIGN_, n*sizeof(T));
+ return (T*)pmem;
+ }
+
+ inline void _dealloc() { free(_data); _data=0; }
+
+ inline void _copy(const T* const src)
+ {
+#pragma omp parallel for LOOPSCHED
+ for (size_t i = 0; i < _N+_SE-_SS; ++i)
+ _data[i] = src[i];
+ }
+
+ inline void _clear()
+ {
+#pragma omp parallel for LOOPSCHED
+ for (size_t i = 0; i < _N+_SE-_SS; ++i)
+ _data[i] = 0.0;
+ }
+
+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); }
+ virtual ~LightVector() { _dealloc(); }
+
+ static const int SS = _SS;
+ static const int SE = _SE;
+ using DataType = T;
+
+ LightVector& operator=(const LightVector& rhs)
+ {
+ if (this != &rhs)
+ {
+ assert(_N+_SE-_SS == rhs._N+_SE-_SS);
+ _copy(rhs._data);
+ }
+ return *this;
+ }
+
+ LightVector& operator=(const T& c)
+ {
+#pragma omp parallel for LOOPSCHED
+ for (size_t i = 0; i < _N+_SE-_SS; ++i)
+ _data[i] = c;
+ return *this;
+ }
+
+ LightVector& operator*=(const T& c)
+ {
+#pragma omp parallel for LOOPSCHED
+ for (size_t i = 0; i < _N+_SE-_SS; ++i)
+ _data[i] *= c;
+ return *this;
+ }
+
+ LightVector& operator+=(const LightVector& c)
+ {
+#pragma omp parallel for LOOPSCHED
+ for (size_t i = 0; i < _N+_SE-_SS; ++i)
+ _data[i] += c._data[i];
+ return *this;
+ }
+
+ inline size_t size() const { return _N; }
+ inline size_t size_all() const { return _N+_SE-_SS; }
+ inline void clear() { _clear(); }
+ inline T& operator[](const size_t i) { assert(i<_N+_SE-_SS); return _data[i]; }
+ inline const T& operator[](const size_t i) const { assert(i<_N+_SE-_SS); return _data[i]; }
+ inline T& operator()(const int ix) { assert((ix-_SS >= 0)&&(ix-_SS < static_cast<int>(_N)+_SE-_SS)); return _data[ix-_SS]; }
+ inline const T& operator()(const int ix) const { assert((ix-_SS >= 0)&&(ix-_SS < static_cast<int>(_N)+_SE-_SS)); return _data[ix-_SS]; }
+ inline T* data() { return _data; }
+ inline const T * const data() const { return _data; }
+};
+
+template <typename T>
+class StateVector;
+
+template <typename T>
+inline std::ostream& operator<<(std::ostream& out, StateVector<T> const& SV);
+
+template <typename T>
+class StateVector
+{
+public:
+ StateVector(size_t const N=1) : m_state(N), m_N(N) {}
+ StateVector(StateVector const& rhs) : m_state(rhs.m_N), m_N(rhs.m_N) { _copy(rhs); }
+ ~StateVector() {}
+
+ using DataType = T;
+
+ // accessors
+ inline T& operator[](const int i)
+ {
+ assert(0 <= i);
+ assert(i < m_N);
+ return m_state[i];
+ }
+
+ inline T operator[](const int i) const
+ {
+ assert(0 <= i);
+ assert(i < m_N);
+ return m_state[i];
+ }
+
+ inline Real* const data() { return &m_state[0].s[0]; }
+ inline Real const * const data() const { return &m_state[0].s[0]; }
+
+ inline size_t size() const { return m_N; }
+
+ // operators
+ StateVector& operator=(StateVector const& rhs)
+ {
+ if (this != &rhs) _copy(rhs);
+ return *this;
+ }
+ StateVector& operator=(double const rhs)
+ {
+ for (size_t i = 0; i < m_state.size(); ++i)
+ m_state[i] = rhs;
+ return *this;
+ }
+
+ StateVector& operator+=(StateVector const& rhs)
+ {
+ for (size_t j=0; j < m_state.size(); ++j)
+ m_state[j] += rhs.m_state[j];
+ return *this;
+ }
+
+ StateVector& operator-=(StateVector const& rhs)
+ {
+ for (size_t j=0; j < m_state.size(); ++j)
+ m_state[j] -= rhs.m_state[j];
+ return *this;
+ }
+
+ StateVector& operator*=(StateVector const& rhs)
+ {
+ for (size_t j=0; j < m_state.size(); ++j)
+ m_state[j] *= rhs.m_state[j];
+ return *this;
+ }
+
+ StateVector& operator/=(StateVector const& rhs)
+ {
+ for (size_t j=0; j < m_state.size(); ++j)
+ m_state[j] /= rhs.m_state[j];
+ return *this;
+ }
+
+ inline StateVector operator+(StateVector const& rhs) const
+ {
+ StateVector thisCopy(*this);
+ return thisCopy += rhs;
+ }
+
+ inline StateVector operator-(StateVector const& rhs) const
+ {
+ StateVector thisCopy(*this);
+ return thisCopy -= rhs;
+ }
+
+ inline StateVector operator*(StateVector const& rhs) const
+ {
+ StateVector thisCopy(*this);
+ return thisCopy *= rhs;
+ }
+
+ inline StateVector operator/(StateVector const& rhs) const
+ {
+ StateVector thisCopy(*this);
+ return thisCopy /= rhs;
+ }
+
+ inline StateVector& operator*=(Real const c)
+ {
+ for (size_t i = 0; i < m_state.size(); ++i)
+ m_state[i] *= c;
+ return *this;
+ }
+
+ inline StateVector operator*(Real const c) const
+ {
+ StateVector thisCopy(*this);
+ return thisCopy *= c;
+ }
+
+ friend std::ostream& operator<< <>(std::ostream& out, StateVector const& SV);
+
+private:
+ LightVector<T> m_state;
+ size_t m_N;
+
+ inline void _copy(StateVector const& rhs)
+ {
+ assert(m_N == rhs.m_N);
+ m_N = rhs.m_N;
+ for(size_t j=0; j < m_state.size(); ++j)
+ m_state[j] = rhs.m_state[j];
+ }
+};
+
+template <typename T>
+inline std::ostream& operator<<(std::ostream& out, StateVector<T> const& SV)
+{
+ for (size_t j = 0; j < SV.m_state.size(); ++j)
+ for (size_t i = 0; i < T::SIZE; ++i)
+ out << '\t' << SV.m_state[j].s[i];
+ return out;
+}
+
+template <typename T>
+inline StateVector<T> operator*(Real const c, StateVector<T> const& SV)
+{
+ StateVector<T> SVcopy(SV);
+ return SVcopy *= c;
+}
+
+#endif /* COMMON_H_4JKNXGVM */
