import ChemPlugin.*;  // import ChemPlugin and CpiLink class

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

class Advection1 {

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

      System.out.println("Model advection-dispersion in one dimension");
      System.out.println("");

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

      System.out.print("Please enter fluid velocity in m/yr: ");
      Scanner input = new Scanner(System.in);
      double veloc_in = input.nextDouble();
      input.nextLine();

      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.0; // years
      double next_output = 0.0; //years

      // Open output file and write instance positions on the first line
      String file = "Advection.txt";
      FileWriter fileWriter = new FileWriter(file);
      BufferedWriter bufferedWriter =  new BufferedWriter(fileWriter);
      bufferedWriter.write("years");
      for (int i=0; i<nx; i++){
         bufferedWriter.write(String.format("\t%8.4f", (i+0.5) * deltax));
      }
      bufferedWriter.newLine();

      // Configure and initialize the inlet and interior instances
      ChemPlugin cp_inlet = new ChemPlugin();
      cp_inlet.Config("Na+= 1 mmol/kg; Cl- = 1 mmol/kg");
      if (cp_inlet.Initialize() != 0)
         System.out.println("Inlet failed to initialize");

      // Interior instances
      ChemPlugin cp[] =  new ChemPlugin[nx];
      for (int i=0; i<nx; i++)
         cp[i] = new ChemPlugin();
      cp[0].Console("stdout");
      cp[0].Config("pluses = banner");

      String cmd = "Na+ = 0.001 mmol/kg; Cl- = 0.001 mmol/kg; " +
         "volume = "  + (deltax * deltay * deltaz) +
         " m3; porosity = " + porosity + "; time end = " + time_end + " year";

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

      next_output = write_line(bufferedWriter, 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");
      }

      // free outlet link
      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){
               bufferedWriter.close();
               String pause = input.nextLine();
               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].AdvanceChemical() != 0){
               System.exit(1);
            }
         }
         next_output = write_line(bufferedWriter, cp, nx, delta_years, next_output);
      } 

   }

   public static double write_line(BufferedWriter bufferedWriter, ChemPlugin cp[],
         int nx, double gap, double then ) throws IOException
   {
      double now = cp[0].Report1("Time", "years");
      if ((then - now) <= gap / 1e4) {
         bufferedWriter.write(String.format("%8.4f",now));
         for (int i =0;i<nx;i++){
            bufferedWriter.write(String.format("\t%8.4f", cp[i].Report1("concentration Na+", "mmol/kg")));
         }
         then += gap;
         bufferedWriter.newLine();
      }
      return then;

   }
}