orderVerification.cpp (4926B)
1 /* File: orderVerification.cpp */ 2 /* Date: Sun Feb 7 22:48:30 2016 */ 3 /* Author: Fabian Wermelinger */ 4 /* Tag: Verify order of time stepper */ 5 /* Copyright 2016 ETH Zurich. All Rights Reserved. */ 6 #include <ODETB/GnuplotDump.h> 7 #include <ODETB/TimeStepper/KernelBase.h> 8 #include <ODETB/TimeStepper/TimeStepper.h> 9 #include <ODETB/common.h> 10 11 #include <cassert> 12 #include <iostream> 13 #include <cmath> 14 #include <vector> 15 #include <string> 16 #include <utility> 17 #include <fstream> 18 19 using namespace std; 20 21 constexpr Real DECAY = -1.0; 22 23 template <typename Tinput, typename Trhs=Tinput> 24 class Dahlquist : public KernelBase<Tinput,Trhs> 25 { 26 GnuplotDump m_gp; 27 28 protected: 29 const Real m_a; 30 31 public: 32 Dahlquist(const Real a=1.0) : m_gp("out", GnuplotDump::ASCII), m_a(a) {} 33 ~Dahlquist() {} 34 35 void compute(const Tinput &U, Trhs &rhs, const Real, void const *const) 36 { 37 rhs = m_a*U; 38 } 39 40 void write(const size_t step, 41 const Real t, 42 const Real dt, 43 const Tinput &U, 44 const void *const) 45 { 46 m_gp.write(step, 47 t, 48 dt, 49 (const Real *)U.data(), 50 U.size() * Tinput::DataType::SIZE); 51 } 52 }; 53 54 template <typename Tinput, typename Trhs=Tinput> 55 class DahlquistLSRK3 : public Dahlquist<Tinput,Trhs> 56 { 57 using Dahlquist<Tinput, Trhs>::m_a; 58 59 public: 60 DahlquistLSRK3(const Real a=1.0) : Dahlquist<Tinput,Trhs>(a) {} 61 ~DahlquistLSRK3() {} 62 63 void compute(const Tinput &U, 64 Trhs &rhs, 65 const Real, 66 void const *const data = nullptr) 67 { 68 assert(data != nullptr); 69 const Real *const LSRKData = static_cast<const Real *>(data); 70 const Real dt = LSRKData[0]; 71 const Real A = LSRKData[1]; 72 rhs = A*rhs + dt*m_a*U; 73 } 74 }; 75 76 // types we use 77 typedef State<Real,1> State1; 78 using vec_t = LightVector<State1>; 79 using baseKernel_t = KernelBase<vec_t>*; 80 using baseStepper_t = StepperBase<vec_t>*; 81 82 inline Real exact(const Real t, const Real c) { return exp(c*t); } 83 84 struct Error 85 { 86 Error() : dt(0), L1(0), L2(0), Linf(0) {} 87 Real dt, L1, L2, Linf; 88 }; 89 90 vector<Error> verify(baseStepper_t stepper, int const stages, StepperSettings& S, const Real dt=1.0, const Real scale=0.5) 91 { 92 vector<Error> error(stages); 93 S.dt = dt; 94 for (int s = 0; s < stages; ++s) 95 { 96 size_t const N = static_cast<size_t>(S.tFinal/S.dt)+1; 97 S.t = 0.0; 98 stepper->U() = 1.0; 99 Real &numerical = stepper->U()[0][0]; 100 for (size_t i = 0; i < N; ++i) 101 { 102 stepper->step(); 103 const Real absErr = abs(numerical - exact(S.t, DECAY)); 104 error[s].L1 += absErr; 105 error[s].L2 += absErr*absErr; 106 error[s].Linf = (error[s].Linf < absErr ? absErr : error[s].Linf); 107 } 108 error[s].L1 /= N; 109 error[s].L2 = sqrt(error[s].L2/N); 110 error[s].dt = S.dt; 111 S.dt *= scale; 112 } 113 return error; 114 } 115 116 int main(int argc, const char** argv) 117 { 118 // set up solution vector 119 vec_t U(1); 120 121 // kernel instances 122 baseKernel_t dahlquist = new Dahlquist<vec_t>(DECAY); 123 baseKernel_t dahlquistLSRK3 = new DahlquistLSRK3<vec_t>(DECAY); 124 125 // setup time steppers 126 ArgumentParser parser(argc, argv); 127 StepperSettings S(parser); 128 129 vector<pair<pair<string,Real>, baseStepper_t> > timeStepperList; 130 timeStepperList.push_back(make_pair(make_pair("Euler",1.0), new Euler<vec_t>(U, S, dahlquist))); 131 timeStepperList.push_back(make_pair(make_pair("RK4",4.0), new RK4<vec_t>(U, S, dahlquist))); 132 timeStepperList.push_back(make_pair(make_pair("LSRK3",3.0), new LSRK3<vec_t>(U, S, dahlquistLSRK3))); 133 // timeStepperList.push_back(make_pair(make_pair("RKF45",5.0), new RKF45<vec_t>(U, S, dahlquist))); 134 // timeStepperList.push_back(make_pair(make_pair("RKV56",6.0), new RKV56<vec_t>(U, S, dahlquist))); 135 #ifdef _USE_SUNDIALS_ 136 timeStepperList.push_back(make_pair(make_pair("BDF",5.0), new BDF<vec_t>(U, S, dahlquist))); 137 #endif /* _USE_SUNDIALS_ */ 138 139 for (size_t i = 0; i < timeStepperList.size(); ++i) 140 { 141 vector<Error> error = verify(timeStepperList[i].second, 7, S); 142 ofstream save(timeStepperList[i].first.first + ".dat"); 143 save.setf(std::ios::scientific, std::ios::floatfield); 144 auto expected = [&](const Real x){ return error[0].L1*pow(x, timeStepperList[i].first.second); }; 145 save << "# dt\t1/dt\tL1\tL2\tLinf\texpected\n"; 146 for (size_t j = 0; j < error.size(); ++j) 147 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; 148 save.close(); 149 } 150 151 for (size_t i = 0; i < timeStepperList.size(); ++i) 152 delete timeStepperList[i].second; 153 delete dahlquist; 154 delete dahlquistLSRK3; 155 156 return 0; 157 }