#include <iostream>
#include "ChemPlugin.h"


int main(int argc, char** argv) {
   std::cout << "Model a flow-through reactor" << std::endl << std::endl;

   // Configure and initialize the inlet fluid.
   ChemPlugin cp_inlet("stdout");
   const char *cmds = "Ca++ = 1 mmol/kg; HCO3- = 1 mmol/kg;"
                "pH = 1; balance on Cl-";
   cp_inlet.Config(cmds);
   cp_inlet.Initialize();

   // Configure and initialize the stirred reactor.
   ChemPlugin cp("stdout");
   cmds = "swap Calcite for Ca++; Calcite = 0.03 free m3; Cl- = 1 mmol/kg;"
          "swap CO2(g) for H+; fugacity CO2(g) = 1; balance on HCO3-;"
          "volume = 1 m3; fix f CO2(g); delxi = 0.01; pluses = banner";
   cp.Config(cmds);
   cp.Initialize(1.0, "day");

   // Link reactor to inlet and free outlet; set rate of flow.
   CpiLink link1 = cp.Link(cp_inlet);
   CpiLink link2 = cp.Outlet();

   link1.FlowRate(10.0, "m3/day");
   link2.FlowRate(-10.0, "m3/day");

   // Time marching loop.
   cp.PlotHeader("FlowThrough.gtp", "char");
   cp.PlotBlock();
   while (true) {
      double deltat = cp.ReportTimeStep();
      if (cp.AdvanceTimeStep(deltat)) break;
      if (cp.AdvanceTransport()) break;
      if (cp.AdvanceChemical()) break;
      cp.PlotBlock();
   }

   // Any keystroke closes the console.
   std::cin.get();
   return 0;
}