Inspect the Robertson integration test for accuracy

[K20shores]

We pulled many of our analytical tests from this website. One of them is the Robertson.

During #576 @sjsprecious asked me to try to reduce the tolerance for test_analytical_surface_rxn, but that made me think of the tolerance we had on robertson which was 0.1. This seemed like a large tolerance.

So, I mimiced the fortran driver code in c++ with micm to test a range of solver tolerances and we found that for some reason we plateau at a relative error of 1e-1.

However, this is probably fine since the concentrations of A and B are essentially zero at the end of the simulation and so the large percent difference doesn't really affect anything. However, the lack of effect of the different tolerances mathematically was worrying.

If in the driver code, at the beginning of each time step we reset the initial concentrations and solve, we get a different distribution of relative errors that's more expected.

We need to figure out if we are actually testing what we really think we are in the Robertson analytical test or not. The branch that contains the robertson test is many_solver_tolerances. In case that gets deleted, here's the modified robertson test function

template<class BuilderPolicy>
void test_analytical_robertson(BuilderPolicy& builder, double tolerance = 1e-8)
{
  /*
   * A -> B, k1 = 0.04
   * B + B -> C + B, k2 = 3e7
   * B + C -> A + C, k3 = 1e4
   *
   * this problem is described in
   * Hairer, E., Wanner, G., 1996. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2nd
   * edition. ed. Springer, Berlin ; New York. Page 3
   *
   * solutions are provided here
   * https://www.unige.ch/~hairer/testset/testset.html
   */

  auto a = micm::Species("A");
  auto b = micm::Species("B");
  auto c = micm::Species("C");

  micm::Phase gas_phase{ std::vector<micm::Species>{ a, b, c } };

  micm::Process r1 = micm::Process::Create()
                         .SetReactants({ a })
                         .SetProducts({ Yields(b, 1) })
                         .SetRateConstant(micm::UserDefinedRateConstant({ .label_ = "r1" }))
                         .SetPhase(gas_phase);

  micm::Process r2 = micm::Process::Create()
                         .SetReactants({ b, b })
                         .SetProducts({ Yields(b, 1), Yields(c, 1) })
                         .SetRateConstant(micm::UserDefinedRateConstant({ .label_ = "r2" }))
                         .SetPhase(gas_phase);

  micm::Process r3 = micm::Process::Create()
                         .SetReactants({ b, c })
                         .SetProducts({ Yields(a, 1), Yields(c, 1) })
                         .SetRateConstant(micm::UserDefinedRateConstant({ .label_ = "r3" }))
                         .SetPhase(gas_phase);

  auto processes = std::vector<micm::Process>{ r1, r2, r3 };
  auto solver =
      builder.SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase })).SetReactions(processes).Build();
  
  // tolerances taken from https://www.unige.ch/~hairer/testset/stiff/rober/driver_rodas.f
  // i don't know why these particular tolerances are used, they just are
  int min_tolerance_power = 2;
  int max_tolerance_power = 10;
  int tolerance_factor_power = 4;
  int num_tolerance_loops = (max_tolerance_power - min_tolerance_power) * tolerance_factor_power + 1;
  double start_tolerance = std::pow(0.1, min_tolerance_power);
  double tolerance_factor = std::pow(0.1, 1.0 / tolerance_factor_power);

  double temperature = 272.5;
  double pressure = 101253.3;
  double air_density = 1e6;

  auto state = solver.GetState();

  double k1 = 0.04;
  double k2 = 3e7;
  double k3 = 1e4;

  state.SetCustomRateParameter("r1", k1);
  state.SetCustomRateParameter("r2", k2);
  state.SetCustomRateParameter("r3", k3);

  constexpr size_t N = 12;

  std::vector<std::vector<double>> model_concentrations(N + 1, std::vector<double>(3));
  std::vector<std::vector<double>> analytical_concentrations(N + 1, std::vector<double>(3));

  analytical_concentrations = { { 1, 0, 0 },
                                { 0.9664597373330035E+00, 0.3074626578578675E-04, 0.3350951640121071E-01 },
                                { 0.8413699238414729E+00, 0.1623390937990473E-04, 0.1586138422491472E+00 },
                                { 0.6172348823960878E+00, 0.6153591274639123E-05, 0.3827589640126376E+00 },
                                { 0.3368745306607069E+00, 0.2013702318261393E-05, 0.6631234556369748E+00 },
                                { 0.1073004285378040E+00, 0.4800166972571660E-06, 0.8926990914454987E+00 },
                                { 0.1786592114209946E-01, 0.7274751468436319E-07, 0.9821340061103859E+00 },
                                { 0.2031483924973415E-02, 0.8142277783356159E-08, 0.9979685079327488E+00 },
                                { 0.2076093439016395E-03, 0.8306077485067610E-09, 0.9997923898254906E+00 },
                                { 0.2082417512179460E-04, 0.8329841429908955E-10, 0.9999791757415798E+00 },
                                { 0.2083229471647004E-05, 0.8332935037760723E-11, 0.9999979167621954E+00 },
                                { 0.2083328471883087E-06, 0.8333315602809495E-12, 0.9999997916663195E+00 },
                                { 0.2083340149701284E-07, 0.8333360770334744E-13, 0.9999999791665152E+00 } };

  state.conditions_[0].temperature_ = temperature;
  state.conditions_[0].pressure_ = pressure;
  state.conditions_[0].air_density_ = air_density;

  for (int toleranceLoop = 1; toleranceLoop <= num_tolerance_loops; ++toleranceLoop) {
    double relative_tolerance = start_tolerance;
    double absolute_tolerance = 1.0e-6 * relative_tolerance;

    solver.solver_.parameters_.relative_tolerance_ = relative_tolerance;
    solver.solver_.parameters_.absolute_tolerance_ = { absolute_tolerance, absolute_tolerance, absolute_tolerance };

    for(auto& conc : model_concentrations)
    {
      conc = { 0, 0, 0 };
    }
    model_concentrations[0] = { 1, 0, 0 };
      
    state.variables_[0] = model_concentrations[0];

    // print the tolerance and number loop
    std::cout << "Loop: " << toleranceLoop << " Relative Tolerance: " << relative_tolerance << " Absolute Tolerance: " << absolute_tolerance << std::endl;

    double time_step = 1.0;
    std::vector<double> times;
    times.push_back(0);
    for (size_t i_time = 0; i_time < N; ++i_time)
    {
      state.variables_[0] = model_concentrations[0];
      double solve_time = time_step + i_time * time_step;
      times.push_back(solve_time);
      solver.CalculateRateConstants(state);
      // Model results
      double actual_solve = 0;
      while (actual_solve < time_step)
      {
        auto result = solver.Solve(time_step - actual_solve, state);
        actual_solve += result.final_time_;
      }
      model_concentrations[i_time + 1] = state.variables_[0];
      time_step *= 10;
    }

    std::vector<std::string> header = { "time", "A", "B", "C" };
    // set the name of the output to include the  tolerance loop
    std::string output_name = "robertson_model_concentrations_" + std::to_string(toleranceLoop) + ".csv";
    writeCSV(output_name, header, model_concentrations, times);
    writeCSV("robertson_analytical_concentrations.csv", header, analytical_concentrations, times);

    start_tolerance *= tolerance_factor;
  }

  // auto map = state.variable_map_;

  // size_t _a = map.at("A");
  // size_t _b = map.at("B");
  // size_t _c = map.at("C");

  // for (size_t i = 0; i < model_concentrations.size(); ++i)
  // {
  //   EXPECT_NEAR(model_concentrations[i][_a], analytical_concentrations[i][0], tolerance)
  //       << "Arrays differ at index (" << i << ", " << 0 << ")";
  //   EXPECT_NEAR(model_concentrations[i][_b], analytical_concentrations[i][1], tolerance)
  //       << "Arrays differ at index (" << i << ", " << 1 << ")";
  //   EXPECT_NEAR(model_concentrations[i][_c], analytical_concentrations[i][2], tolerance)
  //       << "Arrays differ at index (" << i << ", " << 2 << ")";
  // }
}

[K20shores]

Possibly related to #625 and #662