/*============================================================================
 * This function is called at the end of each time step, and has a very
 *  general purpose
 *  (i.e. anything that does not have another dedicated user function)
 *============================================================================*/

/* Code_Saturne version 5.0.5 */

/*
  This file is part of Code_Saturne, a general-purpose CFD tool.

  Copyright (C) 1998-2017 EDF S.A.

  This program is free software; you can redistribute it and/or modify it under
  the terms of the GNU General Public License as published by the Free Software
  Foundation; either version 2 of the License, or (at your option) any later
  version.

  This program is distributed in the hope that it will be useful, but WITHOUT
  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
  details.

  You should have received a copy of the GNU General Public License along with
  this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
  Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/

/*----------------------------------------------------------------------------*/

#include "cs_defs.h"

/*----------------------------------------------------------------------------
 * Standard C library headers
 *----------------------------------------------------------------------------*/

#include <assert.h>
#include <math.h>
#include <stdio.h>

#if defined(HAVE_MPI)
#include <mpi.h>
#endif

/*----------------------------------------------------------------------------
 * PLE library headers
 *----------------------------------------------------------------------------*/

#include <ple_coupling.h>

/*----------------------------------------------------------------------------
 *  Local headers
 *----------------------------------------------------------------------------*/

#include "bft_mem.h"
#include "bft_error.h"
#include "bft_printf.h"

#include "cs_base.h"
#include "cs_field.h"
#include "cs_field_pointer.h"
#include "cs_field_operator.h"
#include "cs_math.h"
#include "cs_mesh.h"
#include "cs_mesh_quantities.h"
#include "cs_halo.h"
#include "cs_halo_perio.h"
#include "cs_log.h"
#include "cs_parall.h"
#include "cs_parameters.h"
#include "cs_physical_constants.h"
#include "cs_post_util.h"
#include "cs_rotation.h"
#include "cs_time_moment.h"
#include "cs_time_step.h"
#include "cs_turbomachinery.h"
#include "cs_selector.h"

#include "cs_post.h"

/*----------------------------------------------------------------------------
 *  Header for the current file
 *----------------------------------------------------------------------------*/

#include "cs_prototypes.h"
#include "cs_base.h"
#include "cs_mesh.h"
#include "cs_mesh_quantities.h"
#include "cs_mesh_bad_cells.h"
#include "cs_probe.h"
#include "cs_volume_zone.h"
#include "cs_domain.h"
/*----------------------------------------------------------------------------*/

BEGIN_C_DECLS


/*============================================================================
 * User function definitions
 *============================================================================*/

static void
_findpt_r(const cs_real_3_t     coords,
          cs_lnum_t            *node,
          cs_lnum_t            *rank_1,
		  cs_real_t 			dis2mn)
{
  cs_real_t d[3];
  cs_lnum_t n_cells = cs_glob_mesh-> n_cells;
  cs_real_3_t *cell_cen = (cs_real_3_t *)cs_glob_mesh_quantities->cell_cen;
  *node = (int)(n_cells/2+1);

  for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
      for (int i = 0; i < 3; i++){
        d[i] = coords[i] - cell_cen[cell_id][i];
	  }
	  cs_real_t dis2 = cs_math_3_square_norm(d);
	  if (dis2 < dis2mn) {
		*node = cell_id;
		dis2mn = dis2;
		}
	}
	cs_parall_min_id_rank_r(node, rank_1, dis2mn);
}

void
cs_user_extra_operations(void)
{
  int const s = 5150;
  static double mean_phi_u_arr_t[5150];
  static double mean_phi_v_arr_t[5150];
  static double mean_phi_w_arr_t[5150];
  static double mean_phi_t_arr_t[5150];
  static int ranks[4200000];
  static int cell_ids[4200000];
  
  /* Mesh-related variables */

  const cs_lnum_t n_b_faces = cs_glob_mesh->n_b_faces;
  const cs_real_3_t  *restrict cell_cen
    = (const cs_real_3_t *restrict)cs_glob_mesh_quantities->cell_cen;

  /* 0. Initialization
     ================= */

  /* Local constants defintion */
  int rank;
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  const cs_real_t _PI = acos(-1.);
  int nt_maxx = (int)cs_glob_time_step->nt_max;
  int nt_curr = (int)cs_glob_time_step->nt_cur;
  /* Example 2: extraction of a velocity profile in cylindrical corordinates
     ======================================================================= */

  if (fmod(nt_curr, 10)==0 && nt_curr!=0){
	cs_lnum_t n_cells = cs_glob_mesh-> n_cells;
    cs_real_3_t *vel = (cs_real_3_t *)CS_F_(u)->val;
	cs_real_t *t = (cs_real_t *)CS_F_(t)->val;
	double u, v, w, t1;
	double sum_phi_u, sum_phi_v, sum_phi_w, sum_phi_t;
	double mean_phi_u, mean_phi_v, mean_phi_w, mean_phi_t;
	double u_fluct_rms, v_fluct_rms, w_fluct_rms;
	double u_fluct_var, v_fluct_var, w_fluct_var;
	double u_fluct_rms_phi, v_fluct_rms_phi, w_fluct_rms_phi;
	double u_fluct_var_phi, v_fluct_var_phi, w_fluct_var_phi;
	double t_fluct_rms, t_fluct_var, t_fluct_rms_phi, t_fluct_var_phi;
	static double u_fluct_rms_arr[5150], v_fluct_rms_arr[5150], w_fluct_rms_arr[5150];
	static double u_fluct_var_arr[5150], v_fluct_var_arr[5150], w_fluct_var_arr[5150];
	static double t_fluct_rms_arr[5150], t_fluct_var_arr[5150];
	static double razd_arr[5150], k_arr[5150];
	FILE *f5 = NULL;
	char filename5[100];
	sprintf(filename5, "coords.txt");
	int j=0;
	int k = 0;
	int a1 = (int)rank;
    double radius = 0.008001;
    double axicoo = 0.03;
	double razd;
	double ux, vy, theta;
	double radius1[51] = {0.006004364905, 0.006013531205, 0.00602361415, 0.0060347054, 0.00604690575, 0.00606032599, 0.00607508839, 0.0060913272, 0.0061091912, 0.00612884, 0.0061504525, 0.0061742275, 0.00620038, 0.0062291475, 0.006260792, 0.006295601, 0.0063389, 0.00637601, 0.0064223405, 0.0064733045, 0.00655, 0.00665, 0.00675, 0.00685, 0.00695, 0.00705, 0.00715, 0.00725, 0.00735, 0.00745, 0.007526695, 0.00757766, 0.00762399, 0.007668885, 0.0077044, 0.007739205, 0.00777085, 0.00779962, 0.007825775, 0.007825775, 0.00784955, 0.00787116, 0.007891075, 0.00790894, 0.00792491, 0.00793967, 0.00795309, 0.00796529, 0.007976385, 0.00798647, 0.007995635};
    cs_lnum_t npoint = 450;
	cs_lnum_t cell_id;
	int stevec_celic = 0; 
	static double z_arr[5150];
	int st=0;
	double dr = 0.0000025;
	for (int jot = 0; jot<51; jot+=1){
		double radius_curr = radius1[jot];
		for (double z_dir = 0.000; z_dir<axicoo; z_dir+=0.0003){
			sum_phi_u = 0.0;
			sum_phi_v = 0.0;
			sum_phi_w = 0.0;
			sum_phi_t = 0.0;
			mean_phi_u = 0.0;
			mean_phi_v = 0.0;
			mean_phi_w = 0.0;
			u_fluct_rms = 0;
			v_fluct_rms = 0;
			w_fluct_rms = 0;
			t_fluct_rms = 0;
			razd = 0;
			double zzz = 0;
	   	    k = 0;
			int cell_id1 = -999;
			int rank_id1 = -999;
			for (cs_lnum_t point_id = 0; point_id < npoint; point_id++) {
			  if(nt_curr==10){
				  cs_real_3_t xyz = {radius_curr*cos((double)point_id/npoint*2.*_PI),
									 radius_curr*sin((double)point_id/npoint*2.*_PI),
																			  z_dir};
				  cell_id = 0;
				  a1 = (int)rank;
				  cs_real_t dis2mn = 100;
				  _findpt_r(xyz, &cell_id, &a1, dis2mn);
				  cell_ids[stevec_celic] = cell_id;
				  ranks[stevec_celic] = a1;
				  stevec_celic++;
			  }
			  else{
				cell_id = (cs_lnum_t)cell_ids[stevec_celic];
				a1 = (cs_lnum_t)ranks[stevec_celic];
				stevec_celic++;
			  }
			  if ((int)cell_id != (int)cell_id1 || (int)rank_id1 != (int)a1) {
				if((int)rank==(int)a1){
					cell_id1 = cell_id;
					rank_id1 = a1;
					k+=1;
					ux = vel[cell_id][0];
					vy = vel[cell_id][1];
					w = vel[cell_id][2];
					theta = (double)point_id/npoint*2.*_PI;
					u =ux*cos(theta) + vy*sin(theta);
					v =-ux*sin(theta) + vy*cos(theta);
					t1 = t[cell_id];
					sum_phi_u += (double)u;
					sum_phi_v += (double)v;
					sum_phi_w += (double)w;
					sum_phi_t += (double)t1;
					u_fluct_rms += pow((double)u, 2);
					v_fluct_rms += pow((double)v, 2);
					w_fluct_rms += pow((double)w, 2);
					t_fluct_rms += pow((double)t1, 2);
					razd += sqrt((double)cell_cen[cell_id][0]*(double)cell_cen[cell_id][0] + (double)cell_cen[cell_id][1]*(double)cell_cen[cell_id][1]);
					zzz+= (double)cell_cen[cell_id][2];
				}
			  } 
			}/* end of loop on point_id */
			cs_parall_sum(1, CS_DOUBLE, &u_fluct_rms);
			cs_parall_sum(1, CS_DOUBLE, &v_fluct_rms);
			cs_parall_sum(1, CS_DOUBLE, &w_fluct_rms);
			cs_parall_sum(1, CS_DOUBLE, &t_fluct_rms);
			cs_parall_sum(1, CS_DOUBLE, &sum_phi_u);
			cs_parall_sum(1, CS_DOUBLE, &sum_phi_v);
			cs_parall_sum(1, CS_DOUBLE, &sum_phi_w);
			cs_parall_sum(1, CS_DOUBLE, &sum_phi_t);
			cs_parall_sum(1, CS_DOUBLE, &razd);
			cs_parall_sum(1, CS_INT32, &k);
			cs_parall_sum(1, CS_DOUBLE, &zzz);
			mean_phi_u = sum_phi_u/k;
			mean_phi_v = sum_phi_v/k;
			mean_phi_w = sum_phi_w/k;
			mean_phi_t = sum_phi_t/k;
			u_fluct_rms_phi= (double)u_fluct_rms/k;
			v_fluct_rms_phi = (double)v_fluct_rms/k;
			w_fluct_rms_phi = (double)w_fluct_rms/k;
			t_fluct_rms_phi = (double)t_fluct_rms/k;
			u_fluct_var_phi = (double)u_fluct_rms_phi;
			v_fluct_var_phi = (double)v_fluct_rms_phi;
			w_fluct_var_phi = (double)w_fluct_rms_phi;
			t_fluct_var_phi = (double)t_fluct_rms_phi;
			u_fluct_rms_phi= (double)u_fluct_rms_phi - pow((double)mean_phi_u, 2);
			v_fluct_rms_phi= (double)v_fluct_rms_phi - pow((double)mean_phi_v, 2);
			w_fluct_rms_phi= (double)w_fluct_rms_phi - pow((double)mean_phi_w, 2);
			t_fluct_rms_phi = (double)t_fluct_rms_phi - pow((double)mean_phi_t, 2);
			if(nt_curr==10 && rank<1){
				mean_phi_u_arr_t[st] = (double)mean_phi_u;
				mean_phi_v_arr_t[st] = (double)mean_phi_v;
				mean_phi_w_arr_t[st] = (double)mean_phi_w;
				mean_phi_t_arr_t[st] = (double)mean_phi_t;
				u_fluct_rms_arr[st] = u_fluct_rms_phi;
				v_fluct_rms_arr[st] = v_fluct_rms_phi;
				w_fluct_rms_arr[st] = w_fluct_rms_phi;
				t_fluct_rms_arr[st] = t_fluct_rms_phi;
				u_fluct_var_arr[st] = u_fluct_var_phi;
				v_fluct_var_arr[st] = v_fluct_var_phi;
				w_fluct_var_arr[st] = w_fluct_var_phi;
				t_fluct_var_arr[st] = t_fluct_var_phi;
				razd_arr[st] = (double)razd/k;
				k_arr[st] = k;
				z_arr[st] = (double)zzz/k;
			}
			else if(rank<1){
				mean_phi_u_arr_t[st] += (double)mean_phi_u;
				mean_phi_v_arr_t[st] += (double)mean_phi_v;
				mean_phi_w_arr_t[st] += (double)mean_phi_w;
				mean_phi_t_arr_t[st] += (double)mean_phi_t;
				u_fluct_rms_arr[st] += u_fluct_rms_phi;
				v_fluct_rms_arr[st] += v_fluct_rms_phi;
				w_fluct_rms_arr[st] += w_fluct_rms_phi;
				t_fluct_rms_arr[st] += t_fluct_rms_phi;
				u_fluct_var_arr[st] += u_fluct_var_phi;
				v_fluct_var_arr[st] += v_fluct_var_phi;
				w_fluct_var_arr[st] += w_fluct_var_phi;
				t_fluct_var_arr[st] += t_fluct_var_phi;
			}
			st+=1;
		}
	}
   	if(nt_curr==10 && rank<1){
		FILE *f2 = NULL;
		char filename2[100];
		sprintf(filename2,"razd_%d.txt", (int)nt_curr);
		f2 = fopen(filename2,"w");
		for (int j = 0; j<5100; j++){
			fprintf(f2, "%e \t %e \t %e \n", razd_arr[j], k_arr[j], z_arr[j]);
		}
		fclose(f2);
	}
   	if(fmod(nt_curr, 1000)==0 && rank<1){
		FILE *f1 = NULL;
		char filename1[100];
		sprintf(filename1,"povprecja_%d.txt", (int)nt_curr);
		f1 = fopen(filename1,"w");
		for (int j = 0; j<5100; j++){
			fprintf(f1, "%e \t %e \t %e \t %e \t %e \t %e  \t %e \t %e \t %e \t %e \t %e \t %e \n", mean_phi_u_arr_t[j], mean_phi_v_arr_t[j], mean_phi_w_arr_t[j], u_fluct_rms_arr[j], v_fluct_rms_arr[j], w_fluct_rms_arr[j], u_fluct_var_arr[j], v_fluct_var_arr[j], w_fluct_var_arr[j], mean_phi_t_arr_t[j], t_fluct_rms_arr[j], t_fluct_var_arr[j]);
		}
		fclose(f1);
	}
  }
}
END_C_DECLS