use strict;
use warnings;

use ChemPlugin;

sub write_line
{
   my ($file, $cp, $nx, $gap, $then)  = @_; 
#open(my $file, '>', $filename) || die("Could not open file '$filename' $!");
   my $now = $cp->[0]->Report1('Time', 'years');
   if($then-$now <= $gap/1e4){
      printf $file "%8.4f", $now;
      for (my $i=0; $i<$nx; ++$i){
         printf $file "\t%8.4f", $cp->[$i]->Report1("concentration Na+", "mmol/kg");
      }
      print $file "\n";
      $then += $gap;
   }
   return $then;
}

# disable stdout buffering
$| = 1;

print "Model advection-dispersion in one dimension\n\n";

# Simulation parameters
my $nx = 400; # number of instances along x
my $length = 100; # m
my $deltax = $length / $nx; 
my $deltay = 1.0;
my $deltaz = 1.0;
my $porosity = 0.25;

print "Please enter fluid velocity in m/yr: ";
my $veloc_in = <STDIN>;
my $velocity = $veloc_in / 31557600;

my $flow = $deltay * $deltaz * $porosity * $velocity; #m3/s

my $diffcoef = 1e-10; # m2/s
my $dispersivity = 1.0; # m
my $dispcoef = $velocity * $dispersivity + $diffcoef; # m2/s
my $trans = $deltay * $deltaz * $porosity * $dispcoef / $deltax; # m3/s

my $time_end = 10.0; # years
my $delta_years = $time_end / 5.0; # years
my $next_output = 0;

# Open output file and write instance positions on the first line
my $filename = 'Advection.txt';
open(my $file, '>', $filename) || die("Could not open file '$filename' $!");
print $file 'years';
for (my $i=0; $i<$nx; $i++){
   printf $file "\t%8.4f", (($i+0.5) * $deltax);
}
print $file "\n";

# Configue and initialize the inlet and interior instances
my $cp_inlet = new ChemPlugin();
$cp_inlet->Config('Na+ = 1 mmol/kg; Cl- = 1 mmol/kg');
if ($cp_inlet->Initialize()){
   print 'Inlet failed to initialize';
   exit;
}

# Array of interior instances
my @cp = ();
for (my $i=0; $i<$nx; $i++){
   push(@cp, new ChemPlugin());
}
$cp[0]->Console('stdout');
$cp[0]->Config('pluses = banner');

my $cmd = "Na+ = 0.001 mmol/kg; Cl- = 0.001 mmol/kg; ".
"volume = ". ($deltax * $deltay * $deltaz) .
" m3; porosity = " . $porosity .
"; time end = " . $time_end . " year";

for (my $i=0; $i<$nx; $i++){
   $cp[$i]->Config($cmd);
   if($cp[$i]->Initialize()){
      print "Instance $i failed to initialize";
      exit;
   }
}

$next_output = write_line($file, \@cp, $nx, $delta_years, $next_output);

# Link the instances
my $link = $cp[0]->Link($cp_inlet);
$link->Transmissivity($trans, 'm3/s');
$link->FlowRate($flow, 'm3/s');

for(my $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');


while(1){
   my $deltat = 1e99;
   for(my $i=0; $i<$nx; $i++){
      my $time_step = $cp[$i]->ReportTimeStep();
      $deltat = ($deltat < $time_step ? $deltat : $time_step);
   }
   for(my $i=0; $i<$nx; $i++){
      if($cp[$i]->AdvanceTimeStep($deltat)) {
	     close $file;
         my $pause = <STDIN>;
         exit;
      }
   }
   for(my $i=0; $i<$nx; $i++){
      if($cp[$i]->AdvanceTransport()) {
         exit;
      }
   }
   for(my $i=0; $i<$nx; $i++){
      if($cp[$i]->AdvanceChemical()) {
         exit;
      }
   }
   $next_output = write_line($file, \@cp, $nx, $delta_years, $next_output);
}
