commit 7e780512a4d3f7ab709394ffcaa96d61434bba24
parent 1bc8a2d5cbaf162ce4f10aa90342b7bbc3008f5e
Author: Fabian Wermelinger <fabianw@mavt.ethz.ch>
Date: Thu, 22 Sep 2016 20:07:25 +0200
added order verification test for time steppers
Diffstat:
3 files changed, 200 insertions(+), 0 deletions(-)
diff --git a/TimeStepper/orderVerification/.gitignore b/TimeStepper/orderVerification/.gitignore
@@ -0,0 +1,3 @@
+orderVerification
+*.dat
+*.dSYM
diff --git a/TimeStepper/orderVerification/Makefile b/TimeStepper/orderVerification/Makefile
@@ -0,0 +1,53 @@
+CC = g++
+
+debug ?= 0
+omp ?= 0
+prec ?= float
+sse ?= 0
+align ?= 16
+sundials ?= 0
+
+LIBS =
+CXXFLAGS = -std=c++11 -g -D_ALIGN_=$(align)
+CXXFLAGS += $(extra)
+
+ifeq "$(debug)" "0"
+ CXXFLAGS += -O3 -DNDEBUG
+ CXXFLAGS += -fno-expensive-optimizations -falign-functions=16
+else
+ CXXFLAGS += -O0
+endif
+ifeq "$(omp)" "1"
+ CXXFLAGS += -fopenmp -D_OMP_
+endif
+ifeq "$(prec)" "float"
+ CXXFLAGS += -D_FLOAT_PRECISION_
+ ifeq "$(sundials)" "1"
+ CXXFLAGS += -DSUNDIALS_SINGLE_PRECISION
+ endif
+else
+ CXXFLAGS += -D_DOUBLE_PRECISION_
+ ifeq "$(sundials)" "1"
+ CXXFLAGS += -DSUNDIALS_DOUBLE_PRECISION
+ endif
+endif
+ifeq "$(sse)" "1"
+ CXXFLAGS += -msse -msse2
+endif
+
+CXXFLAGS += -I../../include/ -I../../include/TimeStepper -I.
+ifeq "$(sundials)" "1"
+ CXXFLAGS += -D_USE_SUNDIALS_
+ CXXFLAGS += -I../thirdParty/sundials/build/include
+ LIBS += -L../thirdParty/sundials/build/lib -lsundials_cvode -lsundials_nvecserial
+endif
+
+HDR = $(wildcard *.h)
+CXXFLAGS += -I.. -I../..
+.PHONY := clean
+
+orderVerification: orderVerification.cpp $(HDR)
+ $(CC) $(CXXFLAGS) orderVerification.cpp -o orderVerification $(LIBS)
+
+clean:
+ rm -f *~ orderVerification
diff --git a/TimeStepper/orderVerification/orderVerification.cpp b/TimeStepper/orderVerification/orderVerification.cpp
@@ -0,0 +1,144 @@
+/* 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 double DECAY = -1.0;
+
+template <typename Tinput, typename Trhs>
+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>
+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
+using vec_t = StateVector<State1>;
+using baseKernel_t = KernelBase<vec_t,vec_t>*;
+using baseStepper_t = StepperBase<vec_t,vec_t>*;
+
+inline double exact(const double t, const double c) { return exp(c*t); }
+
+struct Error
+{
+ Error() : dt(0), L1(0), L2(0), Linf(0) {}
+ double 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 double 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,vec_t>(DECAY);
+ baseKernel_t dahlquistLSRK3 = new DahlquistLSRK3<vec_t,vec_t>(DECAY);
+
+ // setup time steppers
+ ArgumentParser parser(argc, argv);
+ StepperSettings S(parser);
+
+ vector<pair<pair<string,double>, baseStepper_t> > timeStepperList;
+ timeStepperList.push_back(make_pair(make_pair("Euler",1.0), new Euler<vec_t, vec_t>(U, S, dahlquist)));
+ timeStepperList.push_back(make_pair(make_pair("RK4",4.0), new RK4<vec_t, vec_t>(U, S, dahlquist)));
+ timeStepperList.push_back(make_pair(make_pair("LSRK3",3.0), new LSRK3<vec_t, vec_t>(U, S, dahlquistLSRK3)));
+ // timeStepperList.push_back(make_pair(make_pair("RKF45",5.0), new RKF45<vec_t, vec_t>(U, S, dahlquist)));
+ // timeStepperList.push_back(make_pair(make_pair("RKV56",6.0), new RKV56<vec_t, vec_t>(U, S, dahlquist)));
+#ifdef _USE_SUNDIALS_
+ timeStepperList.push_back(make_pair(make_pair("BDF",5.0), new BDF<vec_t, 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 double 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;
+}
