import ChemPlugin.ChemPlugin;
import ChemPlugin.CpiLink;

import java.util.Scanner;
import java.io.*;


public class RTM1 {

   public static void main(String[] args) throws IOException {

      System.out.println("Model reactive transport in one dimension");
      System.out.println("");
      Scanner input = new Scanner(System.in);
      System.out.print("Enter RTM input script: ");
      String inputFile = input.nextLine();

      int nx = 100; // number of instances along x
      double length = 100; // m
      double deltax = length / nx; // m
      double deltay = 1.0, deltaz = 1.0;

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

      // Create the inlet an d interior instances
      ChemPlugin cp_inlet = new ChemPlugin();
      ChemPlugin cp[] = new ChemPlugin[nx];
      for (int i=0;i<nx;i++)
         cp[i] = new ChemPlugin();
      cp_inlet.Console("stdout");
      cp[0].Console("stdout");
      cp[0].Config("pluses = banner");

      String cmd = "volume = " + (deltax * deltay * deltaz) + " m3;" 
         + "time end = " + time_end + " years";

      for(int i=0; i<nx; i++)
         cp[i].Config(cmd);

      //	  		System.out.println("Reading from file " + inputFile + "\n");
      String line = null;
      int scope = 0;

      FileReader fileReader = null;
      fileReader = new FileReader(inputFile);
      BufferedReader bufferedReader = new BufferedReader(fileReader);

      while((line = bufferedReader.readLine()) != null){
         if (line.equals("go"))
            break;
         else if (line.equals("scope inlet"))
            scope = 1;
         else if (line.equals("scope initial"))
            scope = 2;
         else if (scope == 1){
            cp_inlet.Config(line);
            //			    	System.out.println("in inlet");
         }

         else if (scope == 2)
            for(int i=0; i<nx; i++)
               cp[i].Config(line);
      }
      bufferedReader.close();

      if (cp_inlet.Initialize() != 0){
         System.out.println("Inlet failed to initialize");
         System.exit(1);
      }

      cp_inlet.PrintOutput("RTM.txt", "Inlet fluid");

      for (int i=0; i<nx; i++){
         if (cp[i].Initialize() != 0){
            System.out.println("Instance" + i + "failed to initialize");
            System.exit(1);
         }
      }

      next_output = write_results(cp, nx, delta_years, next_output);

      // Query user for velocity; calculate flow rate and transmissivities

      System.out.println("Please enter fluid velocity in m/yr: ");
      double veloc_in = input.nextDouble(); 
      input.close();
      double velocity = veloc_in / 31557600;
      double porosity = cp[0].Report1("porosity");
      double flow = deltay * deltaz * porosity * velocity;

      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 tcond = 2.0; // W/m/K
      double ttrans = deltay * deltaz * tcond / deltax; 

      // Link the instances 
      CpiLink link = cp[0].Link(cp_inlet);
      link.Transmissivity(trans, "m3/s");
      link.HeatTrans(ttrans, "W/K");
      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.HeatTrans(ttrans, "W/K");
         link.FlowRate(flow, "m3/s");
      }

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

      double deltat;
      while (true) {
         deltat = 1e99;
         for(int i=0; i<nx; i++)
            deltat = Math.min(deltat, cp[i].ReportTimeStep());

         for(int i=0; i<nx; i++){
            if (cp[i].AdvanceTimeStep(deltat) != 0){
               System.exit(0);
            }
         }

         for(int i=0; i<nx; i++){
            if (cp[i].AdvanceTransport() != 0){
               System.exit(1);
            }
         }

         for(int i=0; i<nx; i++){
            if (cp[i].AdvanceHeatTransport() != 0){
               System.exit(1);
            }
         }

         for(int i=0; i<nx; i++){
            if (cp[i].AdvanceChemical() != 0){
               System.exit(1);
            }
         }
         next_output = write_results(cp, nx, delta_years, next_output);
      } 
   }

   public static double write_results(ChemPlugin[] cp, int nx, double gap, double then){

      double now = cp[0].Report1("Time", "years");
      if ((then - now) < gap / 1e4) {
         for (int i =0;i<nx;i++){
            String label = String.format( "Instance %2d", i);
            cp[i].PrintOutput("RTM.txt", label);
         }
         then += gap;
      }
      return then;
   }

}