/***************************************************************************
 *   Copyright (C) 2014 by Jordi Farjas				   *
 ***************************************************************************/

// implicit method
// 1D radial coordinate
// substrate located on the left side

#include <iostream>
#include <cstdlib>
#include <conio.h>
#include <time.h>
#include <math.h>
#include "Combustion.h"

using namespace std;

int main(int argc, char *argv[])
{
 FILE *PF,*PF2,*PF3;
 char File_Name[_MAX_PATH];
 char File_Path[_MAX_PATH];
 errno_t err;
 // simulation parameters
 double *Temp,*alfa;			// present temperatures, K, and degree of transformation
 double *Temp_a,*alfa_a;		// temperatures, K, and degree of transformation of previous iteration
 double *Temp_n,*alfa_n;		// new evaluated temperatures, K, and degree of transformation

 double d_Avg_alpha;		// time derivative of the averaged degree of transformation
 double Avg_alpha_a;		// averaged degree of transformation of previous iteration

 int ii;		// integer related to the grid points
 int j,k,kk,ret;
 double tmp,Sum;
 double ErInc_ter;				// Maximum relative error allowed in one iteration
 time_t start, finish;
 double elapsed_time;
 double Err_Rel_Comp;
 bool aceleration;			// if true we increase the gamma value so the time step also inceases
 
 if(argc<2)
  {
   Help();
   return 0;
  }
 time( &start );

  // Input of system data
 Sample System;	// class that contains the system parameters

 if(System.ImputData(argv[1])==0) return 0;
 
 // extraction of the filename path and creation of the log file .log
 ret=SeparaNomExt(argv[1],File_Name);
 if(ret==0)
	{
		printf("\n\n\n\t\t\tWrong file name: %s\n\t\t\tProgram ended",argv[1]);
	 _getch();
	 return 0;
	}
 TreuPath(argv[1],File_Path);
 strcat_s(File_Name,".log");
 err=fopen_s(&PF2,File_Name,"wt");
 if(err)
  {
   printf("\n\n\n\t\t\tUnable to create the file: %s\n\t\t\tProgram ended",File_Name);
   _getch();
   return 0;
  }

 // Creation of the file that contains the evolution of the averaged degree of transformation
 strcpy_s(File_Name,File_Path);
 strcat_s(File_Name,"Avg_alpha.dat");
 err=fopen_s(&PF3,File_Name,"wt");
 if(err)
  {
   printf("\n\n\n\t\t\tUnable to create the file: %s\n\t\t\tProgram ended",File_Name);
   fprintf(PF2,"\n\n\n\t\t\tUnable to create the file: %s\n\t\t\tProgram ended",File_Name);
   fclose(PF2);
   return 0;
  }
 fputs("time[s] furnace_temp[K] dalpha/dt[s-1] alpha\n",PF3);

 // system data initialization
 if(System.Data_Initialization()==0)
	{
	 fclose(PF2);
	 fclose(PF3);
	  _getch();
	 return 0;
	}

 // creaction of the file that contains the evolution of the position of the propagating front

 if(System.Reaction_exists)
	{
	 strcpy_s(File_Name,File_Path);
	 strcat_s(File_Name,"PropFont.dat");
	 err=fopen_s(&PF,File_Name,"wt");
	 if(err)
	  {
	   printf("\n\n\n\t\t\tUnable to create the file: %s\n\t\t\tProgram ended",File_Name);
       fprintf(PF2,"\n\n\n\t\t\tUnable to create the file: %s\n\t\t\tProgram ended",File_Name);
       fclose(PF2);
	   return 0;
	  }
	 fputs("time[s] furnace_temp[K] alpha z_front[m] front_width[m]\n",PF);
	}

 // memory allocation
 Temp= new double [System.n_elem_tot]; 
 Temp_a= new double [System.n_elem_tot]; 
 Temp_n= new double [System.n_elem_tot]; 
 alfa= new double [System.n_elem_tot]; 
 alfa_a= new double [System.n_elem_sample]; 
 alfa_n= new double [System.n_elem_sample]; 
 if(Temp==NULL || Temp_a==NULL || Temp_n==NULL || alfa==NULL || alfa_a==NULL || alfa_n==NULL)
    {
     puts("\n\n\n\t\t\tUnable to allocate the memory\n\t\t\tProgram ended");
	 fputs("\n\n\n\t\t\tUnable to allocate the memory\n\t\t\tProgram ended",PF2);
	 fclose(PF2);
	 fclose(PF3);
	 if(System.Reaction_exists) fclose(PF);
	  _getch();
	 return 0;
    }

// List the system information into the screen and into the log file
 System.ListData(stdout);
 System.ListData(PF2);

// system state initialization
 for(ii=0;ii<System.n_elem_tot;ii++)
	{
 	 alfa[ii]=0.0;
	 Temp[ii]=System.T_min_adim;
	}
 for(ii=0;ii<System.n_elem_sample;ii++)
	{
 	 alfa[ii]=System.alfa_in;
	 alfa_a[ii]=System.alfa_in;
	 alfa_n[ii]=System.alfa_in;
	}
 
 // ignition zone for the SHS combustion mode
 for(ii=0;ii<System.n_elem_0;ii++)
	{
	 alfa[ii]=1.0;
	 Temp[ii]=System.T_max_adim;
	}

 // saves the initial profile (temperature and degree of transformation as a function of the position)
 System.RegiFront(Temp,alfa,PF2,File_Path);
 
 Avg_alpha_a=0.0;
 Err_Rel_Comp=10.0*System.Err_Rel_Max;
 for(j=0;j<_MAX_ITER_;j++)
	{
	 k=0,kk=0;
  	 // saving the result of the previous iteration
	 for(ii=0;ii<System.n_elem_sample;ii++)  alfa_a[ii]=alfa[ii];
	 for(ii=0;ii<System.n_elem_tot;ii++)  Temp_a[ii]=Temp[ii],Temp_n[ii]=Temp[ii];
	
	 // calculation of the first gess of new values of alfa and temperature
	 System.CalculaAlfaIsoterm(alfa_n,alfa_a,Temp_a);
	 System.CalculaImp(alfa_n,alfa_a,Temp_n,Temp_a);
	 aceleration=false;
	 do // iteration for the calculation of the new values of alfa and temperature
		{
		 System.CalculaAlfaNoIsoterm(alfa_n,alfa_a,Temp_n,Temp_a);
		 System.CalculaImp(alfa_n,alfa_a,Temp_n,Temp_a);
	
	 	 // calculation of the iteration error
		 System.Avg_alpha=0.0;
		 for(ii=0;ii<System.n_elem_sample;ii++)	
			{
			 System.Avg_alpha+=alfa_n[ii];
			 alfa[ii]=alfa_n[ii];
			}
		 for(ErInc_ter=0.0,Sum=0.0,ii=0;ii<System.n_elem_tot;ii++)	
			{
			 tmp=Temp_n[ii]-Temp[ii];
			 ErInc_ter+=tmp*tmp;
			 Sum+=Temp_n[ii];
			 Temp[ii]=Temp_n[ii];
			}
		 if(ErInc_ter<0.0)
			{
			 printf("\n\t\tWarning, negative relaive error %lg\n",ErInc_ter);
			 ErInc_ter=0.0;
			}
		 ErInc_ter=sqrt(ErInc_ter)/Sum;
		 System.Avg_alpha/=(double)((System.n_elem_sample));
		 k++;
		 kk++;
		 if(k>6)	// too many iterations, gamma should be reduced to improve convergence
			{
			 System.gamma/=2.0;
 			 System.gamma_subs/=2.0;
			 System.Inc_t=System.gamma*System.Inc_r*System.Inc_r;
			 for(ii=0;ii<System.n_elem_sample;ii++)  alfa[ii]=alfa_a[ii];
			 for(ii=0;ii<System.n_elem_tot;ii++)  Temp[ii]=Temp_a[ii],Temp_n[ii]=Temp_a[ii];
			 System.CalculaAlfaIsoterm(alfa_n,alfa_a,Temp_a);
			 System.CalculaImp(alfa_n,alfa_a,Temp_n,Temp_a);
			 k=0;
			 ErInc_ter=10*System.Err_Rel_Max;
			}
		 if(kk>MAXIT)
			{
			 printf("\n\t\tWarning, no convergence at the time step number %d\n",j);
			 ErInc_ter=0.0;
			}

		}while(ErInc_ter>System.Err_Rel_Max);
	 if(k<4 && System.gamma<System.gamma_max) aceleration=true;		// gamma may be increased

	 // calculation of the time dereivative of the average degree of transformation
	 d_Avg_alpha=(System.Avg_alpha-Avg_alpha_a)/System.Inc_t/System.f_time;	
	 Avg_alpha_a=System.Avg_alpha;		// saving the value of the average degree of transformation
	 System.elap_time_adim+=System.Inc_t;	// new value of dimensionless time
	 System.elap_time=System.elap_time_adim*System.f_time;	// new value of time

	 // saves the temperature and degree of transformation profiles 
	 if(System.Int_Fronts)
		{
		 if(System.Avg_alpha>System.Val_Front) System.RegiFront(Temp,alfa,PF2,File_Path);
		 if(System.Avg_alpha>System.Val_N_Point_r)
			{ // saves the temporal evolution of the degree of transformation
			 fprintf(PF3,"%17.11lg\t%17.11lg\t%lg\t%lg\n",
				 System.elap_time,System.T_furnace*System.f_Temp,d_Avg_alpha,System.Avg_alpha);
			 System.Val_N_Point_r+=System.Inc_N_Point_r;
			}
		}
	 else
		{
		 if(System.elap_time>System.Val_Front) System.RegiFront(Temp,alfa,PF2,File_Path);
		 if(System.elap_time>System.Val_N_Point_r)
			{ // saves the temporal evolution of the degree of transformation
			 fprintf(PF3,"%17.11lg\t%17.11lg\t%lg\t%lg\n",
				 System.elap_time,System.T_furnace*System.f_Temp,d_Avg_alpha,System.Avg_alpha);
			 System.Val_N_Point_r+=System.Inc_N_Point_r;
			}
		}

	 // information of the evolution of the calculation (on the screen)
	 if(j%System.N_Point_Calc==0) printf("%d %lg %lg %lg %d %lg\n",
				 j+1,System.elap_time,System.T_furnace*System.f_Temp,System.Avg_alpha,k,ErInc_ter);

	 // Locates the position of the propagation front if any
	// Prop_Rgt=true searches a front propagating rightwards, false leftwards
	 if(j%System.N_Point_Calc==0) System.RegPosFront(PF,alfa,true);

	 // check for calculation completition
	 if(System.elap_time>System.Tot_time) j=_MAX_ITER_+1;
	 if(System.Avg_alpha>=__FRAC_MAX__) j=_MAX_ITER_+1;
  
	 // temperature increase
	 System.T_furnace+=System.Vel_esc_adim*System.Inc_t;
	 if(aceleration) 	// fast convergence, gamma is increased
		{
		 System.gamma*=1.1;
		 System.gamma_subs*=1.1;
		 System.Inc_t=System.gamma*System.Inc_r*System.Inc_r;
		}
	}

// saves the final profile (temperature and degree of transformation as a function of the position)

 System.RegiFront(Temp,alfa,PF2,File_Path);

 fprintf(PF3,"%17.11lg\t%17.11lg\t%lg\t%lg\n",
	 System.elap_time,System.T_furnace*System.f_Temp,d_Avg_alpha,System.Avg_alpha);
 fclose(PF3);
 
 time( &finish );
 elapsed_time = difftime( finish, start );
 ii=(int)floor(elapsed_time/3600.0);
 ret=(int)floor((elapsed_time-3600.0*ii)/60.0);
 k=(int)(elapsed_time-60.0*ret-3600.0*ii);
 fprintf(PF2,"\n\tThe program lasts: %d h, %d m i %d  s\n",ii,ret,k);
 printf("\n\tThe program lasts: %d h, %d m i %d  s\n",ii,ret,k);

 fclose(PF3);
 if(System.Reaction_exists) fclose(PF);
 fclose(PF2);
 delete [] alfa_n;
 delete [] alfa_a;
 delete [] alfa;
 delete [] Temp_n;
 delete [] Temp_a;
 delete [] Temp;
 
 puts("\n\t\t\tExecution ended normally\n");
 puts("\n\tPlease press a key to end the program.");
 _getch();
 return EXIT_SUCCESS;
}
