added sensBrute.cpp, README, and examples

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 <absolute tolerance> 
+	-r <relative tolerance> 
+	-f <BVP residue tolerance>
+	-T <initial temperature in K>
+	-P <initial pressure in atm>
+	-m <mechanism file (cti or xml)>
+	-c <mole fraction composition>
+	-t <integration time>
+	-s <enable sensitivity analysis>
+
+Execute sensBVP -h for all available options.
+
+Try the cases in the examples directory.

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

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

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))

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
+

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);