initial import

apps/bvp/AxiStagnBVP.cpp

@@ -0,0 +1,144 @@
+/// @file AxiStagnBVP.cpp
+
+#include <cantera/Cantera.h>
+#include "AxiStagnBVP.h"
+
+AxiStagnBVP::AxiStagnBVP(int nsp, int np, double L) : 
+    BVP::BoundaryValueProblem(nsp+4, 
+        np, 0.0, L) {
+
+    // specify the component bounds, error tolerances, and names.
+        BVP::Component u;
+        u.lower = -200.0;
+        u.upper = 200.0;
+        u.rtol = 1.0e-8;
+        u.atol = 1.0e-15;
+        u.name = "u";
+        setComponent(0, u);  // the axial velocity will be component 0
+
+        BVP::Component V;
+        V.lower = -1.0e8;
+        V.upper = 1.0e8;
+        V.rtol = 1.0e-8;
+        V.atol = 1.0e-15;
+        V.name = "V";
+        setComponent(1, V);  // the radial velocity will be component 1
+
+        BVP::Component T;
+        T.lower = 200.0;
+        T.upper = 1.0e9;
+        T.rtol = 1.0e-8;
+        T.atol = 1.0e-15;
+        T.name = "T";
+        setComponent(2, T);  // the temperature will be component 2
+
+        BVP::Component lambda;
+        lambda.lower = -1.0e20;
+        lambda.upper = 1.0e20;
+        lambda.rtol = 1.0e-8;
+        lambda.atol = 1.0e-15;
+        lambda.name = "Lambda";
+        setComponent(3, lambda);  // the pressure-gradient eigenvalue will be 
+                                  //component 3
+        BVP::Component Y;
+        Y.lower = -1.0e-5;
+        Y.upper = 1.0e2;
+        Y.rtol = 1.0e-8;
+        Y.atol = 1.0e-15;
+        for (k = 0; k < nsp; k++) {
+            Y.name = thermo->speciesName(k);
+            setComponent(k+4, Y);  
+        }
+}
+
+
+// destructor
+AxiStagnBVP::~AxiStagnBVP() {}
+
+// specify guesses for the initial values. These can be anything
+// that leads to a converged solution.
+doublereal initialValue(int n, int j) { 
+    switch (n) {
+    case 0:
+        return m_u0;
+    case 1:
+        return m_u0/m_L;
+    case 2:
+        return m_Tinf;
+    case 4:
+        return 1.0;
+    default:
+        return 0.0;
+    }
+}
+
+// Specify the residual. This is where the ODE system and boundary 
+// conditions are specified. The solver will attempt to find a solution
+// x so that this function returns 0 for all n and j.
+doublereal AxiStagnFlow::residual(doublereal* x, int n, int j) {
+    
+    // if n = 0, return the residual for the continuity equation
+    if (n == 0) {
+        if (isRight(j)) {
+            return -rho_u(x,j);  // force u to zero at the right
+        }
+        else {
+            return -(rho_u(x, j+1) - rho_u(x,j))/m_dz[j]
+                -(density(j+1)*V(x,j+1) + density(j)*V(x,j));
+        }
+    }
+
+    else if (n == 1) {
+
+        // if n = 1, then return the residual for radial momentum
+        if (isLeft(j)) {
+            return V(x,j);
+        }
+        else if (isRight(j)) {
+            return V(x,j);             // force V to zero at the wall
+        }
+        else {
+            return (shear(x,j) - lambda(x,j) - rho_u(x,j)*dVdz(x,j) 
+                - m_rho[j]*V(x,j)*V(x,j))/m_rho[j]
+                - rdt*(V(x,j) - V_prev(j));
+        }
+
+    }
+
+    else if (n == 2) {
+        if (isLeft(j)) {
+            return T(x,j) - m_Tinlet;
+        }
+        else if (isRight(j)) {
+            return T(x,j) - m_Tsurf;
+        }
+        else {
+            setGas(x,j);
+
+            // heat release term
+            const vector_fp& h_RT = m_thermo->enthalpy_RT_ref();
+            const vector_fp& cp_R = m_thermo->cp_R_ref();
+
+            sum = 0.0;
+            sum2 = 0.0;
+            doublereal flxk;
+            for (k = 0; k < m_nsp; k++) {
+                flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j));
+                sum += wdot(k,j)*h_RT[k];
+                sum2 += flxk*cp_R[k]/m_wt[k];
+            }
+            sum *= GasConstant * T(x,j);
+            dtdzj = dTdz(x,j);
+            sum2 *= GasConstant * dtdzj;
+
+            rsd =    - m_cp[j]*rho_u(x,j)*dtdzj 
+                     - divHeatFlux(x,j) - sum - sum2;
+            rsd /= (m_rho[j]*m_cp[j]);
+
+            rsd -= rdt*(T(x,j) - T_prev(j));
+        }            
+    }
+
+}
+
+

apps/bvp/AxiStagnBVP.h

@@ -0,0 +1,131 @@
+/// @file AxiStagnBVP.h
+
+#include <cantera/Cantera.h>
+#include "BoundaryValueProblem.h"
+
+
+/**
+ * This class solves 
+ */
+class AxiStagnBVP : public BVP::BoundaryValueProblem {
+
+public:
+
+    AxiStagnBVP(int nsp, int np, double L) : BVP::BoundaryValueProblem(nsp+4, 
+        np, 0.0, L) {
+
+        // specify the component bounds, error tolerances, and names.
+        BVP::Component u;
+        u.lower = -200.0;
+        u.upper = 200.0;
+        u.rtol = 1.0e-8;
+        u.atol = 1.0e-15;
+        u.name = "u";
+        setComponent(0, u);  // the axial velocity will be component 0
+
+        BVP::Component V;
+        V.lower = -1.0e8;
+        V.upper = 1.0e8;
+        V.rtol = 1.0e-8;
+        V.atol = 1.0e-15;
+        V.name = "V";
+        setComponent(1, V);  // the radial velocity will be component 1
+
+        BVP::Component T;
+        T.lower = 200.0;
+        T.upper = 1.0e9;
+        T.rtol = 1.0e-8;
+        T.atol = 1.0e-15;
+        T.name = "T";
+        setComponent(2, T);  // the temperature will be component 2
+
+        BVP::Component lambda;
+        lambda.lower = -1.0e20;
+        lambda.upper = 1.0e20;
+        lambda.rtol = 1.0e-8;
+        lambda.atol = 1.0e-15;
+        lambda.name = "Lambda";
+        setComponent(3, lambda);  // the pressure-gradient eigenvalue will be 
+                                  //component 3
+        BVP::Component Y;
+        Y.lower = -1.0e-5;
+        Y.upper = 1.0e2;
+        Y.rtol = 1.0e-8;
+        Y.atol = 1.0e-15;
+        for (k = 0; k < nsp; k++) {
+            Y.name = thermo->speciesName(k);
+            setComponent(k+4, Y);  
+        }
+    }
+
+
+    // destructor
+    virtual ~AxiStagnBVP() {}
+
+    // specify guesses for the initial values. These can be anything
+    // that leads to a converged solution.
+    virtual doublereal initialValue(int n, int j) { 
+        switch (n) {
+        case 0:
+            return 0.1*z(j);
+        case 1:
+            return 0.5*z(j);
+        default:
+            return 0.0;
+        }
+    }
+
+    // Specify the residual. This is where the ODE system and boundary 
+    // conditions are specified. The solver will attempt to find a solution
+    // x so that this function returns 0 for all n and j.
+    virtual doublereal residual(doublereal* x, int n, int j) {
+
+        // if n = 0, return the residual for the first ODE
+        if (n == 0) {
+            if (isLeft(j))  // here we specify zeta(0) = 0
+                return zeta(x,j);
+            else
+                // this implements d(zeta)/dz = u
+                return (zeta(x,j) - zeta(x,j-1))/(z(j)-z(j-1)) - u(x,j);
+        }
+        // if n = 1, then return the residual for the second ODE
+        else {
+            if (isLeft(j))   // here we specify u(0) = 0
+                return u(x,j);
+            else if (isRight(j)) // and here we specify u(L) = 1
+                return u(x,j) - 1.0;
+            else
+                // this implements the 2nd ODE
+                return cdif2(x,1,j) + 0.5*zeta(x,j)*centralFirstDeriv(x,1,j);
+        }
+    }
+
+
+private:
+
+    // for convenience only. Note that the compiler will inline these.
+    double zeta(double* x, int j) { return value(x,0,j); }
+    double u(double* x, int j) { return value(x,1,j); }
+
+};
+
+
+int main() {
+    try {
+ 
+       // Specify a problem on (0,10), with an initial uniform grid of
+       // 6 points.
+        AxiStagnBVP eqs(6, 10.0);
+        // Solve the equations, refining the grid as needed, and print lots of diagnostic output (loglevel = 4)
+        eqs.solve(4);
+        // write the solution to a CSV file.
+        eqs.writeCSV();
+        return 0;
+    }
+    catch (CanteraError) {
+        showErrors(cerr);
+        return -1;
+    }
+}
+
+

apps/bvp/BoundaryValueProblem.h

@@ -8,6 +8,7 @@
 #include <cantera/Cantera.h>
 #include <cantera/onedim.h>
 using namespace Cantera;
+using namespace std;
 
 /// Namespace for the boundary value problem package.
 namespace BVP {

apps/bvp/stagnation.mak

@@ -0,0 +1,86 @@
+#!/bin/sh
+
+# This Makefile builds a C++ application that uses Cantera.  By
+# default, the main program file is 'demo.cpp,' which prints out some
+# properties of a reacting gas mixture. 
+
+# To build program 'demo', simply type 'make', or 'make -f <this
+# file>' if this file is named something other than 'Makefile.'  
+
+# Once you have verified that the demo runs, edit this file to replace
+# object file 'demo.o' with your own object file or files. 
+
+
+#------------------------  edit this block ---------------------------------
+
+# the name of the executable program to be created
+PROG_NAME = stagnation.x
+
+# the object files to be linked together. 
+OBJS = AxiStagnBVP.o
+
+# additional flags to be passed to the linker. If your program
+# requires other external libraries, put them here
+LINK_OPTIONS = -L/usr/local/lib -framework Accelerate 
+
+#---------------------------------------------------------------------------
+# You probably don't need to edit anything below.
+
+# the C++ compiler
+CXX = g++
+
+# C++ compile flags
+CXX_FLAGS = -O3 -Wall -fPIC 
+
+# external libraries
+EXT_LIBS =  -luser -loneD -lzeroD -lequil -lkinetics -ltransport -lthermo -lctnumerics -lcvode -lctbase  -lctmath -ltpx -lctf2c -lconverters -lctcxx
+
+# Ending C++ linking libraries
+LCXX_END_LIBS = -lm
+
+# the directory where the Cantera libraries are located
+CANTERA_LIBDIR=/Applications/Cantera/lib
+
+# the directory where Cantera include files may be found.
+CANTERA_INCDIR=/Applications/Cantera/include
+
+# flags passed to the C++ compiler/linker for the linking step
+LCXXFLAGS = -L$(CANTERA_LIBDIR)  -O3 -Wall -fPIC
+
+# how to compile C++ source files to object files
+.cpp.o:
+	$(CXX) -c $< -I$(CANTERA_INCDIR) $(CXX_FLAGS)
+
+PROGRAM = $(PROG_NAME)$(EXE_EXT)
+
+DEPENDS = $(OBJS:.o=.d)
+
+all: $(PROGRAM)
+
+$(PROGRAM): $(OBJS)
+	$(CXX) -o $(PROGRAM) $(OBJS) $(LCXXFLAGS)\
+        $(CANTERA_LIBS) $(LINK_OPTIONS) $(EXT_LIBS) \
+          $(LCXX_END_LIBS)
+
+%.d:
+	g++ -MM -I$(CANTERA_INCDIR) $*.cpp > $*.d
+
+clean:
+	$(RM) $(OBJS) $(PROGRAM)
+
+depends: $(DEPENDS)
+	cat *.d > .depends
+	$(RM) $(DEPENDS) 
+
+TAGS: 
+	etags *.h *.cpp
+
+ifeq ($(wildcard .depends), .depends)
+include .depends
+endif
+
+
+
+
+
+