commit 48cca4770f31efa3c771e35d217aee00e20b2f14
Author: vyaas <mutantbex@gmail.com>
Date:   Sun May 10 22:40:41 2020 -0500

    first

diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..8232996
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,91 @@
+#compiler	=icpc
+#CANTERA_DIR	=/opt/scientific/cantera-2.3_intel_mkl
+#IDA_DIR		=/opt/scientific/sundials-3.1.1_intel_mkl
+#EXE		=lagrangianCombustion2
+#DESTDIR		=~/bin
+
+compiler	=g++
+CANTERA_DIR	=/opt/scientific/cantera-2.4_gnu_blas
+IDA_DIR		=/opt/scientific/sundials-3.1.1_gnu_lapack
+EXE		=lagrangianCombustion
+DESTDIR		=~/bin
+
+CANTERA_INCLUDES=-I$(CANTERA_DIR)/include
+GSL_INCLUDES	=-I/usr/include/gsl
+GSL_LIBS 	=-L/usr/lib -lgsl -lgslcblas
+IDA_INCLUDES	=-I$(IDA_DIR)/include
+SUN_INCLUDES	=-I$(IDA_DIR)/include
+IDA_LIBS	=-L$(IDA_DIR)/lib -lsundials_nvecopenmp \
+		 -lsundials_ida -lsundials_sunlinsollapackband
+CANTERA_LIBS 	=-L$(CANTERA_DIR)/lib -lcantera_shared
+RPATH		=-Wl,-rpath=$(IDA_DIR)/lib,-rpath=$(CANTERA_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
+
+OBJS = parse.o UserData.o gridRoutines.o residue.o solution.o main.o 
+
+all: $(EXE)
+
+# pull in dependency info for *existing* .o files
+-include $(OBJS:.o=.d)
+
+parse.o: parse.cpp parse.h parse.hpp
+	$(CPP) $(CPPFLAGS) -c parse.cpp -o parse.o
+	$(CPP) -MM parse.cpp > parse.d
+
+solution.o: solution.cpp solution.h 
+	$(CPP) $(CPPFLAGS) $(CANTERA_INCLUDES) $(IDA_INCLUDES) \
+		-c solution.cpp -o solution.o
+	$(CPP) -MM $(CANTERA_INCLUDES) $(IDA_INCLUDES) \
+		solution.cpp > solution.d
+
+residue.o: residue.cpp residue.h UserData.h gridRoutines.h
+	$(CPP) $(CPPFLAGS) $(CANTERA_INCLUDES) \
+		$(IDA_INCLUDES) $(GSL_INCLUDES) -c residue.cpp \
+		-o residue.o -fopenmp
+	$(CPP) -MM $(CANTERA_INCLUDES) \
+		$(IDA_INCLUDES) $(GSL_INCLUDES) \
+		residue.cpp > residue.d
+
+gridRoutines.o: gridRoutines.cpp gridRoutines.h parse.h
+	$(CPP) $(CPPFLAGS) -c gridRoutines.cpp -o gridRoutines.o
+	$(CPP) -MM gridRoutines.cpp > gridRoutines.d
+
+UserData.o: UserData.cpp UserData.h gridRoutines.h parse.h 
+	$(CPP) $(CPPFLAGS) $(CANTERA_INCLUDES) \
+		-c  UserData.cpp -o UserData.o
+	$(CPP) -MM $(CANTERA_INCLUDES) \
+		UserData.cpp > UserData.d
+
+main.o: main.cpp UserData.h solution.h residue.h macros.h
+	$(CPP) $(CPPFLAGS) $(CANTERA_INCLUDES) $(IDA_INCLUDES) \
+		-c  main.cpp -o main.o
+	$(CPP) -MM $(CANTERA_INCLUDES) $(IDA_INCLUDES) \
+		main.cpp > main.d
+
+$(EXE): $(OBJS)
+	$(CPP) $(CPPFLAGS) $(IDA_INCLUDES) \
+		$(CANTERA_INCLUDES) $(GSL_INCLUDES) $(RPATH) \
+		$(OBJS) -o $(EXE) -fopenmp $(IDA_LIBS)  \
+		$(CANTERA_LIBS) $(GSL_LIBS)
+
+.PHONY: install
+install:
+	cp $(EXE) $(DESTDIR)/$(EXE) 
+
+.PHONY: clean
+clean:
+	rm -f $(EXE) $(OBJS) *.d
+
diff --git a/UserData.cpp b/UserData.cpp
new file mode 100644
index 0000000..de04232
--- /dev/null
+++ b/UserData.cpp
@@ -0,0 +1,409 @@
+#include "UserData.h"
+#include "parse.h"
+
+void freeUserData(UserData data){
+	if(data!=NULL){
+		if(data->trmix!=NULL){
+			delete data->trmix;
+			printf("Transport Deleted!\n");
+
+		}
+		if(data->gas!=NULL){
+			delete data->gas;
+			printf("Gas Deleted!\n");
+
+		}
+		if(data->adaptiveGrid){
+			if(data->grid->xOld!=NULL){
+				delete[] data->grid->xOld;
+				printf("old grid array Deleted!\n");
+			}
+			if(data->grid->x!=NULL){
+				delete[] data->grid->x;
+				printf("current grid array Deleted!\n");
+			}
+			if(data->grid!=NULL){
+				free(data->grid);
+				printf("grid object Freed!\n");
+			}
+		}
+		else{
+			if(data->uniformGrid!=NULL){
+				delete[] data->uniformGrid;
+				printf("uniformGrid deleted!\n");
+			}
+		}
+		if(data->innerMassFractions!=NULL){
+			delete[] data->innerMassFractions;
+			printf("innerMassFractions array Deleted!\n");
+		}
+		if(data->output!=NULL){
+  			fclose(data->output);
+			printf("Output File Cleared from Memory!\n");
+		}
+		if(data->gridOutput!=NULL){
+  			fclose(data->gridOutput);
+			printf("Grid Output File Cleared from Memory!\n");
+		}
+		//if(data->ratesOutput!=NULL){
+  		//	fclose(data->ratesOutput);
+		//	printf("Rates Output File Cleared from Memory!\n");
+		//}
+		if(data->globalOutput!=NULL){
+  			fclose(data->globalOutput);
+			printf("Global Output File Cleared from Memory!\n");
+		}
+
+	}
+  	free(data);             /* Free the user data */
+	printf("\n\n");
+}
+
+UserData allocateUserData(FILE *input){
+
+	UserData data;
+  	data = (UserData) malloc(sizeof *data);
+
+	if(!data){
+		printf("Allocation Failed!\n");
+		return(NULL);
+	}
+	setSaneDefaults(data);
+
+	int ier;
+
+	ier=parseNumber<size_t>(input, "basePts"	, MAXBUFLEN, &data->npts);
+	if(ier==-1 || data->npts<=0){
+		printf("Enter non-zero basePts!\n");
+		return(NULL);
+	}
+
+	ier=parseNumber<double>(input, "domainLength", MAXBUFLEN, &data->domainLength);
+	if(ier==-1 || data->domainLength<=0.0e0){
+		printf("domainLength error!\n");
+		return(NULL);
+	}
+
+	ier=parseNumber<int>(input, "constantPressure"	, MAXBUFLEN, &data->constantPressure);
+	if(ier==-1 || (data->constantPressure!=0 && data->constantPressure!=1)){
+		printf("constantPressure error!\n");
+		return(NULL);
+	}
+
+	ier=parseNumber<int>(input, "problemType"	, MAXBUFLEN, &data->problemType);
+	if(ier==-1 || (data->problemType!=0 
+		       && data->problemType!=1
+			&& data->problemType!=2
+			&& data->problemType!=3)){
+		printf("include valid problemType!\n");
+		printf("0: premixed combustion with NO equilibrated ignition kernel\n");
+		printf("1: premixed combustion WITH equilibrated ignition kernel\n");
+		printf("2: arbitrary initial conditions\n");
+		printf("3: Restart\n");
+		return(NULL);
+	}
+
+	ier=parseNumber<double>(input, "dPdt"		, MAXBUFLEN, &data->dPdt);
+
+	ier=parseNumber<int>   (input, "reflectProblem"	, MAXBUFLEN, &data->reflectProblem);
+	if(data->reflectProblem!=0 && data->reflectProblem!=1){
+		printf("Invalid entry for reflectProblem! Can be only 1 or 0.\n");
+		return(NULL);
+	}
+	
+	ier=parseNumber<double>(input, "mdot"		, MAXBUFLEN, &data->mdot);
+	
+	ier=parseNumber<double>(input, "initialTemperature", MAXBUFLEN,
+			&data->initialTemperature);
+	if(ier==-1 || data->initialTemperature<=0.0e0){
+		printf("Enter positive initialTemperature in K!\n");
+		return(NULL);
+	}
+
+
+	ier=parseNumber<double>(input, "initialPressure", MAXBUFLEN,
+			&data->initialPressure);
+	if(ier==-1 || data->initialTemperature<=0.0e0){
+		printf("Enter positive initialPressure in atm!\n");
+		return(NULL);
+	}
+
+	ier=parseNumber<int>   (input, "metric"	, MAXBUFLEN, &data->metric);
+	if(data->metric!=0 && data->metric!=1 && data->metric!=2){
+		printf("Invalid entry for metric!\n");
+		printf("0: Cartesian\n");
+		printf("1: Cylindrical\n");
+		printf("2: Spherical\n");
+		return(NULL);
+	}
+
+	ier=parseNumber<double>(input, "QDot", MAXBUFLEN, &data->maxQDot);
+	if(ier==-1 && data->problemType==0){
+		printf("Need to specify QDot for problemType 0!\n");
+		return(NULL);
+	}
+
+	ier=parseNumber<double>(input, "kernelSize", MAXBUFLEN, &data->kernelSize);
+	if(ier==-1 && data->problemType==0){
+		printf("Need to specify kernelSize for problemType 0!\n");
+		return(NULL);
+	}
+
+	ier=parseNumber<double>(input, "ignTime", MAXBUFLEN, &data->ignTime);
+	if(ier==-1 && data->problemType==0){
+		printf("Need to specify ignTime for problemType 0!\n");
+		return(NULL);
+	}
+
+	ier=parseNumber<double>(input, "mixingWidth", MAXBUFLEN,
+			&data->mixingWidth);
+	if(ier==-1){
+		printf("Need to specify mixingWidth!\n");
+		return(NULL);
+	}
+	
+	ier=parseNumber<double>(input, "shift", MAXBUFLEN, &data->shift);
+
+	ier=parseNumber<double>(input, "firstRadius", MAXBUFLEN, &data->firstRadius);
+	
+	ier=parseNumber<double>(input, "wallTemperature", MAXBUFLEN, &data->wallTemperature);
+
+	ier=parseNumber<int>   (input, "dirichletInner"	, MAXBUFLEN,
+			&data->dirichletInner);
+	if(data->dirichletInner!=0 && data->dirichletInner!=1){
+		printf("dirichletInner can either be 0 or 1!\n");
+		return(NULL);
+	}
+
+	ier=parseNumber<int>   (input, "dirichletOuter"	, MAXBUFLEN,
+			&data->dirichletOuter);
+	if(data->dirichletOuter!=0 && data->dirichletOuter!=1){
+		printf("dirichletOuter can either be 0 or 1!\n");
+		return(NULL);
+	}
+
+	ier=parseNumber<int>   (input, "adaptiveGrid"	, MAXBUFLEN,
+			&data->adaptiveGrid);
+	if(ier==-1 || (data->adaptiveGrid!=0 && data->adaptiveGrid!=1)){
+		printf("specify adaptiveGrid as 0 or 1!\n");
+		return(NULL);
+	}
+
+	ier=parseNumber<int>   (input, "moveGrid"	, MAXBUFLEN,
+			&data->moveGrid);
+	if(ier==-1 || (data->moveGrid!=0 && data->moveGrid!=1)){
+		printf("specify moveGrid as 0 or 1!\n");
+		return(NULL);
+	}
+
+	ier=parseNumber<double>   (input, "isotherm"	, MAXBUFLEN,
+			&data->isotherm);
+	if(ier==-1){
+		printf("specify temperature of isotherm!\n");
+		return(NULL);
+	}
+
+	ier=parseNumber<double>(input, "gridOffset", MAXBUFLEN, &data->gridOffset);
+
+	ier=parseNumber<int>   (input, "nSaves"	, MAXBUFLEN, &data->nSaves);
+	if(data->nSaves<0 ){
+		printf("nSaves must be greater than 0!\n");
+		return(NULL);
+	}
+	
+	ier=parseNumber<int>   (input, "writeRates"	, MAXBUFLEN,
+			&data->writeRates);
+	if(data->writeRates!=0 && data->writeRates!=1){
+		printf("writeRates must either be 0 or 1!\n");
+		return(NULL);
+	}
+
+	ier=parseNumber<int>   (input, "writeEveryRegrid", MAXBUFLEN,
+			&data->writeEveryRegrid);
+	if(data->writeEveryRegrid!=0 && data->writeEveryRegrid!=1){
+		printf("writeEveryRegrid must either be 0 or 1!\n");
+		return(NULL);
+	}
+
+	ier=parseNumber<double>   (input, "writeDeltaT", MAXBUFLEN,
+			&data->writeDeltaT);
+
+	ier=parseNumber<int>   (input, "setConstraints"	, MAXBUFLEN,
+			&data->setConstraints);
+	if(data->setConstraints!=0 && data->setConstraints!=1){
+		printf("setConstraints must either be 0 or 1!\n");
+		return(NULL);
+	}
+
+	ier=parseNumber<int>   (input, "suppressAlg"	, MAXBUFLEN,
+			&data->suppressAlg);
+	if(data->setConstraints!=0 && data->suppressAlg!=1){
+		printf("suppressAlg must either be 0 or 1!\n");
+		return(NULL);
+	}
+
+	ier=parseNumber<int>   (input, "dryRun"	, MAXBUFLEN,
+			&data->dryRun);
+	if(data->dryRun!=0 && data->dryRun!=1){
+		printf("dryRun must either be 0 or 1!\n");
+		return(NULL);
+	}
+
+	ier=parseNumber<double>   (input, "finalTime"	, MAXBUFLEN,
+			&data->finalTime);
+
+	ier=parseNumber<double>   (input, "relativeTolerance"	, MAXBUFLEN,
+			&data->relativeTolerance);
+	ier=parseNumber<double>   (input, "radiusTolerance"	, MAXBUFLEN,
+			&data->radiusTolerance);
+	ier=parseNumber<double>   (input, "temperatureTolerance", MAXBUFLEN,
+			&data->temperatureTolerance);
+	ier=parseNumber<double>   (input, "pressureTolerance", MAXBUFLEN,
+			&data->pressureTolerance);
+	ier=parseNumber<double>   (input, "massFractionTolerance", MAXBUFLEN,
+			&data->massFractionTolerance);
+	ier=parseNumber<double>   (input, "bathGasTolerance", MAXBUFLEN,
+			&data->bathGasTolerance);
+
+	char chem[MAXBUFLEN],mix[MAXBUFLEN],tran[MAXBUFLEN];
+
+	ier=parseNumber<char>(input, "chemistryFile"	, MAXBUFLEN, chem);
+	if(ier==-1){
+		printf("Enter chemistryFile!\n");
+		return(NULL);
+	}else{
+		try{
+			data->gas = new Cantera::IdealGasMix(chem);
+			data->nsp=data->gas->nSpecies(); //assign no: of species
+
+    		} catch (Cantera::CanteraError& err) {
+		    printf("Error:\n");
+		    printf("%s\n",err.what());
+		    return(NULL);
+    		}
+	}
+
+	ier=parseNumber<char>(input, "transportModel", MAXBUFLEN, tran);
+	if(ier==-1){
+		printf("Enter transportModel!\n");
+		return(NULL);
+	}else{
+		try{
+  			data->trmix = Cantera::newTransportMgr(tran, data->gas);
+		}catch (Cantera::CanteraError& err) {
+			printf("Error:\n");
+			printf("%s\n",err.what());
+			return(NULL);
+    		}
+	}
+
+	ier=parseNumber<char>(input, "mixtureComposition", MAXBUFLEN, mix);
+	if(ier==-1){
+		printf("Enter mixtureComposition!\n");
+		return(NULL);
+	}else{
+		if(data->gas!=NULL){
+			try{
+				data->gas->setState_TPX(data->initialTemperature,
+							data->initialPressure*Cantera::OneAtm,
+							mix);
+			}catch (Cantera::CanteraError& err) {
+		    		printf("Error:\n");
+		    		printf("%s\n",err.what());
+		    		return(NULL);
+    			}
+		}
+	}
+
+	ier=parseNumber<int>   (input, "nThreads"	, MAXBUFLEN, &data->nThreads);
+	if(data->nThreads<0 ){
+		printf("nThreads must be greater than 0!\n");
+		return(NULL);
+	}
+
+	data->nr=0;
+	//data->np=data->nr+1;
+	data->nt=data->nr+1;
+	data->ny=data->nt+1;
+	data->np=data->ny+data->nsp;
+				
+	data->nvar=data->nsp+3;		 //assign no: of variables (R,T,P,nsp species)
+
+	if(!data->adaptiveGrid){
+		data->uniformGrid = new double [data->npts];
+		data->neq=data->nvar*data->npts;
+	}
+	else{
+		data->grid=(UserGrid) malloc(sizeof *data->grid);
+		ier=getGridSettings(input,data->grid);
+		if(ier==-1)return(NULL);
+
+		ier=initializeGrid(data->grid);
+		if(ier==-1)return(NULL);
+
+		ier=reGrid(data->grid, data->grid->position);
+		if(ier==-1)return(NULL);
+
+		data->npts=data->grid->nPts;
+		data->neq=data->nvar*data->npts;
+	}
+
+	data->output=fopen("output.dat","w");
+	data->globalOutput=fopen("globalOutput.dat","w");
+	data->gridOutput=fopen("grid.dat","w");
+	//data->ratesOutput=fopen("rates.dat","w");
+	
+	data->innerMassFractions = new double [data->nsp];
+
+	return(data);
+}
+
+void setSaneDefaults(UserData data){
+	data->domainLength=1.0e-02;
+	data->constantPressure=1;
+	data->problemType=1;
+	data->dPdt=0.0e0;
+	data->reflectProblem=0;
+	data->mdot=0.0;
+	data->initialTemperature=300.0;
+	data->initialPressure=1.0;
+	data->metric=0;
+	data->ignTime=1e-09;
+	data->maxQDot=0.0e0;
+	data->maxTemperature=300.0e0;
+	data->mixingWidth=1e-04;
+	data->shift=3.0e-03;
+	data->firstRadius=1e-04;
+	data->wallTemperature=298.0e0;
+	data->dirichletInner=0;
+	data->dirichletOuter=0;
+	data->adaptiveGrid=0;
+	data->moveGrid=0;
+	data->gridOffset=0.0e0;
+	data->isotherm=1000.0;
+	data->nSaves=30;
+	data->writeRates=0;
+	data->writeEveryRegrid=0;
+	data->writeDeltaT=1e-04;
+	data->relativeTolerance=1e-06;
+	data->radiusTolerance=1e-08;
+	data->temperatureTolerance=1e-06;
+	data->pressureTolerance=1e-06;
+	data->massFractionTolerance=1e-09;
+	data->bathGasTolerance=1e-06;
+	data->finalTime=1e-02;
+	data->tNow=0.0e0;
+	data->setConstraints=0;
+	data->suppressAlg=1;
+	data->regrid=0;
+	data->gridDirection=1;
+	data->dryRun=0;
+	data->nThreads=1;
+
+	data->flamePosition[0]=0.0e0;
+	data->flamePosition[1]=0.0e0;
+	data->flameTime[0]=0.0e0;
+	data->flameTime[1]=0.0e0;
+	data->nTimeSteps=0;
+}
+
diff --git a/UserData.h b/UserData.h
new file mode 100644
index 0000000..7578c4c
--- /dev/null
+++ b/UserData.h
@@ -0,0 +1,170 @@
+#ifndef CANTERA_DEF
+#define CANTERA_DEF
+#include <cantera/IdealGasMix.h>
+#include <cantera/transport.h>
+#endif
+
+#include "gridRoutines.h"
+
+#ifndef USER_DEF
+#define USER_DEF
+typedef struct UserDataTag{
+	/*An ideal gas object from Cantera. Contains thermodynamic and kinetic
+	 * info of all species.*/
+	Cantera::IdealGasMix* gas;
+	/*A Transport object from Cantera. Contains all transport info of all
+	 * species.*/
+	Cantera::Transport* trmix;
+	/*Length of the domain (in meters):*/			  
+	double domainLength;	  
+	/*Mass of gas in domain (in kg):*/			  
+	double mass;		  
+	/*Parameter that indicates the symmetry of the problem;*/			  
+	/*metric=0:Planar*/			  
+	/*metric=1:Cylindrical*/			  
+	/*metric=2:Spherical*/			  
+	int metric;		  
+	/*No: of species:*/
+	size_t nsp;	  
+	/*No: of equations:*/
+	size_t neq;
+	/*No: of variables:*/
+	size_t nvar;
+	/*Pointer indices (see "macros.h" for aliases that use these):*/
+	/*Pointer index for temperature:*/
+	size_t nt;	  
+	/*Pointer index for species:*/
+	size_t ny;
+	/*Pointer index for spatial coordinate:*/
+	size_t nr;
+	/*Pointer index for pressure:*/
+	size_t np;
+	/*Species index of bath gas:*/			  
+	size_t k_bath;		  
+	/*Species index of oxidizer:*/			  
+	size_t k_oxidizer;		  
+	size_t k_OH;
+	size_t k_HO2;
+	/*User-defined mass flux (kg/m^2/s):*/			  
+	double mdot;
+	/*Flag to solve isobaric/isochoric problem;*/
+	/*constantPressure=1: isobaric*/
+	/*constantPressure=0: isochoric*/
+	int constantPressure;
+	/*User-defined dPdt (Pa/s), activates when problem is "isobaric":*/			  
+	double dPdt;	
+	/*Initial temperature of the gas (K):*/
+	double initialTemperature;
+	/*Initial Pressure of the gas (atm):*/
+	double initialPressure;
+	/*Classification of problem type;*/
+	/*problemType=0: Mixture is premixed and spatially uniform initially.
+	 * In order for mixture to ignite, an external heat source (finite
+	 * maxQDot) must be used.*/
+	/*problemType=1: Mixture is premixed but spatially non-uniform
+	 * initially. Equilibrium products are contained within a hot kernel of
+	 * size given by "shift" and a mixing length scale given by
+	 * "mixingWidth".*/
+	/*problemType=2: User specified initial condition. Use file
+	 * "initialCondition.dat".*/
+	int problemType;
+	/*Maximum External heat source (K/s):*/
+	double maxQDot;
+	/*Ignition kernel size:*/
+	double kernelSize;
+
+	double maxTemperature;
+	/*Maximum time for which the external heat source is applied (s):*/
+	double ignTime;
+	/*Vector of Mass Fractions used to impose Robin Boundary Condition for
+	 * species at the domain origin:*/
+	double* innerMassFractions;
+	/*Value of temperature to be used if Dirichlet Boundary Conditions are
+	 * imposed for temperature:*/
+	double innerTemperature;
+	double wallTemperature;
+	/*Isotherm chosen to find the location of a "burning" front (K):*/
+	double isotherm;
+	/*Interval of time integration:*/
+	double finalTime;
+	/*Current time:*/
+	double tNow;
+	/*Flag to reflect initial conditions across center of the domain:*/
+	int reflectProblem;
+	/*Parameters for initial conditions in setting up profiles:
+	increasing function of x: g=0.5*(erf(x-3*w-shift)/w)+1)
+	decreasing function of x: f=1-g*/
+	double mixingWidth;
+	double shift;
+	double firstRadius;
+	/*Flag to run program without time-integration i.e. simply layout the
+	 * initial conditions and quit:*/
+	int dryRun;
+	/*Relative Tolerance:*/
+	double relativeTolerance;
+	/*Absolute Tolerance for spatial coordinate:*/
+	double radiusTolerance;
+	/*Absolute Tolerance for Temperature:*/
+	double temperatureTolerance;
+	/*Absolute Tolerance for Pressure:*/
+	double pressureTolerance;
+	/*Absolute Tolerance for Mass Fractions:*/
+	double massFractionTolerance;
+	/*Absolute Tolerance for bath gas mass fraction:*/
+	double bathGasTolerance;
+	/*Flag to set constraints on Mass fractions so they don't acquire
+	 * negative values:*/
+	int setConstraints;
+	/*Flag to suppress error checking on algebraic variables:*/
+	int suppressAlg;
+	/*Number of time-steps elapsed before saving the solution:*/
+	int nSaves;		  
+	/*Flag to set write for every regrid:*/
+	int writeEveryRegrid;
+	double writeDeltaT;
+	/*Solution output file:*/
+	FILE* output;
+	/*Flag to write the rates (ydot) of solution components into the
+	 * "ratesOutput" file:*/
+	int writeRates;
+	/*Grid output file:*/
+	FILE* gridOutput;
+	///*Rate of change (ydot) output file (see "writeRates"):*/
+	//FILE* ratesOutput;
+	/*Global properties (mdot, radius of flame, etc.) output file:*/
+	FILE* globalOutput;
+
+	/*Flag to adapt grid:*/
+	int adaptiveGrid;
+	/*Flag to move grid:*/
+	int moveGrid;
+	/*Flag to initiate regrid:*/
+	int regrid;
+	/*Integer that specifies grid movement direction:
+	 * gridDirection = -1: Move Left
+	 * gridDirection = +1: Move Right*/ 
+	int gridDirection;
+	/*Total number of points for grid:*/
+	size_t npts;
+
+	double gridOffset;
+
+	UserGrid grid;
+	double* uniformGrid;
+
+	int dirichletInner,dirichletOuter;
+
+	int nThreads;
+	double clockStart;
+	
+	/*These arrays are used to compute dr/dt, which in turn is used to
+	 * compute the flame speed S_u:*/
+	double flamePosition[2];
+	double flameTime[2];
+	size_t nTimeSteps;
+} *UserData;
+UserData allocateUserData(FILE *input);
+void setSaneDefaults(UserData data);
+void freeUserData(UserData data);
+#endif
+
diff --git a/gridRoutines.cpp b/gridRoutines.cpp
new file mode 100644
index 0000000..f4a9065
--- /dev/null
+++ b/gridRoutines.cpp
@@ -0,0 +1,244 @@
+#include "gridRoutines.h"
+#include <stdio.h>
+inline double l(const double* x, 
+	 	const double* a,
+	 	const double* w,
+		const double* fac,
+		const int* refineLeft){
+	if(*refineLeft==0){
+		return(tanh(-*a*(*x+*w*100.0e0)));
+	}else{
+		return(tanh(-*a*(*x-*w*(*fac))));
+	}
+}
+
+inline double r(const double* x, 
+		const double* a, 
+		const double* w,
+		const double* fac,
+		const int* refineRight){
+	if(*refineRight==0){
+		return(tanh(*a*(*x-(1.0e0+*w*100.0e0))));
+	}else{
+		return(tanh(*a*(*x-(1.0e0-*w*(*fac)))));
+	}
+}
+
+inline double f(const double* x,
+		const double* a,
+		const double* c,
+		const double* w){
+	return(tanh(-*a*(*x-(*c+*w)))
+		+tanh(-*a*((*x-1.0e0)-(*c+*w)))
+		+tanh(-*a*((*x+1.0e0)-(*c+*w))));
+}
+
+inline double g(const double* x, 
+		const double* a, 
+		const double* c, 
+		const double* w){
+	return(tanh(*a*(*x-(*c-*w)))
+	       +tanh(*a*((*x-1.0e0)-(*c-*w)))
+	       +tanh(*a*((*x+1.0e0)-(*c-*w))));
+}
+
+inline double rho(const double* x,
+	   	  const double* a,
+	   	  const double* c,
+	   	  const double* w,
+	   	  const double* mag,
+	   	  const double* leftFac,
+	   	  const double* rightFac,
+		  const int* refineLeft,
+		  const int* refineRight){
+
+	return(((2.0e0+f(x,a,c,w)
+		 +g(x,a,c,w)
+		 +l(x,a,w,leftFac,refineLeft)
+		 +r(x,a,w,rightFac,refineRight))*0.5e0)
+		  *(*mag-1.0e0)+1.0e0);
+}
+
+size_t maxPoints(const size_t basePts,
+	   	 const double* a,
+	   	 const double* w,
+	   	 const double* mag,
+	   	 const double* leftFac,
+	   	 const double* rightFac,
+	   	 const int* refineLeft,
+	   	 const int* refineRight){
+	double dx=1.0e0/((double)(basePts)-1.0e0);
+	double y=0.0e0;
+	size_t i=0;
+	double r=0.0e0;
+	double t=0.5e0;
+	while(y<=1.0e0){
+		r=rho(&y,a,&t,w,mag,leftFac,rightFac,refineLeft,refineRight);
+		y=y+(dx/r);
+		i++;
+	}
+	return(i);
+}
+
+void fillGrid(const size_t* basePts,
+	      const size_t* nPts,
+	      const double* a,
+	      const double* c,
+	      const double* w,
+	      const double* mag,
+	      const double* leftFac,
+	      const double* rightFac,
+	      const int* refineLeft,
+	      const int* refineRight,
+	      double x[]){
+
+	FILE* out;out=fopen("tmp.dat","w");
+
+	double y=0.0e0;
+	double r=0.0e0;
+	double dx=1.0e0/((double)(*basePts)-1.0e0);
+	for(size_t j=0;j<*nPts;j++){
+		r=rho(&y,a,c,w,mag,leftFac,rightFac,refineLeft,refineRight);
+		fprintf(out, "%15.15e\n",dx/r);
+		y=y+(dx/r);
+	}
+	fclose(out);
+	
+	double dxp[*nPts-1];
+	for (size_t j = 0; j < *nPts; j++) {
+		dxp[j]=0.0e0;
+	}
+
+	FILE* tmp;tmp=fopen("tmp.dat","r");
+	char buf[MAXBUFLEN];
+	size_t i=0;
+	while (fgets(buf,MAXBUFLEN, tmp)!=NULL){
+		sscanf(buf, "%lf", &y);
+		dxp[i]=y;
+		i++;
+	}
+	fclose(tmp);
+
+	double sum=0.0e0;
+	double err=0.0e0;
+	double fix=0.0e0;
+	for(size_t j=0;j<*nPts-1;j++){
+		sum+=dxp[j];
+	}
+	err=1.0e0-sum;
+	printf("sum before correction: %15.6e\n",sum);
+	printf("err before correction: %15.6e\n",err);
+	fix=err/((double)(*nPts));
+	sum=0.0e0;
+	for(size_t j=0;j<*nPts-1;j++){
+		dxp[j]+=fix;
+		sum+=dxp[j];
+	}
+	err=1.0e0-sum;
+	printf("sum after correction:%15.6e\n",sum);
+	printf("err after correction: %15.6e\n",err);
+
+	x[0]=0.0e0;
+	for(size_t j=0;j<*nPts-1;j++){
+		x[j+1]=x[j]+dxp[j];
+	}
+	x[*nPts-1]=1.0e0;
+
+}
+
+double safePosition(double c, double w){
+	if(c<w){
+		return(w);
+	}
+	else if(c>1.0e0-w){
+		return(1.0e0-w);
+	}
+	else{
+		return(c);
+	}
+}
+
+int reGrid(UserGrid grid, double position){
+
+	printf("before regrid: %ld\n", grid->nPts);
+	double xx[grid->nPts];
+
+	fillGrid(&grid->basePts,
+		 &grid->nPts,
+		 &grid->a,
+		 &position,
+		 &grid->w,
+		 &grid->mag,
+		 &grid->leftFac,
+		 &grid->rightFac,
+		 &grid->refineLeft,
+		 &grid->refineRight,
+		 xx);
+
+	for (size_t i = 0; i < grid->nPts; i++) {
+		grid->x[i]=xx[i];
+	}
+	return(0);
+}
+
+void storeGrid(const double* x, double *y, const size_t nPts){
+	for(size_t i=0;i<nPts;i++){
+		y[i]=x[i];
+	}
+}
+
+int initializeGrid(UserGrid grid){
+	grid->nPts=maxPoints(grid->basePts,
+			     &grid->a,
+			     &grid->w,
+			     &grid->mag,
+			     &grid->leftFac,
+			     &grid->rightFac,
+			     &grid->refineLeft,
+			     &grid->refineRight);
+	printf("nPts: %ld\n",grid->nPts);
+	grid->leastMove=grid->w;
+
+	grid->x    = new double [grid->nPts];
+	grid->xOld = new double [grid->nPts];
+	for (size_t i = 0; i < grid->nPts; i++) {
+		grid->x[i]=0.0e0;
+		grid->xOld[i]=0.0e0;
+	}
+	return(0);
+}
+
+int getGridSettings(FILE *input, UserGrid grid){
+
+	int ier=0;
+
+	ier=parseNumber<size_t>(input, "basePts"	, MAXBUFLEN, &grid->basePts);
+	if(ier==-1)return(-1);
+
+	ier=parseNumber<double>(input, "gridDensitySlope", MAXBUFLEN, &grid->a);
+	if(ier==-1)return(-1);
+
+	ier=parseNumber<double>(input, "fineGridHalfWidth", MAXBUFLEN, &grid->w);
+	if(ier==-1)return(-1);
+
+	ier=parseNumber<double>(input, "gridRefinement", MAXBUFLEN, &grid->mag);
+	if(ier==-1)return(-1);
+
+	ier=parseNumber<double>(input, "leftRefineFactor", MAXBUFLEN, &grid->leftFac);
+	if(ier==-1)return(-1);
+
+	ier=parseNumber<double>(input, "rightRefineFactor", MAXBUFLEN, &grid->rightFac);
+	if(ier==-1)return(-1);
+
+	ier=parseNumber<int>(input, "refineLeft"	, MAXBUFLEN, &grid->refineLeft);
+	if(ier==-1)return(-1);
+
+	ier=parseNumber<int>(input, "refineRight"	, MAXBUFLEN, &grid->refineRight);
+	if(ier==-1)return(-1);
+
+	ier=parseNumber<double>(input, "position"	, MAXBUFLEN, &grid->position);
+	if(ier==-1)return(-1);
+
+	return(0);
+}
+
diff --git a/gridRoutines.h b/gridRoutines.h
new file mode 100644
index 0000000..06231ac
--- /dev/null
+++ b/gridRoutines.h
@@ -0,0 +1,82 @@
+#include "parse.h"
+
+#ifndef GSL_DEF
+#define GSL_DEF
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_spline.h>
+#endif
+
+#ifndef GRID_DEF
+#define GRID_DEF
+
+typedef struct gridTag{
+
+	size_t basePts;
+	size_t nPts;
+	double position;
+	double leastMove;
+	double a;
+	double w;
+	double mag;
+	double leftFac;
+	double rightFac;
+	int refineLeft;
+	int refineRight;
+	double* x;
+	double* xOld;
+
+} *UserGrid;
+
+inline double l(const double* x, 
+	 	const double* a,
+	 	const double* w,
+		const double* fac,
+		const int* refineLeft);
+
+inline double r(const double* x, 
+		const double* a, 
+		const double* w,
+		const double* fac,
+		const int* refineRight);
+
+inline double f(const double* x,
+		const double* a,
+		const double* c,
+		const double* w);
+
+inline double g(const double* x, 
+		const double* a, 
+		const double* c, 
+		const double* w);
+
+inline double rho(const double* x,
+	   	  const double* a,
+	   	  const double* c,
+	   	  const double* w,
+	   	  const double* mag,
+	   	  const double* leftFac,
+	   	  const double* rightFac,
+		  const int* refineLeft,
+		  const int* refineRight);
+
+size_t maxPoints(const size_t basePts,
+	   	 const double* a,
+	   	 const double* w,
+	   	 const double* mag,
+	   	 const double* leftFac,
+	   	 const double* rightFac,
+	   	 const int* refineLeft,
+	   	 const int* refineRight);
+
+
+double safePosition(double c, double w);
+
+int reGrid(UserGrid grid, double position);
+
+void storeGrid(const double* x, double *y, const size_t nPts);
+
+int initializeGrid(UserGrid grid);
+
+int getGridSettings(FILE *input, UserGrid grid);
+
+#endif
diff --git a/input.dat b/input.dat
new file mode 100644
index 0000000..a139184
--- /dev/null
+++ b/input.dat
@@ -0,0 +1,50 @@
+#Problem constants:
+domainLength=10.16e-02		#in meters
+constantPressure=0		#1=isobaric, 0=isochoric
+dPdt=0.0e0			#Units: Pa/s
+problemType=1			#0=heatsource to ignite(premix),1=premix, 2=custom initial condition
+reflectProblem=0		#If 1, the initial conditions are reflected across the center of the domain. Useful for spherically imploding flames for instance.
+mdot=0.0
+initialTemperature=300.0	#in Kelvin
+initialPressure=1.0		#in atmospheres
+chemistryFile=./ch4_1step.xml	#File containing thermodynamic, kinetic, and transport info
+transportModel=Mix		#Transport model
+mixtureComposition=CH4:0.45,O2:1.0,N2:3.76	#Mole Fractions!
+metric=2			#0=planar,1=cylindrical,2=spherical
+ignTime=5e-07			#Time over which energy is supplied in seconds
+QDot=1e+14			#Power supplied in units of K/s
+kernelSize=0e-03
+mixingWidth=2e-03		#mixing layer width in m (for diffusion problem)
+shift=4e-03		#shift of mixing layer into domain (for diffusion problem)
+dirichletInner=0		#Boundary conditions for scalars (T,Y)
+dirichletOuter=1
+wallTemperature=300.0
+#Grid settings:
+adaptiveGrid=1
+moveGrid=1
+gridOffset=0
+basePts=32
+gridDensitySlope=70.0e0
+fineGridHalfWidth=8e-03
+gridRefinement=1000.0e0;
+leftRefineFactor=0.01
+rightRefineFactor=0.01
+refineLeft=1
+refineRight=1
+position=0.0
+isotherm=1000.0
+#Solver settings:
+nSaves=100			#Save solution every nSaves timesteps
+writeEveryRegrid=1
+writeRates=0
+relativeTolerance=1e-06
+radiusTolerance=1e-09		#absolute tolerance for radius
+temperatureTolerance=1e-09
+pressureTolerance=1e-09
+massFractionTolerance=1e-09
+bathGasTolerance=1e-09
+finalTime=1.0e0			#in seconds
+setConstraints=0
+suppressAlg=0
+nThreads=1
+dryRun=0
diff --git a/macros.h b/macros.h
new file mode 100644
index 0000000..dcf36db
--- /dev/null
+++ b/macros.h
@@ -0,0 +1,63 @@
+#ifndef MACRO_DEF
+#define MACRO_DEF
+//#define SUNDIALS_DOUBLE_PRECISION 1
+//#define SUNDIALS_SINGLE_PRECISION 1
+
+#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 FOUR  RCONST(4.0)
+#define TEN   RCONST(10.0)
+
+/* 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 psi(i)   psidata[i-1]
+
+#define T(i)   ydata[((i-1)*data->nvar)+data->nt]
+#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 Tdot(i)   ydotdata[((i-1)*data->nvar)+data->nt]
+#define Ydot(i,k) ydotdata[((i-1)*data->nvar)+data->ny+k-1]
+#define Rdot(i)   ydotdata[((i-1)*data->nvar)+data->nr]
+#define Pdot(i)   ydotdata[((i-1)*data->nvar)+data->np]
+
+#define Tres(i)   resdata[((i-1)*data->nvar)+data->nt]
+#define Yres(i,k) resdata[((i-1)*data->nvar)+data->ny+k-1]
+#define Rres(i)   resdata[((i-1)*data->nvar)+data->nr]
+#define Pres(i)   resdata[((i-1)*data->nvar)+data->np]
+
+#define Tid(i)   iddata[((i-1)*data->nvar)+data->nt]
+#define Yid(i,k)   iddata[((i-1)*data->nvar)+data->ny+k-1]
+#define Rid(i)   iddata[((i-1)*data->nvar)+data->nr]
+#define Pid(i)   iddata[((i-1)*data->nvar)+data->np]
+
+#define Yav(i) Yav[i-1]
+#define YAvg(i) YAvg[i-1]
+#define YVmhalf(i) YVmhalf[i-1]
+#define YVphalf(i) YVphalf[i-1]
+#define X(i) X[i-1]
+#define Xp(i) Xp[i-1]
+#define Xgradhalf(i) Xgradhalf[i-1]
+#define XLeft(i) XLeft[i-1]
+#define XRight(i) XRight[i-1]
+#define gradX(i) gradX[i-1]
+#define wdot(i) wdot[i-1]
+#define enthalpy(i) enthalpy[i-1]
+#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]
+#define atolP(i)   atolvdata[((i-1)*data->nvar)+data->np]
+
+
+#define constraintsY(i,k)   constraintsdata[((i-1)*data->nvar)+data->ny+k-1]
+
+#endif
diff --git a/main.cpp b/main.cpp
new file mode 100644
index 0000000..f9dc320
--- /dev/null
+++ b/main.cpp
@@ -0,0 +1,338 @@
+/*
+ _____ ___  ____   ____ 
+|_   _/ _ \|  _ \ / ___|
+  | || | | | |_) | |    
+  | || |_| |  _ <| |___ 
+  |_| \___/|_| \_\\____|
+                        
+*/
+
+#include "UserData.h"
+#include "solution.h"
+#include "residue.h"
+#include "macros.h"
+#include "timing.h"
+
+#include <ida/ida.h>
+#include <ida/ida_direct.h>
+#include <sunmatrix/sunmatrix_band.h>
+#include <sunlinsol/sunlinsol_lapackband.h>
+//#include <ida/ida_band.h>
+
+static int check_flag(void *flagvalue, 
+		      const char *funcname, 
+		      int opt);
+
+void freeAtLast(void* mem, N_Vector *y, 
+		N_Vector *ydot, 
+		N_Vector *res, 
+		N_Vector *id, 
+		N_Vector *atolv,
+		N_Vector *constraints,UserData data);
+
+int main(){
+
+  	FILE *input;input=fopen("input.dat","r");
+  	UserData data;data=NULL;data=allocateUserData(input);
+	fclose(input);
+	data->clockStart=get_wall_time();
+
+	if(data==NULL){
+		printf("check input file!\n");
+  		freeUserData(data);
+		return(-1);
+	}
+	long int ier,mu,ml,count,netf,ncfn,njevals,nrevals;
+	realtype tNow,*atolvdata,*constraintsdata,finalTime,tolsfac;
+
+  	N_Vector y,ydot,id,res,atolv,constraints;
+  	y=ydot=id=res=atolv=constraints=NULL;
+ 	ier=allocateSolution(data->neq,data->nThreads,&y,&ydot,&id,&res,&atolv,&constraints);
+	ier=setAlgebraicVariables(&id,data);
+	ier=setInitialCondition(&y,&ydot,data);
+	if(ier==-1)return(-1);
+	tNow=data->tNow;
+	finalTime=data->finalTime;
+
+	double* ydata;
+	double* ydotdata;
+  	ydata    = N_VGetArrayPointer_OpenMP(y);
+  	ydotdata = N_VGetArrayPointer_OpenMP(ydot);
+
+	void *mem;mem=NULL;mem = IDACreate();
+  	ier = IDASetUserData(mem, data);
+  	ier = IDASetId(mem, id);
+  	ier = IDAInit(mem, residue, tNow, y, ydot);
+
+        // Atol array
+        atolvdata = N_VGetArrayPointer_OpenMP(atolv);
+	for (size_t i = 1; i <=data->npts; i++) {
+             atolT(i) = data->temperatureTolerance;
+             for (size_t k = 1; k <=data->nsp; k++) {
+		     	if(k!=data->k_bath){
+                		atolY(i,k) = data->massFractionTolerance;
+			}
+			else{
+                		atolY(i,k) = data->bathGasTolerance;
+			}
+
+             }
+             atolR(i) = data->radiusTolerance;
+             atolP(i) = data->pressureTolerance;
+	}
+  	ier = IDASVtolerances(mem, data->relativeTolerance, atolv);
+
+  	mu = 2*data->nvar; ml = mu;
+	SUNMatrix A; A=NULL;
+	A=SUNBandMatrix(data->neq,mu,ml,mu+ml);
+	SUNLinearSolver LS; LS=NULL;
+	//LS=SUNBandLinearSolver(y,A);
+	LS=SUNLapackBand(y,A);
+	ier=IDADlsSetLinearSolver(mem,LS,A);
+  	//ier = IDABand(mem, data->neq, mu, ml);
+
+        constraintsdata = N_VGetArrayPointer_OpenMP(constraints);
+	if(data->setConstraints){
+		for (size_t i = 1; i <=data->npts; i++) {
+             		for (size_t k = 1; k <=data->nsp; k++) {
+             		   constraintsY(i,k) = ONE;
+             		}
+		}
+  		ier=IDASetConstraints(mem, constraints);
+	}
+
+	ier = IDASetSuppressAlg(mem, data->suppressAlg);
+  	if(check_flag(&ier, "IDASetSuppressAlg", 1)) return(1);
+
+	//ier= IDASetMaxNumStepsIC(mem, 1);
+	//ier= IDASetMaxNumJacsIC(mem,8);
+	//ier= IDASetMaxNumItersIC(mem,100);
+	//ier= IDASetMaxBacksIC(mem,2000);
+	//ier = IDASetLineSearchOffIC(mem,SUNTRUE);
+	
+	if(!data->dryRun){
+		printf("Calculating Initial Conditions:\n");
+		printf("Cross your fingers...\n");
+  		ier = IDACalcIC(mem, IDA_YA_YDP_INIT, 1e-05*finalTime);
+		//If at first IDACalcIC doesn't succeed, try, try, try again:
+		for (int i = 0; i < 10; i++) {
+  			ier = IDACalcIC(mem, IDA_YA_YDP_INIT, (1e-01+pow(10,i)*finalTime));
+			if(ier==0){
+				break;
+			}
+		}
+		//...or else cry again :(
+  		if(check_flag(&ier, "IDACalcIC", 1)){
+			freeAtLast(mem,&y,&ydot,&res,&id,&atolv,&constraints,data);
+		       	return(-1);
+		}else{
+			printf("Initial (Consistent) Conditions Calculated!\n");
+		}
+		ier = IDASetInitStep(mem,1e-12);
+	}
+
+	printSpaceTimeHeader(data);
+	printGlobalHeader(data);
+	printSpaceTimeOutput(tNow, &y, data->output, data);
+	printSpaceTimeOutput(tNow, &y, data->gridOutput, data);
+
+	if(!data->dryRun){
+		count=0;
+		double dt=1e-08;
+		double t1=0.0e0;
+		double xOld=0.0e0;
+		double x=0.0e0;
+		double dx=0.0e0;
+		double dxMin=1.0e0;
+		double dxRatio=dx/dxMin;
+		double writeTime=tNow;
+		int move=0;
+		int kcur=0;
+		if(data->adaptiveGrid){
+			dxMin=data->grid->leastMove;
+			//xOld=maxCurvPosition(ydata, data->nt, data->nvar,
+			//		     data->grid->x, data->npts);
+			xOld=isothermPosition(ydata, data->isotherm, data->nt, 
+	                		   data->nvar, data->grid->x, data->npts);
+		}
+		while (tNow<=finalTime) {
+
+			t1=tNow;
+  			ier = IDASolve(mem, finalTime, &tNow, y, ydot, IDA_ONE_STEP);
+  			if(check_flag(&ier, "IDASolve", 1)){
+				freeAtLast(mem,&y,&ydot,&res,&id,&atolv,&constraints,data);
+		       		return(-1);
+			}
+			dt=tNow-t1;
+
+  			ier = IDAGetCurrentOrder(mem, &kcur);
+
+			if(data->adaptiveGrid==1 && data->moveGrid==1){
+				//x=maxCurvPosition(ydata, data->nt, data->nvar,
+				//		  data->grid->x, data->npts);
+
+				x=isothermPosition(ydata, data->isotherm, data->nt, 
+	                			   data->nvar, data->grid->x, data->npts);
+
+				//x=maxGradPosition(ydata, data->nt, data->nvar,
+				//		  data->grid->x, data->npts);
+				dx=x-xOld;
+
+				if(dx*dxMin>0.0e0){
+					move=1;
+				}else{
+					move=-1;
+				}
+
+				//if(fabs(dx)>=dxMin && x+(double)(move)*0.5e0*dxMin<=1.0e0){
+				dxRatio=fabs(dx/dxMin);
+				if(dxRatio>=1.0e0 && dxRatio<=2.0e0){
+					printf("Regridding!\n");
+
+					data->regrid=1;
+					//printSpaceTimeOutput(tNow, &y, data->gridOutput, data);
+
+					ier=reGrid(data->grid, x+(double)(move)*0.5e0*dxMin);
+					if(ier==-1){
+						freeAtLast(mem,&y,&ydot,&res,&id,&atolv,&constraints,data);
+						return(-1);
+					}
+
+					updateSolution(ydata, ydotdata, data->nvar,
+					               data->grid->xOld,data->grid->x,data->npts);
+					storeGrid(data->grid->x,data->grid->xOld,data->npts);
+					xOld=x;
+
+					printf("Regrid Complete! Restarting Problem at %15.6e s\n",tNow);
+					ier = IDAReInit(mem, tNow, y, ydot);
+  					if(check_flag(&ier, "IDAReInit", 1)){
+						freeAtLast(mem,&y,&ydot,&res,&id,&atolv,&constraints,data);
+		       				return(-1);
+					}
+					ier = IDASetInitStep(mem,1e-01*dt);
+					printf("Reinitialized! Calculating Initial Conditions:\n");
+					printf("Cross your fingers...\n");
+  					ier = IDACalcIC(mem, IDA_YA_YDP_INIT, tNow+1e-01*dt);
+  					if(check_flag(&ier, "IDACalcIC", 1)){
+  						ier = IDACalcIC(mem, IDA_YA_YDP_INIT, tNow+1e+01*dt);
+					}
+					//Every once in a while, for reasons
+					//that befuddle this mathematically
+					//lesser author, IDACalcIC fails. Here,
+					//I desperately try to make it converge
+					//again by sampling increasingly larger
+					//time-steps:
+					for (int i = 0; i < 10; i++) {
+  						ier = IDACalcIC(mem, IDA_YA_YDP_INIT, tNow+(1e-01+pow(10,i)*dt));
+						if(ier==0){
+							break;
+						}
+					}
+					//Failure :( Back to the drawing board:
+  					if(check_flag(&ier, "IDACalcIC", 1)){
+						freeAtLast(mem,&y,&ydot,&res,&id,&atolv,&constraints,data);
+		       				return(-1);
+					}
+					printf("Initial (Consistent) Conditions Calculated!\n");
+					printf("------------------------------------------\n\n");
+					//if(data->writeEveryRegrid){
+
+				}
+			}
+
+			if(count%data->nSaves==0 && !data->writeEveryRegrid){
+				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;
+				data->flamePosition[1]=0.0e0;
+				data->flameTime[0]=tNow;
+				data->flameTime[1]=tNow;
+			}
+			if(kcur>=2 ){
+			      printGlobalVariables(tNow,&y,&ydot,data);
+			}
+
+			if(tNow>writeTime){
+				writeTime+=data->writeDeltaT;
+				printSpaceTimeOutput(tNow, &y, data->output, data);
+				FILE* fp;
+				fp=fopen("restart.bin","w");
+				writeRestart(tNow, &y, &ydot, fp, data);
+				fclose(fp);
+			}
+
+  			ier = IDAGetNumErrTestFails(mem, &netf);
+  			ier = IDAGetNumNonlinSolvConvFails(mem, &ncfn);
+  			ier = IDADlsGetNumJacEvals(mem, &njevals);
+  			ier = IDADlsGetNumResEvals(mem, &nrevals);
+  			printf("etf = %ld ,"
+			       "nlsf= %ld ," 
+			       "J=%ld ,"
+			       "R=%ld ,"
+			       "o=%d ,",netf, ncfn, njevals, nrevals, kcur);
+			printf("Time=%15.6e s,",tNow);
+			printf("frac: %15.6e\n",dxRatio);
+			count++;
+			data->nTimeSteps=count;
+		}
+	}
+
+	SUNLinSolFree(LS);
+	SUNMatDestroy(A);
+	freeAtLast(mem,&y,&ydot,&res,&id,&atolv,&constraints,data);
+
+	return(0);
+}
+
+void freeAtLast(void* mem, 
+		N_Vector *y, 
+		N_Vector *ydot, 
+		N_Vector *res, 
+		N_Vector *id, 
+		N_Vector *atolv,
+		N_Vector *constraints,UserData data){
+
+  	IDAFree(&mem);
+	freeSolution(y,ydot,res,id,atolv,constraints);
+  	freeUserData(data);
+}
+
+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);
+}
diff --git a/parse.cpp b/parse.cpp
new file mode 100644
index 0000000..f872ca5
--- /dev/null
+++ b/parse.cpp
@@ -0,0 +1,20 @@
+#include "parse.h"
+void getFromString (const char* buf, int* n){
+	*n=atoi(buf);
+	printf("%d\n",*n);
+}
+
+void getFromString (const char* buf, size_t* n){
+	*n=(size_t)(atoi(buf));
+	printf("%lu\n",*n);
+}
+
+void getFromString (const char* buf, double* n){
+	*n=(double)(atof(buf));
+	printf("%15.6e\n",*n);
+}
+
+void getFromString (const char* buf, char* n){
+	sscanf(buf,"%s",n);
+	printf("%s\n",n);
+}
diff --git a/parse.h b/parse.h
new file mode 100644
index 0000000..45f3a2c
--- /dev/null
+++ b/parse.h
@@ -0,0 +1,31 @@
+#ifndef PRINT_DEF
+#define PRINT_DEF
+#include <string.h>	//for strings
+#include <stdio.h>	//for printf,scanf
+#include <stdlib.h>	//for atoi, atof
+#endif
+
+#ifndef PARSE_DEF
+#define PARSE_DEF
+#define MAXBUFLEN 200
+
+void getFromString (const char* buf, int* n);
+void getFromString (const char* buf, size_t* n);
+void getFromString (const char* buf, double* n);
+void getFromString (const char* buf, char* n);
+
+int parseString(FILE* input, const char* keyword, const size_t bufLen, char* n);
+
+template<typename T>
+int parseNumber(FILE* input, const char* keyword, const size_t bufLen, T* n);
+
+template<typename T>
+int parseArray(FILE* input, const char* keyword, const size_t bufLen,
+		const size_t arrLen, T y[]);
+
+
+#include "parse.hpp"
+
+#endif
+
+
diff --git a/parse.hpp b/parse.hpp
new file mode 100644
index 0000000..634db96
--- /dev/null
+++ b/parse.hpp
@@ -0,0 +1,76 @@
+template<typename T>
+int parseNumber(FILE* input, const char* keyword, const size_t bufLen, T* n){
+
+	char buf[bufLen];
+	char buf1[bufLen];
+	char comment[1];
+	char *ret;
+
+	while (fgets(buf,bufLen, input)!=NULL){
+		comment[0]=buf[0];
+		if(strncmp(comment,"#",1)==0){
+			//printf("Comment!:%s\n",buf);
+		}
+		else{
+			ret=strtok(buf,"=");
+			if(strcmp(ret,keyword)==0){
+				/*offset buf by keyword size + 1 for the "="*/
+				strncpy (buf1, buf+strlen(keyword)+1, bufLen);	
+       				printf("%30s: ",keyword);
+				getFromString(buf1,n);
+				rewind(input);
+				return(0);
+			}
+		}
+	}
+	rewind(input);
+	return(-1);
+}
+
+template<typename T>
+int parseArray(FILE* input, const char* keyword, const size_t bufLen,
+		const size_t arrLen, T y[]){
+
+	char buf[bufLen];
+	char buf1[bufLen];
+	char comment[1];
+	char *ret;
+
+	while (fgets(buf,bufLen, input)!=NULL){
+		comment[0]=buf[0];
+		if(strncmp(comment,"#",1)==0){
+			//printf("Comment!:%s\n",buf);
+		}
+		else{
+			ret=strtok(buf,"=");
+
+			if(strcmp(ret,keyword)==0){
+				/*offset buf by keyword size + 1 for the "="*/
+				strncpy (buf1, buf+strlen(keyword)+1, bufLen);	
+       				printf("%30s:\n",keyword);
+				ret=strtok(buf1,",");
+				size_t j=0;
+				while(ret!=NULL){
+					if(j<arrLen){
+						//y[j]=atof(ret);
+						getFromString(ret,&y[j]);
+					}
+   		  			ret=strtok(NULL,",");
+					j++;
+				}
+				rewind(input);
+				if(j!=arrLen){
+					printf("Check no: of values entered for %s\n",keyword);
+					printf("%lu values required!\n",arrLen);
+					return(-1);
+				}
+				else{
+					return(0);
+				}
+			}
+		}
+	}
+	rewind(input);
+	return(-1);
+}
+
diff --git a/residue.cpp b/residue.cpp
new file mode 100644
index 0000000..57e20cd
--- /dev/null
+++ b/residue.cpp
@@ -0,0 +1,1303 @@
+#ifndef GSL_DEF
+#define GSL_DEF
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_spline.h>
+#endif
+#include "residue.h"
+#include "macros.h"
+#include <cmath>
+#include <stdio.h>
+#include "timing.hpp"
+
+double maxGradPosition(const double* y, const size_t nt, 
+	               const size_t nvar, const double* x, size_t nPts){
+	double maxGradT=0.0e0;
+	double gradT=0.0e0;
+	double pos=0.0e0;
+	size_t j,jm;
+	for (size_t i = 1; i <nPts; i++) {
+		j=i*nvar+nt;
+		jm=(i-1)*nvar+nt;
+		gradT=fabs((y[j]-y[jm])/(x[i]-x[i-1]));
+		if (gradT>=maxGradT) {
+			maxGradT=gradT;
+			pos=x[i];
+		}
+	}
+	return(pos);
+}
+
+int maxGradIndex(const double* y, const size_t nt, 
+	         const size_t nvar, const double* x, size_t nPts){
+	double maxGradT=0.0e0;
+	double gradT=0.0e0;
+	int pos=0.0e0;
+	size_t j,jm;
+	for (size_t i = 1; i <nPts; i++) {
+		j=i*nvar+nt;
+		jm=(i-1)*nvar+nt;
+		gradT=fabs((y[j]-y[jm])/(x[i]-x[i-1]));
+		if (gradT>=maxGradT) {
+			maxGradT=gradT;
+			pos=i;
+		}
+	}
+	return(pos);
+}
+
+double maxCurvPosition(const double* y, const size_t nt, 
+	               const size_t nvar, const double* x, size_t nPts){
+	double maxCurvT=0.0e0;
+	double gradTp=0.0e0;
+	double gradTm=0.0e0;
+	double curvT=0.0e0;
+	double dx=0.0e0;
+	double pos=0.0e0;
+	size_t j,jm,jp;
+	for (size_t i = 1; i <nPts-1; i++) {
+		j=i*nvar+nt;
+		jm=(i-1)*nvar+nt;
+		jp=(i+1)*nvar+nt;
+		gradTp=fabs((y[jp]-y[j])/(x[i+1]-x[i]));
+		gradTm=fabs((y[j]-y[jm])/(x[i]-x[i-1]));
+		dx=0.5e0*((x[i]+x[i+1])-(x[i-1]+x[i]));
+		curvT=(gradTp-gradTm)/dx;
+		if (curvT>=maxCurvT) {
+			maxCurvT=curvT;
+			pos=x[i];
+		}
+	}
+	return(pos);
+}
+
+int maxCurvIndex(const double* y, const size_t nt, 
+	               const size_t nvar, const double* x, size_t nPts){
+	double maxCurvT=0.0e0;
+	double gradTp=0.0e0;
+	double gradTm=0.0e0;
+	double curvT=0.0e0;
+	double dx=0.0e0;
+	int pos=0;
+	size_t j,jm,jp;
+	for (size_t i = 1; i <nPts-1; i++) {
+		j=i*nvar+nt;
+		jm=(i-1)*nvar+nt;
+		jp=(i+1)*nvar+nt;
+		gradTp=fabs((y[jp]-y[j])/(x[i+1]-x[i]));
+		gradTm=fabs((y[j]-y[jm])/(x[i]-x[i-1]));
+		dx=0.5e0*((x[i]+x[i+1])-(x[i-1]+x[i]));
+		curvT=(gradTp-gradTm)/dx;
+		if (curvT>=maxCurvT) {
+			maxCurvT=curvT;
+			pos=i;
+		}
+	}
+	return(pos);
+}
+
+double isothermPosition(const double* y, const double T, const size_t nt, 
+	                const size_t nvar, const double* x, const size_t nPts){
+	double pos=x[nPts-1];
+	size_t j;
+	for (size_t i = 1; i <nPts; i++) {
+		j=i*nvar+nt;
+		if (y[j]<=T) {
+			pos=x[i];
+			break;
+		}
+	}
+	return(pos);
+}
+
+void updateSolution(double* y, double* ydot, const size_t nvar,
+		    const double xOld[],const double xNew[],const size_t nPts){
+
+	double ytemp[nPts],ydottemp[nPts];
+
+	gsl_interp_accel* acc;
+	gsl_spline* spline;
+	acc = gsl_interp_accel_alloc();
+	spline = gsl_spline_alloc(gsl_interp_steffen, nPts);
+
+	gsl_interp_accel* accdot;
+	gsl_spline* splinedot;
+	accdot    = gsl_interp_accel_alloc();
+	splinedot = gsl_spline_alloc(gsl_interp_steffen, nPts);
+
+	for (size_t j = 0; j < nvar; j++) {
+
+		for (size_t i = 0; i < nPts; i++) {
+			ytemp[i]=y[j+i*nvar];
+			ydottemp[i]=ydot[j+i*nvar];
+		}
+
+		gsl_spline_init(spline,xOld,ytemp,nPts);
+		gsl_spline_init(splinedot,xOld,ydottemp,nPts);
+
+		for (size_t i = 0; i < nPts; i++) {
+
+		      	y[j+i*nvar]=gsl_spline_eval(spline,xNew[i],acc);
+		      	ydot[j+i*nvar]=gsl_spline_eval(splinedot,xNew[i],accdot);
+		}
+	}
+
+	//Exploring "fixing" boundary conditions:
+	//for (size_t j = 1; j < nvar; j++) {
+	//	//printf("%15.6e\t%15.6e\n", y[j],y[j+nvar]);
+	//	y[j]=y[j+nvar];
+	//	//y[j+(nPts-1)*nvar]=y[j+(nPts-2)*nvar];
+	//	//ydot[j+nvar]=ydot[j];
+	//}
+	//y[0]=0.0e0;
+	
+	gsl_interp_accel_free(acc);
+	gsl_spline_free(spline);
+
+	gsl_interp_accel_free(accdot);
+	gsl_spline_free(splinedot);
+}
+
+//Locate bath gas:
+size_t BathGasIndex(UserData data){
+	size_t index=0;
+	double max1;
+	double max=data->gas->massFraction(0);
+	for(size_t k=1;k<data->nsp;k++){
+		max1=data->gas->massFraction(k);
+		if(max1>=max){
+			max=max1;
+			index=k;
+		}
+	}
+	return(index+1);
+}
+
+//Locate Oxidizer:
+size_t oxidizerIndex(UserData data){
+	size_t index=0;
+	for(size_t k=1;k<data->nsp;k++){
+		if(data->gas->speciesName(k-1)=="O2"){
+			index=k;
+		}
+	}
+	return(index);
+}
+
+//Locate OH:
+size_t OHIndex(UserData data){
+	size_t index=0;
+	for(size_t k=1;k<data->nsp;k++){
+		if(data->gas->speciesName(k-1)=="OH"){
+			index=k;
+		}
+	}
+	return(index);
+}
+
+//Locate HO2:
+size_t HO2Index(UserData data){
+	size_t index=0;
+	for(size_t k=1;k<data->nsp;k++){
+		if(data->gas->speciesName(k-1)=="HO2"){
+			index=k;
+		}
+	}
+	return(index);
+}
+
+int setAlgebraicVariables(N_Vector* id, UserData data){
+	double *iddata;
+  	N_VConst(ONE, *id);
+  	iddata   = N_VGetArrayPointer_OpenMP(*id);
+	data->k_bath=BathGasIndex(data);
+	data->k_oxidizer=oxidizerIndex(data);
+	data->k_OH=OHIndex(data);
+	data->k_HO2=HO2Index(data);
+	printf("Oxidizer index: %lu\n",data->k_oxidizer);
+	printf("Bath gas index:%lu\n",data->k_bath);
+	for (size_t i = 1; i <=data->npts; i++) {
+		/*Algebraic variables: indicated by ZERO.*/
+		Rid(i)=ZERO;
+		Pid(i)=ZERO;
+		Yid(i,data->k_bath)=ZERO;
+	}	
+	Rid(1)=ONE;
+	Tid(1)=ZERO;
+	Tid(data->npts)=ZERO;
+	if(data->constantPressure){
+		Pid(data->npts)=ONE;
+	}
+	else{
+		Pid(data->npts)=ZERO;
+		Rid(data->npts)=ONE;
+	}
+	for (size_t k = 1; k <=data->nsp; k++) {
+		Yid(1,k)=ZERO;
+		Yid(data->npts,k)=ZERO;
+	}
+	return(0);
+}
+
+inline double calc_area(double x,int* i){
+	switch (*i) {
+		case 0:
+			return(ONE);
+		case 1:
+			return(x);
+		case 2:
+			return(x*x);
+		default:
+			return(ONE);
+	}
+}
+
+void readInitialCondition(FILE* input, double* ydata, const size_t nvar, const size_t nr, const size_t nPts){
+
+	FILE* output;output=fopen("test.dat","w");
+
+	size_t bufLen=10000;
+	size_t nRows=0;
+	size_t nColumns=nvar;
+
+	char buf[bufLen];
+	char buf1[bufLen];
+	char comment[1];
+	char *ret;
+
+	while (fgets(buf,bufLen, input)!=NULL){
+		comment[0]=buf[0];
+		if(strncmp(comment,"#",1)!=0){
+			nRows++;
+		}
+	}
+	rewind(input);
+
+	printf("nRows: %ld\n", nRows);
+	double y[nRows*nColumns];
+
+	size_t i=0;
+	while (fgets(buf,bufLen, input)!=NULL){
+		comment[0]=buf[0];
+		if(strncmp(comment,"#",1)==0){
+		}
+		else{
+			ret=strtok(buf,"\t");
+			size_t j=0;
+			y[i*nColumns+j]=(double)(atof(ret));
+			j++;
+			while(ret!=NULL){
+				ret=strtok(NULL,"\t");
+				if(j<nColumns){
+					y[i*nColumns+j]=(double)(atof(ret));
+				}
+				j++;
+			}
+			i++;
+		}
+	}
+
+	for (i = 0; i < nRows; i++) {
+		for (size_t j = 0; j < nColumns; j++) {
+			fprintf(output, "%15.6e\t",y[i*nColumns+j]);
+		}
+		fprintf(output, "\n");
+	}
+	fclose(output);
+
+	double xOld[nRows],xNew[nPts],ytemp[nPts];
+
+	for (size_t j = 0; j < nRows; j++) {
+		xOld[j]=y[j*nColumns+nr];
+	}
+
+	double dx=xOld[nRows-1]/((double)(nPts)-1.0e0);
+	for (size_t j = 0; j < nPts; j++) {
+		xNew[j]=(double)j*dx;
+	}
+
+	gsl_interp_accel* acc;
+	gsl_spline* spline;
+	acc = gsl_interp_accel_alloc();
+	spline = gsl_spline_alloc(gsl_interp_steffen, nRows);
+
+	for (size_t j = 0; j < nColumns; j++) {
+
+		for (size_t k = 0; k < nRows; k++) {
+			ytemp[k]=y[j+k*nColumns];
+		}
+
+		gsl_spline_init(spline,xOld,ytemp,nRows);
+
+		for (size_t k = 1; k < nPts-1; k++) {
+		      	ydata[j+k*nColumns]=gsl_spline_eval(spline,xNew[k],acc);
+		}
+	}
+	
+	gsl_interp_accel_free(acc);
+	gsl_spline_free(spline);
+
+}
+
+double systemMass(double* ydata,UserData data){
+	double mass=0.0e0;
+	double rho;
+	for (size_t i = 2; i <=data->npts; i++) {
+  		data->gas->setState_TPY(T(i), P(i), &Y(i,1));   
+		rho=data->gas->density();                                           
+		//psi(i)=psi(i-1)+rho*(R(i)-R(i-1))*calc_area(HALF*(R(i)+R(i-1)),&m);
+		mass+=rho*(R(i)-R(i-1))*calc_area(R(i),&data->metric);
+	}
+	return(mass);
+}
+
+
+int initializePsiGrid(double* ydata, double* psidata, UserData data){
+
+	double rho;
+	/*Create a psi grid that corresponds CONSISTENTLY to the spatial grid
+	 * "R" created above. Note that the Lagrangian variable psi has  units
+	 * of kg. */
+	psi(1)=ZERO;
+	for (size_t i = 2; i <=data->npts; i++) {
+  		data->gas->setState_TPY(T(i), P(i), &Y(i,1));   
+		rho=data->gas->density();                                           
+		//psi(i)=psi(i-1)+rho*(R(i)-R(i-1))*calc_area(HALF*(R(i)+R(i-1)),&data->metric);
+		psi(i)=psi(i-1)+rho*(R(i)-R(i-1))*calc_area(R(i),&data->metric);
+	}
+
+	/*The mass of the entire system is the value of psi at the last grid
+	 * point. Normalize psi by this mass so that it varies from zero to
+	 * one. This makes psi dimensionless. So the mass needs to be
+	 * multiplied back in the approporiate places in the governing
+	 * equations so that units match.*/
+	data->mass=psi(data->npts);
+	for (size_t i = 1; i <=data->npts; i++) {
+		psi(i)=psi(i)/data->mass;
+	}
+	return(0);
+}
+
+
+int setInitialCondition(N_Vector* y, 
+	    		N_Vector* ydot, 
+            		UserData data){
+
+	double* ydata;
+	double* ydotdata;
+	double* psidata;
+	double* innerMassFractionsData;
+	double f=ZERO;
+	double g=ZERO;
+
+	double perturb,rho;
+	double epsilon=ZERO;
+	int m,ier;
+  	ydata    = N_VGetArrayPointer_OpenMP(*y);
+  	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);
+	innerMassFractionsData =  data->innerMassFractions;
+
+	if(data->adaptiveGrid){
+  		psidata  = data->grid->xOld;
+	}
+	else{
+  		psidata  = data->uniformGrid;
+	}
+
+	m=data->metric;
+
+	data->innerTemperature=data->initialTemperature;
+	for (size_t k = 1; k <=data->nsp; k++) {
+  		innerMassFractionsData[k-1]=data->gas->massFraction(k-1);
+	}
+
+	//Define Grid:
+	//double Rmin=1e-03*data->domainLength;
+	double Rmin=0.0e0;
+	double dR=(data->domainLength-Rmin)/((double)(data->npts)-1.0e0);
+	double dv=(pow(data->domainLength,1+data->metric)-pow(data->firstRadius*data->domainLength,1+data->metric))/((double)(data->npts)-1.0e0);
+	for (size_t i = 1; i <=data->npts; i++) {
+		if(data->metric==0){
+			R(i)=Rmin+(double)((i-1)*dR);
+		}else{
+			if(i==1){
+				R(i)=ZERO;
+			}else if(i==2){
+				R(i)=data->firstRadius*data->domainLength;
+			}else{
+				R(i)=pow(pow(R(i-1),1+data->metric)+dv,1.0/((double)(1+data->metric)));
+			}
+		}
+		T(i)=data->initialTemperature;
+		for (size_t k = 1; k <=data->nsp; k++) {
+  	        	Y(i,k)=data->gas->massFraction(k-1);	//Indexing different in Cantera
+		}
+  	        P(i)=data->initialPressure*Cantera::OneAtm;
+	}	
+	R(data->npts)=data->domainLength;
+
+	double Tmax;
+	double Tmin=data->initialTemperature;
+	double w=data->mixingWidth;
+	double YN2=ZERO;
+	double YO2=ZERO;
+	double YFuel,YOxidizer,sum;
+
+	//if(data->problemType==0){
+	//	data->gas->equilibrate("HP");
+	//	data->maxTemperature=data->gas->temperature();
+	//}
+	//else if(data->problemType==1){
+	//	/*Premixed Combustion: Equilibrium products comprise ignition
+	//	 * kernel at t=0. The width of the kernel is "mixingWidth"
+	//	 * shifted by "shift" from the center.*/
+	//	data->gas->equilibrate("HP");
+	//	Tmax=data->gas->temperature();
+	//	for (size_t i = 1; i <=data->npts; i++) {
+	//		g=HALF*(tanh((R(i)-data->shift)/w)+ONE);	//increasing function of x
+	//		f=ONE-g;					//decreasing function of x
+	//		T(i)=(Tmax-Tmin)*f+Tmin;
+	//		for (size_t k = 1; k <=data->nsp; k++) {
+	//			Y(i,k)=(data->gas->massFraction(k-1)-Y(i,k))*f+Y(i,k);
+	//		}
+	//	}
+	//	if(data->dirichletOuter){
+	//		T(data->npts)=data->wallTemperature;
+	//	}
+	//}
+	//else if(data->problemType==2){
+	//	FILE* input;
+	//	if(input=fopen("initialCondition.dat","r")){
+	//		readInitialCondition(input, ydata, data->nvar, data->nr, data->npts);
+	//		fclose(input);
+	//	}
+	//	else{
+	//		printf("file initialCondition.dat not found!\n");
+	//		return(-1);
+	//	}
+	//}
+
+	initializePsiGrid(ydata,psidata,data);
+
+	if(data->adaptiveGrid){
+	//	if(data->problemType!=0){
+	//		data->grid->position=maxGradPosition(ydata, data->nt, data->nvar,
+	//						 data->grid->xOld, data->npts);
+	//		//data->grid->position=maxCurvPosition(ydata, data->nt, data->nvar,
+	//		//				 data->grid->xOld, data->npts);
+	//	}
+	//	else{
+	//	}
+		if(data->problemType!=3){
+			data->grid->position=0.0e0;
+			double x=data->grid->position+data->gridOffset*data->grid->leastMove;
+			printf("New grid center:%15.6e\n",x);
+			ier=reGrid(data->grid, x);
+			if(ier==-1)return(-1);
+			updateSolution(ydata, ydotdata, data->nvar,
+			               data->grid->xOld,data->grid->x,data->npts);
+			storeGrid(data->grid->x,data->grid->xOld,data->npts);
+		}
+	}
+	else{
+		double psiNew[data->npts];
+		double dpsi=1.0e0/((double)(data->npts)-1.0e0);
+		for (size_t i = 0; i < data->npts; i++) {
+			psiNew[i]=(double)(i)*dpsi;
+		}
+		printf("Last point:%15.6e\n",psiNew[data->npts-1]);
+		updateSolution(ydata, ydotdata, data->nvar,
+		               data->uniformGrid,psiNew,data->npts);
+		storeGrid(psiNew,data->uniformGrid,data->npts);
+	}
+
+	if(data->problemType==0){
+		data->gas->equilibrate("HP");
+		data->maxTemperature=data->gas->temperature();
+	}
+	else if(data->problemType==1){
+		/*Premixed Combustion: Equilibrium products comprise ignition
+		 * kernel at t=0. The width of the kernel is "mixingWidth"
+		 * shifted by "shift" from the center.*/
+		data->gas->equilibrate("HP");
+		Tmax=data->gas->temperature();
+		for (size_t i = 1; i <=data->npts; i++) {
+			g=HALF*(tanh((R(i)-data->shift)/w)+ONE);	//increasing function of x
+			f=ONE-g;					//decreasing function of x
+			T(i)=(Tmax-Tmin)*f+Tmin;
+			for (size_t k = 1; k <=data->nsp; k++) {
+				Y(i,k)=(data->gas->massFraction(k-1)-Y(i,k))*f+Y(i,k);
+			}
+		}
+		if(data->dirichletOuter){
+			T(data->npts)=data->wallTemperature;
+		}
+	}
+	else if(data->problemType==2){
+		FILE* input;
+		if(input=fopen("initialCondition.dat","r")){
+			readInitialCondition(input, ydata, data->nvar, data->nr, data->npts);
+			fclose(input);
+		}
+		else{
+			printf("file initialCondition.dat not found!\n");
+			return(-1);
+		}
+	}
+	else if(data->problemType==3){
+		FILE* input;
+		if(input=fopen("restart.bin","r")){
+			readRestart(y, ydot, input, data);
+			fclose(input);
+			printf("Restart solution loaded!\n");
+			printf("Problem starting at t=%15.6e\n",data->tNow);
+			return(0);
+		}
+		else{
+			printf("file restart.bin not found!\n");
+			return(-1);
+		}
+	}
+
+	if(data->reflectProblem){
+		double temp;
+		int j=1;
+		while (data->npts+1-2*j>=0) {
+			temp=T(j);
+			T(j)=T(data->npts+1-j);
+			T(data->npts+1-j)=temp;
+			for (size_t k = 1; k <=data->nsp; k++) {
+				temp=Y(j,k);
+				Y(j,k)=Y(data->npts+1-j,k);
+				Y(data->npts+1-j,k)=temp;
+			}
+			j=j+1;
+		}
+	}
+
+	/*Floor small values to zero*/
+	for (size_t i = 1; i <=data->npts; i++) {
+		for (size_t k = 1; k <=data->nsp; k++) {
+			if(fabs(Y(i,k))<=data->massFractionTolerance){
+				Y(i,k)=0.0e0;
+			}
+		}
+	}
+
+	//Set grid to location of maximum curvature:
+	if(data->adaptiveGrid){
+		data->grid->position=maxCurvPosition(ydata, data->nt, data->nvar,
+				  		data->grid->x, data->npts);
+		ier=reGrid(data->grid, data->grid->position);
+		updateSolution(ydata, ydotdata, data->nvar,
+		               data->grid->xOld,data->grid->x,data->npts);
+		storeGrid(data->grid->x,data->grid->xOld,data->npts);
+	}
+
+	/*Ensure consistent boundary conditions*/	
+	T(1)=T(2);
+	for (size_t k = 1; k <=data->nsp; k++) {
+		Y(1,k)=Y(2,k);
+		Y(data->npts,k)=Y(data->npts-1,k);
+	}
+
+
+	return(0);
+}
+
+inline double Qdot(double* t, 
+	    	   double* x, 
+	    	   double* ignTime, 
+	    	   double* kernelSize, 
+	    	   double* maxQdot){
+	double qdot;
+	if(*x<=*kernelSize){
+		if((*t)<=(*ignTime)){
+			qdot=(*maxQdot);
+		}
+		else{
+			qdot=0.0e0;
+		}
+	}else{
+		qdot=0.0e0;
+	}
+	return(qdot);
+}
+
+inline void setGas(UserData data, 
+	    double *ydata, 
+	    size_t gridPoint){
+  	data->gas->setTemperature(T(gridPoint)); 	        
+  	data->gas->setMassFractions_NoNorm(&Y(gridPoint,1)); 	
+  	data->gas->setPressure(P(gridPoint));
+}
+
+void getTransport(UserData data, 
+		  double *ydata, 
+		  size_t gridPoint,
+		  double *rho,
+		  double *lambda,
+		  double YV[]){
+
+	double YAvg[data->nsp],
+		 XLeft[data->nsp],
+		 XRight[data->nsp],
+		 gradX[data->nsp];
+
+	setGas(data,ydata,gridPoint);
+	data->gas->getMoleFractions(XLeft);                                     
+	setGas(data,ydata,gridPoint+1);
+	data->gas->getMoleFractions(XRight);                                     
+                                                                             
+	for (size_t k = 1; k <=data->nsp; k++) {                                          
+		YAvg(k)=HALF*(Y(gridPoint,k)+
+			      Y(gridPoint+1,k));                                 
+		gradX(k)=(XRight(k)-XLeft(k))/
+			 (R(gridPoint+1)-R(gridPoint));                       
+	}                                                                    
+	double TAvg = HALF*(T(gridPoint)+T(gridPoint+1));
+	double gradT=(T(gridPoint+1)-T(gridPoint))/
+		       (R(gridPoint+1)-R(gridPoint));                                  
+
+
+  	data->gas->setTemperature(TAvg); 	        
+  	data->gas->setMassFractions_NoNorm(YAvg); 	
+  	data->gas->setPressure(P(gridPoint));
+
+	*rho=data->gas->density();                                       
+	*lambda=data->trmix->thermalConductivity();                      
+	data->trmix->getSpeciesFluxes(1,&gradT,data->nsp,                     
+				      gradX,data->nsp,YV);                
+	//setGas(data,ydata,gridPoint);
+}
+
+int residue(double t, 
+	    N_Vector y, 
+	    N_Vector ydot, 
+	    N_Vector res, 
+            void *user_data){
+
+	/*Declare and fetch nvectors and user data:*/
+
+  	double *ydata, *ydotdata, *resdata, *psidata, *innerMassFractionsData;
+  	
+	UserData data;
+	data = (UserData)user_data;
+	size_t npts=data->npts;
+	size_t nsp=data->nsp;
+	size_t k_bath = data->k_bath;
+
+	ydata 	= N_VGetArrayPointer_OpenMP(y); 
+	ydotdata= N_VGetArrayPointer_OpenMP(ydot); 
+	resdata = N_VGetArrayPointer_OpenMP(res);
+	if(data->adaptiveGrid==1){
+		psidata = data->grid->x;
+	}else{
+		psidata = data->uniformGrid;
+	}
+
+	innerMassFractionsData =  data->innerMassFractions;
+
+	/* Grid stencil:*/
+
+	/*-------|---------*---------|---------*---------|-------*/
+	/*-------|---------*---------|---------*---------|-------*/
+	/*-------|---------*---------|---------*---------|-------*/
+	/*-------m-------mhalf-------j-------phalf-------p-------*/
+	/*-------|---------*---------|---------*---------|-------*/
+	/*-------|---------*---------|---------*---------|-------*/
+	/*-------|<=======dxm=======>|<=======dxp=======>|-------*/
+	/*-------|---------*<======dxav=======>*---------|-------*/
+	/*-------|<================dxpm=================>|-------*/
+
+	/* Various variables defined for book-keeping and storing previously
+	 * calculated values:
+	 * rho		: densities at points  m, mhalf, j, p, and phalf.
+	 * area		: the matric at points m, mhalf, j, p, and phalf.
+	 * m 		: exponent that determines geometry;
+	 * lambda	: thermal conductivities at mhalf and phalf.
+	 * mdot		: mass flow rate at m, j, and p.
+	 * X		: mole fractions at j and p.
+	 * YV		: diffusion fluxes at mhalf and phalf.
+	 * Tgrad	: temperature gradient at mhalf and phalf.
+	 * Tav		: average temperature between two points.
+	 * Pav		: average pressure between two points.
+	 * Yav		: average mass fractions between two points.
+	 * Xgradhalf	: mole fraction gradient at j.
+	 * Cpb		: mass based bulk specific heat.
+	 * tranTerm	: transient terms.
+	 * advTerm	: advection terms.
+	 * diffTerm	: diffusion terms.
+	 * srcTerm	: source terms.
+	 */
+
+	double rhomhalf,  rhom, lambdamhalf, YVmhalf[nsp], 
+		 rho, 
+		 rhophalf, lambdaphalf, YVphalf[nsp],
+		 Cpb, Cvb,  Cp[nsp],   wdot[nsp],   enthalpy[nsp], energy[nsp],
+		 tranTerm, diffTerm,    srcTerm, advTerm,
+		 area,areamhalf,areaphalf,aream,areamhalfsq,areaphalfsq;
+	
+	/*Aliases for difference coefficients:*/
+	double cendfm, cendfc, cendfp;
+	/*Aliases for various grid spacings:*/
+	double dpsip, dpsiav, dpsipm, dpsim, dpsimm;
+	dpsip=dpsiav=dpsipm=dpsim=dpsimm=ONE;
+	double mass, mdotIn;
+	double sum, sum1, sum2, sum3;
+
+	size_t j,k;
+	int m;
+	m=data->metric;					//Unitless
+	mass=data->mass;				//Units: kg
+	mdotIn=data->mdot*calc_area(R(npts),&m);	//Units: kg/s
+
+//	/*evaluate properties at j=1*************************/
+	setGas(data,ydata,1);
+	rhom=data->gas->density();                                       
+        Cpb=data->gas->cp_mass();       	//J/kg/K
+        Cvb=data->gas->cv_mass();       	//J/kg/K
+	aream= calc_area(R(1),&m);	    
+
+	/*******************************************************************/
+	/*Calculate values at j=2's m and mhalf*****************************/
+
+	getTransport(data, ydata, 1, &rhomhalf,&lambdamhalf,YVmhalf);
+	areamhalf= calc_area(HALF*(R(1)+R(2)),&m);                                 
+	areamhalfsq= areamhalf*areamhalf;                                   
+                                                                             
+	/*******************************************************************/
+
+	/*Fill up res with left side (center) boundary conditions:**********/
+	/*We impose zero fluxes at the center:*/                        
+
+	/*Mass:*/                                                            
+	Rres(1)=Rdot(1);
+
+	/*Energy:*/                                                     
+
+	if(data->dirichletInner){
+		//Tres(1)=Tdot(1);
+		Tres(1)=T(1)-data->innerTemperature;
+	}
+	else{
+		Tres(1)=T(2)-T(1);
+		//Tres(1)=Tdot(1) - (Pdot(1)/(rhom*Cpb))
+		//	+(double)(data->metric+1)*(rhomhalf*lambdamhalf*areamhalfsq*(T(2)-T(1))/psi(2)-psi(1));
+	}
+                                                                        
+	/*Species:*/                                                    
+	sum=ZERO;
+	for (k = 1; k <=nsp; k++) {                                        
+		if(k!=k_bath){
+			if(fabs(mdotIn)>1e-14){
+		      		Yres(1,k)=innerMassFractionsData[k-1]-
+					  Y(1,k)-
+					  (YVmhalf(k)*areamhalf)/mdotIn;
+			}
+			else{
+		      		//Yres(1,k)=Y(1,k)-innerMassFractionsData[k-1];
+		      		Yres(1,k)=Y(2,k)-Y(1,k);
+				/*The below flux boundary condition makes the
+				 * problem more prone to diverge. How does one
+				 * fix this?*/
+		      		//Yres(1,k)=YVmhalf(k);
+			}
+		      sum=sum+Y(1,k);
+		}
+	}                                                                 
+	Yres(1,k_bath)=ONE-sum-Y(1,k_bath);                                     
+
+                                                                             
+	/*Pressure:*/                                                        
+	Pres(1)=P(2)-P(1);
+
+	/*Fill up res with governing equations at inner points:*************/
+	for (j = 2; j < npts; j++) {                                         
+                                                                             
+		/*evaluate various mesh differences*///                      
+        	dpsip =        (psi(j+1) - psi(j)  )*mass;                            
+        	dpsim =        (psi(j)   - psi(j-1))*mass;                            
+        	dpsiav =  HALF*(psi(j+1) - psi(j-1))*mass;                            
+        	dpsipm =       (psi(j+1) - psi(j-1))*mass;                            
+		/***********************************///                      
+                                                                             
+		/*evaluate various central difference coefficients*/       
+        	cendfm = - dpsip / (dpsim*dpsipm);                               
+        	cendfc =   (dpsip-dpsim) / (dpsip*dpsim);                         
+        	cendfp =   dpsim / (dpsip*dpsipm);                               
+		/**************************************************/
+                                                                   
+                                                                  
+		/*evaluate properties at j*************************/
+		setGas(data,ydata,j);
+		rho=data->gas->density();		//kg/m^3
+        	Cpb=data->gas->cp_mass();       	//J/kg/K
+        	Cvb=data->gas->cv_mass();       	//J/kg/K
+		data->gas->getNetProductionRates(wdot); //kmol/m^3
+		data->gas->getEnthalpy_RT(enthalpy);	//unitless
+		data->gas->getCp_R(Cp);			//unitless
+		area = calc_area(R(j),&m);	        //m^2
+                                                                           
+		/*evaluate properties at p*************************/
+		getTransport(data, ydata, j, &rhophalf,&lambdaphalf,YVphalf);
+		areaphalf= calc_area(HALF*(R(j)+R(j+1)),&m);	    
+		areaphalfsq= areaphalf*areaphalf;		    
+		/**************************************************///     
+
+		/*Mass:*/
+		/* ∂r/∂ψ = 1/ρA */
+		Rres(j)=((R(j)-R(j-1))/dpsim)-(TWO/(rhom*aream+rho*area));
+
+		/*Energy:*/
+		/* ∂T/∂t = - ṁ(∂T/∂ψ) 
+		 * 	   + (1/cₚ)(∂/∂ψ)(λρA²∂T/∂ψ) 
+		 * 	   - (A/cₚ) ∑ YᵢVᵢcₚᵢ(∂T/∂ψ) 
+		 * 	   - (1/ρcₚ)∑ ώᵢhᵢ 
+		 * 	   + (1/ρcₚ)(∂P/∂t) */
+		/*Notes: 
+		 * λ has units J/m/s/K.
+		 * YᵢVᵢ has units kg/m^2/s.
+		 * 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.
+		 * cₚᵢ has units J/kg/K, so we must multiply the specific heat
+		 * defined above (getCp_R) by the gas constant (J/kmol/K) and
+		 * divide by the molecular weight (kg/kmol) to get the right
+		 * units.
+		 * */
+
+		//enthalpy formulation:
+		tranTerm = Tdot(j) - (Pdot(j)/(rho*Cpb));
+		sum=ZERO;
+		sum1=ZERO;
+		for (k = 1; k <=nsp; k++) {
+			sum=sum+wdot(k)*enthalpy(k);
+			sum1=sum1+(Cp(k)/data->gas->molecularWeight(k-1))*rho
+			     *HALF*(YVmhalf(k)+YVphalf(k));
+		}
+		sum=sum*Cantera::GasConstant*T(j);
+		sum1=sum1*Cantera::GasConstant;
+		diffTerm  =(( (rhophalf*areaphalfsq*lambdaphalf*(T(j+1)-T(j))/dpsip)
+			      -(rhomhalf*areamhalfsq*lambdamhalf*(T(j)-T(j-1))/dpsim) )
+			     /(dpsiav*Cpb) )
+			  -(sum1*area*(cendfp*T(j+1)
+			              +cendfc*T(j)
+			              +cendfm*T(j-1))/Cpb);
+		srcTerm   = (sum-Qdot(&t,&R(j),&data->ignTime,&data->kernelSize,&data->maxQDot))/(rho*Cpb);
+		advTerm   = (mdotIn*(T(j)-T(j-1))/dpsim);
+		Tres(j)= tranTerm
+			+advTerm 
+			-diffTerm
+			+srcTerm;
+
+	//	//energy formulation:
+	//	tranTerm = Tdot(j);
+	//	sum=ZERO;
+	//	sum1=ZERO;
+	//	sum2=ZERO;
+	//	sum3=ZERO;
+	//	for (k = 1; k <=nsp; k++) {
+	//		energy(k)=enthalpy(k)-ONE;
+	//		sum=sum+wdot(k)*energy(k);
+	//		sum1=sum1+(Cp(k)/data->gas->molecularWeight(k-1))*rho
+	//		     *HALF*(YVmhalf(k)+YVphalf(k));
+	//		sum2=sum2+(YVmhalf(k)/data->gas->molecularWeight(k-1));
+	//		sum3=sum3+(YVphalf(k)/data->gas->molecularWeight(k-1));
+	//	}
+	//	sum=sum*Cantera::GasConstant*T(j);
+	//	sum1=sum1*Cantera::GasConstant;
+	//	diffTerm  =(( (rhophalf*areaphalfsq*lambdaphalf*(T(j+1)-T(j))/dpsip)
+	//		      -(rhomhalf*areamhalfsq*lambdamhalf*(T(j)-T(j-1))/dpsim) )
+	//		     /(dpsiav*Cvb) )
+	//		  -(sum1*area*(cendfp*T(j+1)
+	//		              +cendfc*T(j)
+	//		              +cendfm*T(j-1))/Cvb);
+	//	srcTerm   = (sum-Qdot(&t,&R(j),&data->ignTime,&data->kernelSize,&data->maxQDot))/(rho*Cvb);
+	//	advTerm   = (mdotIn*(T(j)-T(j-1))/dpsim);
+	//	advTerm = advTerm + (Cantera::GasConstant*T(j)*area/Cvb)*((sum3-sum2)/dpsiav);
+	//	Tres(j)= tranTerm
+	//		+advTerm 
+	//		-diffTerm
+	//		+srcTerm;
+
+		/*Species:*/
+		/* ∂Yᵢ/∂t = - ṁ(∂Yᵢ/∂ψ) 
+		 * 	    - (∂/∂ψ)(AYᵢVᵢ) 
+		 * 	    + (ώᵢWᵢ/ρ)  */
+
+		sum=ZERO;
+		for (k = 1; k <=nsp; k++) {
+			if(k!=k_bath){
+				tranTerm  = Ydot(j,k);
+				diffTerm  = (YVphalf(k)*areaphalf
+					    -YVmhalf(k)*areamhalf)/dpsiav;
+				srcTerm   = wdot(k)
+			     		    *(data->gas->molecularWeight(k-1))/rho;
+				advTerm   = (mdotIn*(Y(j,k)-Y(j-1,k))/dpsim);
+
+				Yres(j,k)= tranTerm 
+					  +advTerm
+					  +diffTerm
+					  -srcTerm;
+
+		      		sum=sum+Y(j,k);
+			}
+		}
+		Yres(j,k_bath)=ONE-sum-Y(j,k_bath);
+
+		/*Pressure:*/
+		Pres(j) = P(j+1)-P(j);
+
+		/*Assign values evaluated at p and phalf to m       
+		 * and mhalf to save some cpu cost:****************/
+		areamhalf=areaphalf;                                
+		areamhalfsq=areaphalfsq;                            
+		aream=area;
+		rhom=rho;
+		rhomhalf=rhophalf;                                  
+		lambdamhalf=lambdaphalf;                            
+		for (k = 1; k <=nsp; k++) {                         
+			YVmhalf(k)=YVphalf(k);                      
+		}                                                   
+		/**************************************************/
+		                                                    
+	}                                                           
+	/*******************************************************************///
+
+
+	/*Fill up res with right side (wall) boundary conditions:***********/
+	/*We impose zero fluxes at the wall:*/                               
+	setGas(data,ydata,npts);
+	rho=data->gas->density();                               
+	area = calc_area(R(npts),&m);				   
+
+	/*Mass:*/                                                          
+	dpsim=(psi(npts)-psi(npts-1))*mass;
+	Rres(npts)=((R(npts)-R(npts-1))/dpsim)-(TWO/(rhom*aream+rho*area));
+
+	/*Energy:*/                                             
+	if(data->dirichletOuter){
+		Tres(npts)=T(npts)-data->wallTemperature;
+	}
+	else{
+		Tres(npts)=T(npts)-T(npts-1);
+	}
+                                                                
+	/*Species:*/                                            
+	sum=ZERO;
+	for (k = 1; k <=nsp; k++) {                             
+		if(k!=k_bath){
+		    	//Yres(npts,k)=Y(npts,k)-Y(npts-1,k);
+		    	Yres(npts,k)=YVmhalf(k);         
+		    	sum=sum+Y(npts,k);
+		}
+	}                                                       
+	Yres(npts,k_bath)=ONE-sum-Y(npts,k_bath);             
+                                                                
+	
+	/*Pressure:*/                                                      
+	if(data->constantPressure){
+		Pres(npts)=Pdot(npts)-data->dPdt;                                     
+	}
+	else{
+		Pres(npts)=R(npts)-data->domainLength;                     
+		//Pres(npts)=Rdot(npts);                     
+	}
+
+
+	/*******************************************************************/
+
+	//for (j = 1; j <=npts; j++) {
+	//	//for (k = 1; k <=nsp; k++) {
+	//	//	Yres(j,k)=Ydot(j,k);
+	//	//}
+	//	//Tres(j)=Tdot(j);
+	//}
+
+	return(0);
+
+}
+
+void printSpaceTimeHeader(UserData data)
+{
+	fprintf((data->output), "%15s\t","#1");
+	for (size_t k = 1; k <=data->nvar+2; k++) {
+		fprintf((data->output), "%15lu\t",k+1);
+	}
+	fprintf((data->output), "\n");
+
+	fprintf((data->output), "%15s\t%15s\t%15s\t","#psi","time(s)","dpsi");
+	fprintf((data->output), "%15s\t%15s\t","radius(m)","Temp(K)");
+	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\n","Pressure(Pa)");
+}
+
+void printSpaceTimeOutput(double t, N_Vector* y, FILE* output, UserData data)
+{
+	double *ydata,*psidata;
+  	ydata    = N_VGetArrayPointer_OpenMP(*y);
+
+	if(data->adaptiveGrid){
+  		psidata  = data->grid->x;
+	}else{
+  		psidata  = data->uniformGrid;
+	}
+
+	for (size_t i = 0; i < data->npts; i++) {
+		fprintf(output, "%15.6e\t%15.6e\t",psi(i+1),t);
+		if(i==0){
+			fprintf(output, "%15.6e\t",psi(2)-psi(1));
+		}
+		else{
+			fprintf(output, "%15.6e\t",psi(i+1)-psi(i));
+		}
+		for (size_t j = 0; j < data->nvar; j++) {
+			fprintf(output, "%15.9e\t",ydata[j+i*data->nvar]);
+		}
+		fprintf(output, "\n");
+	}
+	fprintf(output, "\n\n");
+}
+
+void writeRestart(double t, N_Vector* y, N_Vector* ydot, FILE* output, UserData data){
+	double *ydata,*psidata, *ydotdata;
+	ydata = N_VGetArrayPointer_OpenMP(*y);
+	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);
+	if(data->adaptiveGrid){
+  		psidata  = data->grid->x;
+	}else{
+  		psidata  = data->uniformGrid;
+	}
+ 	fwrite(&t, sizeof(t), 1, output);		//write time
+ 	fwrite(psidata, data->npts*sizeof(psidata), 1, output);	//write grid
+ 	fwrite(ydata, data->neq*sizeof(ydata), 1, output);	//write solution
+ 	fwrite(ydotdata, data->neq*sizeof(ydotdata), 1, output);	//write solutiondot
+}
+
+void readRestart(N_Vector* y, N_Vector* ydot, FILE* input, UserData data){
+	double *ydata,*psidata, *ydotdata;
+	double t;
+	if(data->adaptiveGrid){
+  		psidata  = data->grid->x;
+	}else{
+  		psidata  = data->uniformGrid;
+	}
+	ydata = N_VGetArrayPointer_OpenMP(*y);
+	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);
+ 	fread(&t, sizeof(t), 1, input);
+	data->tNow=t;
+ 	fread(psidata, data->npts*sizeof(psidata), 1, input);
+ 	fread(ydata, data->neq*sizeof(ydata), 1, input);
+ 	fread(ydotdata, data->neq*sizeof(ydotdata), 1, input);
+	if(data->adaptiveGrid){
+		storeGrid(data->grid->x,data->grid->xOld,data->npts);
+	}
+}
+
+void printGlobalHeader(UserData data)
+{
+	fprintf((data->globalOutput), "%8s\t","#Time(s)");
+	//fprintf((data->globalOutput), "%15s","S_u(exp*)(m/s)");
+	fprintf((data->globalOutput), "%15s","bFlux(kg/m^2/s)");
+	fprintf((data->globalOutput), "%16s","  IsothermPos(m)");
+	fprintf((data->globalOutput), "%15s\t","Pressure(Pa)");
+	fprintf((data->globalOutput), "%15s\t","Pdot(Pa/s)");
+	fprintf((data->globalOutput), "%15s\t","gamma");
+	fprintf((data->globalOutput), "%15s\t","S_u(m/s)");
+	fprintf((data->globalOutput), "%15s\t","Tu(K)");
+	fprintf((data->globalOutput), "\n");
+}
+//
+//void printSpaceTimeRates(double t, N_Vector ydot, UserData data)
+//{
+//	double *ydotdata,*psidata;
+//  	ydotdata    = N_VGetArrayPointer_OpenMP(ydot);
+//  	psidata    = N_VGetArrayPointer_OpenMP(data->grid);
+//	for (int i = 0; i < data->npts; i++) {
+//		fprintf((data->ratesOutput), "%15.6e\t%15.6e\t",psi(i+1),t);
+//		for (int j = 0; j < data->nvar; j++) {
+//			fprintf((data->ratesOutput), "%15.6e\t",ydotdata[j+i*data->nvar]);
+//		}
+//		fprintf((data->ratesOutput), "\n");
+//	}
+//	fprintf((data->ratesOutput), "\n\n");
+//}
+//
+void printGlobalVariables(double t, N_Vector* y, N_Vector* ydot, UserData data)
+{
+	double *ydata,*ydotdata, *innerMassFractionsData, *psidata;
+	innerMassFractionsData =  data->innerMassFractions;
+
+	if(data->adaptiveGrid){
+  		psidata  = data->grid->x;
+	}else{
+  		psidata  = data->uniformGrid;
+	}
+
+	double TAvg, RAvg, YAvg, psiAvg;
+  	ydata    = N_VGetArrayPointer_OpenMP(*y);
+  	ydotdata = N_VGetArrayPointer_OpenMP(*ydot);
+	TAvg=data->isotherm;
+	double sum=ZERO;
+	double dpsim,area,aream,drdt;
+	double Cpb,Cvb,gamma,rho,flameArea,Tu;
+
+
+	/*Find the isotherm chosen by the user*/
+	size_t j=1;
+	size_t jj=1;
+	size_t jjj=1;
+	double wdot[data->nsp];
+	double wdotMax=0.0e0;
+	double advTerm=0.0e0;
+	psiAvg=0.0e0;
+	if(T(data->npts)>T(1)){
+		while (T(j)<TAvg) {
+			j=j+1;
+		}
+		YAvg=innerMassFractionsData[data->k_oxidizer-1]-Y(data->npts,data->k_oxidizer);
+		while (fabs((T(jj+1)-T(jj))/T(jj))>1e-08) {
+			jj=jj+1;
+		}
+
+		setGas(data,ydata,jj);
+		Tu=T(jj);
+		rho=data->gas->density();
+        	Cpb=data->gas->cp_mass();       	//J/kg/K
+        	Cvb=data->gas->cv_mass();       	//J/kg/K
+		gamma=Cpb/Cvb;
+		
+	}
+	else{
+		while (T(j)>TAvg) {
+			j=j+1;
+		}
+		YAvg=innerMassFractionsData[data->k_oxidizer-1]-Y(1,data->k_oxidizer);
+		while (fabs((T(data->npts-jj-1)-T(data->npts-jj))/T(data->npts-jj))>1e-08) {
+			jj=jj+1;
+		}
+
+		setGas(data,ydata,data->npts-jj);
+		Tu=T(data->npts-jj);
+		rho=data->gas->density();
+        	Cpb=data->gas->cp_mass();       	//J/kg/K
+        	Cvb=data->gas->cv_mass();       	//J/kg/K
+		gamma=Cpb/Cvb;
+	}
+	if(T(j)<TAvg){
+		RAvg=((R(j+1)-R(j))/(T(j+1)-T(j)))*(TAvg-T(j))+R(j);
+	}
+	else{
+		RAvg=((R(j)-R(j-1))/(T(j)-T(j-1)))*(TAvg-T(j-1))+R(j-1);
+	}
+
+	////Experimental burning speed calculation:	
+	//int nMax=0;
+	////nMax=maxCurvIndex(ydata, data->nt, data->nvar,
+	////     	data->grid->x, data->npts);
+	//nMax=maxGradIndex(ydata, data->nt, data->nvar,
+	//     	data->grid->x, data->npts);
+	//advTerm=(T(nMax)-T(nMax-1))/(data->mass*(psi(nMax)-psi(nMax-1)));
+	//aream=calc_area(R(nMax),&data->metric);	    
+	////setGas(data,ydata,nMax);
+	////rho=data->gas->density();
+	//psiAvg=-Tdot(nMax)/(rho*aream*advTerm);
+	////if(t>data->ignTime){
+	////	for(size_t n=2;n<data->npts;n++){
+	////		setGas(data,ydata,n);
+	////		data->gas->getNetProductionRates(wdot); //kmol/m^3
+	////		advTerm=(T(n)-T(n-1))/(data->mass*(psi(n)-psi(n-1)));
+	////		if(fabs(wdot[data->k_oxidizer-1])>=wdotMax){
+	////			aream=calc_area(R(n),&data->metric);	    
+	////			psiAvg=-Tdot(n)/(rho*aream*advTerm);
+	////			wdotMax=fabs(wdot[data->k_oxidizer-1]);
+	////		}
+	////	}
+	////}
+	////else{
+	////	psiAvg=0.0e0;
+	////}
+
+
+	//drdt=(RAvg-data->flamePosition[1])/(t-data->flameTime[1]);
+	//data->flamePosition[0]=data->flamePosition[1];
+	//data->flamePosition[1]=RAvg;
+	//data->flameTime[0]=data->flameTime[1];
+	//data->flameTime[1]=t;
+	//flameArea=calc_area(RAvg,&data->metric);
+
+	/*Use the Trapezoidal rule to calculate the mass burning rate based on
+	 * the consumption of O2*/
+	aream= calc_area(R(1)+1e-03*data->domainLength,&data->metric);	    
+	for (j = 2; j <data->npts; j++) {
+        	dpsim=(psi(j)-psi(j-1))*data->mass;                            
+		area= calc_area(R(j),&data->metric);	    
+		sum=sum+HALF*dpsim*((Ydot(j-1,data->k_oxidizer)/aream)
+				   +(Ydot(j,data->k_oxidizer)/area));
+		aream=area;
+	}
+
+	//double maxOH,maxHO2;
+	//maxOH=0.0e0;
+	//maxHO2=0.0e0;
+	//for(j=1;j<data->npts;j++){
+	//	if(Y(j,data->k_OH)>maxOH){
+	//		maxOH=Y(j,data->k_OH);
+	//	}
+	//}
+	//for(j=1;j<data->npts;j++){
+	//	if(Y(j,data->k_HO2)>maxHO2){
+	//		maxHO2=Y(j,data->k_HO2);
+	//	}
+	//}
+
+	fprintf((data->globalOutput), "%15.6e\t",t);
+	//fprintf((data->globalOutput), "%15.6e\t",psiAvg);
+	fprintf((data->globalOutput), "%15.6e\t",fabs(sum)/YAvg);
+	fprintf((data->globalOutput), "%15.6e\t",RAvg);
+	fprintf((data->globalOutput), "%15.6e\t",P(data->npts));
+	fprintf((data->globalOutput), "%15.6e\t",Pdot(data->npts));
+	fprintf((data->globalOutput), "%15.6e\t",gamma);
+	fprintf((data->globalOutput), "%15.6e\t",fabs(sum)/(YAvg*rho));
+	fprintf((data->globalOutput), "%15.6e\t",Tu);
+	fprintf((data->globalOutput), "\n");
+}
+
+//
+//void printSpaceTimeOutputInterpolated(double t, N_Vector y, UserData data)
+//{
+//	double *ydata,*psidata;
+//  	ydata    = N_VGetArrayPointer_OpenMP(y);
+//  	psidata    = N_VGetArrayPointer_OpenMP(data->grid);
+//	for (int i = 0; i < data->npts; i++) {
+//		fprintf((data->gridOutput), "%15.6e\t%15.6e\t",psi(i+1),t);
+//		for (int j = 0; j < data->nvar; j++) {
+//			fprintf((data->gridOutput), "%15.6e\t",ydata[j+i*data->nvar]);
+//		}
+//		fprintf((data->gridOutput), "\n");
+//	}
+//	fprintf((data->gridOutput), "\n\n");
+//}
+//
+//
+////void repairSolution(N_Vector y, N_Vector ydot, UserData data){
+////	int npts=data->npts;
+////  	double *ydata;
+////  	double *ydotdata;
+////	ydata    = N_VGetArrayPointer_OpenMP(y); 
+////	ydotdata    = N_VGetArrayPointer_OpenMP(ydot); 
+////
+////	T(2)=T(1);
+////	T(npts-1)=T(npts);
+////	Tdot(2)=Tdot(1);
+////	Tdot(npts-1)=Tdot(npts);
+////	for (int k = 1; k <=data->nsp; k++) {
+////		    Y(2,k)=Y(1,k);
+////		    Y(npts-1,k)=Y(npts,k);
+////
+////		    Ydot(2,k)=Ydot(1,k);
+////		    Ydot(npts-1,k)=Ydot(npts,k);
+////	}
+////}
diff --git a/residue.h b/residue.h
new file mode 100644
index 0000000..6222fc9
--- /dev/null
+++ b/residue.h
@@ -0,0 +1,90 @@
+#ifndef SUNDIALS_DEF
+#define SUNDIALS_DEF
+#include <sundials/sundials_types.h>
+#include <nvector/nvector_openmp.h>
+#endif
+
+#ifndef PRINT_DEF
+#define PRINT_DEF
+#include <string.h>	//for strings
+#include <stdio.h>	//for printf,scanf
+#include <stdlib.h>	//for atoi, atof
+#endif
+
+#ifndef CANTERA_DEF
+#define CANTERA_DEF
+#include <cantera/IdealGasMix.h>
+#include <cantera/transport.h>
+#endif
+
+#include "UserData.h"
+
+
+double maxGradPosition(const double* y, const size_t nt, 
+	               const size_t nvar, const double* x, size_t nPts);
+
+int maxGradIndex(const double* y, const size_t nt, 
+	         const size_t nvar, const double* x, size_t nPts);
+
+double maxCurvPosition(const double* y, const size_t nt, 
+	               const size_t nvar, const double* x, size_t nPts);
+
+int maxCurvIndex(const double* y, const size_t nt, 
+	               const size_t nvar, const double* x, size_t nPts);
+
+double isothermPosition(const double* y, const double T, const size_t nt, 
+	                const size_t nvar, const double* x, const size_t nPts);
+
+int setAlgebraicVariables(N_Vector *id,UserData data);
+
+inline double calc_area(double x,int* i);
+
+void updateSolution(double* y, double* ydot, const size_t nvar,
+		    const double xOld[],const double xNew[],const size_t npts);
+
+void readInitialCondition(FILE* input, double* ydata, const size_t nvar, const size_t nr, const size_t nPts);
+
+double systemMass(double* ydata, UserData data);
+
+int initializePsiGrid(double* ydata, double* psidata, UserData data);
+
+int setInitialCondition(N_Vector* y, 
+	    		N_Vector* ydot, 
+            		UserData data);
+
+
+inline void setGas(UserData data, double *ydata, size_t gridPoint);
+
+void getTransport(UserData data, 
+		  double *ydata, 
+		  size_t gridPoint,
+		  double *rho,
+		  double *lambda,
+		  double *YV);
+
+int residue(double t, 
+	    N_Vector y, 
+	    N_Vector ydot, 
+	    N_Vector res, 
+            void *user_data);
+
+void trackFlameOH(N_Vector y,UserData data);
+void trackFlame(N_Vector y,UserData data);
+size_t BathGasIndex(UserData data);
+size_t oxidizerIndex(UserData data);
+
+inline double Qdot(double* t, 
+	           double* x, 
+	           double* ignTime, 
+	           double* kernelSize, 
+	           double* maxQdot);
+
+void printSpaceTimeHeader(UserData data);
+void printSpaceTimeOutput(double t, N_Vector* y, FILE* output, UserData data);
+void printSpaceTimeRates(double t, N_Vector ydot, UserData data);
+void printGlobalHeader(UserData data);
+void printGlobalVariables(double t, N_Vector* y, N_Vector* ydot, UserData data);
+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);
diff --git a/solution.cpp b/solution.cpp
new file mode 100644
index 0000000..9e13b06
--- /dev/null
+++ b/solution.cpp
@@ -0,0 +1,90 @@
+#include "solution.h"
+int allocateSolution(size_t neq, 
+		     int nThreads,
+		     N_Vector *y, 
+		     N_Vector *ydot, 
+		     N_Vector *res, 
+		     N_Vector *id, 
+		     N_Vector *atolv,
+		     N_Vector *constraints){
+
+
+  	*y = N_VNew_OpenMP(neq,nThreads);
+  	if(*y==NULL){
+		printf("Allocation Failed!\n");
+		return(1);
+	}
+  	N_VConst((realtype)(0.0e0), *y);
+
+  	*ydot = N_VNew_OpenMP(neq,nThreads);
+  	if(*ydot==NULL){
+		printf("Allocation Failed!\n");
+		return(1);
+	}
+  	N_VConst((realtype)(0.0e0), *ydot);
+
+  	*res = N_VNew_OpenMP(neq,nThreads);
+  	if(*res==NULL){
+		printf("Allocation Failed!\n");
+		return(1);
+	}
+  	N_VConst((realtype)(0.0e0), *res);
+
+  	*id = N_VNew_OpenMP(neq,nThreads);
+  	if(*id==NULL){
+		printf("Allocation Failed!\n");
+		return(1);
+	}
+  	N_VConst((realtype)(1.0e0), *id);
+
+  	*atolv = N_VNew_OpenMP(neq,nThreads);
+  	if(*atolv==NULL){
+		printf("Allocation Failed!\n");
+		return(1);
+	}
+  	N_VConst((realtype)(0.0e0), *atolv);
+
+  	*constraints = N_VNew_OpenMP(neq,nThreads);
+  	if(*constraints==NULL){
+		printf("Allocation Failed!\n");
+		return(1);
+	}
+  	N_VConst((realtype)(0.0e0), *constraints);
+
+	return(0);
+}
+
+void freeSolution(N_Vector *y, 
+		  N_Vector *ydot, 
+		  N_Vector *res, 
+		  N_Vector *id, 
+		  N_Vector *atolv,
+		  N_Vector *constraints){
+
+  	if(*y!=NULL){
+		N_VDestroy_OpenMP(*y);
+		printf("y Destroyed!\n");
+	}
+	if(*ydot!=NULL){
+  		N_VDestroy_OpenMP(*ydot);
+		printf("ydot Destroyed!\n");
+	}
+	if(*res!=NULL){
+  		N_VDestroy_OpenMP(*res);
+		printf("res Destroyed!\n");
+	}
+	if(*id!=NULL){
+  		N_VDestroy_OpenMP(*id);
+		printf("id Destroyed!\n");
+	}
+	if(*atolv!=NULL){
+  		N_VDestroy_OpenMP(*atolv);
+		printf("atolv Destroyed!\n");
+	}
+
+	if(*constraints!=NULL){
+  		N_VDestroy_OpenMP(*constraints);
+		printf("constraints Destroyed!\n");
+	}
+	printf("\n\n");
+}
diff --git a/solution.h b/solution.h
new file mode 100644
index 0000000..7f220c3
--- /dev/null
+++ b/solution.h
@@ -0,0 +1,15 @@
+#ifndef PRINT_DEF
+#define PRINT_DEF
+#include <string.h>	//for strings
+#include <stdio.h>	//for printf,scanf
+#include <stdlib.h>	//for atoi, atof
+#endif
+
+#ifndef SUNDIALS_DEF
+#define SUNDIALS_DEF
+#include <sundials/sundials_types.h>
+#include <nvector/nvector_openmp.h>
+#endif
+
+int allocateSolution(size_t neq, int nThreads, N_Vector *y, N_Vector *ydot, N_Vector *res, N_Vector *id, N_Vector *atolv, N_Vector *constraints);
+void freeSolution(N_Vector *y, N_Vector *ydot, N_Vector *res, N_Vector *id, N_Vector *atolv, N_Vector *constraints);
diff --git a/timing.h b/timing.h
new file mode 100644
index 0000000..b65b3d3
--- /dev/null
+++ b/timing.h
@@ -0,0 +1,2 @@
+double get_wall_time();
+double get_cpu_time();
diff --git a/timing.hpp b/timing.hpp
new file mode 100644
index 0000000..e07e3eb
--- /dev/null
+++ b/timing.hpp
@@ -0,0 +1,45 @@
+//  Windows
+#ifdef _WIN32
+#include <Windows.h>
+double get_wall_time(){
+    LARGE_INTEGER time,freq;
+    if (!QueryPerformanceFrequency(&freq)){
+        //  Handle error
+        return 0;
+    }
+    if (!QueryPerformanceCounter(&time)){
+        //  Handle error
+        return 0;
+    }
+    return (double)time.QuadPart / freq.QuadPart;
+}
+double get_cpu_time(){
+    FILETIME a,b,c,d;
+    if (GetProcessTimes(GetCurrentProcess(),&a,&b,&c,&d) != 0){
+        //  Returns total user time.
+        //  Can be tweaked to include kernel times as well.
+        return
+            (double)(d.dwLowDateTime |
+            ((unsigned long long)d.dwHighDateTime << 32)) * 0.0000001;
+    }else{
+        //  Handle error
+        return 0;
+    }
+}
+
+//  Posix/Linux
+#else
+#include <time.h>
+#include <sys/time.h>
+double get_wall_time(){
+    struct timeval time;
+    if (gettimeofday(&time,NULL)){
+        //  Handle error
+        return 0;
+    }
+    return (double)time.tv_sec + (double)time.tv_usec * .000001;
+}
+double get_cpu_time(){
+    return (double)clock() / CLOCKS_PER_SEC;
+}
+#endif