From f9d512661967e2050b55cc6bad44a3451134b94b Mon Sep 17 00:00:00 2001 From: Weiye Wang Date: Thu, 6 Apr 2023 11:13:48 -0700 Subject: [PATCH] Characteristic time scales --- UserData.cpp | 9 ++++++ UserData.h | 10 +++++++ macros.h | 8 ++++++ main.cpp | 11 ++++++++ residue.cpp | 77 ++++++++++++++++++++++++++++++++++++++++++++++++++++ residue.h | 4 +++ 6 files changed, 119 insertions(+) diff --git a/UserData.cpp b/UserData.cpp index 8d2dc60..c0aa148 100644 --- a/UserData.cpp +++ b/UserData.cpp @@ -53,6 +53,10 @@ void freeUserData(UserData data){ fclose(data->globalOutput); printf("Global Output File Cleared from Memory!\n"); } + if(data->timescaleOutput!=NULL){ + fclose(data->timescaleOutput); + printf("Characteristic Timescale Output File Cleared from Memory!\n"); + } } free(data); /* Free the user data */ @@ -383,9 +387,14 @@ UserData allocateUserData(FILE *input){ data->output=fopen("output.dat","w"); data->globalOutput=fopen("globalOutput.dat","w"); data->gridOutput=fopen("grid.dat","w"); + data->timescaleOutput=fopen("timeScale.dat","w") ; //data->ratesOutput=fopen("rates.dat","w"); data->innerMassFractions = new double [data->nsp]; + //data->wdot_mole = new double [data->nsp] ; + //data->wdot_mass = new double [data->nsp] ; + data->time_scale = new double [(data->npts) * (data->nsp)] ; + //data->MW = new double [data->nsp] ; return(data); } diff --git a/UserData.h b/UserData.h index c7b13a3..da595f3 100644 --- a/UserData.h +++ b/UserData.h @@ -184,6 +184,16 @@ typedef struct UserDataTag{ double flamePosition[2]; double flameTime[2]; size_t nTimeSteps; + + /*Following arrays are used to compute the characteristic time scale of + species*/ + //double* wdot_mole ; + //double* wdot_mass ; + double* time_scale ; + //double* MW; + FILE* timescaleOutput; + + } *UserData; UserData allocateUserData(FILE *input); void setSaneDefaults(UserData data); diff --git a/macros.h b/macros.h index b70d477..3c827b8 100644 --- a/macros.h +++ b/macros.h @@ -66,4 +66,12 @@ #define constraintsY(i,k) constraintsdata[((i-1)*data->nvar)+data->ny+k-1] #define constraintsR(i) constraintsdata[((i-1)*data->nvar)+data->nr] +/*Following marcos are defined to calculate the characteristic time-scale */ +#define wdot_mole(i) wdot_mole[i-1] +#define wdot_mass(i) wdot_mass[i-1] +#define MW(i) MW[i-1] +#define time_scale(i,k) time_scale[(i-1)*data->nsp+k-1] +#define concentra(i) concentra[i-1] + + #endif diff --git a/main.cpp b/main.cpp index 4f3b809..6f101de 100644 --- a/main.cpp +++ b/main.cpp @@ -175,8 +175,11 @@ int main(){ printSpaceTimeHeader(data); printGlobalHeader(data); + printTimescaleHeader(data); printSpaceTimeOutput(tNow, &y, data->output, data); printSpaceTimeOutput(tNow, &y, data->gridOutput, data); + printTimescaleOutput(tNow, &y, data->timescaleOutput,data); + if(!data->dryRun){ count=0; @@ -304,6 +307,14 @@ int main(){ //} } + if(count%data->nSaves==0 ){ + printTimescaleOutput(tNow,&y, data->timescaleOutput,data); + //printSpaceTimeOutput(tNow, &y, data->output, data); + //if(data->writeRates){ + // printSpaceTimeRates(tNow, ydot, data); + //} + } + /*Print global variables only if time-step is of high order!*/ if(data->nTimeSteps==0){ data->flamePosition[0]=0.0e0; diff --git a/residue.cpp b/residue.cpp index 1500a1a..237ff30 100644 --- a/residue.cpp +++ b/residue.cpp @@ -1444,3 +1444,80 @@ void printGlobalVariables(double t, N_Vector* y, N_Vector* ydot, UserData data) //// Ydot(npts-1,k)=Ydot(npts,k); //// } ////} + + +/*Following functions are added to derive the characteristic time scale of species*/ +void getTimescale(UserData data, N_Vector* y){ + size_t i, k, nsp,npts ; + nsp = data->nsp ; + npts = data->npts ; + double rho, wdot_mole[nsp], wdot_mass[nsp],MW[nsp],concentra[nsp]; + //double time_scale[npts*nsp] ; + + double* ydata ; + ydata = N_VGetArrayPointer_OpenMP(*y); //access the data stored in N_Vector* y + + for(i=1; i<= npts;i++){ + setGas(data,ydata,i); //set the gas state at each grid point + rho = data->gas->density(); //get the averaged density at each grid point, Unit:kg/m^3 + data->gas->getNetProductionRates(wdot_mole) ; //Unit:kmol/(m^3 * s) + + for(k=1; k<= nsp; k++){ + MW(k) = data->gas->molecularWeight(k-1) ; //Unit:kg/kmol + } + + for(k=1; k<= nsp; k++){ + wdot_mass(k) = wdot_mole(k) * MW(k) ; //Unit:kg/(m^3 * s) + } + + for(k=1; k<= nsp; k++){ + concentra(k) = Y(i,k) * rho ; //Unit:kg/m^3 + } + + for(k=1;k<= nsp;k++){ + data->time_scale(i,k) = concentra(k)/(wdot_mass(k)+ 1.00e-16) ; + } + + + } + +} + + +void printTimescaleHeader(UserData data) +{ + fprintf((data->timescaleOutput), "%15s\t","#1"); + for (size_t k = 1; k <=data->nsp+1; k++) { + fprintf((data->timescaleOutput), "%15lu\t",k+1); + } + fprintf((data->timescaleOutput), "%15lu\n",data->nsp+3); + + + fprintf((data->timescaleOutput), "%15s\t%15s\t%15s\t","#time","radius","Temp(K)"); + //fprintf((data->output), "%15s\t%15s\t","radius(m)","Temp(K)"); + for (size_t k = 1; k <=(data->nsp); k++) { + fprintf((data->timescaleOutput), "%15s\t",data->gas->speciesName(k-1).c_str()); + } + //fprintf((data->output), "%15s\t","Pressure(Pa)"); + //fprintf((data->timescaleOutput), "%15s\n",data->gas->speciesName(data->nsp-1).c_str()); + //fprintf((data->output), "%15s\n","Mdot (kg/s)"); + fprintf((data->timescaleOutput), "\n"); +} + +void printTimescaleOutput(double t,N_Vector* y,FILE* output,UserData data) +{ + double* ydata; + ydata = N_VGetArrayPointer_OpenMP(*y); + + for(size_t i=1 ; i<= data->npts; i++){ + fprintf(output, "%15.9e\t%15.9e\t",t,R(i)); + + for(size_t k=1; k<=data->nsp;k++){ + fprintf(output, "%15.9e\t",data->time_scale(i,k)); + } + fprintf(output, "\n"); + } + fprintf(output, "\n"); + + +} \ No newline at end of file diff --git a/residue.h b/residue.h index 6222fc9..4b168b7 100644 --- a/residue.h +++ b/residue.h @@ -88,3 +88,7 @@ void printSpaceTimeOutputInterpolated(double t, N_Vector y, UserData data); void writeRestart(double t, N_Vector* y, N_Vector* ydot, FILE* output, UserData data); void readRestart(N_Vector* y, N_Vector* ydot, FILE* input, UserData data); + +void getTimescale(UserData data, N_Vector* y); +void printTimescaleHeader(UserData data); +void printTimescaleOutput(double t,N_Vector* y,FILE* output,UserData data); \ No newline at end of file