From 48cca4770f31efa3c771e35d217aee00e20b2f14 Mon Sep 17 00:00:00 2001
From: vyaas <mutantbex@gmail.com>
Date: Sun, 10 May 2020 22:40:41 -0500
Subject: [PATCH] first

---
 Makefile         |   91 ++++
 UserData.cpp     |  409 +++++++++++++++
 UserData.h       |  170 ++++++
 gridRoutines.cpp |  244 +++++++++
 gridRoutines.h   |   82 +++
 input.dat        |   50 ++
 macros.h         |   63 +++
 main.cpp         |  338 ++++++++++++
 parse.cpp        |   20 +
 parse.h          |   31 ++
 parse.hpp        |   76 +++
 residue.cpp      | 1303 ++++++++++++++++++++++++++++++++++++++++++++++
 residue.h        |   90 ++++
 solution.cpp     |   90 ++++
 solution.h       |   15 +
 timing.h         |    2 +
 timing.hpp       |   45 ++
 17 files changed, 3119 insertions(+)
 create mode 100644 Makefile
 create mode 100644 UserData.cpp
 create mode 100644 UserData.h
 create mode 100644 gridRoutines.cpp
 create mode 100644 gridRoutines.h
 create mode 100644 input.dat
 create mode 100644 macros.h
 create mode 100644 main.cpp
 create mode 100644 parse.cpp
 create mode 100644 parse.h
 create mode 100644 parse.hpp
 create mode 100644 residue.cpp
 create mode 100644 residue.h
 create mode 100644 solution.cpp
 create mode 100644 solution.h
 create mode 100644 timing.h
 create mode 100644 timing.hpp

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