#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>
#include <algorithm>
#include "ChemPlugin.h"

int exit_client(int status)
{
   std::cin.get();
   return status;
}

void write_line(std::ofstream& f, ChemPlugin *cp, int nx,
                double gap, double& then)
{
   double now = cp[0].Report1("Time", "years");
   if ((then - now) < gap / 1e4) {
      f << now;
      for (int i=0; i<nx; i++)
         f << "\t" << cp[i].Report1("concentration Na+", "mmol/kg");
      f << std::endl;
      then += gap;
   }
}

int main(int argc, char** argv) {
   std::cout << "Model advection-dispersion in one dimension"
             << std::endl << std::endl;

   // Simulation parameters.
   int nx = 400; // number of instances along x
   double length = 100; // m
   double deltax = length / nx; // m
   double deltay = 1.0, deltaz = 1.0; // m
   double porosity = 0.25;  // volume fraction

   double veloc_in; // m/yr
   std::cout << "Please enter fluid velocity in m/yr: ";
   std::cin >> veloc_in;
   std::cin.ignore();

   double velocity = veloc_in / 31557600.;  // m/s
   double flow = deltay * deltaz * porosity * velocity; // m3/s

   double diffcoef = 1e-10;  // m2/s
   double dispersivity = 1.0; // m
   double dispcoef = velocity * dispersivity + diffcoef; // m2/s
   double trans = deltay * deltaz * porosity * dispcoef / deltax;  // m3/s

   double time_end = 10.0;  // years
   double delta_years = time_end / 5; // years
   double next_output = 0.0;  // years

   // Open output file and write instance positions on the first line.
   std::ofstream f;
   f.open("Advection.txt");
   if (f.is_open()) {
      f << "years";
      for (int i=0; i<nx; i++)
         f << "\t" << (i+0.5) * deltax;
      f << std::endl;
   }
   else {
      std::cout << "Failed to open output file" << std::endl;
      return exit_client(-1);
   }

   // Configure and initialize the inlet and interior instances.
   ChemPlugin cp_inlet;
   cp_inlet.Config("Na+ = 1 mmol/kg; Cl- = 1 mmol/kg");
   if (cp_inlet.Initialize()) {
      std::cout << "Inlet failed to initialize" << std::endl;
      return exit_client(-1);
   }

   ChemPlugin *cp = new ChemPlugin[nx];
   cp[0].Console("stdout");
   cp[0].Config("pluses = banner");

   std::string cmd = "Na+ = 0.001 mmol/kg; Cl- = 0.001 mmol/kg; "
                     "volume = " + std::to_string(deltax * deltay * deltaz) +
                     " m3; porosity = " + std::to_string(porosity) +
                     "; time end = " + std::to_string(time_end) + " years";
   
   for (int i=0; i<nx; i++) {
      cp[i].Config(cmd);
      if (cp[i].Initialize()) {
         std::cout << "Instance " << i << " failed to initialize" << std::endl;
         return exit_client(-1);
      }
   }
   write_line(f, cp, nx, delta_years, next_output);

   // Link the instances.
   CpiLink link = cp[0].Link(cp_inlet);
   link.Transmissivity(trans, "m3/s");
   link.FlowRate(flow, "m3/s");

   for (int i=1; i<nx; i++) {
      link = cp[i].Link(cp[i-1]);
      link.Transmissivity(trans, "m3/s");
      link.FlowRate(flow, "m3/s");
   }

   link = cp[nx-1].Link();
   link.FlowRate(-flow, "m3/s");

   // Time marching loop.
   while (true) {
      double deltat = 1e99; 
      for (int i=0; i<nx; i++)
         deltat = std::min(deltat, cp[i].ReportTimeStep());
      for (int i=0; i<nx; i++)
         if (cp[i].AdvanceTimeStep(deltat)) return exit_client(0);
      for (int i=0; i<nx; i++)
         if (cp[i].AdvanceTransport()) return exit_client(-1);
      for (int i=0; i<nx; i++)
         if (cp[i].AdvanceChemical()) return exit_client(-1);

      write_line(f, cp, nx, delta_years, next_output);
   }

   // Never gets here.
   return 0;
}