use strict;
use warnings;

use ChemPlugin;

sub write_results
{
   my ($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){
      for (my $i=0; $i<$nx; ++$i){
         my $label = sprintf("Instance %2d",$i);
         $cp->[$i]->PrintOutput('RTM.txt', $label);
      }
      $then += $gap;
   }
   return $then;
}

# disable output buffering
$| = 1;

print "Model reactive transport in one dimension\n\n";

# Query user for input file
print 'Enter RTM input script: ';
my $inputFile = <STDIN>;
chomp $inputFile;

#simulation parameters
my $nx = 100; # number of instances along x
my $length = 100; #m
my $deltax = $length / $nx;
my $deltay = 1.0;
my $deltaz = 1.0;

my $time_end = 10.0; #years
my $delta_years = $time_end / 5;
my $next_output = 0;

# Create the inlet and interior instances
my $cp_inlet = new ChemPlugin();
my @cp = ();
for (my $i=0; $i<$nx; $i++){
   push(@cp, new ChemPlugin());
}
$cp_inlet->Console('stdout');
$cp[0]->Console('stdout');
$cp[0]->Config('pluses = banner');

my $cmd = "volume = ". ($deltax * $deltay * $deltaz) . " m3; ".
"time end = " . $time_end . " year";

for (my $i=0; $i<$nx; $i++){
   $cp[$i]->Config($cmd);
}

my $scope = 0;
# Configure inlet, initial domain
open(my $file, '<',  $inputFile) || die("Could not open file '$inputFile' $!");
while (my $line = <$file>){
   if($line =~ /go/) {
      last;
   }
   elsif($line =~ /scope inlet/) {
      $scope = 1;
   }
   elsif($line =~ /scope initial/) {
      $scope = 2;
   }
   elsif($scope == 1){
      $cp_inlet->Config($line);
   }
   elsif($scope == 2){
      for (my $i=0; $i<$nx; $i++){
         $cp[$i]->Config($line);
      }
   }
}

if ($cp_inlet->Initialize()){
   print 'Inlet failed to initialize';
   exit;
}
$cp_inlet->PrintOutput('RTM.txt', 'Inlet fluid');

for (my $i=0; $i<$nx; $i++){
   if($cp[$i]->Initialize()){
      print "Instance $i failed to initialize";
      exit;
   }
}

$next_output = write_results(\@cp, $nx, $delta_years, $next_output);

# Query user for velocity; calculate flow rate and transmissivities
print "Please enter fluid velocity in m/yr: \n";
my $veloc_in = <STDIN>;

my $velocity = $veloc_in / 31557600;
my $porosity = $cp[0]->Report1("porosity");
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 $tcond = 2.0; #W/m/K
my $ttrans = $deltay * $deltaz * $tcond / $deltax;

# Link the instances
my $link = $cp[0]->Link($cp_inlet);
$link->Transmissivity($trans, 'm3/s');
$link->HeatTrans($ttrans, 'W/K');
$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->HeatTrans($ttrans, 'W/K');
   $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)) {
         exit;
      }
   }
   for(my $i=0; $i<$nx; $i++){
      if($cp[$i]->AdvanceTransport()) {
         exit;
      }
   }
   for(my $i=0; $i<$nx; $i++){
      if($cp[$i]->AdvanceHeatTransport()) {
         exit;
      }
   }
   for(my $i=0; $i<$nx; $i++){
      if($cp[$i]->AdvanceChemical()) {
         exit;
      }
   }
   $next_output = write_results(\@cp, $nx, $delta_years, $next_output);
}
