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

// Alternating direction implicit method.
// cylindrical coordinates
// Sourrounding vessel at the bottom and right sides

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

using namespace std;

#define __TEMP_BUSC_MAX__  // if defined in bmp files colors are graded according present maximum and minimum temperatures
						   // if not defined colors are graded according absolute maximum and minimum temperatures	

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_2,**alfa_2;		// temperatures, K, and degree of transformation at half 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,jj;						// integers related to the radial and axial 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
 double Temp_Max,Temp_Min;	// amximum and minimum temperatures for color in bmp files

 
 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;
 
 // exttarction 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
 // first column time in second
 // second and third column radial and radial width of the front, meters
 // second and tird column axial and axial width of the front, meters

 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 r_front[m] r_width[m] z_front[m] z_width[m]\n",PF);
	}

 // memory allocation
 Temp= new double * [System.n_r_tot]; 
 Temp_a= new double * [System.n_r_tot];
 Temp_2= new double * [System.n_r_tot];
 Temp_n= new double * [System.n_r_tot]; 
 alfa= new double * [System.n_r_tot]; 
 alfa_a= new double * [System.n_r_sample]; 
 alfa_2= new double * [System.n_r_sample]; 
 alfa_n= new double * [System.n_r_sample]; 
 if(Temp==NULL || Temp_a==NULL || Temp_2==NULL || Temp_n==NULL ||alfa==NULL || alfa_a==NULL || alfa_n==NULL || alfa_2==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;
    }

 for(ii=0;ii<System.n_r_tot;ii++)
    {
     Temp[ii] = new double [System.n_z_tot];
	 Temp_a[ii] = new double [System.n_z_tot];
	 Temp_2[ii] = new double [System.n_z_tot];
 	 Temp_n[ii] = new double [System.n_z_tot];
	 alfa[ii] = new double [System.n_z_tot];
     if(Temp[ii]==NULL || Temp_a[ii]==NULL || Temp_2[ii]==NULL || Temp_n[ii]==NULL || alfa[ii]==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;
        }
    }
 for(ii=0;ii<System.n_r_sample;ii++)
    {
	 alfa_a[ii] = new double [System.n_z_sample];
	 alfa_2[ii] = new double [System.n_z_sample];
 	 alfa_n[ii] = new double [System.n_z_sample];
     if(alfa_a[ii]==NULL || alfa_2[ii]==NULL || alfa_n[ii]==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_r_tot;ii++)
 	 for(jj=0;jj<System.n_z_tot;jj++)
		{
		 alfa[ii][jj]=0.0;
		 Temp[ii][jj]=System.T_min_adim;
		}
 for(ii=0;ii<System.n_r_sample;ii++)
	{
 	 for(jj=0;jj<System.n_z_sample;jj++)
		{
		 alfa[ii][jj]=System.alfa_in;
		 alfa_a[ii][jj]=System.alfa_in;
		 alfa_n[ii][jj]=System.alfa_in;
		}
	}

 // ignition zone for the SHS combustion mode
 for(ii=0;ii<System.n_r_0;ii++)
	{
	 for(jj=0;jj<System.n_z_0;jj++)
		{
		 alfa[ii][jj]=1.0;
		 Temp[ii][jj]=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);
 System.RegBMP(alfa,"alpha",File_Path,1.0,System.alfa_in);
#ifdef __TEMP_BUSC_MAX__
 Temp_Max=Temp[0][0];
 Temp_Min=Temp[0][0];
 for(ii=0;ii<System.n_r_sample;ii++)
		 for(jj=0;jj<System.n_z_sample;jj++)
			{
			 if(Temp[ii][jj]>Temp_Max) Temp_Max=Temp[ii][jj];
 			 if(Temp[ii][jj]<Temp_Min) Temp_Min=Temp[ii][jj];
			}
 System.RegBMP(Temp,"Temp",File_Path,Temp_Max,Temp_Min);
#else
 System.RegBMP(Temp,"Temp",File_Path,System.T_max_adim,System.T_min_adim);
#endif
 
 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_r_sample;ii++)
		 for(jj=0;jj<System.n_z_sample;jj++) alfa_a[ii][jj]=alfa[ii][jj];
	 for(ii=0;ii<System.n_r_tot;ii++)
		 for(jj=0;jj<System.n_z_tot;jj++)
			 Temp_a[ii][jj]=Temp[ii][jj],Temp_n[ii][jj]=Temp[ii][jj];
	
	 // calculation of the first gess of new values of alfa i temperatura
	 System.CalculaAlfaIsoterm(alfa_2,alfa_a,Temp_a);
	 System.CalculaImp_r(alfa_2,alfa_a,Temp_2,Temp_a);
	 System.CalculaAlfaIsoterm(alfa_n,alfa_2,Temp_2);
	 System.CalculaImp_z(alfa_n,alfa_2,Temp_n,Temp_2);
	 aceleration=false;
	 do // iteration for the calculation of the new values of alfa i temperatura
		{
		 System.CalculaAlfaNoIsoterm(alfa_2,alfa_a,Temp_2,Temp_a);
		 System.CalculaImp_r(alfa_2,alfa_a,Temp_2,Temp_a);
	
		 System.CalculaAlfaNoIsoterm(alfa_n,alfa_2,Temp_n,Temp_2);
		 System.CalculaImp_z(alfa_n,alfa_2,Temp_n,Temp_2);
	
	 	 // calculation of the iteration error
		 System.Avg_alpha=0.0;
		 for(ii=0;ii<System.n_r_sample;ii++)	
			 for(jj=0;jj<System.n_z_sample;jj++)
				{
				 tmp=(1.0+2.0*ii);
				 System.Avg_alpha+=tmp*alfa_n[ii][jj];
				 alfa[ii][jj]=alfa_n[ii][jj];
				}
		 for(ErInc_ter=0.0,Sum=0.0,ii=0;ii<System.n_r_tot;ii++)	
			 for(jj=0;jj<System.n_z_tot;jj++)
				{
				 tmp=Temp_n[ii][jj]-Temp[ii][jj];
				 ErInc_ter+=tmp*tmp;
				 Sum+=Temp_n[ii][jj];
				 Temp[ii][jj]=Temp_n[ii][jj];
				}
		 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_r_sample)*(System.n_r_sample)*(System.n_z_sample));
		 k++;
		 kk++;
		 if(k>6)	// too many iterations, gamma should be reduced to improve convergence
			{
			 System.gamma/=2.0;
			 System.gamma_2/=2.0;
 			 System.gamma_subs/=2.0;
			 System.gamma_subs_2/=2.0;
			 System.Inc_t_2=System.gamma_2*System.Inc_s*System.Inc_s;
			 System.Inc_t=2.0*System.Inc_t_2;
			 for(ii=0;ii<System.n_r_sample;ii++)	
				for(jj=0;jj<System.n_z_sample;jj++)
					 alfa[ii][jj]=alfa_a[ii][jj];
			 for(ii=0;ii<System.n_r_tot;ii++)	
				 for(jj=0;jj<System.n_z_tot;jj++)
					{
					 Temp[ii][jj]=Temp_a[ii][jj];
 					 Temp_n[ii][jj]=Temp_a[ii][jj];
					}
			 System.CalculaAlfaIsoterm(alfa_2,alfa_a,Temp_a);
			 System.CalculaImp_r(alfa_2,alfa_a,Temp_2,Temp_a);
			 System.CalculaAlfaIsoterm(alfa_n,alfa_2,Temp_2);
			 System.CalculaImp_z(alfa_n,alfa_2,Temp_n,Temp_2);
			 k=0;
			 ErInc_ter=10*System.Err_Rel_Max;
			}
		 if(kk>MAXIT)
			{
			 fprintf(PF2,"\n\t\tWarning, no convergence at the time step number %d\n",j);
			 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_bmp)
			 {
			  System.RegBMP(alfa,"alpha",File_Path,1.0,0.0);
#ifdef __TEMP_BUSC_MAX__
			  Temp_Max=Temp[0][0];
			  Temp_Min=Temp[0][0];
			  for(ii=0;ii<System.n_r_sample;ii++)
				 for(jj=0;jj<System.n_z_sample;jj++)
					{
					 if(Temp[ii][jj]>Temp_Max) Temp_Max=Temp[ii][jj];
					 if(Temp[ii][jj]<Temp_Min) Temp_Min=Temp[ii][jj];
					}
			 System.RegBMP(Temp,"Temp",File_Path,Temp_Max,Temp_Min);
#else
			 System.RegBMP(Temp,"Temp",File_Path,System.T_max_adim,System.T_min_adim);
#endif
 			 System.k_front_bmp++;	
			 System.Val_Front_bmp+=System.Inc_Front_bmp;	
			}
		 if(System.Avg_alpha>System.Val_Front_prof) System.RegiFront(Temp,alfa,PF2,File_Path);
		 if(System.Avg_alpha>System.Val_N_Point_s)
			{
			 fprintf(PF3,"%lg %lg %lg %lg %lg\n",
					System.elap_time,System.T_furnace*System.f_Temp,d_Avg_alpha,System.Avg_alpha);
			 System.Val_N_Point_s+=System.Inc_N_Point_s;
			}
		}
	 else
		{
		 if(System.elap_time>System.Val_Front_bmp)
			 {
			  System.RegBMP(alfa,"alpha",File_Path,1.0,0.0);
#ifdef __TEMP_BUSC_MAX__
			  Temp_Max=Temp[0][0];
			  Temp_Min=Temp[0][0];
			  for(ii=0;ii<System.n_r_sample;ii++)
				 for(jj=0;jj<System.n_z_sample;jj++)
					{
					 if(Temp[ii][jj]>Temp_Max) Temp_Max=Temp[ii][jj];
		 			 if(Temp[ii][jj]<Temp_Min) Temp_Min=Temp[ii][jj];
					}
			 System.RegBMP(Temp,"Temp",File_Path,Temp_Max,Temp_Min);
#else
			 System.RegBMP(Temp,"Temp",File_Path,System.T_max_adim,System.T_min_adim);
#endif
			 System.k_front_bmp++;	
			 System.Val_Front_bmp+=System.Inc_Front_bmp;	
			}
		 if(System.elap_time>System.Val_Front_prof) System.RegiFront(Temp,alfa,PF2,File_Path);
		 if(System.elap_time>System.Val_N_Point_s)
			{
			 fprintf(PF3,"%lg %lg %lg %lg %lg\n",
					System.elap_time,System.T_furnace*System.f_Temp,d_Avg_alpha,System.Avg_alpha);
			 System.Val_N_Point_s+=System.Inc_N_Point_s;
			}
		}	

	 // screen information of the evolution of the calculation
	 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
	 if(j%System.N_Point_Calc==0) System.RegPosFront(PF,alfa,true,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_2*=1.1;
		 System.gamma_subs*=1.1;
 		 System.gamma_subs_2*=1.1;
		 System.Inc_t_2=System.gamma_2*System.Inc_s*System.Inc_s;
 		 System.Inc_t=2.0*System.Inc_t_2;
		}
	}

 // saves the final profile (temperature and degree of transformation as a function of the position)
 System.RegBMP(alfa,"alpha",File_Path,1.0,0.0);
#ifdef __TEMP_BUSC_MAX__
 Temp_Max=Temp[0][0];
 Temp_Min=Temp[0][0];
 for(ii=0;ii<System.n_r_sample;ii++)
		 for(jj=0;jj<System.n_z_sample;jj++)
			{
			 if(Temp[ii][jj]>Temp_Max) Temp_Max=Temp[ii][jj];
 			 if(Temp[ii][jj]<Temp_Min) Temp_Min=Temp[ii][jj];
			}
 System.RegBMP(Temp,"Temp",File_Path,Temp_Max,Temp_Min);
#else
 System.RegBMP(Temp,"Temp",File_Path,System.T_max_adim,System.T_min_adim);
#endif

 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);

 if(System.Reaction_exists) fclose(PF);
 fclose(PF2);
 for(ii=0;ii<System.n_r_sample;ii++)
    {
	 delete [] alfa_n[ii];
	 delete [] alfa_2[ii];
	 delete [] alfa_a[ii];
    }
 for(ii=0;ii<System.n_r_tot;ii++)
    {
	 delete [] alfa[ii];
	 delete [] Temp_n[ii];
	 delete [] Temp_2[ii];
	 delete [] Temp_a[ii];
	 delete [] Temp[ii];
    }
 delete [] alfa_n;
 delete [] alfa_2;
 delete [] alfa_a;
 delete [] alfa;
 delete [] Temp_n;
 delete [] Temp_2;
 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;
}
