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("temperature", "C");
      }
      print $file "\n";
      $then += $gap;
   }
   return $then;
}

# disable stdout buffering
$| = 1;

print "Model heat conduction in one dimension\n\n";

# 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 $porosity = 0.25;

my $tcond = 2.0; # W/m/k
my $trans = $deltay * $deltaz * $tcond / $deltax;

my $time_end = 15.0; # years
my $delta_years = $time_end / 3.0; # years
my $next_output = 0;

# Open output file and write instance positions on the first line
my $filename = 'HeatConduction.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";

my @cp = ();
for (my $i=0; $i<$nx; $i++){
   push(@cp, new ChemPlugin());
}
$cp[0]->Console('stdout');
$cp[0]->Config('pluses = banner');

# Configue and initialize the instances
my $cmd = "span 20 C to 100 C; Na+ = 0.001 mmol/kg; Cl- = 0.001 mmol/kg; volume = ". 
($deltax * $deltay * $deltaz) .
" m3; porosity = " . $porosity .
"; time end = " . $time_end . " years;";

for (my $i=0; $i<$nx; $i++){
   $cp[$i]->Config($cmd);
   if($i<$nx/2) {
      $cp[$i]->Config('T = 100 C');
   }
   else {
      $cp[$i]->Config('T = 20 C');
   }
   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
for(my $i=1;$i<$nx;$i++){
   my $link = $cp[$i]->Link($cp[$i-1]);
   $link->HeatTrans($trans, 'W/K');
}

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;
         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_line($file, \@cp, $nx, $delta_years, $next_output);
}
