diff --git a/README.md b/README.md new file mode 100644 index 0000000..9fbe56c --- /dev/null +++ b/README.md @@ -0,0 +1,27 @@ +This is the public repository for the code sensBVP. Its use is in +calculating sensitivity coefficients of ignition delay time to rate +coefficients. In sensBVP, this is carried out by transforming the +initial value problem to a boundary value problem, reducing the time to +calculate the sensitivity coefficients by at least an order of +magnitude. + +In order to install sensBVP, modify the various paths in the Makefile. +Make sure that SUNDIALS (v3.1.1 or higher), Cantera (2.3 or higher), and +GSL are installed. Run make all, followed by make install in the source +directory. + +After installing, the executable sensBVP can be executed as follows: + +sensBVP -a + -r + -f + -T + -P + -m + -c + -t + -s + +Execute sensBVP -h for all available options. + +Try the cases in the examples directory. diff --git a/examples/BVP/methane_BVP.sh b/examples/BVP/methane_BVP.sh new file mode 100755 index 0000000..ec59ad3 --- /dev/null +++ b/examples/BVP/methane_BVP.sh @@ -0,0 +1 @@ +sensBVP -r 1e-06 -a 1e-12 -f 1e-07 -T 900.0 -P 10.0 -m ./mech-FFCM1.cti -c CH4:0.0499002,O2:0.199601,N2:0.750499 -t 1 -s diff --git a/examples/BruteForce/methane_Brute_Force.sh b/examples/BruteForce/methane_Brute_Force.sh new file mode 100755 index 0000000..c94de5b --- /dev/null +++ b/examples/BruteForce/methane_Brute_Force.sh @@ -0,0 +1 @@ +sensBrute -r 1e-06 -a 1e-12 -T 900.0 -P 10.0 -m ./mech-FFCM1.cti -c CH4:0.0499002,O2:0.199601,N2:0.750499 -t 1 -s diff --git a/examples/JacobianBasedSensitivity/jacVarSens.py b/examples/JacobianBasedSensitivity/jacVarSens.py new file mode 100644 index 0000000..b200606 --- /dev/null +++ b/examples/JacobianBasedSensitivity/jacVarSens.py @@ -0,0 +1,64 @@ +""" +Constant-pressure, adiabatic kinetics simulation with sensitivity analysis +""" + +import sys +import time +import numpy as np + +start = time.time() +import cantera as ct + +gas = ct.Solution('./mech-FFCM1.cti') +temp = 900.0 +pres = 10.0*ct.one_atm + +gas.TPX = temp, pres, 'CH4:0.0499002,O2:0.199601,N2:0.750499' +r = ct.IdealGasConstPressureReactor(gas, name='R1') +sim = ct.ReactorNet([r]) + +# enable sensitivity with respect to all reactions +for i in range(0,gas.n_reactions): + r.add_sensitivity_reaction(i) + +# set the tolerances for the solution and for the sensitivity coefficients +sim.rtol = 1.0e-6 +sim.atol = 1.0e-12 +sim.rtol_sensitivity = 1.0e-6 +sim.atol_sensitivity = 1.0e-12 + +#states = ct.SolutionArray(gas, extra=['t','s2']) +ignitionStateFound=False +Told=temp +TIgn=1.115367e+03 #K +tEnd=1.0 #s +tauIgn=0.0 +nPts=1000 +dt=tEnd/(float(nPts)-1.0) +sens=np.zeros(gas.n_reactions) +out=open("ignitionSensitivities.dat","w") + +for t in np.arange(0, tEnd, dt): + + TOld=r.T + sim.advance(t) + TCurrent=r.T + for i in range(0,gas.n_reactions): + sens[i] = sim.sensitivity('temperature', i) + + if(ignitionStateFound==False): + if(r.T>=TIgn): + print("Ignition state found!\n") + print("T=%15.6e K\n"%(TCurrent)) + print("tau=%15.6e s\n"%(t)) + ignitionStateFound=True + tauIgn=t + dTdtau=(TCurrent-TOld)/dt + for i in range(0,gas.n_reactions): + sens[i]=sens[i]*(-1.0e0/dTdtau)*(TCurrent/tauIgn) + out.write("%15d\t%15.6e\n"%(i,sens[i])) + break + +out.close() +end = time.time() +print("Elapsed time: %15.6e s\n"%(end-start)) diff --git a/src/Makefile b/src/Makefile new file mode 100644 index 0000000..37f6ff0 --- /dev/null +++ b/src/Makefile @@ -0,0 +1,61 @@ +compiler =g++ +CANTERA_DIR =/opt/scientific/cantera-2.4_gnu_blas +CVODE_DIR =/opt/scientific/sundials-3.1.1_gnu_blas +KINSOL_DIR =/opt/scientific/sundials-3.1.1_gnu_blas +BVPEXE =sensBVP +BRUTEEXE =sensBrute +DESTDIR =~/bin + +CANTERA_INCLUDES=-I$(CANTERA_DIR)/include +CVODE_INCLUDES =-I$(CVODE_DIR)/include +KINSOL_INCLUDES =-I$(KINSOL_DIR)/include +CVODE_LIBS =-L$(CVODE_DIR)/lib -lsundials_nvecserial -lsundials_cvode +KINSOL_LIBS =-L$(KINSOL_DIR)/lib -lsundials_kinsol +CANTERA_LIBS =-L$(CANTERA_DIR)/lib -lcantera_shared +GSL_INCLUDES =-I/usr/include/gsl +GSL_LIBS =-L/usr/lib -lgsl -lgslcblas +RPATH =-Wl,-rpath=$(CVODE_DIR)/lib,-rpath=$(KINSOL_DIR)/lib +RM=rm -f + +compiler?=g++ + +ifeq ($(compiler),g++) + CPPFLAGS= -Wall -O3 + CPP=g++ +endif + +ifeq ($(compiler),icpc) + export GXX_INCLUDE=/usr/lib/gcc/x86_64-pc-linux-gnu/7.4.1/include/c++ + CPPFLAGS= -Wall -O3 -gxx-name=/usr/bin/g++-7 -std=c++11 + CPP=icpc +endif + +all: $(BVPEXE) $(BRUTEEXE) + +sensBVP.o: sensBVP.cpp + $(CPP) $(CPPFLAGS) $(CANTERA_INCLUDES) $(CVODE_INCLUDES) \ + $(KINSOL_INCLUDES) $(GSL_INCLUDES) \ + -c sensBVP.cpp -o sensBVP.o + +sensBrute.o: sensBrute.cpp + $(CPP) $(CPPFLAGS) $(CANTERA_INCLUDES) $(CVODE_INCLUDES) \ + $(GSL_INCLUDES) \ + -c sensBrute.cpp -o sensBrute.o + +$(BVPEXE): sensBVP.o + $(CPP) $(CPPFLAGS) \ + sensBVP.o -o $(BVPEXE) $(RPATH) $(CVODE_LIBS) \ + $(KINSOL_LIBS) $(CANTERA_LIBS) $(GSL_LIBS) + +$(BRUTEEXE): sensBrute.o + $(CPP) $(CPPFLAGS) \ + sensBrute.o -o $(BRUTEEXE) $(RPATH) $(CVODE_LIBS) \ + $(CANTERA_LIBS) $(GSL_LIBS) + +.PHONY: install +install: + cp $(BVPEXE) $(BRUTEEXE) $(DESTDIR) + +clean: + rm -f $(BVPEXE) $(BRUTEEXE) *.o *.d + diff --git a/src/sensBVP.cpp b/src/sensBVP.cpp index ba4c66a..ed025e5 100644 --- a/src/sensBVP.cpp +++ b/src/sensBVP.cpp @@ -1135,30 +1135,12 @@ static int residueKS(N_Vector yp, N_Vector res, void *user_data) static void lookupDpdt(UserData data, realtype t, realtype *P, realtype *dPdt) { - realtype *tData, *PData, *dPdtData; - tData = N_VGetArrayPointer_Serial(data->tArray); - PData = N_VGetArrayPointer_Serial(data->PArray); + realtype *dPdtData; dPdtData = N_VGetArrayPointer_Serial(data->dPdtArray); - int jguess=data->jguess; - int k=0; - int safel=0; - int nOrder=4; - realtype dy; int npts=data->tArraySize; if(t<=data->t1){ - //k=hunt(t, tData, npts, jguess); - ////printf("Desired value at %d: %15.6e\t%15.6e\t%15.6e\n",k,t,P(k),dPdt(k)); - //safel=IMIN(IMAX(k-(nOrder)/2,0),data->tArraySize+1-nOrder-1); - //data->jguess=safel; - ////printf("safel %d: \n",safel); - ////printf("safel val %15.6e: \n",*(tData+safel)); - //polint(tData+safel-1, PData+safel-1, nOrder+1, t, P, &dy); - ////printf("P error %15.6e: \n",dy); - //polint(tData+safel-1, dPdtData+safel-1, nOrder+1, t, dPdt, &dy); - ////printf("dPdt error %15.6e: \n",dy); - //// *P=gsl_spline_eval(data->spline,t,data->acc); *dPdt=gsl_spline_eval(data->splinedot,t,data->accdot); diff --git a/src/sensBrute.cpp b/src/sensBrute.cpp new file mode 100644 index 0000000..439aa3a --- /dev/null +++ b/src/sensBrute.cpp @@ -0,0 +1,1182 @@ +#include +#include +#include +#include +#include +#include + +/*Cantera include files*/ +#include + +/*Sundials include files*/ +#include /* prototypes for CVODE fcts., consts. */ +#include /* access to serial N_Vector */ +#include /* access to dense SUNMatrix */ +#include /* access to dense SUNLinearSolver */ +#include /* access to CVDls interface */ +#include /* defs. of realtype, sunindextype */ +#include + +#include +#include + +static int imaxarg1,imaxarg2; +#define IMAX(a,b) (imaxarg1=(a),imaxarg2=(b),(imaxarg1) > (imaxarg2) ?\ + (imaxarg1) : (imaxarg2)) + +static int iminarg1,iminarg2; +#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\ + (iminarg1) : (iminarg2)) + +/*quick-sort related routines*/ +static unsigned long *lvector(long nl, long nh); +static void free_lvector(unsigned long *v, long nl, long nh); +static void sort(unsigned long n, double arr[], int ind[]); + +#define ZERO RCONST(0.0) +#define HALF RCONST(0.5) +#define ONE RCONST(1.0) +#define TWO RCONST(2.0) +#define THREE RCONST(3.0) +#define TEN RCONST(10.0) +#define BUFSIZE 1000 +#define EPSILON RCONST(0.1) + +/* In order to keep begin the index numbers from 1 instead of 0, we define + * macros here. Also, we define macros to ease referencing various variables in + * the sundials nvector. + */ + +#define Y0(k) NV_Ith_S(data->Y0,k-1) +#define x(i) NV_Ith_S(data->x,i-1) + +#define t(i) NV_Ith_S(data->tArray,i) +#define P(i) NV_Ith_S(data->PArray,i) +#define dPdt(i) NV_Ith_S(data->dPdtArray,i) + +#define T(i) NV_Ith_S(y,((i-1)*data->nvar)+data->nt) +#define Y(i,k) NV_Ith_S(y,((i-1)*data->nvar)+data->ny+k-1) + +#define Tdot(i) NV_Ith_S(ydot,((i-1)*data->nvar)+data->nt) +#define Ydot(i,k) NV_Ith_S(ydot,((i-1)*data->nvar)+data->ny+k-1) + +#define Pp(i) NV_Ith_S(data->Pp,i-1) +#define dPdtp(i) NV_Ith_S(data->dPdtp,i-1) + +#define Tres(i) NV_Ith_S(res,((i-1)*(data->nvar+1))+data->nt) +#define Yres(i,k) NV_Ith_S(res,((i-1)*(data->nvar+1))+data->ny+k-1) +#define taures(i) NV_Ith_S(res,((i-1)*(data->nvar+1))+data->ntau) + +#define wdot(i) wdot[i-1] +#define enthalpy(i) enthalpy[i-1] + + +typedef struct { + /* This structure contains all information relevant to evaluating the + * residue. + */ + Cantera::IdealGasMix* gas;//Ideal gas object containing thermodynamic and + //kinetic data + + realtype P; //Pressure (in Pascals) + realtype T0; //Initial Temperature (in K) + N_Vector Y0; //Initial Mass Fractions + realtype rho0; //initial density + realtype TIgn; //Ignition Temperature (in K) + realtype tauIgn; //Ignition delay time (in s) + N_Vector x; //The grid + N_Vector tArray; + N_Vector PArray; + N_Vector dPdtArray; + N_Vector Pp; + N_Vector dPdtp; + realtype dPdt; + int tArraySize; + bool constantVolume; //boolean to enable/disable constant volume + //version of the ignition delay problem + bool imposedPdt; //boolean to enable/disable manual entry for P and dPdt + int jguess; //saved index to accelerate lookup via "hunt" + bool writeSpeciesSensitivities; //boolean to enable/disable writing + //of species sensitivities + bool IVPSuccess; + int nsp,neq,npts,nvar; //key quantities required in for loops: + //nsp-> no: of chemical species + //neq-> number of governing equations + //npts-> no: of grid points + //nvar->no: of variables + int nt,ny,ntau; //array indices for various variables: + //nt-> temperature (0) + //ny-> species mass fractions + //ntau-> ignition delay time + int KSneq; //No: of equations in the BVP formulation of + //the ignition delay problem + realtype rel; //the relative perturbation factor used in + //computing difference quotients + realtype atol; + realtype rtol; + + realtype t1; //Maximum time to run simulation + int nreac; + int pert_index; //index of reaction whose rate is perturbed + bool sensitivityAnalysis;// boolean to activate perturbations for + //sensitivity analysis + bool sens; //if true, perform sensitivity analysis + bool secondOrder; //if true, use second order finite differences + FILE *output; //output file for temperature and species profiles + FILE *ignSensOutput; //output file for ignition sensitivities + FILE *speSensOutput; //output file for species sensitivities + FILE *dPdtInput; //output file for species sensitivities + + gsl_interp_accel* acc; + gsl_interp_accel* accdot; + gsl_spline* spline; + gsl_spline* splinedot; +} *UserData; + + +/* Set the initial conditions using this subroutine: */ +static int setInitialProfile(UserData data, N_Vector y); + +/* Evaluate the residue for the IVP using this subroutine: */ +static int residue(realtype t, N_Vector y, N_Vector ydot, void *user_data); + +static int parseInput(UserData data, int argc, char *argv[]); + +static int hunt(realtype x, realtype xx[], int n, int jguess); + +static void polint(realtype *xdata, realtype *f, int n, realtype x,realtype *y, realtype *dy); + +static void lookupDpdt(UserData data, realtype t, realtype *P, realtype *dPdt); + +/* Subroutine that prints the column titles in the output file: */ +static void printInstructions(); + +/* Subroutine that prints the column titles in the output file: */ +static void printHeader(UserData data); +static void printSensitivitiesHeader(UserData data); + +/* Subroutine that prints the output of the IVP into the output file contained + * in the UserData structure: */ +static void printOutput(realtype t, N_Vector y, UserData data); + +/* Subroutine that prints the sensitivities into the output file (ignSensOutput) + * contained in the UserData structure: */ +static void printIgnSensitivities(realtype sensCoeffs[], int indices[], UserData data); + +static int parsedPdt(UserData data); + +/* Print solver statistics for the IVP: */ +static void PrintFinalStats(void *cvode_mem); + +/* Subroutine that reports failure if memory hasn't been successfully allocated + * to various objects: */ +static int check_flag(void *flagvalue, const char *funcname, int opt); + +int main(int argc, char *argv[]) +{ + clock_t start, end; + double cpu_time_used; + start = clock(); + + int ier; //error flag + + UserData data; //User defined data + data = NULL; + + /* Create and load problem data block. */ + data = (UserData) malloc(sizeof *data); + ier=parseInput(data,argc,argv); + if(ier==-1)return(-1); + + { + /* Solver Memory:*/ + void *mem; + mem=NULL; + + /* Save the initial mass fractions: */ + data->Y0=N_VNew_Serial(data->nsp); + for(int k=1;k<=data->nsp;k++){ + Y0(k)=data->gas->massFraction(k-1); + } + + /*Assign variable indices:*/ + data->nt=0; + data->ny=data->nt+1; + data->ntau=data->ny+data->nsp; + + /*Get and save the no: of reactions:*/ + data->nreac=data->gas->nReactions(); + + /* Solution vector of the IVP: */ + N_Vector y; + y = NULL; + y = N_VNew_Serial(data->neq); + ier=check_flag((void *)y, "N_VNew_Serial", 0); + + realtype t0, tret; //times + t0=tret=ZERO; + + /*Set the initial values:*/ + setInitialProfile(data, y); + /*Print out the column names and the initial values:*/ + printHeader(data); + printOutput(0.0e0,y,data); + + /* Create a CVode solver object for solution of the IVP: */ + mem = CVodeCreate(CV_BDF,CV_NEWTON); + ier=check_flag((void *)mem, "CVodeCreate", 0); + + /*Associate the user data with the solver object: */ + ier = CVodeSetUserData(mem, data); + ier=check_flag(&ier, "CVodeSetUserData", 1); + + /* Initialize the solver object by connecting it to the solution vector + * y: */ + ier = CVodeInit(mem, residue, t0, y); + ier=check_flag(&ier, "CVodeInit", 1); + + /*Set the tolerances: */ + ier = CVodeSStolerances(mem, data->rtol, data->atol); + ier=check_flag(&ier, "IDASStolerances", 1); + + /* Create dense SUNMatrix for use in linear solves: */ + SUNMatrix A; + A = SUNDenseMatrix(data->neq, data->neq); + ier=check_flag((void *)A, "SUNDenseMatrix", 0); + + /* Create dense SUNLinearSolver object for use by CVode: */ + SUNLinearSolver LS; + LS = SUNDenseLinearSolver(y, A); + ier=check_flag((void *)LS, "SUNDenseLinearSolver", 0); + + /* Call CVDlsSetLinearSolver to attach the matrix and linear solver to + * CVode: */ + ier = CVDlsSetLinearSolver(mem, LS, A); + ier=check_flag(&ier, "CVDlsSetLinearSolver", 1); + + /* Use a difference quotient Jacobian: */ + ier = CVDlsSetJacFn(mem, NULL); + ier=check_flag(&ier, "CVDlsSetJacFn", 1); + + /*Begin IVP solution; Save solutions in the temporary solution vector; + * stop problem when the temperature reaches a certain value (like 400 + * K) corresponding to the time for ignition :*/ + int i=1; + bool delayFound=false; + while (tret<=data->t1) { + if(T(1)>=data->TIgn && !delayFound){ + printf("Ignition Delay: %15.6es\n", tret); + data->tauIgn=tret; //Save the ignition delay time + delayFound=true; + } + ier = CVode(mem, data->t1, y, &tret, CV_ONE_STEP); + if(check_flag(&ier, "CVode", 1)) { + data->IVPSuccess=false; + break; + } + else{ + data->IVPSuccess=true; + } + printOutput(tret,y,data); + } + + if(data->IVPSuccess && data->sens){ + /*Enable sensitivity analysis:*/ + data->sensitivityAnalysis=true; + data->rel=EPSILON; + realtype oneOverRel=ONE/(data->rel); + + /*Create an array to store the logarithmic sensitivity + * coefficients:*/ + realtype sensCoeffs[data->nreac]; + + /*Create an array to store the indices of the reactions (needed + * for sorting the sensitivity coefficients):*/ + int indices[data->nreac]; + double tauIgnPerturbedForward=0.0e0; + double tauIgnPerturbedBackward=0.0e0; + + for(int j=0;jnreac;j++){ + + /*Save the index of the reaction whose collision + * frequency (A in the Arrhenius law k=A*exp(-Eₐ/RT)) + * is to be perturbed:*/ + data->pert_index=j; + + /*Perturb forward:*/ + data->rel=EPSILON; + /*Set the initial values:*/ + setInitialProfile(data, y); + + /* Initialize the solver object by connecting it to the solution vector + * y: */ + ier = CVodeReInit(mem, t0, y); + ier=check_flag(&ier, "CVodeReInit", 1); + + /*Rerun the ignition delay problem!*/ + printf("Solving forward perturbed problem %d:\n",j); + tret=0.0e0; + delayFound=false; + while (tret<=data->t1) { + if(T(1)>=data->TIgn && !delayFound){ + printf("\tIgnition Delay(forward): %15.6es\n", tret); + tauIgnPerturbedForward=tret; //Save the ignition delay time + delayFound=true; + /*No point continuing solution*/ + break; + } + ier = CVode(mem, data->t1, y, &tret, CV_ONE_STEP); + if(check_flag(&ier, "CVode", 1)) { + data->IVPSuccess=false; + break; + } + else{ + data->IVPSuccess=true; + } + //printOutput(tret,y,data); + } + if(!delayFound)tauIgnPerturbedForward=1e+99; + + if(data->secondOrder){ + /*Perturb backward:*/ + data->rel=-EPSILON; + /*Set the initial values:*/ + setInitialProfile(data, y); + + /* Initialize the solver object by connecting it to the solution vector + * y: */ + ier = CVodeReInit(mem, t0, y); + ier=check_flag(&ier, "CVodeReInit", 1); + + /*Rerun the ignition delay problem!*/ + printf("Solving backward perturbed problem %d:\n",j); + tret=0.0e0; + delayFound=false; + while (tret<=data->t1) { + if(T(1)>=data->TIgn && !delayFound){ + printf("\tIgnition Delay(backward): %15.6es\n", tret); + tauIgnPerturbedBackward=tret; //Save the ignition delay time + delayFound=true; + /*No point continuing solution*/ + break; + } + ier = CVode(mem, data->t1, y, &tret, CV_ONE_STEP); + if(check_flag(&ier, "CVode", 1)) { + data->IVPSuccess=false; + break; + } + else{ + data->IVPSuccess=true; + } + //printOutput(tret,y,data); + } + if(!delayFound)tauIgnPerturbedBackward=1e+99; + + /*Take the finite difference quotient as an + * approximation to the logarithmic sensitivity + * coefficient:*/ + sensCoeffs[j]=oneOverRel*(tauIgnPerturbedForward-tauIgnPerturbedBackward)/(2.0e0*data->tauIgn); + }else{ + sensCoeffs[j]=oneOverRel*(tauIgnPerturbedForward-data->tauIgn)/(data->tauIgn); + } + + indices[j]=j; + printf("\n"); + } + + /*Sort the sensitivities in ascending order. Note the advancing + * of the beginning indices of the sensCoeffs and indices + * arrays. This is due to the Numerical recipes convention for + * array indexing used in sort subroutine:*/ + //sort(data->nreac,sensCoeffs-1,indices-1); + + /*Print out the sensitivities:*/ + printIgnSensitivities(sensCoeffs,indices,data); + } + + /* Print remaining counters and free memory. */ + PrintFinalStats(mem); + CVodeFree(&mem); + N_VDestroy_Serial(y); + } + + /*Free memory and delete all the vectors and user data:*/ + if(data->Y0!=NULL){ + N_VDestroy_Serial(data->Y0); + printf("Initial Mass fraction Vector Deleted!\n"); + } + if(data->x!=NULL){ + N_VDestroy_Serial(data->x); + printf("Grid for BVP deleted!\n"); + } + if(data->Pp!=NULL){ + N_VDestroy_Serial(data->Pp); + printf("P array deleted!\n"); + } + if(data->dPdtp!=NULL){ + N_VDestroy_Serial(data->dPdtp); + printf("dPdt array deleted!\n"); + } + if(data->gas!=NULL){ + delete data->gas; + printf("Gas Deleted!\n"); + + } + if(data->imposedPdt){ + if(data->acc!=NULL){ + gsl_interp_accel_free(data->acc); + } + if(data->accdot!=NULL){ + gsl_interp_accel_free(data->accdot); + } + if(data->spline!=NULL){ + gsl_spline_free(data->spline); + } + if(data->splinedot!=NULL){ + gsl_spline_free(data->splinedot); + } + } + if(data!=NULL){ + /* Free the user data */ + free(data); + printf("User data structure deleted!\n"); + } + + end = clock(); + cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC; + printf("Elapsed cpu time: %15.6e s\n", cpu_time_used); + return(0); +} + +static int parseInput(UserData data, int argc, char *argv[]){ + + /*Set defaults here:*/ + /*****************************************************/ + /*Relative tolerance for IVP:*/ + data->rtol = RCONST(1e-14); + /*Absolute tolerance for IVP:*/ + data->atol = RCONST(1e-14); + /*Final Time:*/ + data->t1 = RCONST(10.0); + /*Solve constant pressure problem:*/ + data->constantVolume=false; + /*Do not impose a P vs t curve:*/ + data->imposedPdt=false; + /*Disable writing species sensitivities:*/ + data->writeSpeciesSensitivities=false; + /*Disable sensitivity analysis:*/ + data->sens=false; + /*Index to start with for pressure lookup (see hunt):*/ + data->jguess=0; + /*Set rate of change of pressure to zero:*/ + data->dPdt=ZERO; + /*Disable sensitivity analysis:*/ + data->sensitivityAnalysis=false; + /*Use first order finite differences:*/ + data->secondOrder=false; + /*Find the relative perturbation constant:*/ + //data->rel=SUNRsqrt(UNIT_ROUNDOFF); + data->rel=0.1; + /*Set flags that indicate success of various stages:*/ + data->IVPSuccess=false; + /*****************************************************/ + + int ier; + int opt; + char mech[BUFSIZE+1]; + char comp[BUFSIZE+1]; + bool enteredT0, enteredP, enteredMech, enteredComp; + enteredT0=enteredP=enteredMech=enteredComp=false; + while((opt=getopt(argc,argv,"a:r:T:P:m:c:t:vsd2")) != -1){ + switch(opt){ + case 'a': + data->atol=RCONST(atof(optarg)); + break; + case 'r': + data->rtol=RCONST(atof(optarg)); + break; + case 'T': + data->T0=RCONST(atof(optarg)); + enteredT0=true; + break; + case 'P': + data->P=RCONST(atof(optarg))*Cantera::OneAtm; + enteredP=true; + break; + case 'm': + snprintf(mech,BUFSIZE,"%s",optarg); + enteredMech=true; + break; + case 'c': + snprintf(comp,BUFSIZE,"%s",optarg); + enteredComp=true; + break; + case 't': + data->t1=RCONST(atof(optarg)); + enteredT0=true; + break; + case 'v': + data->constantVolume=true; + break; + case 'S': + data->writeSpeciesSensitivities=true; + break; + case 's': + data->sens=true; + break; + case 'd': + data->imposedPdt=true; + break; + case '2': + data->secondOrder=true; + break; + default: + printInstructions(); + return(-1); + } + } + if(!enteredT0){ + printf("Not specified initial Temperature! Exiting...\n"); + printInstructions(); + return(-1); + } + else if(!enteredP){ + printf("Not specified Pressure! Exiting...\n"); + printInstructions(); + return(-1); + } + else if(!enteredMech){ + printf("Not specified Mechanism! Exiting...\n"); + printInstructions(); + return(-1); + } + else if(!enteredComp){ + printf("Not specified Composition! Exiting...\n"); + printInstructions(); + return(-1); + } + else{ + if(data->imposedPdt){ + ier=parsedPdt(data); + if(ier==-1){ + return(-1); + } + } + + printf("Required inputs provided:\n"); + printf("\natol: %15.6e\n", data->atol); + printf("\nrtol: %15.6e\n", data->rtol); + printf("\nT0 : %15.6e K\n", data->T0); + printf("\nP : %15.6e Pa\n", data->P); + printf("\nmech: %s\n", mech); + printf("\ncomp: %s\n", comp); + printf("\nt1 : %15.6e s\n", data->t1); + printf("\nconstantVolume : %d\n", data->constantVolume); + printf("\nimposedPdt : %d\n", data->imposedPdt); + printf("\nwriteSpeciesSensitivities : %d\n", + data->writeSpeciesSensitivities); + printf("Proceeding...\n\n"); + } + + /*Define Gas here:*/ + data->gas = new Cantera::IdealGasMix(mech); + + /* Set the initial state of the gas: */ + data->gas->setState_TPX(data->T0, + data->P, + comp); + + /* Create output file for the solution of the IVP: */ + data->output=fopen("output.dat","w"); + + /* Create output file for the ignition sensitivities (post BVP): */ + data->ignSensOutput=fopen("ignitionSensitivities.dat","w"); + + /* Create output file for the species sensitivities (post BVP): */ + data->speSensOutput=fopen("speciesSensitivities.dat","w"); + + data->rho0=data->gas->density(); + if(data->constantVolume){ + data->gas->equilibrate("UV"); + } + else{ + data->gas->equilibrate("HP"); + } + data->TIgn=data->T0+(data->gas->temperature() + -data->T0)*RCONST(0.20); + printf("Ignition Temperature: %15.6e K\n",data->TIgn); + data->gas->setState_TPX(data->T0, + data->P, + comp); + + /*Get and save the no: of species:*/ + data->nsp=data->gas->nSpecies(); + /*set the no: of variables for the IVP:*/ + data->nvar=data->nsp+1; + /*set the no: of equations for the IVP:*/ + data->neq=data->nvar; + + return(0); +} + +static int setInitialProfile(UserData data, N_Vector y) +{ + /*This routine sets the initial temperature and mass fractions for the + * solution of the initial value problem:*/ + N_VConst(ZERO, y); + + T(1)=data->T0; + for (int k = 1; k <=data->nsp; k++) { + Y(1,k)=Y0(k); + } + + return(0); + +} + +static int residue(realtype t, N_Vector y, N_Vector ydot, void *user_data) +{ + + /*This is the sub-routine that computes the right hand side F(y) in the + * problem ∂y/∂t = F(y). + * + * The energy conservation equation is + * ∂T/∂t = (1/ρcₚ)(dP/dt -∑ ώᵢhᵢ) + * + * T has units of K + * hᵢ has units J/kmol, so we must multiply the enthalpy + * defined above (getEnthalpy_RT) by T (K) and the gas constant + * (J/kmol/K) to get the right units. + * ρ has units of kg/m³ + * cₚ has units J/kg/K + * ώᵢ has units of kmol/m³/s + * dP/dt has units of Pa/s + * + * In the case of a constant volume problem, the energy conservation + * equation can be re-written as + * ∂T/∂t = -(1/ρcᵥ)∑ ώᵢeᵢ + * eᵢ has units J/kmol + * cᵥ has units J/kg/K + * + * + * The species conservation equation is + * ∂Yᵢ/∂t = ώᵢWᵢ/ρ + * Wᵢ has units of kg/kmol + */ + + /*Assign data structure type to user_data:*/ + UserData data; + data=(UserData)user_data; + + /*Get no: of species:*/ + int nsp=data->nsp; + + /*Create temporary arrays to store enthalpy (non-dimensional) and ώᵢ:*/ + realtype wdot[nsp]; + realtype enthalpy[nsp]; + + try { + /*If sensitivity analysis has been enabled, perturb the collision + * frequency by using cantera's setMultiplier function:*/ + if(data->sensitivityAnalysis){ + data->gas->setMultiplier(data->pert_index,ONE+data->rel); + //printf("%d\t%15.9e\n",data->pert_index,ONE+data->rel); + } + + /*Set the gas conditions:*/ + if(data->constantVolume){ + data->gas->setMassFractions_NoNorm(&Y(1,1)); + data->gas->setTemperature(T(1)); + data->gas->setDensity(data->rho0); + data->P=data->gas->pressure(); + } + else if(data->imposedPdt){ + data->gas->setMassFractions_NoNorm(&Y(1,1)); + data->gas->setTemperature(T(1)); + realtype P,dPdt; + P=dPdt=0.0e0; + lookupDpdt(data, t, &P, &dPdt); + data->P=P*Cantera::OneAtm; + data->dPdt=dPdt*Cantera::OneAtm; + data->gas->setPressure(data->P); + } + else{ + data->gas->setMassFractions_NoNorm(&Y(1,1)); + data->gas->setTemperature(T(1)); + data->gas->setPressure(data->P); + } + + realtype rho=data->gas->density(); //get ρ + realtype cp=data->gas->cp_mass(); //get cₚ + data->gas->getNetProductionRates(wdot); //get ώᵢ + data->gas->getEnthalpy_RT(enthalpy); //get hᵢ/RT + + if(data->constantVolume){ + for(int k=1;k<=nsp;k++){ + enthalpy(k)=enthalpy(k)-ONE; + } + cp=data->gas->cv_mass(); + } + + realtype sum=ZERO; + for(int k=1;k<=nsp;k++){ + Ydot(1,k)=wdot(k)*(data->gas->molecularWeight(k-1))/rho; + sum=sum+wdot(k)*enthalpy(k)*T(1); + } + sum=sum*Cantera::GasConstant; + Tdot(1)=(data->dPdt-sum)/(rho*cp); + + /*If sensitivity analysis has been enabled, *un*perturb the collision + * frequency:*/ + if(data->sensitivityAnalysis){ + data->gas->setMultiplier(data->pert_index,ONE); + } + + } catch (Cantera::CanteraError& err) { + printf("Error:\n"); + printf("%s\n",err.what()); + return(-1); + } + + return(0); +} + +static void lookupDpdt(UserData data, realtype t, realtype *P, realtype *dPdt) +{ + + realtype *tData, *PData, *dPdtData; + tData = N_VGetArrayPointer_Serial(data->tArray); + PData = N_VGetArrayPointer_Serial(data->PArray); + dPdtData = N_VGetArrayPointer_Serial(data->dPdtArray); + int jguess=data->jguess; + int k=0; + int safel=0; + int nOrder=4; + + realtype dy; + int npts=data->tArraySize; + + if(t<=data->t1){ + *P=gsl_spline_eval(data->spline,t,data->acc); + *dPdt=gsl_spline_eval(data->splinedot,t,data->accdot); + } + else{ + *P=P(npts-1); + *dPdt=dPdt(npts-1); + } +} + + +static int hunt(realtype x, realtype xx[], int n, int jguess) +{ + int jlo=jguess; + int jm,jhi,inc; + int ascnd; + + ascnd=(xx[n-1] >= xx[0]); + if (jlo <= 0 || jlo > n-1) { + jlo=0; + jhi=n; + } else { + inc=1; + if ((x >= xx[jlo]) == ascnd) { + if (jlo == n-1) return(jlo); + jhi=(jlo)+1; + while ((x >= xx[jhi]) == ascnd) { + jlo=jhi; + inc += inc; + jhi=(jlo)+inc; + if (jhi > n-1) { + jhi=n; + break; + } + } + } else { + if (jlo == 1) { + jlo=0; + return (jlo); + } + jhi=(jlo)--; + while ((x < xx[jlo]) == ascnd) { + jhi=(jlo); + inc <<= 1; + if (inc >= jhi) { + jlo=0; + break; + } + else jlo=jhi-inc; + } + } + } + while (jhi-(jlo) != 1) { + jm=(jhi+(jlo)) >> 1; + if ((x >= xx[jm]) == ascnd) + jlo=jm; + else + jhi=jm; + } + if (x == xx[n-1]) jlo=n-2; + if (x == xx[0]) jlo=1; + + return(jlo); +} + +static void polint(realtype *xdata, realtype *f, int n, realtype x, realtype *y, realtype *dy){ + int i,m,ns=1; + realtype den,dif,dift,ho,hp,w; + realtype c[n+1],d[n+1]; + dif=fabs(x-xdata[1]); + for (i=1;i<=n;i++) { + if ( (dift=fabs(x-xdata[i])) < dif) { + ns=i; + dif=dift; + } + c[i]=f[i]; + d[i]=f[i]; + } + *y=f[ns--]; + for (m=1;m=l;i--) { + if (arr[i] <= a) break; + arr[i+1]=arr[i]; + ind[i+1]=ind[i]; + } + arr[i+1]=a; + ind[i+1]=b; + } + if (jstack == 0) break; + ir=istack[jstack--]; + l=istack[jstack--]; + } else { + k=(l+ir) >> 1; + SWAP(arr[k],arr[l+1]) + SWAP(ind[k],ind[l+1]) + if (arr[l] > arr[ir]) { + SWAP(arr[l],arr[ir]) + SWAP(ind[l],ind[ir]) + } + if (arr[l+1] > arr[ir]) { + SWAP(arr[l+1],arr[ir]) + SWAP(ind[l+1],ind[ir]) + } + if (arr[l] > arr[l+1]) { + SWAP(arr[l],arr[l+1]) + SWAP(ind[l],ind[l+1]) + } + i=l+1; + j=ir; + a=arr[l+1]; + b=ind[l+1]; + for (;;) { + do i++; while (arr[i] < a); + do j--; while (arr[j] > a); + if (j < i) break; + SWAP(arr[i],arr[j]); + SWAP(ind[i],ind[j]); + } + arr[l+1]=arr[j]; + ind[l+1]=ind[j]; + arr[j]=a; + ind[j]=b; + jstack += 2; + if (jstack > NSTACK) printf("NSTACK too small in sort."); + if (ir-i+1 >= j-l) { + istack[jstack]=ir; + istack[jstack-1]=i; + ir=j-1; + } else { + istack[jstack]=j-1; + istack[jstack-1]=l; + l=i; + } + } + } + free_lvector(istack,1,NSTACK); +} +#undef M +#undef NSTACK +#undef SWAP + +static void printInstructions(){ + printf("\nInputs either incomplete or wrong!\n"); + printf("\n-a \n"); + printf("\n-r \n"); + printf("\n-T \n"); + printf("\n-P \n"); + printf("\n-m \n"); + printf("\n-c \n"); + printf("\n-t \n"); + printf("\n-v :Enables constant volume simulation\n"); + printf("\n-s :Enables ignition delay sensitivity output\n"); + printf("\n-S :Enables species mass fraction sensitivity output (works only with -s)\n"); + printf("\n-d :Enables manual dPdt entryspecies sensitivity output\n"); + printf("\n-2 :Enables 2nd order accurate finite differences\n"); + printf("\nexample: "); + printf("-T 1200.0 -P 1.0 -m gri30.cti"); + printf(" -c H2:1.0,N2:3.76,O2:1.0"); + printf(" -t 1e-03\n"); +} + +static void printHeader(UserData data) +{ + fprintf((data->output), "%15s\t","#1"); + for (int k = 1; k <=data->nvar; k++) { + fprintf((data->output), "%15d\t",k+1); + } + fprintf((data->output), "\n"); + + fprintf((data->output), "%15s\t%15s\t%15s\t", + "#time(s)","Temp(K)","Pressure(Pa)"); + for (int k = 1; k <=data->nsp; k++) { + fprintf((data->output), "%15s\t", + data->gas->speciesName(k-1).c_str()); + } + fprintf((data->output), "\n"); + +} + +static void printSensitivitiesHeader(UserData data) +{ + fprintf((data->speSensOutput), "%15s\t","#1"); + for (int k = 1; k <=data->nvar; k++) { + fprintf((data->speSensOutput), "%15d\t",k+1); + } + fprintf((data->speSensOutput), "\n"); + + fprintf((data->speSensOutput), "%15s\t%15s\t", + "#time(s)","Temp(K)"); + for (int k = 1; k <=data->nsp; k++) { + fprintf((data->speSensOutput), "%15s\t", + data->gas->speciesName(k-1).c_str()); + } + fprintf((data->speSensOutput), "\n"); + +} + +static void printOutput(realtype t, N_Vector y, UserData data) +{ + fprintf((data->output), "%15.6e\t%15.6e\t%15.6e\t",t,T(1),data->P); + for (int k = 1; k <=data->nsp; k++) { + fprintf((data->output), "%15.6e\t",Y(1,k)); + } + fprintf((data->output), "\n"); +} + +static void printIgnSensitivities(realtype sensCoeffs[], int indices[], UserData data) +{ + char buf[100]; + std::string line; + for(int j=0;jnreac;j++){ + line=data->gas->reactionString(indices[j]); + sprintf(buf,"%s",line.c_str()); + fprintf((data->ignSensOutput), "%d\t\"%35s\"\t%15.6e\n",indices[j],buf,sensCoeffs[j]); + } + //for(int j=0;jnreac;j++){ + // fprintf((data->ignSensOutput), "%d\t%15.6e\n",indices[j],sensCoeffs[j]); + //} +} + +static int parsedPdt(UserData data){ + /* Open file containing P and dPdt as functions of time: */ + data->dPdtInput=fopen("dPdt.dat","r"); + if(data->dPdtInput==NULL){ + printf("The dPdt.dat file wasn't found!\n"); + return(-1); + } + char buf[1000]; + char comment[1]; + char *ret; + int i=0; + int j=0; + while (fgets(buf,100, data->dPdtInput)!=NULL){ + comment[0]=buf[0]; + if(strncmp(comment,"#",1)!=0 && + strncmp(comment,"\n",1)!=0){ + i++; + } + } + printf("No: of rows with numbers: %d\n",i); + rewind(data->dPdtInput); + int nPts=i; + data->tArraySize=nPts; + i=0; + + data->tArray = N_VNew_Serial(nPts); + data->PArray = N_VNew_Serial(nPts); + data->dPdtArray = N_VNew_Serial(nPts); + + while (fgets(buf,100, data->dPdtInput)!=NULL){ + comment[0]=buf[0]; + if(strncmp(comment,"#",1)==0 || + strncmp(comment,"\n",1)==0){ + printf("Comment! Skip Line!\n"); + } + else{ + j=0; + ret=strtok(buf,", \t"); + t(i)=atof(ret); + j++; + while(ret!=NULL){ + ret=strtok(NULL,", \t"); + if(j==1){ + P(i)=atof(ret); + } + else if(j==2){ + dPdt(i)=atof(ret); + } + j++; + } + i++; + } + } + fclose(data->dPdtInput); + //for(int k=0;kt1=t(nPts-1); + data->P=P(0)*Cantera::OneAtm; + + /*check the polynomial interpolation (testing only)*/ + //realtype t,P,dPdt; + //t=1e-03; + //P=dPdt=0.0e0; + //lookupDpdt(data, t, &P, &dPdt); + //printf("Interpolated P value %15.6e. \n",P); + //printf("Interpolated dPdt value %15.6e. \n",dPdt); + // + + //GSL additions here: + data->acc=gsl_interp_accel_alloc(); + data->spline=gsl_spline_alloc(gsl_interp_steffen,data->tArraySize); + + data->accdot=gsl_interp_accel_alloc(); + data->splinedot=gsl_spline_alloc(gsl_interp_steffen,data->tArraySize); + + double* Pdata; + double* dPdtdata; + double* tdata; + + tdata=N_VGetArrayPointer_Serial(data->tArray); + Pdata=N_VGetArrayPointer_Serial(data->PArray); + dPdtdata=N_VGetArrayPointer_Serial(data->dPdtArray); + + gsl_spline_init(data->spline,tdata,Pdata,data->tArraySize); + gsl_spline_init(data->splinedot,tdata,dPdtdata,data->tArraySize); + + return(0); +} + +static void PrintFinalStats(void *cvode_mem) +{ + long int nst, nfe, nsetups, nje, nfeLS, nni, ncfn, netf, nge; + int flag; + + flag = CVodeGetNumSteps(cvode_mem, &nst); + check_flag(&flag, "CVodeGetNumSteps", 1); + flag = CVodeGetNumRhsEvals(cvode_mem, &nfe); + check_flag(&flag, "CVodeGetNumRhsEvals", 1); + flag = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups); + check_flag(&flag, "CVodeGetNumLinSolvSetups", 1); + flag = CVodeGetNumErrTestFails(cvode_mem, &netf); + check_flag(&flag, "CVodeGetNumErrTestFails", 1); + flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni); + check_flag(&flag, "CVodeGetNumNonlinSolvIters", 1); + flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn); + check_flag(&flag, "CVodeGetNumNonlinSolvConvFails", 1); + + flag = CVDlsGetNumJacEvals(cvode_mem, &nje); + check_flag(&flag, "CVDlsGetNumJacEvals", 1); + flag = CVDlsGetNumRhsEvals(cvode_mem, &nfeLS); + check_flag(&flag, "CVDlsGetNumRhsEvals", 1); + + flag = CVodeGetNumGEvals(cvode_mem, &nge); + check_flag(&flag, "CVodeGetNumGEvals", 1); + + printf("\nFinal CVode Statistics:\n"); + printf("nst = %-6ld nfe = %-6ld nsetups = %-6ld nfeLS = %-6ld nje = %ld\n", + nst, nfe, nsetups, nfeLS, nje); + printf("nni = %-6ld ncfn = %-6ld netf = %-6ld nge = %ld\n \n", + nni, ncfn, netf, nge); +} + +static int check_flag(void *flagvalue, const char *funcname, int opt) +{ + int *errflag; + + /* Check if SUNDIALS function returned NULL pointer - no memory allocated */ + if (opt == 0 && flagvalue == NULL) { + fprintf(stderr, + "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n", + funcname); + return(1); + } + + /* Check if flag < 0 */ + else if (opt == 1) { + errflag = (int *) flagvalue; + if (*errflag < 0) { + fprintf(stderr, + "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n", + funcname, *errflag); + return(1); + } + } + + /* Check if function returned NULL pointer - no memory allocated */ + else if (opt == 2 && flagvalue == NULL) { + fprintf(stderr, + "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n", + funcname); + return(1); + } + + return(0); +}