diff --git a/.gitignore b/.gitignore
index 1193b2d..23da90b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -18,3 +18,9 @@ lnInclude
/laminarReactingFoam/Make/*
!/laminarReactingFoam/Make/file
!/laminarReactingFoam/Make/option
+
+/testApp/Make/*
+!/testApp/Make/file
+!/testApp/Make/option
+
+canteraTest
diff --git a/testApp/Make/files b/testApp/Make/files
new file mode 100644
index 0000000..33e8507
--- /dev/null
+++ b/testApp/Make/files
@@ -0,0 +1,3 @@
+chemFoam.C
+
+EXE = canteraTest
diff --git a/testApp/Make/options b/testApp/Make/options
new file mode 100644
index 0000000..d08ab34
--- /dev/null
+++ b/testApp/Make/options
@@ -0,0 +1,15 @@
+EXE_INC = \
+ -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+ -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+ -I$(LIB_SRC)/ODE/lnInclude\
+ -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+ -lchemistryModel \
+ -lfiniteVolume \
+ -lmeshTools \
+ -lcantera
diff --git a/testApp/YEqn.H b/testApp/YEqn.H
new file mode 100644
index 0000000..edb4072
--- /dev/null
+++ b/testApp/YEqn.H
@@ -0,0 +1,12 @@
+{
+ forAll(Y, specieI)
+ {
+ volScalarField& Yi = Y[specieI];
+
+ solve
+ (
+ fvm::ddt(rho, Yi) - chemistry.RR(specieI),
+ mesh.solver("Yi")
+ );
+ }
+}
\ No newline at end of file
diff --git a/testApp/chemFoam.C b/testApp/chemFoam.C
new file mode 100644
index 0000000..a68d7d0
--- /dev/null
+++ b/testApp/chemFoam.C
@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see .
+
+Application
+ chemFoam
+
+Description
+ Solver for chemistry problems, designed for use on single cell cases to
+ provide comparison against other chemistry solvers, that uses a single cell
+ mesh, and fields created from the initial conditions.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "psiReactionThermo.H"
+#include "psiChemistryModel.H"
+#include "chemistrySolver.H"
+#include "OFstream.H"
+#include "thermoPhysicsTypes.H"
+#include "basicMultiComponentMixture.H"
+#include "cellModeller.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+ argList::noParallel();
+
+ #define CREATE_MESH createSingleCellMesh.H
+ #define NO_CONTROL
+ #include "postProcess.H"
+
+ #include "setRootCase.H"
+ #include "createTime.H"
+ #include "createSingleCellMesh.H"
+ #include "createFields.H"
+ #include "createFieldRefs.H"
+ #include "readInitialConditions.H"
+ #include "createControls.H"
+
+ // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ Info<< "\nStarting time loop\n" << endl;
+
+ while (runTime.run())
+ {
+ #include "readControls.H"
+
+ #include "setDeltaT.H"
+
+ runTime++;
+ Info<< "Time = " << runTime.timeName() << nl << endl;
+
+ #include "solveChemistry.H"
+ #include "YEqn.H"
+ #include "hEqn.H"
+ #include "pEqn.H"
+
+ #include "output.H"
+
+ Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+ << " ClockTime = " << runTime.elapsedClockTime() << " s"
+ << nl << endl;
+ }
+
+ Info << "Number of steps = " << runTime.timeIndex() << endl;
+ Info << "End" << nl << endl;
+
+ return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/testApp/createBaseFields.H b/testApp/createBaseFields.H
new file mode 100644
index 0000000..0762f27
--- /dev/null
+++ b/testApp/createBaseFields.H
@@ -0,0 +1,57 @@
+// write base thermo fields - not registered since will be re-read by
+// thermo package
+
+Info<< "Creating base fields for time " << runTime.timeName() << endl;
+{
+ volScalarField Ydefault
+ (
+ IOobject
+ (
+ "Ydefault",
+ runTime.timeName(),
+ mesh,
+ IOobject::READ_IF_PRESENT,
+ IOobject::NO_WRITE,
+ false
+ ),
+ mesh,
+ dimensionedScalar("Ydefault", dimless, 1)
+ );
+
+ Ydefault.write();
+
+ volScalarField p
+ (
+ IOobject
+ (
+ "p",
+ runTime.timeName(),
+ mesh,
+ IOobject::READ_IF_PRESENT,
+ IOobject::NO_WRITE,
+ false
+ ),
+ mesh,
+ dimensionedScalar("p", dimPressure, p0)
+ );
+
+ p.write();
+
+ volScalarField T
+ (
+ IOobject
+ (
+ "T",
+ runTime.timeName(),
+ mesh,
+ IOobject::READ_IF_PRESENT,
+ IOobject::NO_WRITE,
+ false
+ ),
+ mesh,
+ dimensionedScalar("T", dimTemperature, T0)
+ );
+
+ T.write();
+}
+
diff --git a/testApp/createControls.H b/testApp/createControls.H
new file mode 100644
index 0000000..6aeb2ff
--- /dev/null
+++ b/testApp/createControls.H
@@ -0,0 +1,2 @@
+Switch adjustTimeStep(runTime.controlDict().lookup("adjustTimeStep"));
+scalar maxDeltaT(readScalar(runTime.controlDict().lookup("maxDeltaT")));
diff --git a/testApp/createFieldRefs.H b/testApp/createFieldRefs.H
new file mode 100644
index 0000000..ef34180
--- /dev/null
+++ b/testApp/createFieldRefs.H
@@ -0,0 +1,4 @@
+scalar dtChem = refCast(chemistry).deltaTChem()[0];
+basicMultiComponentMixture& composition = thermo.composition();
+PtrList& Y = composition.Y();
+volScalarField& p = thermo.p();
diff --git a/testApp/createFields.H b/testApp/createFields.H
new file mode 100644
index 0000000..62be3ce
--- /dev/null
+++ b/testApp/createFields.H
@@ -0,0 +1,84 @@
+ if (mesh.nCells() != 1)
+ {
+ FatalErrorInFunction
+ << "Solver only applicable to single cell cases"
+ << exit(FatalError);
+ }
+
+ Info<< "Reading initial conditions.\n" << endl;
+ IOdictionary initialConditions
+ (
+ IOobject
+ (
+ "initialConditions",
+ runTime.constant(),
+ runTime,
+ IOobject::MUST_READ_IF_MODIFIED,
+ IOobject::NO_WRITE
+ )
+ );
+
+ scalar p0 = readScalar(initialConditions.lookup("p"));
+ scalar T0 = readScalar(initialConditions.lookup("T"));
+
+ #include "createBaseFields.H"
+
+ Info<< nl << "Reading thermophysicalProperties" << endl;
+ autoPtr pChemistry(psiChemistryModel::New(mesh));
+
+ psiChemistryModel& chemistry = pChemistry();
+
+ psiReactionThermo& thermo = chemistry.thermo();
+ thermo.validate(args.executable(), "h");
+
+ volScalarField rho
+ (
+ IOobject
+ (
+ "rho",
+ runTime.timeName(),
+ runTime,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ thermo.rho()
+ );
+
+ volScalarField Rspecific
+ (
+ IOobject
+ (
+ "Rspecific",
+ runTime.timeName(),
+ runTime,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar
+ (
+ "zero",
+ dimensionSet(dimEnergy/dimMass/dimTemperature),
+ 0.0
+ )
+ );
+
+ volVectorField U
+ (
+ IOobject
+ (
+ "U",
+ runTime.timeName(),
+ runTime,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ mesh,
+ dimensionedVector("zero", dimVelocity, Zero)
+ );
+
+ #include "createPhi.H"
+
+ OFstream post(args.path()/"chemFoam.out");
+ post<< "# Time" << token::TAB << "Temperature [K]" << token::TAB
+ << "Pressure [Pa]" << endl;
diff --git a/testApp/createSingleCellMesh.H b/testApp/createSingleCellMesh.H
new file mode 100644
index 0000000..9b9d0f3
--- /dev/null
+++ b/testApp/createSingleCellMesh.H
@@ -0,0 +1,46 @@
+Info<< "Constructing single cell mesh" << nl << endl;
+
+labelList owner(6, label(0));
+labelList neighbour(0);
+
+pointField points(8);
+points[0] = vector(0, 0, 0);
+points[1] = vector(1, 0, 0);
+points[2] = vector(1, 1, 0);
+points[3] = vector(0, 1, 0);
+points[4] = vector(0, 0, 1);
+points[5] = vector(1, 0, 1);
+points[6] = vector(1, 1, 1);
+points[7] = vector(0, 1, 1);
+
+const cellModel& hexa = *(cellModeller::lookup("hex"));
+faceList faces = hexa.modelFaces();
+
+fvMesh mesh
+(
+ IOobject
+ (
+ fvMesh::defaultRegion,
+ runTime.timeName(),
+ runTime,
+ IOobject::READ_IF_PRESENT
+ ),
+ xferMove>(points),
+ faces.xfer(),
+ owner.xfer(),
+ neighbour.xfer()
+);
+
+List patches(1);
+
+patches[0] = new emptyPolyPatch
+(
+ "boundary",
+ 6,
+ 0,
+ 0,
+ mesh.boundaryMesh(),
+ emptyPolyPatch::typeName
+);
+
+mesh.addFvPatches(patches);
diff --git a/testApp/hEqn.H b/testApp/hEqn.H
new file mode 100644
index 0000000..4dd593e
--- /dev/null
+++ b/testApp/hEqn.H
@@ -0,0 +1,14 @@
+{
+ volScalarField& h = thermo.he();
+
+ if (constProp == "volume")
+ {
+ h[0] = u0 + p[0]/rho[0] + integratedHeat;
+ }
+ else
+ {
+ h[0] = h0 + integratedHeat;
+ }
+
+ thermo.correct();
+}
diff --git a/testApp/output.H b/testApp/output.H
new file mode 100644
index 0000000..5620fe6
--- /dev/null
+++ b/testApp/output.H
@@ -0,0 +1,11 @@
+ runTime.write();
+
+ Info<< "Sh = " << Sh
+ << ", T = " << thermo.T()[0]
+ << ", p = " << thermo.p()[0]
+ << ", " << Y[0].name() << " = " << Y[0][0]
+ << endl;
+
+ post<< runTime.value() << token::TAB << thermo.T()[0] << token::TAB
+ << thermo.p()[0] << endl;
+
diff --git a/testApp/pEqn.H b/testApp/pEqn.H
new file mode 100644
index 0000000..a919ff1
--- /dev/null
+++ b/testApp/pEqn.H
@@ -0,0 +1,16 @@
+{
+ rho = thermo.rho();
+ if (constProp == "volume")
+ {
+ scalar invW = 0.0;
+ forAll(Y, i)
+ {
+ invW += Y[i][0]/specieData[i].W();
+ }
+
+ Rspecific[0] = 1000.0*constant::physicoChemical::R.value()*invW;
+
+ p[0] = rho0*Rspecific[0]*thermo.T()[0];
+ rho[0] = rho0;
+ }
+}
diff --git a/testApp/readControls.H b/testApp/readControls.H
new file mode 100644
index 0000000..c56206d
--- /dev/null
+++ b/testApp/readControls.H
@@ -0,0 +1,3 @@
+runTime.controlDict().lookup("adjustTimeStep") >> adjustTimeStep;
+
+maxDeltaT = readScalar(runTime.controlDict().lookup("maxDeltaT"));
diff --git a/testApp/readInitialConditions.H b/testApp/readInitialConditions.H
new file mode 100644
index 0000000..8e340a7
--- /dev/null
+++ b/testApp/readInitialConditions.H
@@ -0,0 +1,108 @@
+ word constProp(initialConditions.lookup("constantProperty"));
+ if ((constProp != "pressure") && (constProp != "volume"))
+ {
+ FatalError << "in initialConditions, unknown constantProperty type "
+ << constProp << nl << " Valid types are: pressure volume."
+ << abort(FatalError);
+ }
+
+ word fractionBasis(initialConditions.lookup("fractionBasis"));
+ if ((fractionBasis != "mass") && (fractionBasis != "mole"))
+ {
+ FatalError << "in initialConditions, unknown fractionBasis type " << nl
+ << "Valid types are: mass or mole."
+ << fractionBasis << abort(FatalError);
+ }
+
+ label nSpecie = Y.size();
+ PtrList specieData(Y.size());
+ forAll(specieData, i)
+ {
+ specieData.set
+ (
+ i,
+ new gasHThermoPhysics
+ (
+ dynamic_cast&>
+ (thermo).speciesData()[i]
+ )
+ );
+ }
+
+ scalarList Y0(nSpecie, 0.0);
+ scalarList X0(nSpecie, 0.0);
+
+ dictionary fractions(initialConditions.subDict("fractions"));
+ if (fractionBasis == "mole")
+ {
+ forAll(Y, i)
+ {
+ const word& name = Y[i].name();
+ if (fractions.found(name))
+ {
+ X0[i] = readScalar(fractions.lookup(name));
+ }
+ }
+
+ scalar mw = 0.0;
+ const scalar mTot = sum(X0);
+ forAll(Y, i)
+ {
+ X0[i] /= mTot;
+ mw += specieData[i].W()*X0[i];
+ }
+
+ forAll(Y, i)
+ {
+ Y0[i] = X0[i]*specieData[i].W()/mw;
+ }
+ }
+ else // mass fraction
+ {
+ forAll(Y, i)
+ {
+ const word& name = Y[i].name();
+ if (fractions.found(name))
+ {
+ Y0[i] = readScalar(fractions.lookup(name));
+ }
+ }
+
+ scalar invW = 0.0;
+ const scalar mTot = sum(Y0);
+ forAll(Y, i)
+ {
+ Y0[i] /= mTot;
+ invW += Y0[i]/specieData[i].W();
+ }
+ const scalar mw = 1.0/invW;
+
+ forAll(Y, i)
+ {
+ X0[i] = Y0[i]*mw/specieData[i].W();
+ }
+ }
+
+ scalar h0 = 0.0;
+ forAll(Y, i)
+ {
+ Y[i] = Y0[i];
+ h0 += Y0[i]*specieData[i].Hs(p[0], T0);
+ }
+
+ thermo.he() = dimensionedScalar("h", dimEnergy/dimMass, h0);
+ thermo.correct();
+
+ rho = thermo.rho();
+ scalar rho0 = rho[0];
+ scalar u0 = h0 - p0/rho0;
+ scalar R0 = p0/(rho0*T0);
+ Rspecific[0] = R0;
+
+ scalar integratedHeat = 0.0;
+
+ Info << constProp << " will be held constant." << nl
+ << " p = " << p[0] << " [Pa]" << nl
+ << " T = " << thermo.T()[0] << " [K] " << nl
+ << " rho = " << rho[0] << " [kg/m3]" << nl
+ << endl;
diff --git a/testApp/setDeltaT.H b/testApp/setDeltaT.H
new file mode 100644
index 0000000..46d9f7b
--- /dev/null
+++ b/testApp/setDeltaT.H
@@ -0,0 +1,6 @@
+if (adjustTimeStep)
+{
+ runTime.setDeltaT(min(dtChem, maxDeltaT));
+ Info<< "deltaT = " << runTime.deltaT().value() << endl;
+}
+
diff --git a/testApp/solveChemistry.H b/testApp/solveChemistry.H
new file mode 100644
index 0000000..cb4258c
--- /dev/null
+++ b/testApp/solveChemistry.H
@@ -0,0 +1,3 @@
+ dtChem = chemistry.solve(runTime.deltaT().value());
+ scalar Sh = chemistry.Sh()()[0]/rho[0];
+ integratedHeat += Sh*runTime.deltaT().value();