From 6faf03b05fd6df473841accf5726c18c774e032e Mon Sep 17 00:00:00 2001
From: Weiye Wang <weiyewan@usc.edu>
Date: Wed, 12 Oct 2022 10:59:46 -0700
Subject: [PATCH] Add the HRR

---
 UserData.cpp |  1 +
 UserData.h   |  2 ++
 macros.h     |  2 ++
 residue.cpp  | 21 +++++++++++++++++++++
 4 files changed, 26 insertions(+)

diff --git a/UserData.cpp b/UserData.cpp
index de04232..b91d4cb 100644
--- a/UserData.cpp
+++ b/UserData.cpp
@@ -354,6 +354,7 @@ UserData allocateUserData(FILE *input){
 	//data->ratesOutput=fopen("rates.dat","w");
 	
 	data->innerMassFractions = new double [data->nsp];
+	data->HRRdata = new double [data->npts];
 
 	return(data);
 }
diff --git a/UserData.h b/UserData.h
index 7578c4c..f2b73e0 100644
--- a/UserData.h
+++ b/UserData.h
@@ -162,6 +162,8 @@ typedef struct UserDataTag{
 	double flamePosition[2];
 	double flameTime[2];
 	size_t nTimeSteps;
+	/*heat release rate (HRR) data */
+	double* HRRdata;
 } *UserData;
 UserData allocateUserData(FILE *input);
 void setSaneDefaults(UserData data);
diff --git a/macros.h b/macros.h
index dcf36db..35cad0d 100644
--- a/macros.h
+++ b/macros.h
@@ -21,6 +21,7 @@
 #define Y(i,k) ydata[((i-1)*data->nvar)+data->ny+k-1]
 #define R(i)   ydata[((i-1)*data->nvar)+data->nr]
 #define P(i)   ydata[((i-1)*data->nvar)+data->np]
+#define HRRdata(i)  HRRdata[i-1]
 
 #define Tdot(i)   ydotdata[((i-1)*data->nvar)+data->nt]
 #define Ydot(i,k) ydotdata[((i-1)*data->nvar)+data->ny+k-1]
@@ -52,6 +53,7 @@
 #define energy(i) energy[i-1]
 #define Cp(i) Cp[i-1]
 
+
 #define atolT(i)   atolvdata[((i-1)*data->nvar)+data->nt]
 #define atolY(i,k)   atolvdata[((i-1)*data->nvar)+data->ny+k-1]
 #define atolR(i)   atolvdata[((i-1)*data->nvar)+data->nr]
diff --git a/residue.cpp b/residue.cpp
index 6e205b2..a3b5dbb 100644
--- a/residue.cpp
+++ b/residue.cpp
@@ -741,6 +741,11 @@ int residue(double t,
 	double cendfm, cendfc, cendfp;
 	/*Aliases for various grid spacings:*/
 	double dpsip, dpsiav, dpsipm, dpsim, dpsimm;
+	/*define the heat release rate related parameters*/
+	double Tsp=298.0; 
+	double HRR = 0 ; 
+	double Hf = 0 ;
+
 	dpsip=dpsiav=dpsipm=dpsim=dpsimm=ONE;
 	double mass, mdotIn;
 	double sum, sum1, sum2, sum3;
@@ -751,6 +756,11 @@ int residue(double t,
 	mass=data->mass;				//Units: kg
 	mdotIn=data->mdot*calc_area(R(npts),&m);	//Units: kg/s
 
+	/*Initialize the HRR data*/ 
+	for (j=1; j<= npts ; j++) {
+		HRRdata(j) = 0 ;
+	}
+
 //	/*evaluate properties at j=1*************************/
 	setGas(data,ydata,1);
 	rhom=data->gas->density();                                       
@@ -888,6 +898,14 @@ int residue(double t,
 			+advTerm 
 			-diffTerm
 			+srcTerm;
+		
+		
+		/*Calculate the Heat Release Rate */
+		for (k = 1; k <=nsp; k++) {
+			Hf = enthalpy(k) - Cp(k) * (T(j)-Tsp) ; 
+			HRR = wdot(k) * Hf ; 
+			HRRdata(j) = HRR ; 
+		}
 
 	//	//energy formulation:
 	//	tranTerm = Tdot(j);
@@ -1030,7 +1048,9 @@ void printSpaceTimeHeader(UserData data)
 	for (size_t k = 1; k <=data->nsp; k++) {
 		fprintf((data->output), "%15s\t",data->gas->speciesName(k-1).c_str());
 	}
+	fprintf((data->output), "%15s\t","HRR()");
 	fprintf((data->output), "%15s\n","Pressure(Pa)");
+
 }
 
 void printSpaceTimeOutput(double t, N_Vector* y, FILE* output, UserData data)
@@ -1055,6 +1075,7 @@ void printSpaceTimeOutput(double t, N_Vector* y, FILE* output, UserData data)
 		for (size_t j = 0; j < data->nvar; j++) {
 			fprintf(output, "%15.9e\t",ydata[j+i*data->nvar]);
 		}
+		fprintf(output,"%15.6e\t",HRRdata(i+1));
 		fprintf(output, "\n");
 	}
 	fprintf(output, "\n\n");