/* --------------------------------------------
* File Name : global.h
* Purpose : contains all the parameters, reference values and boundary conditions
* Output files : -
* Creation Date : 19/03/2020
* Created by : Victoria Hamtiaux
 ----------------------------------------------*/

#ifndef GLOBAL_H_
#define GLOBAL_H_

#include <math.h>
#include <petsc.h>
#include <petscsys.h>
#include <stdio.h>
#include <stdlib.h>
#include "global_evap.h"

/****************
  COEFFICIENTS
****************/
#define PI 3.14159265
#define gama 1.4

/****************
  TIME STABILITY
****************/
const double CFL = 0.0001;           //Coefficient for time stability condition in time step calculation : convection time (Courant–Friedrichs–Lewy condition)
const double Fr = 0.125;          //Coefficient for time stability condition in time step calculation : diffusion time (Fourier number)
const double Tol = 1e-7;          //Tolerance for the Poisson solver (Infinity norm)
const double explicit_diff = 0.;  // 0 or 1: implicit or explicit treatment of the diffusive terms in y direction in the momentum equation for u and w respectively.
const int corrector_step = 1;     // 0 or 1: only predictor or predictor+corrector steps

/****************
     FEATURES
****************/
const int forcing = 0;		        // 0 or 1 : without or with forcing(external pressure gradient applied)
const int c_xy = 0;		            // 0 or 1 or 2 : In 3D or in 2D in x-y direction or in 2D in z-y direction
const int turb_stat = 1;		      // 0 or 1 : without or with computation of turbulence statistics
const int fluid_type = 0;         // 0 or 1 : water or air
const int probes = 0;             // 0 or 1: without or with record of vertical velocity at probes
const int step_probes = 20;       // iterations in between output sampling
const double phi_ref = 1.;        // for no porous media = 1. --- otherwize volume fraction of the fluid phase in [-]
                                  //      -> properties to define in global_por.h
const int evaporation = 0;        // 0 or 1: without or with evaporation at the free-surface
                                  //      -> properties to define in global_evap.h
const int descentFS = 0;          // 0 or 1: without or with descent of FS (only taken into acount if evaporation == 1)
                                  //      -> properties to define in global_evap.h
const int particles = 1;          // 0 or 1: without or with particles
                                  //      -> properties to define in global_particles.h
const int wm_LES = 0;             // 0 or 1: without or with the wall model
const int wm_temp_LES = 0;        // 0 or 1: without or with the temperature (never use wm_temp without wm!)
                                  //      -> properties to define in global_wm.h
const int LES_Dyn = 0;            // 0 or 1 : without or with LES model
// TYPE LES model
// s : smagorinsky
// l : lagrangian
// w : wale
// d : dynamic smago
const char LES_mod = 'l';

/*********************
 REFERENCE PROPERTIES
*********************/
const double l_ref= 2.78e-02;     // in [m] : reference length
const double T_ref= 373.15; 		  // in [K] : reference temperature
const double p_atm= 101325.;      // in [Pa] : atmospheric pressure
const double g = 9.81;            // in [m/s^2] : gravity force

/***************** Interpolation for temperature range  *****************/
//Fits for  373.1 up to 373.2
// Interpolation constants for pure water density :
const double a_rhoD =  8.673709208309e+02 ;      //[kg/m^3]
const double b_rhoD =  1.018730297243e+00 ;      //[kg/m^3 / K]
const double c_rhoD =  -7.664264437836e-04 ;     //[kg/m^3 / K^2]
const double d_rhoD =  -5.679861788799e-06 ;     //[kg/m^3 / K^3]
const double e_rhoD =  5.811569567603e-09 ;      //[kg/m^3 / K^4]
// Interpolation constants for pure water specific heat capacity :
const double a_cpD =  4.496815829050e+03 ;       //[J/kg K]
const double b_cpD =  -8.937322550577e-01 ;      //[J/kg K / K]
const double c_cpD =  -4.230539999120e-03 ;      //[J/kg K / K^2]
const double d_cpD =  1.129738286552e-05 ;       //[J/kg K / K^3]
const double e_cpD =  2.857596308360e-09 ;       //[J/kg K / K^4]
// Interpolation constants for pure water dynamic viscosity :
const double a_muD =  7.047336658260e-05 ;       //[kg/m s]
const double b_muD =  1.888859999534e-07 ;       //[kg/m s / K]
const double c_muD =  5.062610558921e-10 ;       //[kg/m s / K^2]
const double d_muD =  1.356904465002e-12 ;       //[kg/m s / K^3]
const double e_muD =  0.000000000000e+00 ;       //[kg/m s / K^4]
// Interpolation constants for pure water thermal conductivity :
const double a_kD =  -3.641581157026e-01 ;       //[W/m K ]
const double b_kD =  4.178010984685e-03 ;        //[W/m K / K]
const double c_kD =  1.817489777242e-06 ;        //[W/m K / K^2]
const double d_kD =  -2.251132067554e-08 ;       //[W/m K / K^3]
const double e_kD =  2.060340114369e-11 ;        //[W/m K / K^4]
// Interpolation constants for pure water Thermal expansion coef :
const double a_betaD =  -1.301227686814e-03 ;    //[1/K]
const double b_betaD =  3.610534787800e-06 ;     //[1/K / K]
const double c_betaD =  1.084580726202e-08 ;     //[1/K / K^2]
const double d_betaD =  -1.563770657069e-11 ;    //[1/K / K^3]

const double rho_ref= a_rhoD + b_rhoD * T_ref + c_rhoD * PetscSqr(T_ref) + d_rhoD * T_ref*PetscSqr(T_ref) + e_rhoD * PetscSqr(PetscSqr(T_ref));		                    // in [kg/m^3] : reference density : Warning those are values for liquid water at T_ref, it is preferable not to change them (since calculated in code for water)
const double mu_ref= a_muD + b_muD * T_ref + c_muD * PetscSqr(T_ref) + d_muD * T_ref*PetscSqr(T_ref) + e_muD * PetscSqr(PetscSqr(T_ref)) ;                            	// in [kg/m/s] : reference dynamic viscosity : Warning those are values for liquid water at T_ref, it is preferable not to change them (since calculated in code for water)
const double cp_ref= a_cpD + b_cpD * T_ref + c_cpD * PetscSqr(T_ref) + d_cpD * T_ref*PetscSqr(T_ref) + e_cpD * PetscSqr(PetscSqr(T_ref)) ;      // in [J/kg/K] : reference heat capacity : Warning those are values for liquid water at T_ref, it is preferable not to change them (since calculated in code for water)
const double k_ref= a_kD + b_kD * T_ref + c_kD * PetscSqr(T_ref) + d_kD * T_ref*PetscSqr(T_ref) + e_kD * PetscSqr(PetscSqr(T_ref)) ; 	          // in [W/m/K] = [J/s/m/K] : reference thermal conductivity : Warning those are values for liquid water at T_ref, it is preferable not to change them (since calculated in code for water)
const double beta_ref = a_betaD + b_betaD * T_ref + c_betaD * PetscSqr(T_ref) + d_betaD * T_ref*PetscSqr(T_ref) ;                                               // in [1/K] : Expansion coefficient

/*********************
 BOUNDARY CONDITIONS
*********************/

/***************** VELOCITY *****************/
// defining the U boundary condition  at N boundary
// w : wall
// f : freeslip
// o : outflow
const char Bcomb_U_N = 'f';

// defining the U boundary condition  at S boundary
// w : wall
// o : outflow
const char Bcomb_U_S = 'w';

// defining the U boundary condition  at W-E boundary
// w : wall-wall
// p : periodic
// i : inflow-outflow
const char Bcomb_U_EW = 'w';

// defining the U boundary condition  at F-B boundary
// w : wall-wall
// p : periodic
// o : outflow-outflow
const char Bcomb_U_FR = 'w';

/***************** TEMPERATURE *****************/
// n : Neumann
// d : Dirichlet
// p : periodic (for TR ot TW only)
// North :
const char Bcomb_TN = 'd' ;
const double grad_N = 0. * l_ref / T_ref;
const double temperature_N = 373.11658 / T_ref;
const double Tval_N = grad_N * (Bcomb_TN == 'n') + temperature_N* (Bcomb_TN == 'd');

// South :
const char Bcomb_TS = 'd' ;
const double grad_S = 0.;
const double temperature_S = 373.18342 / T_ref;
const double Tval_S = grad_S * (Bcomb_TS == 'n') + temperature_S* (Bcomb_TS == 'd');

// Background :
const char Bcomb_TR = 'n' ;                       //Warning!!! if bcomb_U_FB == 'p',every variable is periodic so temperature and pressure too
const double Tval_R = 0.;

// Foreground :
const char Bcomb_TF = 'n' ;
const double Tval_F = 0.;

// West :
const char Bcomb_TW = 'n' ;                       //Warning!!! if bcomb_U_EW == 'p',every variable is periodic so temperature and pressure too
const double Tval_W = 0.;

// East :
const char Bcomb_TE = 'n' ;
const double Tval_E = 0.;

/**********************************
  COMPUTATION OF VARIOUS VARIABLES
**********************************/
const double alpha_ref = k_ref / rho_ref / cp_ref ;                       // in [m^2/s]: thermal diffusivity
const double nu_ref = mu_ref/rho_ref;							                        // in [m^2/s]: kinematic viscosity
const double dT = (temperature_S - temperature_N >= 0 ? temperature_S - temperature_N : temperature_N - temperature_S)*T_ref; // in [K]: temperature difference between north and south boundaries
const double Ra_l = g*beta_ref*dT*(l_ref*l_ref*l_ref)/nu_ref/alpha_ref ;  // in [-]: Rayleigh number
//const double u_ref = PetscSqrtReal(g * l_ref * beta_ref * dT);            // in [m/s]: free-fall velocity use as
// FAKE
const double u_ref = (g * l_ref * beta_ref * dT);            // in [m/s]: free-fall velocity use as reference velocity
const double t_ref = l_ref / u_ref;                                       // in [s]: free-fall time
const double Re = rho_ref * u_ref * l_ref / mu_ref ;                      // in [-]: Reynolds number
const double Pr = mu_ref * cp_ref / k_ref ;                               // in [-]: Prandtl number
const double Ri = g * l_ref / (u_ref * u_ref);                            // in [-]: Richardson number
const double p_ref = u_ref * u_ref * rho_ref;                            	// in [Pa]: reference pressure

/**************************************
   NON-DIM of PREVIOUS INTERPOLATIONS
**************************************/
// Interpolation constants for pure water density :
const double a_rho = a_rhoD / rho_ref ;
const double b_rho = b_rhoD / rho_ref * T_ref;
const double c_rho = c_rhoD / rho_ref * PetscSqr(T_ref);
const double d_rho = d_rhoD / rho_ref * T_ref*PetscSqr(T_ref);
const double e_rho = e_rhoD / rho_ref * PetscSqr(PetscSqr(T_ref));

// Interpolation constants for pure water specific heat capacity :
const double a_cp = a_cpD / cp_ref ;
const double b_cp = b_cpD / cp_ref * T_ref;
const double c_cp = c_cpD / cp_ref * PetscSqr(T_ref);
const double d_cp = d_cpD / cp_ref * T_ref*PetscSqr(T_ref);
const double e_cp = e_cpD / cp_ref * PetscSqr(PetscSqr(T_ref));
// Interpolation constants for pure water dynamic viscosity :
const double a_mu = a_muD / mu_ref;
const double b_mu = b_muD / mu_ref * T_ref;
const double c_mu = c_muD / mu_ref * PetscSqr(T_ref);
const double d_mu = d_muD / mu_ref * T_ref*PetscSqr(T_ref);
const double e_mu = e_muD / mu_ref * PetscSqr(PetscSqr(T_ref));

// Interpolation constants for pure water thermal conductivity :
const double a_k = a_kD / k_ref;
const double b_k = b_kD / k_ref * T_ref;
const double c_k = c_kD / k_ref * PetscSqr(T_ref);
const double d_k = d_kD / k_ref * T_ref*PetscSqr(T_ref);
const double e_k = e_kD / k_ref * PetscSqr(PetscSqr(T_ref));


#endif /* GLOBAL_H_ */