Refactor Model::calc_rate

model.cpp

@@ -18,19 +18,18 @@ void const Model::printParameters()
   }
 }
 
-void SFOR::calcRate(const dvector &y, dvector &dydt, double t)
+void SFOR::calcRate(const dvector &y, dvector &dydt, double t, double T)
 {
     double vol = y[0];
-    double T = y[1];
+    //double T = y[1];
     dvector par = getParameters();
     dydt[0] = par[0] * exp(-par[1]/8314.33/T) * (par[2] - vol);
-    dydt[1] = 0.0;
 }
 
-void C2SM::calcRate(const dvector &y, dvector &dydt, double t)
+void C2SM::calcRate(const dvector &y, dvector &dydt, double t, double T)
 {
     const dvector par = getParameters();
-    double T = y[2];
+    //double T = y[2];
     //double v = y[0];
 
     // fraction of solid

model.hpp

@@ -53,7 +53,7 @@ class Model
   size_t getNParameters(){return parametersDefault.size();}
   dvector const getParameters(){return parameters;}
   dvector const getParametersDefault(){return parametersDefault;}
-  virtual void calcRate(const dvector &y, dvector &dydt, double t)=0;
+  virtual void calcRate(const dvector &y, dvector &dydt, double t, double T)=0;
   std::string const getName(){return name;}
   parVector const getParametersNames(){return parametersNames;}
   void const printParameters();
@@ -71,7 +71,7 @@ class SFOR: public Model
   SFOR(): Model(sforParDefault, sforParDefault, "SFOR", sforParNames,
     sforInitState){};
   ~SFOR(){}
-  virtual void calcRate(const dvector &y, dvector &dydt, double t);
+  virtual void calcRate(const dvector &y, dvector &dydt, double t, double T);
 };
 
 //
@@ -85,7 +85,7 @@ class C2SM: public Model
   C2SM(): Model(c2smParDefault, c2smParDefault, "C2SM", c2smParNames,
     c2smInitState){};
   ~C2SM(){};
-  virtual void calcRate(const dvector &y, dvector &dydt, double t);
+  virtual void calcRate(const dvector &y, dvector &dydt, double t, double T);
 };
 
 

plotSolution.py

@@ -0,0 +1,17 @@
+import matplotlib.pyplot as plt
+import pandas as pd
+import argparse
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser()
+
+    parser.add_argument('solution', action='store')
+
+    argument = parser.parse_args()
+
+    df = pd.read_csv(argument.solution)
+
+    fig, ax = plt.subplots()
+    ax.plot('time', 'y0', data=df)
+
+    plt.show()

runner.cpp

@@ -43,9 +43,14 @@ void Runner::solve(double t, double dt)
   times = solution.m_times;
 }
 
-void Runner::dydt(const dvector y, dvector dydt, double t)
+void Runner::dydt(const dvector &y, dvector &dydt, double t)
 {
-    model->calcRate(y, dydt, t);
+    size_t n = dydt.size();
+    //double T = dydt.back();
+    //double T = y[n-1];
+    double T = y.back();
+    model->calcRate(y, dydt, t, T);
+    dydt[n-1] = 0.0; // const temperature
 }
 
 

runner.hpp

@@ -14,6 +14,9 @@
 #include <iostream>
 #include <string>
 
+// include gtest for testing purpose
+#include "gtest/gtest_prod.h"
+
 void printv(const dvector &y, const double t);
 
 struct push_back_state_and_time
@@ -37,7 +40,12 @@ class Runner
   Model * model;
   std::vector<double> times;
   std::vector<dvector> states;
-  void dydt(const dvector y, dvector dydt, double t);
+  //
+  //    calculate the derivative
+  //
+  void dydt(const dvector &y, dvector &dydt, double t);
+
+  FRIEND_TEST(RunnerTest2, dydt);
   //std::unique_ptr<Model> model;
 public:
   enum models {sfor, c2sm};

test.cpp

@@ -13,10 +13,11 @@
 // google test
 #include "gtest/gtest.h"
 
-//
 #include "model.hpp"
 #include "runner.hpp"
 
+//#define private public
+
 
 //using namespace test;
 
@@ -96,7 +97,7 @@ TEST_F(SFORTest, rate)
   double v = 0.0;
   dvector y = {v, T};
   dvector dydt(2);
-  model->calcRate(y, dydt, 0.0);
+  model->calcRate(y, dydt, 0.0, T);
   EXPECT_EQ(dydt[0], rate_sfor(parameters, T, v));
 }
 
@@ -107,7 +108,7 @@ TEST_F(C2SMTest, rate)
   double s = 1.0;
   dvector y = {v, s, T};
   dvector dydt(3);
-  model->calcRate(y, dydt, 0.0);
+  model->calcRate(y, dydt, 0.0, T);
   EXPECT_EQ(dydt[0], rate_c2sm(model->getParameters(), T, s));
 }
 
@@ -133,6 +134,7 @@ TEST_F(RunnerTest, parameters)
   EXPECT_EQ(parameters, runner->getParameters());
 }
 
+
 TEST_F(RunnerTest, solve)
 {
   //dvector y0 = {0.0, 1000};
@@ -147,6 +149,26 @@ TEST_F(RunnerTest, solve)
   EXPECT_EQ(states[0].size(), sforInitState.size()+1);
 }
 
+TEST(RunnerTest2, dydt)
+{
+    Runner runner(Runner::sfor);
+    dvector parameters = runner.getParameters();
+    double T = 1000;
+    dvector y = {0.0, T};
+    dvector dydt(y.size());
+    double rate = rate_sfor(parameters, T, 0);
+
+    // check model rate
+    double t=0;
+    runner.model->calcRate(y, dydt, t, T);
+    EXPECT_EQ(rate, dydt[0]);
+
+    //std::cout << "calc dydt" << '\n';
+    runner.dydt(y, dydt, t);
+    std::cout << "calc dydt...ok" << '\n';
+    EXPECT_EQ(rate, dydt[0]);
+}
+
 TEST(PushBackTest, pointer)
 {
   std::vector<dvector> state = {{0.0, 0.0}};