Change main.cpp to run with template

.gitignore

@@ -1,3 +1,7 @@
 build/
 .stdafx.c++.pch
 compiling_error.txt
+.idea
+cmake-build*
+.clang_complete
+solution.csv

input.yml

@@ -5,5 +5,5 @@ model:
 
 operating_conditions:
     T: 1000
-    t: 0.1
+    t: 0.01
     dt: 1e-6

main.cpp

@@ -8,82 +8,80 @@
 
 #include <iostream>
 
-#include <cmath>
-#include <string>
-#include <boost/numeric/odeint.hpp>
+// #include <cmath>
+// #include <string>
+// #include <boost/numeric/odeint.hpp>
 
 // yaml reader
 #include <yaml-cpp/yaml.h>
 
-
 // #include <functional>
-
-
 #include "model.hpp"
 #include "SFOR.hpp"
 #include "C2SM.hpp"
-#include "reactor.hpp"
+#include "reactort.hpp"
 #include <string.h>
 
+//template <typename T>
+//void runReactor(ReactorT<T> reactor)
 
 int main(int argc, char **argv)
 {
-  //dvector parameters = {1e6, 60e6, 0.1};
-  std::cout << "Main program" << std::endl;
-
-
-  //YAML::Node config = YAML::LoadFile(argv[1]);
-  YAML::Node config = YAML::LoadFile(argv[1]);
-
-  Reactor::models model;
-  std::string modelName = config["model"]["type"].as<std::string>();
-  dvector parameters = config["model"]["parameters"].as<dvector>();
-
-
-
-  if(modelName==std::string("SFOR"))
-  {
-      model = Reactor::sfor;
-  }
-  else if(modelName==std::string("C2SM"))
-  {
-      model = Reactor::c2sm;
-  }
-  else
-  {
-      std::cout << "Model " << argv[1] << " not defined" << std::endl;
-      return -1;
-    }
-
-  // set runner with SFOR=0
-  std::cout << "Init Reactor with " << modelName <<" - " << model << std::endl;
-  Reactor run(model);
-
-  std::cout << "Parameters:";
-  //dvector parameters = run.getParameters();
-  //for (auto p: parameters)
-  //{
-  //  std::cout << p << "\t";
-  //}
-  //std::cout << std::endl;
-  run.printParameters();
-
-  // define parameters
-  //dvector y = {0, 1000};
-  double t = config["operating_conditions"]["t"].as<double>();
-  double T = config["operating_conditions"]["T"].as<double>();
-  //double t0 = 0.0;
-  double dt = config["operating_conditions"]["dt"].as<double>();
-
-  // solve
-  std::cout << "Solve" << std::endl;
-  run.solve(t, T, dt, true);
-
-  // get solutions
-  std::vector<double> times = run.getTimes();
-  std::vector<dvector> states = run.getStates();
-
-  run.dump("solution.csv");
-
-  return 0;
+    //dvector parameters = {1e6, 60e6, 0.1};
+    std::cout << "Main program" << std::endl;
+
+    //YAML::Node config = YAML::LoadFile(argv[1]);
+    std::cout << "Read YAML file";
+    YAML::Node config = YAML::LoadFile(argv[1]);
+    std::cout << "...done" << '\n';
+    // read model and parameters
+    std::string modelName = config["model"]["type"].as<std::string>();
+    dvector parameters = config["model"]["parameters"].as<dvector>();
+
+    // // define reactor object
+    // if(modelName==std::string("SFOR"))
+    // {
+    //     //model = Reactor::sfor;
+    //     std::cout << "Create reactor with SFOR" << '\n';
+    //     ReactorT<SFOR> reactor(parameters);
+    //     reactor.printParameters();
+    // }
+    // else if(modelName==std::string("C2SM"))
+    // {
+    //     //model = Reactor::c2sm;
+    //     std::cout << "Create reactor with C2SM" << '\n';
+    //     ReactorT<C2SM> reactor(parameters);
+    // }
+    // else
+    // {
+    //     std::cout << "Model " << argv[1] << " not defined" << std::endl;
+    //     return -1;
+    // }
+    std::cout << "Init Reactor with SFOR" << '\n';
+    ReactorT<SFOR> reactor(parameters);
+
+    // set runner with SFOR=0
+    //std::cout << "Init Reactor with " << modelName <<" - " << model <<
+    //std::endl;
+    //Reactor run(model);
+
+    reactor.printParameters();
+
+    // define parameters
+    //dvector y = {0, 1000};
+    double t = config["operating_conditions"]["t"].as<double>();
+    double T = config["operating_conditions"]["T"].as<double>();
+    //double t0 = 0.0;
+    double dt = config["operating_conditions"]["dt"].as<double>();
+
+    // solve
+    std::cout << "Solve" << std::endl;
+    reactor.solve(t, T, dt, false);
+
+    // get solutions
+    //std::vector<double> times = reactor.getTimes();
+    //std::vector<dvector> states = reactor.getStates();
+    reactor.dump("solution.csv");
+
+    return 0;
 }

plot_solution.py

@@ -0,0 +1,14 @@
+import matplotlib.pyplot as plt
+import pandas as pd
+from dualplot import DualPlot
+
+if __name__ == '__main__':
+    data = pd.read_csv('solution.csv')
+    #plt.plot('time', 'y0', data=data, marker='x')
+    dplot = DualPlot()
+    dplot.plot_left('time', 'y0', data=data, marker='x')
+    dplot.plot_right('time', 'Temperature',  data=data)
+    dplot.set_ylabel_left('Volatile yield')
+    dplot.set_ylabel_right('Temp, K')
+    dplot.set_xlabel('Time')
+    plt.show()

reactort.hpp

@@ -13,6 +13,7 @@
 #include <fstream>
 
 #include "model.hpp"
+#include "csv.hpp"
 
 // include gtest for testing purpose
 #include "gtest/gtest_prod.h"
@@ -48,13 +49,14 @@ class ReactorT
     void dTdt(const dvector &y, dvector &dydt, double t);
 
     dvector getModelParameters();
+    void setModelParameters(const dvector &parameters);
 
     void solve(double t, double temperature=1000, double dt=1e-4, bool verbose=false);
     //void solve(std::vector<std::vector<double>> points, double dt=1e-4,
     //    bool verbose=false);
-    //std::vector<double> getTimes(){return times;}
-    //std::vector<dvector> getStates(){return states;}
-    //void dump(const std::string &csv, std::string sep = ",");
+    std::vector<double> getTimes(){return times;}
+    std::vector<dvector> getStates(){return states;}
+    void dump(const std::string &csv, std::string sep = ",");
     void printParameters();//{model.printParameters();}
 };
 
@@ -69,8 +71,9 @@ ReactorT<T>::ReactorT(): isoThermal(true)
 template <class T>
 ReactorT<T>::ReactorT(const dvector & parameters): isoThermal(true)
 {
-    std::cout << "Construct ReactorT(parameters)" << '\n';
-    model.setParameters(parameters);
+    //std::cout << "Construct ReactorT(parameters)" << '\n';
+    //model.setParameters(parameters);
+    setModelParameters(parameters);
 }
 
 template <class T>
@@ -102,6 +105,12 @@ void ReactorT<T>::dTdt(const dvector &y, dvector &dydt, double t)
         //dydt.emplace(dydt.end()-1, 0.0);
 }
 
+template<class T>
+void ReactorT<T>::setModelParameters(const dvector &parameters)
+{
+    model.setParameters(parameters);
+}
+
 template <class T>
 void ReactorT<T>::solve(
     double t,
@@ -123,10 +132,10 @@ void ReactorT<T>::solve(
 
     //double last=1;
 
-    controlled_stepper_type stepper = make_controlled<error_stepper_type >(1.0e-10 , 1.0e-6 );
+    controlled_stepper_type stepper = make_controlled<error_stepper_type>(1.0e-10 , 1.0e-6 );
     //
     std::for_each(
-        make_adaptive_time_iterator_begin(stepper, fct, y0, 0.0, 0.1, 1e-3),
+        make_adaptive_time_iterator_begin(stepper, fct, y0, 0.0, t, dt),
         make_adaptive_time_iterator_end(stepper, fct, y0),
         [&tsaved, &ysaved, verbose](std::pair< const dvector & , const double & > x )
         {
@@ -135,6 +144,32 @@ void ReactorT<T>::solve(
             tsaved.push_back(x.second);
             ysaved.push_back(x.first);
           });
+
+    times = tsaved;
+    states = ysaved;
+}
+
+template <class T>
+void ReactorT<T>::dump(const std::string &csv, std::string sep)
+{
+    /// Dump results in a file
+    csvfile file(csv);
+    file << "time" ;
+    char buffer[2];
+    for (int j=0; j< states[0].size()-1; ++j)
+    {
+        std::sprintf (buffer, "y%d", j);
+        file << std::string(buffer);
+    }
+    file << "Temperature" << endrow;
+
+    for (int i=0; i < times.size(); ++i)
+    {
+        file << times[i];
+        for (double &status: states[i])
+            file << status;
+        file << endrow;
+    }
 }
 
 #endif /* model_hpp */

test.cpp

@@ -206,9 +206,6 @@ TEST(PushBackTest, parameters)
 }
 
 
-
-
-
 int main(int argc, char **argv)
 {
   ::testing::InitGoogleTest(&argc, argv);