/*********************************************************************
UDF for computing the biot number along a particle trajectory
**********************************************************************/
#include "udf.h"
#include "math.h"
#include "dpm_types.h"
#include "dpm_laws.h"
#include "mem.h"
#include "dpm.h"
#include "materials.h"

DEFINE_INIT(biot_setup,domain)
{
  /* if memory for the particle variable titles has not been
   * allocated yet, do it now */
  	if (NULLP(user_particle_vars)) Init_User_Particle_Vars();
 		/* now set the name and label */
  		strcpy(user_particle_vars[0].name,"biot-number");
  		strcpy(user_particle_vars[0].label,"Biot number");
  		strcpy(user_particle_vars[1].name,"test-value");
  		strcpy(user_particle_vars[1].label,"Test value");
		Message0("\nBiot number UDF: Particle variable titles have been allocated.");
}

/* update the user scalar variables */
DEFINE_DPM_SCALAR_UPDATE(biot_number,cell,thread,initialize,p)
{
  cphase_state_t *c = &(p->cphase);
  if (initialize)
    {
      	/* this is the initialization call, set:
	* P_USER_REAL(p,0) contains the biot number, initialize to 0
	* P_USER_REAL(p,1) contains the test value, initialize to 0*/
      	P_USER_REAL(p,0) = 0.;
      	P_USER_REAL(p,1) = 0.;
	Message0("\nBiot number UDF: Particle variable has been initialized.");
    }
  else
    {
      	
	// Thermal conductivity of particle material (order of magnitude, INSPYRO analysis)
	real lambdaP = 1.0;
	// Sphericity of the particles
	real psi = 0.806;

	// Determine cell in which particle is located
	cell_t ci = P_CELL(p); 		
	// Determine thread in which particle is locted
	Thread *t0 = P_CELL_THREAD(p); 
	
	// Determine relative velocity between particle and fluid for Reynolds number determination
	real v_rel = sqrt(pow((P_VEL(p)[0] - C_U(ci,t0)),2) + pow((P_VEL(p)[1] - C_V(ci,t0)),2) + pow((P_VEL(p)[2] - C_W(ci,t0)),2)); 
	// Determine gas density at particle location
	real rhoF = C_R(ci,t0); 
	// Determine gas viscosity at particle location
	real mueF = C_MU_L(ci,t0);
	// Determine gas thermal conductivity
	real lambdaF = C_K_L(ci,t0);
	// Determine gas specific heat capacity
	real cpF = C_CP(ci,t0);

	// Particle reynolds number
	real ReP = P_DIAM(p)*v_rel*rhoF/mueF;

	// Prandtl number
	real Pr =  mueF*cpF/lambdaF;

	// Nusselt number (Ranz-Marshall correlation)
	real Nu = (2+0.6*pow(ReP,0.5)*pow(Pr,0.3333));

	// Particle volume
	real VP = P_MASS(p)/P_RHO(p);
	
	// Approximate particle surface area, accounting for non-sphericity
	real AP = P_DIAM(p)*P_DIAM(p)*3.1415926/psi;

	// Characteristic lenght for Biot number
	real L = VP/AP;

	/* assign biot number to the first user value */
      	P_USER_REAL(p,0) = (Nu*lambdaF/P_DIAM(p))*L/lambdaP;
      	/* assign value to the second user value */
      	P_USER_REAL(p,1) = Nu;

	Message0("\nBiot number UDF: Particle variable has been modified.");
    }
}

/* write biot number when sorting particles at surfaces */
DEFINE_DPM_OUTPUT(biot_output,header,fp,p,thread,plane)
{
  char name[100];
  if (header)
    {
      if (NNULLP(thread))
	par_fprintf_head(fp,"(%s %d)\n",THREAD_HEAD(thread)->
			 dpm_summary.sort_file_name,11);
      else
	par_fprintf_head(fp,"(%s %d)\n",plane->sort_file_name,11);
      par_fprintf_head(fp,"(%10s %10s %10s %10s %10s %10s %10s"
		       " %10s %10s %10s %10s %s)\n",
		       "X","Y","Z","U","V","W","diameter","T","mass-flow",
		       "time","biot-number","name");
    }
  else
    {
      sprintf(name,"%s:%d",P_INJECTION(p)->name,p->part_id);

      /* add P_INJ_ID(P_INJECTION(p)) and part_id for sorting in parallel */
      par_fprintf(fp,
		  "%d %d ((%10.6g %10.6g %10.6g %10.6g %10.6g %10.6g "
		  "%10.6g %10.6g %10.6g %10.6g %10.6g) %s)\n",
		  P_INJ_ID(P_INJECTION(p)), p->part_id,
		  P_POS(p)[0], P_POS(p)[1], P_POS(p)[2],
		  P_VEL(p)[0], P_VEL(p)[1], P_VEL(p)[2],
		  P_DIAM(p), P_T(p), P_FLOW_RATE(p), P_TIME(p),
		  P_USER_REAL(p,0), name);
    }
}