/*
What's new:
C44 - [Martin Zoubek] first public release
*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "FileIO.h"
#include "Memory.h"
#include "Mesh.h"
#include "Solver.h"
#include "Time.h"
#include "Types.h"

int main(int argc,char **argv){
  FILE *f1,*f2,*f3;
  char *gridfile,*solfile,*sol2file,*convfile,*statfile,*integrate,*fixed,*s2,*localtimestep;
  char s3[1000];
  int i,surfnr,local,init;
  tMesh_Mesh *m;
  double rho,v1,v2,p,kappa,linfty=-1,physicaltime=-1,physicaltimestep=-1,
         limiting;
  tSolver_Boundary *b;
  void (*integrateproc)(tSolver_State *)=0;
  tSolver_State *s;
  double time0,time1, CFL, Re, Pr,Minfty,thetainfty;
  double rhoL,v1L,v2L,pL,rhoR,v1R,v2R,pR;
  BOOL again;
  #define WriteTimeStep	300
  #define Exit0(s,c){\
    printf(s);\
    return c;\
    }
  #define Exit1(s,s1,c){\
    printf(s,s1);\
    return c;\
    }
  #define eq(a,b) (fabs((a)-(b))<1e-20)


  if(argc<2)
    Exit0("Syntax: Solve inifile\n",1);

  if((f1=fopen(argv[1],"rb"))==0)
    Exit1("Error: Can't open file '%s'\n",argv[1],2);
  gridfile=FileIO_ReadString(f1);  //input .grid file
  printf("Source grid file: %s\n",gridfile);
  if((f2=fopen(gridfile,"rb"))==0)
    Exit1("Error: Can't open file '%s'\n",gridfile,2);

  kappa=FileIO_ReadFloat(f1);

  init=FileIO_ReadInteger(f1);  //initial condition (1) or .sol file (2)
  if(init==1 || init==4){  /* initial condition, init=4  Ringleb flow */
    m=Mesh_Read(f2,0);
    rho=FileIO_ReadFloat(f1);
    v1=FileIO_ReadFloat(f1);
    v2=FileIO_ReadFloat(f1);
    p=FileIO_ReadFloat(f1);
    printf("Initial condition: rho %0.10g, v %0.10g,%0.10g, p %0.10g, kappa %0.10g\n",
      rho,v1,v2,p,kappa);
    Solver_SetInitialCondition(m,rho,v1,v2,p,kappa);
    }
  else if(init==3) { /* Sod Channel */ 
    m=Mesh_Read(f2,0);
    rhoL=FileIO_ReadFloat(f1);
    v1L=FileIO_ReadFloat(f1);
    v2L=FileIO_ReadFloat(f1);
    pL=FileIO_ReadFloat(f1);
    printf("Init con: rho %0.10g, v %0.10g,%0.10g, p %0.10g, kappa %0.10g\n",
      rho,v1,v2,p,kappa);
    rhoR=FileIO_ReadFloat(f1);
    v1R=FileIO_ReadFloat(f1);
    v2R=FileIO_ReadFloat(f1);
    pR=FileIO_ReadFloat(f1);
    printf("Init con3: rho %0.10g, v %0.10g,%0.10g, p %0.10g, kappa %0.10g\n",
      rhoR,v1R,v2R,pR,kappa);
     Solver_SetInitialCondition_Sod(m,rhoL,v1L,v2L,pL,rhoR,v1R,v2R,pR,kappa);

  }
  else{  //initial solution file
    solfile=FileIO_ReadString(f1);
    printf("Initial solution file: %s\n",solfile);
    if((f3=fopen(solfile,"rb"))==0)
      Exit1("Error: Can't open file '%s'\n",solfile,2);
    m=Mesh_Read(f2,f3);  
    m->kappa=kappa;
    fclose(f3);
    }
  fclose(f2);

  sol2file=FileIO_ReadString(f1);  //output .sol file
  printf("Output solution file: %s\n",sol2file);
  convfile=FileIO_ReadString(f1);  //convergence file (- to suppress)
  if(strcmp(convfile,"-"))
    printf("Convergence file: %s\n",convfile);
  statfile=FileIO_ReadString(f1);  //statistics file (- to suppress)
  if(strcmp(statfile,"-"))
    printf("Statistics file: %s\n",statfile);
  s2=FileIO_ReadString(f1);  //stopping criterion: L^infty norm (- to suppress)
  if(strcmp(s2,"-")){
    linfty=atof(s2);
    printf("Stopping criterion: L^infty norm: %0.10g\n",linfty);
    }
  s2=FileIO_ReadString(f1);  //stopping criterion: physical time (- to suppress)
  if(strcmp(s2,"-")){
    physicaltime=atof(s2);
    printf("Stopping criterion: Physical time: %0.10g s\n",physicaltime);
    }
  if(eq(linfty,-1)&&eq(physicaltime,-1))
    Exit0("Error: No stopping criterion defined\n",3);
  s2=FileIO_ReadString(f1);  //numbered solution output physical time step (- to suppress)
  if(strcmp(s2,"-")){
    physicaltimestep=atof(s2);
    printf("Numbered solution output physical time step: %0.10g s\n",
      physicaltimestep);
    }
  integrate=FileIO_ReadString(f1);  //integrate proc (euler/ec/rk4)
  if(!strcmp(integrate,"euler")){
    integrateproc=Solver_Integrate_Euler;
    printf("Integrate using Euler formula\n");
    }
  if(!strcmp(integrate,"ec")){
    integrateproc=Solver_Integrate_EulerCauchy;
    printf("Integrate using Euler-Cauchy formula\n");
    }
  if(!strcmp(integrate,"rk4")){
    integrateproc=Solver_Integrate_RK4;
    printf("Integrate using Runge-Kutta 4 formula\n");
    }
  if(!integrateproc)
    Exit1("Error: Unknown integrate formula '%s'\nValid formulae are 'euler', 'ec' and 'rk4'\n",
      integrate,3);
  
  localtimestep=FileIO_ReadString(f1); // local time stepping
  if( !strcmp(localtimestep, "L") ||!strcmp(localtimestep, "local") ||
      !strcmp(localtimestep, "LOCAL") ){
    local = 1;
    printf("Local time stepping\n");
  }
  else{
    printf("Global time stepping\n");
    local = 0;
  }

  limiting=FileIO_ReadFloat(f1); 
  printf("Tolerance for limiting: %7.5g \n",limiting);

  Re = FileIO_ReadFloat(f1);
  Pr = FileIO_ReadFloat(f1);
  if(Re) printf("Reynolds number=  %7.5g,  Prandtl number= %7.5g\n",Re,Pr);
  else Re = 0;
	
  i=FileIO_ReadInteger(f1);  //inlet/outlet count
  printf("Inlet/outlet sides: %i\n",i);
  b=Solver_BoundaryCreate(m);
  while(i){  //for each inlet/outlet: surfnr physical/fixed rho v1 v2 v3 p
    surfnr=FileIO_ReadInteger(f1);
    fixed=FileIO_ReadString(f1);
    if(!strcmp(fixed,"physical")&&!strcmp(fixed,"fixed"))
      Exit1("Error: Unknown inlet/outlet type '%s'\nValid types are 'physical' and 'fixed'\n",
        fixed,3);
    rho=FileIO_ReadFloat(f1);
    v1=FileIO_ReadFloat(f1);
    v2=FileIO_ReadFloat(f1);
    p=FileIO_ReadFloat(f1);
    printf("Side %i: %s, rho %0.10g, v %0.10g,%0.10g, p %0.10g\n",
      surfnr,fixed,rho,v1,v2,p);
    Solver_BoundarySetInletOutlet(b,surfnr,rho,v1,v2,p,
      (BOOL)!strcmp(fixed,"fixed"),m->kappa);
    if(surfnr == 1){
      Minfty=sqrt((v1*v1+v2*v2)*rho/m->kappa/p);
      thetainfty=p/(m->kappa-1)/rho;
    }
  
    i--;
    }
  CFL = FileIO_ReadFloat(f1);
  fclose(f1);

  /*  printf("Before  Curved_Element_Detecting \n"); */
  if(!Curved_Element_Detecting(m))
         Exit0("Error: No file with curved boundary describtion\n",3);
  /*  printf("After Curved_Element_Detecting \n"); */

  s=Solver_Init(m,b,CFL,integrateproc);
  printf("Solver_Init finished\n");

  s->local = local;
  s->limiting = limiting;
  s->number_of_limiting=0;
  s->number_of_nolimiting=0;
  s->Minfty = Minfty;
  s->thetainfty=thetainfty;
  s->Re= Re;
  if(Re) s->kappaRePr = s->kappa/s->Re/Pr;

  if(!eq(physicaltime,-1))
    s->nextt=physicaltime;
  if(!eq(physicaltimestep,-1)){
    sprintf(s3,"%s%0.10g",sol2file,0.0);
    if((f2=fopen(s3,"wb"))==0)
      Exit1("Error: Can't create file '%s'\n",s3,2);
    Mesh_Write(m,0,f2);
    fclose(f2);
    if(physicaltimestep<s->nextt)
      s->nextt=physicaltimestep;
    }
  if(strcmp(convfile,"-")){
    if((f1=fopen(convfile,"wb"))==0)
      Exit1("Error: Can't create file '%s'\n",convfile,2);
    }
  else
    f1=0;

  printf("\n");
  time0=time1=Time_CPUGet();
  i=0;
  do{
    Solver_Solve(s);
    if(i%500 == 0){
      printf("\n");
      printf("    Iter      Time       Reziduum   Mach: minimal maximal    limit \n");
      printf("--------------------------------------------------------------------\n");
    }
    if(i%25 == 0)
      printf(" %6d %12.4f %14.4e  %10.4f %10.4f %5d(%d)\n",
        i,s->t,s->linfty,s->Machmin,s->Machmax,s->number_of_limiting,
	     s->number_of_nolimiting);

    if(f1 && i%10 == 0){
      //fprintf(f1,"%g %0.10g\n",s->t,s->linfty);

      if(s->linfty>0) {
	fprintf(f1,"%i %10.6g %10.6g %10.6g %10.6g %10.6g\n",
		i,log(s->linfty)/log(10),
		s->cL, s->cD, s->cDp,s->cD-s->cDp );
      }
      else
	fprintf(f1,"%i %10.6g\n",i,-30);
    }
    i++;

    if(time1+WriteTimeStep<Time_CPUGet()||s->linfty<=linfty||
      eq(s->t,s->nextt)){
      printf(" %6d %12.4f %14.4e  %10.4f %10.4f \n",
        i,s->t,s->linfty,s->Machmin,s->Machmax);
      if(strcmp(statfile,"-")){
        if((f2=fopen(statfile,"wb"))==0)
          Exit1("Error: Can't create file '%s'\n",statfile,2);
        fprintf(f2,"Iterations: %i\n",i);
        fprintf(f2,"Time: %0.10g s\n",Time_CPUGet()-time0);
        fprintf(f2,"Memory used: %i KB\n",vMemory_MaxAlloc>>10);
        fprintf(f2,"Physical time: %0.10g s\n",s->t);
        fprintf(f2,"L^infty: %0.10g\n",s->linfty);
        fclose(f2);
        }
      Solver_MeshUpdate(s);
      Plot_Walls(s);
      if((f2=fopen(sol2file,"wb"))==0)
        Exit1("Error: Can't create file '%s'\n",sol2file,2);
      Mesh_Write(m,0,f2);
      fclose(f2);
      if(eq(s->t,s->nextt)&&!eq(physicaltimestep,-1)){
        sprintf(s3,"%s%0.10g",sol2file,s->t);
        if((f2=fopen(s3,"wb"))==0)
          Exit1("Error: Can't create file '%s'\n",s3,2);
        Mesh_Write(m,0,f2);
        fclose(f2);
        s->nextt+=physicaltimestep;
        if(!eq(physicaltime,-1)&&s->nextt>physicaltime)
          s->nextt=physicaltime;
        }
      time1=Time_CPUGet();
      }
    again=TRUE;
    

    //printf("Phys time %7.5g %7.5g %i \n",s->t,physicaltime,!eq(linfty,-1));

    if(!eq(linfty,-1))
      if(s->linfty<=linfty)
        again=FALSE;
      else{
        if(!eq(physicaltime,-1))
          if(s->t>=physicaltime)
            again=FALSE;
	  
	
      }
    else
      if(eq(s->t,physicaltime))
        again=FALSE;

    }while(again);
  if(init==3){ /* Sod Channel */
    printf("Computational error\n");
    Sod_Error(s,rhoL,v1L,v2L,pL,rhoR,v1R,v2R,pR);
  }
  if(init==4 || init==5){ /* Ringleb flow */
    printf("Computational error\n");
    Ringleb_Error(s);
  }
  if(f1)
    fclose(f1);
  return 0;
  }

