/** Page containing all diagnostics functions and events of the Idealized case (diagnostics are identical to the Krabbendijke case). Note that the imported output_silces.h can be found [here](http://basilisk.fr/sandbox/vheusinkveld/myfunctions/output/sliceOf3D.h) */

#include "utils.h"
#if dimension == 3
	#include "lambda2.h"
#endif
#include "output_slices.h"

/** Define variables and structures to do: diagnostics, ouput data, output movies. */

struct sDiag dia; 			// Diagnostics

struct sDiag {
	double Ekin;			// Total kinetic energy
	double EkinOld;			// Track changes in kin energy 
	double WdoneOld;		// Track changes in work done 
	double rotVol;			// Diagnosed rotor volum
	double bE;			// Buoyancy energy
	double bEold;			// Track changes in buoyancy energy
	double diss;			// Diagnose dissipation
};

struct sEquiDiag {
    int level;			// Level for which diagnostics should be done
    int ii; 			// Keep track of how many additions are done
    double dtDiag;
    double startDiag;
    double endDiag;
    double dtOutput;
};

struct sOutput {
    double dtDiag;
    double dtVisual;
    double dtSlices;
    double dtProfile;
    double startAve;
    double dtAve;
    char main_dir[12];
    char dir[30];
    char dir_profiles[60];
    char dir_slices[60];
    char dir_equifields[60];
    char dir_strvel[60];
    char dir_refvel[60];
    char dir_diffbins[60];
    char dir_dts[60];
    int sim_i;
};
	
struct sbViewSettings {
    double phi; 			// Phi for 3D bview movie
    double theta;			// Theta for 3D bview movie
    double sphi; 			// Polar angle for sliced image RED
    double stheta;			// Azimuthal angle for sliced image RED
};

/** Initialize structures */
struct sOutput out = {.dtDiag = 1., .dtVisual=1., .dtSlices=5., .dtProfile=60., .main_dir="results", .sim_i=0};

struct sEquiDiag ediag = {.level = 5, .ii = 0, .startDiag = 0., .dtDiag = 0., .dtOutput = 0.};

struct sbViewSettings bvsets = {.phi=0., .theta=0., .sphi=0., .stheta=0.};

event init(i = 0){
	bvsets.phi = 0.;
	bvsets.theta = -M_PI/6.;
	bvsets.sphi = 0.;
	bvsets.stheta = 0.;
}

/** Diagnosing: kinetic energy, diagnosed rotor volume, buoyancy energy, ammount of cells used.*/
event diagnostics (t+=out.dtDiag){
	int n = 0.;
	scalar ekin[]; 		// Kinetic energy field 
	double tempVol = 0.;    // Temp volume 
	double tempEkin = 0.;   // Temp kinetic energy
	double tempDiss = 0.;   // Temp dissipation 
	double maxVel = 0.;     // Maximum velocity in fan
	double bEnergy = 0.;    // Buoyant energy
		
	/** Loop over cells to get diagnostics */ 
	foreach(reduction(+:n) reduction(+:tempVol) reduction(+:tempEkin) 
		reduction(max:maxVel) reduction(+:bEnergy) reduction(+:tempDiss)) {

		tempVol += dv()*fan[];
		if(y + Delta/2. <= rot.y0){
		     bEnergy += dv()*y*(b[] - STRAT(y));
		}
		foreach_dimension() {
			ekin[] += sq(u.x[]);
		}
		maxVel = max(maxVel, sq(ekin[]));
		ekin[] *= 0.5*rho[]*dv();	
		tempEkin += ekin[];
		n++;
	}
	/** Assign values to respective global sturcture vars */ 
	dia.diss = 1.*tempDiss;
	dia.bE = 1.*bEnergy;
	dia.rotVol = 1.*tempVol;
	dia.Ekin = 1.*tempEkin;
	
	if (pid() == 0){
	/** Write away simulation data and case setup for main thread */ 
	char nameOut[90];
	char nameCase[90];
     	snprintf(nameOut, 90, "./%s/output", out.dir);
     	snprintf(nameCase, 90, "./%s/case", out.dir);
	static FILE * fpout = fopen(nameOut, "w");
	static FILE * fpca = fopen(nameCase, "w");

	if(t==0.){
		fprintf(fpca,"L0\tinversion\thubU\tTref\tLambda\txr\tyr\tzr\ttheta\tphit\tr\tW\tP\tcu\trampT\tmaxlvl\tminlvl\teps\n");
		fprintf(fpca, "%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%d\t%d\t%g\n", 
				L0,TREF/gCONST*STRAT(rot.y0), WIND(rot.y0),TREF, Lambda, rot.x0, rot.y0, rot.z0, rot.theta, rot.phit, rot.R, rot.W, rot.P, rot.cu, rot.rampT, maxlevel, minlevel, eps);
		
	        fprintf(stderr,"n\tred\tEkin\tWork\tbE\n");
		fprintf(fpout,"i\tt\tn\tred\tEkin\tWork\tbE\n");
	}
	fprintf(fpout, "%d\t%g\t%d\t%g\t%g\t%g\t%g\n",
		i,t,n,(double)((1<<(maxlevel*3))/n),dia.Ekin,rot.Work, dia.bE);
	
	fprintf(stderr, "n=%d\t%g\t%g\t%g\t%g\n",n,(double)((1<<(maxlevel*dimension))/n),dia.Ekin,rot.Work,dia.bE);
	
	fflush(fpout);
	fflush(fpca);	

	}

	dia.EkinOld = 1.*dia.Ekin;
	dia.WdoneOld = 1.*rot.Work;
	dia.bEold = 1.*dia.bE;
}

#if dimension == 3
/** diagnose velocity of the passing wind machine*/
event refvelocity(t+=1) {
    if(pid() == 0) {
    char nameRefvel[90];
    snprintf(nameRefvel, 90, "%st=%05g", out.dir_refvel, t);
    FILE * fpstr = fopen(nameRefvel, "w");
    fprintf(fpstr, "x,v,vx,vy,vz\n");
        
    double length = 300.;
    int ntot = 300.;
    double vels1[ntot];	

    double xf0 = rot.x0 - length/2*cos(M_PI/6);
    double yf0 = 3;
    double zf0 = rot.z0 + length/2*sin(M_PI/6);

    for(int n = 0; n < ntot; n++) {
	double dist = length*n/ntot;
	double xx = xf0 + dist*cos(-M_PI/6); 
	double yy = yf0; 
	double zz = zf0 + dist*sin(-M_PI/6); 

	double valx1 = interpolate(u.x, xx, yy, zz);
	double valy1 = interpolate(u.y, xx, yy, zz);
	double valz1 = interpolate(u.z, xx, yy, zz);
         
        vels1[n] = sqrt(sq(valx1) + sq(valy1) + sq(valz1));
	fprintf(fpstr, "%g,%g,%g,%g,%g,%g\n", xx, zz, vels1[n], valx1, valy1, valz1);     
    }  
    fclose(fpstr); 
    }
}
/** 'DTS' measurments*/
event dts_meas(t += 1) {
    if(pid() == 0) {
	char nameDtshorS[90];
    	snprintf(nameDtshorS, 90, "%shorS_t=%05g", out.dir_dts, t);
        FILE * fpstrhorS = fopen(nameDtshorS, "w");
        fprintf(fpstrhorS, "x,y,z,b\n");
        
	double lengthhorS = 600.;
        int ntothorS = 600;
	
	double xf0S = rot.x0 - lengthhorS/2*cos(M_PI/6);
	double yf0S = 1;
	double zf0S = rot.z0 + lengthhorS/2*sin(M_PI/6);

	for(int n = 0; n <= ntothorS; n++) {
	    double dist = lengthhorS*n/ntothorS;
	    double xx = xf0S + dist*cos(-M_PI/6); 
	    double yy = yf0S; 
	    double zz = zf0S + dist*sin(-M_PI/6); 
       
	    double valb = interpolate(b, xx, yy, zz);
          
	    fprintf(fpstrhorS, "%g,%g,%g,%g\n", xx, yy, zz, valb);     
	}  
 	fclose(fpstrhorS); 

	char nameDtshorL[90];
    	snprintf(nameDtshorL, 90, "%shorL_t=%05g", out.dir_dts, t);
        FILE * fpstrhorL = fopen(nameDtshorL, "w");
        fprintf(fpstrhorL, "x,y,z,b\n");
        
	double lengthhorL = 600.;
        int ntothorL = 600;
	
	double xf0L = rot.x0 - lengthhorL/2*cos(120*M_PI/180);
	double yf0L = 1;
	double zf0L = rot.z0 + lengthhorL/2*sin(120*M_PI/180);

	for(int n = 0; n <= ntothorL; n++) {
	    double dist = lengthhorL*n/ntothorL;

	    double xx = xf0L + dist*cos(-120*M_PI/180); 
	    double yy = yf0L; 
	    double zz = zf0L + dist*sin(-120*M_PI/180); 
       
	    double valb = interpolate(b, xx, yy, zz);
          
	    fprintf(fpstrhorL, "%g,%g,%g,%g\n", xx, yy, zz, valb);     
	}  
 	fclose(fpstrhorL); 

	char nameDtsver[90];
    	snprintf(nameDtsver, 90, "%sverS_t=%05g", out.dir_dts, t);
        FILE * fpstrver = fopen(nameDtsver, "w");
        fprintf(fpstrver, "x,y,z,b\n");
        
	double lengthver = 20.;
        int ntotver = 160;

	for(int n = 0; n <= ntotver; n++) {
	    double dist = lengthver*n/ntotver;

	    double xx = rot.x0 + 30*cos(-M_PI/6); 
	    double yy = dist; 
	    double zz = rot.z0 + 30*sin(-M_PI/6); 
       
	    double valb = interpolate(b, xx, yy, zz);
          
	    fprintf(fpstrver, "%g,%g,%g,%g\n", xx, yy, zz, valb);     
	}  
 	fclose(fpstrver); 
    }
    
}


#endif

/** Ouputting movies in 2- or 3-D*/

#if dimension == 2
event movies(t += 0.5) {
    vertex scalar omega[]; 	// Vorticity
    scalar lev[];	 	// Grid depth
    scalar ekinRho[]; 		// Kinetic energy

    foreach() {
        omega[] = ((u.y[1,0] - u.y[-1,0]) - (u.x[0,1] - u.x[0,-1]))/(2*Delta); // Curl(u) 
        ekinRho[] = 0.5*rho[]*(sq(u.x[]) + sq(u.y[]));
        lev[] = level;
    }

    boundary ({b, lev, omega, ekinRho});
    output_ppm (b, file = "ppm2mp4 ./results/buoyancy.mp4", n = 1<<maxlevel, linear = true, max=STRAT(L0), min=STRAT(0.));
    output_ppm (ekinRho, file = "ppm2mp4 ./results/ekin.mp4", n = 1<<maxlevel, min = 0, max = 1.*sq(rot.cu));
    output_ppm (omega, file = "ppm2mp4 ./results/vort.mp4", n = 1<<maxlevel, linear = true); 
    output_ppm (lev, file = "ppm2mp4 ./results/grid_depth.mp4", n = 1<<maxlevel, min = minlevel, max = maxlevel);
}
#elif dimension == 3
event movies(t += out.dtVisual) {
    scalar l2[];

    lambda2(u,l2);
    boundary({l2});

    view(fov=25, tx = 0., ty = 0., phi=M_PI/2., theta=M_PI/2., width = 800, height = 800);
    
    translate(-rot.x0,-rot.y0,-rot.z0) {
        box(notics=false);
        isosurface("l2", v=-0.02, color="b", min=STRAT(0.9*roughY0h), max=STRAT(1.1*rot.y0));
	draw_vof("fan", fc = {1,0,0});
    }
   
    translate(-rot.x0,0,-rot.z0){
        squares("u.y", n = {0,1,0}, alpha=2, min=-1, max=1);
    }

    /** Save file with certain fps*/
    char nameVid4[90];
    snprintf(nameVid4, 90, "ppm2mp4 -r %g ./%s/visual_3d_ux.mp4", 10., out.dir);
    save(nameVid4);
    clear();

    view(fov=25, tx = 0., ty = 0., phi=M_PI/2., theta=M_PI/2., width = 800, height = 800);
    
    translate(-rot.x0,-rot.y0,-rot.z0) {
        box(notics=false);
        isosurface("l2", v=-0.02, color="b", min=STRAT(0.9*roughY0h), max=STRAT(1.1*rot.y0));
	draw_vof("fan", fc = {1,0,0});
    }
   
    translate(-rot.x0,0,-rot.z0){
        squares("u.y", n = {0,1,0}, alpha=2, min=-.05, max=.05);
    }

    /** Save file with certain fps*/
    char nameVid3[90];
    snprintf(nameVid3, 90, "ppm2mp4 -r %g ./%s/visual_3d_uy.mp4", 10., out.dir);
    save(nameVid3);
    clear();

    view(fov=25, tx = 0., ty = 0., phi=M_PI/2., theta=M_PI/2., width = 800, height = 800);
    
    translate(-rot.x0,-rot.y0,-rot.z0) {
        box(notics=false);
        isosurface("l2", v=-0.02, color="b", min=STRAT(roughY0h), max=STRAT(2.*rot.y0));
	draw_vof("fan", fc = {1,0,0});
    }
   
    translate(-rot.x0,0,-rot.z0){
        squares("b", n = {0,1,0}, alpha=2, min=0.5*STRAT(2), max=1.5*STRAT(2));
    }

    /** Save file with certain fps*/
    char nameVid1[90];
    snprintf(nameVid1, 90, "ppm2mp4 -r %g ./%s/visual_3d_b2.mp4", 10., out.dir);
    save(nameVid1);
    clear();

    view(fov=25, tx = 0., ty = 0., phi=bvsets.phi, theta=bvsets.theta, width = 1200, height = 1200);
    
    translate(-rot.x0,-rot.y0,-rot.z0) {
        box(notics=false);
        isosurface("l2", v=-0.02, color="b", min=STRAT(0.), max=STRAT(1.5*rot.y0));
	draw_vof("fan", fc = {1,0,0});
    }
    translate(-rot.z0,-rot.y0, -L0){
      	squares("u.x", n = {0,0,1}, alpha=rot.z0, min=-1, max=1);
        cells(n = {0,0,1}, alpha = rot.z0);
    }
    
    translate(0.,-rot.y0,-rot.z0){
        squares("b", n = {1,0,0}, alpha=rot.x0, min=STRAT(0.), max=STRAT(1.5*rot.y0));
    }

    /** Save file with certain fps*/
    char nameVid2[90];
    snprintf(nameVid2, 90, "ppm2mp4 -r %g ./%s/visual_3d_b.mp4", 10., out.dir);
    save(nameVid2);
    clear();
}

/** Take relevant field slices and write away */
event slices(t=out.dtSlices; t+=out.dtSlices) {
    char nameSlice[90];
    coord slice = {1., 0., 1.};
    int res = L0/2;

    for(double yTemp = 0.5; yTemp<=1; yTemp+=0.5) {
	slice.y = yTemp/L0;
   
    	snprintf(nameSlice, 90, "%st=%05gy=%03g", out.dir_slices, t, yTemp);
    	FILE * fpsli = fopen(nameSlice, "w");
    	output_slice(list = (scalar *){b}, fp = fpsli, n = res, linear = true, plane=slice);
    	fclose(fpsli);
    }

    for(double yTemp = 2; yTemp<=4.; yTemp+=2.) {
	slice.y = yTemp/L0;
   
    	snprintf(nameSlice, 90, "%st=%05gy=%03g", out.dir_slices, t, yTemp);
    	FILE * fpsli = fopen(nameSlice, "w");
    	output_slice(list = (scalar *){b}, fp = fpsli, n = res, linear = true, plane=slice);
    	fclose(fpsli);
    }


}
#endif

/** Usefull functions */ 

/** Checks if required folders exists, if not they get created. */
void sim_dir_create(){
   
    sprintf(out.dir, "./%s/%s%02d", out.main_dir, sim_ID, out.sim_i);
    sprintf(out.dir_profiles, "%s/profiles/", out.dir);
    sprintf(out.dir_slices, "%s/slices/", out.dir);
    sprintf(out.dir_equifields, "%s/equifields/", out.dir);
    sprintf(out.dir_strvel, "%s/strvel/", out.dir);
    sprintf(out.dir_refvel, "%s/refvel/", out.dir);
    sprintf(out.dir_diffbins, "%s/diffbins/", out.dir);
    sprintf(out.dir_dts, "%s/dts/", out.dir);

    
    if (pid() == 0){
    struct stat st = {0};
    if (stat(out.main_dir, &st) == -1) {
        mkdir(out.main_dir, 0777);
    }
   
    if (stat(out.dir, &st) == -1) {
        mkdir(out.dir, 0777);
    }
    if (stat(out.dir_slices, &st) == -1) {
        mkdir(out.dir_slices, 0777);
    }
    if (stat(out.dir_profiles, &st) == -1) {
        mkdir(out.dir_profiles, 0777);
    }  
    if (stat(out.dir_equifields, &st) == -1) {
        mkdir(out.dir_equifields, 0777);
    }  
    if (stat(out.dir_strvel, &st) == -1) {
	mkdir(out.dir_strvel, 0777);
    }
    if (stat(out.dir_refvel, &st) == -1) {
	mkdir(out.dir_refvel, 0777);
    }
    if (stat(out.dir_diffbins, &st) == -1) {
	mkdir(out.dir_diffbins, 0777);
    }
    if (stat(out.dir_dts, &st) == -1) {
	mkdir(out.dir_dts, 0777);
    }

    }
}