//This is a hybrid lattice-Boltzmann code for nematic hydrodynamics

//Create folder data!!
//Compile with "gcc -O3 -o prog an_01.c -lm -lpthread -march=native -std=gnu99 -fopenmp -Wall -Wno-unused-variable -Wno-unused-but-set-variable"
//Compile with "gcc -O3 -o prog an_01.c -lm -lpthread -march=native -std=gnu99 -fopenmp -Wall -Wno-unused-variable -Wno-unused-but-set-variable -Wno-format-overflow"

//ActiveDroplet 
//Compile with "gcc -O3 -o prog2 an_01.c diff_eq_system.c -lm -lpthread -march=native -std=gnu99 -fopenmp -Wall -Wno-unused-variable -Wno-unused-but-set-variable -Wno-format-overflow"
//gcc -O3 -o prog2 -g an_01.c diff_eq_system.c -lm -lpthread -march=native -std=gnu99 -fopenmp -Wall -Wno-unused-variable -Wno-unused-but-set-variable -Wno-format-overflow

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <omp.h>
#include <string.h>


int I;	//Number of points in each direction
int J;
int K;
#define NMAX			(I*J*K)		//NMAX is no longer defined as (I+1)*(J+1)*(K+1)

#define ST_HITROSTI		19		//Number of velocity vectors within the model -- in this case 3DQ19 -- 19 velocities		To change the model, one needs to set new values of ST_HITROSTI, new vectors E and also set FEQ, P, and bounce-back

#define STCASKOR		5000000
#define STCASKORINIT		0

#define STCASKORQ		1		//How many steps in Q are performed for each step in LB

#define STKONV			100		//On how many steps the fields are written to a file
#define STKONVV			10000		//On how many steps the fields are written to a file
#define STKONG          100     //On how many steps the defect grad and config are written to a file


#define TWRITEU			1		//Writes velocity
#define TWRITEDEN		1		//Writes density
#define TWRITEDIR		1		//Writes director field
#define TWRITEOP		1       //Writes order parameter field
#define TWRITEQ			1       //Writes Q-tensor field

#define STDVF			100		//On how many steps the defect volume fraction is written to a file
#define TWRITEDVF		1       //Writes defect volume fraction
#define STFE			500		//On how many steps the free energy is written to a file
#define TWRITEFE		1       //Writes free energy
#define TWRITELOC       1       //Writes location of defects
#define TWRITEgrad      1       //Writes defet gradients

// #define STU             300     //Change the initial condition
// #define initialU        1

#define STCOPY          10000      //On how many steps the copy of DATA.bin is written to a file 
#define TWRITECOPY      1      //Writes copy of DATA.bin



#define ZAPISF			1		//Write all necessary functions to a file at the end of the computation
#define BERIF			0		//Read initialisation from the file

#define READINITIAL     0




// #define ZETAchange		70	    //change zeta value
// 
// #define ZETAstart		50	//change zeta value (start)
// #define ZETAend		    80	//change zeta value (end)
// 

int STPROC=64;					//Number of CPU cores for the paralelisation

float **U, **F, **FEQ, **FNEW, **P;				//Allocation of global matrices
float **Q, **QNIC, **QNEW;
float *WKONST;
int **E;
char *LMARK;
float **UNIC;
float **H;
float **FORCE;
float **TAU;

char* QMARK;

// number of surface data points
int Qmax;
float* QNC;

float **gradU;
float **gradU_defect;

char Izpis[] = "droplet";
char imeIzpis[100];


#include "acn_funk_head.h"		//function headers
#include "diff_eq_system.h"
#include "konstante.c"
#include "num_rec.c"			//All needed function from numerical recipes in a single file
#include "acn_utils.c"			//Useful functions -- utility
#include "acn_active_droplet.c"    //defects on shell
#include "acn_funk_lb.c"		//Lattice-Boltzmann part of the code
#include"acn_funk_op.c"
#include "acn_funk_zapis.c"




int main(int argc, char** args){

  ZETA=(float) atof(args[1]);
  I=(int) atoi(args[2]);
  J=(int) atoi(args[3]);
  K=(int) atoi(args[4]);
  
  sprintf(imeIzpis, "%s-%dx%dx%d_%f", Izpis, I, J, K, ZETA);
  
  long Seed=12345; 					//Seed for random initial conditions
  char buf[256];
  int l,m;
  
  U=matrix(1,4,0,NMAX-1);				// velocity field; 4th component is density
  F=matrix(0,ST_HITROSTI-1,0,NMAX-1);		// Lattice Boltzmann scalar distribution functions to calculate velocity U
  FNEW=matrix(0,ST_HITROSTI-1,0,NMAX-1);		// matrix field to temporarily store currently calculated F
  FEQ=matrix(0,ST_HITROSTI-1,0,NMAX-1);		// Lattice Boltzmann scalar distribution functions to calculate velocity U - equilibrium; nedded in colloision operators
  P=matrix(0,ST_HITROSTI-1,0,NMAX-1);		// Lattice Boltzmann scalar distribution functions to calculate velocity U - forcing term; needed in colloision operators
  TAU=matrix(1,9,0,NMAX-1);
  
  Q=matrix(1,6,0,NMAX-1);				//Q-tensor
  QNEW=matrix(1,6,0,NMAX-1);
  QNIC=matrix(1,6,0,NMAX-1);			//Qnic include info about the normal, prefered surface orientation and possibly also surface velocity
  
  WKONST=vector(0,18);
  LMARK=malloc(NMAX*sizeof(char));
  UNIC=matrix(1,3,0,NMAX-1);			//Surface velocity (usually commented out to save time and RAM)
  gradU=matrix(1,3,0,NMAX-1);
  gradU_defect=matrix(1,3,0,NMAX-1);
  
  H=matrix(1,6,0,NMAX-1);				//H tensor -- molecular field
  FORCE=matrix(1,3,0,NMAX-1);				//FORCE = divergence of a total stress tensor -- needed to correct the velocity due to discrete effects
  
  E=imatrix(0,ST_HITROSTI-1,1,3);			// matrix of all LB characteristic velocity lattice vectors
  
  initialiseE();				//Characteristic lattice vectors for the chosen LB model

// if(BERIF){ sprintf(buf, "droplet-200x200x200_0.010000-DATA.bin", imeIzpis); beriDATA(buf);		//Reads F
  
  if(BERIF){ sprintf(buf, "%s-DATA.bin", imeIzpis); beriDATA(buf);		//Reads F
    for(l=0;l<NMAX;l++){								
      for(m=0;m<ST_HITROSTI;m++) FNEW[m][l]=F[m][l];				//Transcribe F to FNEW
      FORCE[1][l]=0; FORCE[2][l]=0; FORCE[3][l]=0; calcF2Ul(l,F);		//Compute U
    }
  }
  else{
    zacPriblQ(Seed);
    zacPriblUF(Seed);				// initial condition for U and F
   
  }
  
  //read initial conditions (janus core)
//   if(READINITIAL) {sprintf(buf,"dim_par_eightL_relax_E01.raw"); beriDROPLETDATA(buf);} 
//   if(READINITIAL) {sprintf(buf,"dim_par_satrn_relax_E01.raw"); beriDROPLETDATA(buf);}


//     ActiveDroplet sphere1;
//     ActiveDroplet* sphere1_ptr = &sphere1;
//     init_droplet(sphere1_ptr);
//     printf("Len after init: %d\n", sphere1.num_points);
//printf("%d\n", sphere1.num_points);
      
  izrQU();					//Main computational loop

  
  
  //Calculates the upper estimates for the Reynolds and Ericksen number
  float umax=0;
  for(l=0;l<NMAX;l++) if(LMARK[l]&LMARK1){
   float delovni=fsqr(U[1][l]) + fsqr(U[2][l]) + fsqr(U[3][l]);
   if(delovni>umax) umax=delovni;
  }
  printf("reynolds: %e\nericksen: %e\n gostota: %f\n", 3*sqrt(umax)*I/(TAUF/DT-0.5f), 1.f/3.f*(TAUF/DT-0.5f)*DENSITYINIT*sqrt(umax)*I, DENSITYINIT);
  

  
  
  //Write the desired variables in .raw files
  //Write the desired variables in .vtk files
  sprintf(buf,"./data/%s_u.vtk",imeIzpis);
  zapisUraw(buf);					// write velocity profile

  sprintf(buf,"./data/%s_den.vtk",imeIzpis);
  zapisDENraw(buf);					// write density profile
  
  sprintf(buf,"./data/%s_dir.vtk",imeIzpis);
  zapisQ2DIRraw(buf);
  
  sprintf(buf,"./data/%s_op.vtk",imeIzpis);
  zapisQ2OPraw(buf);
  
  //sprintf(buf,"./data/%s_u.dat",imeIzpis);		//velocity profile for gnuplot
  //za_profilU(buf,K/2,I/2);
  
//   zapisDVF(buf);
  


  //Write all the variables in a .bin file
  if(ZAPISF){ sprintf(buf, "%s-DATA.bin", imeIzpis); zapisDATA(buf);}
  
  
  //free all matrices	
  free_matrix(U,1,4,0,NMAX-1);			free_matrix(F,0,ST_HITROSTI-1,0,NMAX-1);		free_matrix(FNEW,0,ST_HITROSTI-1,0,NMAX-1);	free_matrix(FEQ,0,ST_HITROSTI-1,0,NMAX-1);			free_matrix(P,0,ST_HITROSTI-1,0,NMAX-1);	free_imatrix(E,0,ST_HITROSTI-1,1,3);			free_matrix(Q,1,6,0,NMAX-1);	free_matrix(QNEW,1,6,0,NMAX-1);
  free_matrix(QNIC,1,6,0,NMAX-1);		free_vector(WKONST,0,18);				free(LMARK);					free_matrix(UNIC,1,3,0,NMAX-1);
  free_matrix(TAU,1,9,0,NMAX-1);
  free_matrix(H,1,6,0,NMAX-1);			free_matrix(FORCE,1,3,0,NMAX-1); free_matrix(gradU,1,3,0,NMAX-1); free_matrix(gradU_defect,1,3,0,NMAX-1);

  return 0;
}