test for n2ls

README.md

@@ -14,19 +14,18 @@ but once you have homebrew installed and setup, here is the short version:
 
 
     brew install git
-    brew install python
-    brew install gfortran
-    pip install numpy
-    pip install scipy
-    brew install freetype
-    pip install matplotlib
+    brew tap Homebrew/python
+    brew install scipy
+
 
 
 ### Cantera
-In addition, this plugin requires [Cantera](https://code.google.com/p/cantera/) 
-and the python wrapper for it. You can [compile cantera from scratch](http://cantera.github.io/docs/sphinx/html/compiling.html), 
-or follow the instructions below for a bit easier route. 
+In addition, this plugin requires [Cantera](http://sourceforge.net/projects/cantera/files/cantera/2.0.2/cantera-2.0.2.tar.gz/download)
+and the python wrapper for it. Because of some backwards incompatible changes made to the most recent 
+versions of Cantera, you need to install version 2.0.2, and not the latest version. 
 
+You can [compile cantera from scratch](http://cantera.github.io/docs/sphinx/html/compiling.html), 
+or follow the instructions below for a bit easier route on windows. 
 
 
 #### Windows
@@ -40,10 +39,47 @@ release version of it!
 https://code.google.com/p/cantera/wiki/WindowsInstallation
 
 #### Mac OS-X
-Assuming you've used homebrew to get OpenMDAO setup, then just use it to install Cantera too! 
 
+You'll need to compile from source here. Unfortunately, installing Cantera on OS-X is a bit of a mess.  
+You can use homebrew to get some of the pre-requisites, but at least one of them needs to be downloaded manually first
+
+first 
+    brew tap homebrew/science
+
+Next go to the [sundials page](http://computation.llnl.gov/casc/sundials/download/download.php) and download 
+sundials-2.5.0.tar.gz. You need to manually put the zip file you just downloaded in 
+
+    /Library/Caches/Homebrew/sundials-2.5.0.tar.gz
+
+Now you're ready to install the cantera pre-reqs: 
+
+    brew install sundials
+    brew install scons
+
+After that, go unzip the [Cantera source code](http://sourceforge.net/projects/cantera/files/cantera/2.0.2/cantera-2.0.2.tar.gz/download) code you downwloaded. GO into that folder and edit a file called cantera.conf with the following lines: 
+
+    CXX = 'llvm-g++'
+    CC = 'llvm-gcc'
+    python_package = 'full'
+    f90_interface = 'n'
+    cxx_flags = '-ftemplate-depth-128 -DGTEST_USE_OWN_TR1_TUPLE=1'
+
+Save that file, then run 
+
+    scons build
+
+If that works, you'll see something like this in the terminal 
+
+    *******************************************************
+    Compilation completed successfully.
+
+    - To run the test suite, type 'scons test'.
+    - To install, type '[sudo] scons install'.
+    *******************************************************
 
-    brew install cantera
+To finish the install, 
+
+   scons install
 
 
 #### Linux

src/pycycle/api.py

@@ -7,3 +7,4 @@
 from heat_exchanger import HeatExchanger
 from cycle_component import CycleComponent
 from flowstation import FlowStation, FlowStationVar
+from thermal import CEA_tp

src/pycycle/test/openmdao_log.txt

@@ -5227,3 +5227,7 @@ Sep 17 13:39:37 E root: Variable 'sub_or_super' must be in ('sub', 'super'), but
 
 *********** BEGIN NEW LOG ************** (2013-10-25 16:27:10.266000) PID=2308
 
+
+
+*********** BEGIN NEW LOG ************** (2014-10-21 10:31:31.399637) PID=52802
+

src/pycycle/test/test_thermal.py

@@ -0,0 +1,57 @@
+import unittest
+
+from openmdao.main.api import set_as_top, Assembly
+from openmdao.util.testutil import assert_rel_error
+
+import numpy as np
+from pycycle.api import CEA_tp
+
+def rel_error(rhs, lhs, tol=0.01):
+    return np.linalg.norm(rhs - lhs) / np.linalg.norm(rhs) < tol
+
+
+class StartTestCase(unittest.TestCase):
+
+    def test_n2ls_4000(self):
+        comp = CEA_tp()
+
+        this_nguess = np.array([ 0.02040748,  0.00231478,  0.01020431,  0.03292581])
+
+        this_ch = np.array([[  2.27222641e-02,   2.50370446e-02,   2.27222641e-02],
+         [  2.50370446e-02,   7.04838323e-02,   4.54456580e-02],
+         [  2.27222641e-02,   4.54456580e-02,   7.59390390e-07]])
+
+        this_rhs = np.array([-0.81078449, -1.59262881, -1.14346367])
+
+        this_mu = np.array([-34.02175706, -50.32264785, -32.60178159])
+
+        comp.nguess = this_nguess
+        comp.T = 4000
+        comp.P = 1.034210
+        comp.run()
+
+        rel_error(comp.rhs, this_rhs, 0.001)
+        rel_error(comp.ch, this_ch, 0.001)
+        rel_error(comp.mu, this_mu, 0.001)
+
+    def test_n2ls_1500(self):
+        comp = CEA_tp()
+
+        this_nguess = np.array([6.58936767e-06,   2.27158983e-02,   6.24254002e-06,   2.27252139e-02])
+
+        this_ch = np.array([[  2.27224876e-02,   4.54383859e-02,   2.27224876e-02],
+                     [  4.54383859e-02,   9.08951525e-02,   4.54508710e-02],
+                     [  2.27224876e-02,   4.54508710e-02,   3.51630007e-06]])
+
+        this_rhs = np.array([-1.40217711, -2.80452075, -1.40240467])
+
+        this_mu = np.array([-43.73879404, -61.71398144, -35.95037481])
+
+        comp.nguess = this_nguess
+        comp.T = 1500
+        comp.P = 1.034210
+        comp.run()
+
+        rel_error(comp.rhs, this_rhs, 0.001)
+        rel_error(comp.ch, this_ch, 0.001)
+        rel_error(comp.mu, this_mu, 0.001)

src/pycycle/thermal.py

@@ -0,0 +1,112 @@
+import numpy as np 
+
+from openmdao.main.api  import Assembly,Component
+from openmdao.lib.datatypes.api import Float, Array
+
+
+R =1.987 #universal gas constant
+
+a=np.array([
+[ #CO
+4.619197250e+05,-1.944704863e+03,5.916714180e+00,-5.664282830e-04,
+1.398814540e-07,-1.787680361e-11,9.620935570e-16,-2.466261084e+03,
+-1.387413108e+01],
+[#CO2
+1.176962419e+05,-1.788791477e+03,8.291523190e+00,-9.223156780e-05,
+4.863676880e-09,-1.891053312e-12,6.330036590e-16,-3.908350590e+04,
+-2.652669281e+01],
+[ #O2
+-1.037939022e+06,2.344830282e+03,1.819732036e+00,1.267847582e-03,
+-2.188067988e-07,2.053719572e-11,-8.193467050e-16,-1.689010929e+04,
+1.738716506e+01
+]])
+wt_mole = np.array([ 28.01, 44.01, 32. ])    
+element_wt = [ 12.01, 16.0 ]
+aij = np.array([ [1,1,0], [1,2,2] ])
+
+
+_num_element = 2
+_num_react = 3
+
+#pre-computed constants used in calculations
+_aij_prod = np.empty((_num_element,_num_element, _num_react))
+for i in range( 0, _num_element ):
+    for j in range( 0, _num_element ):
+        _aij_prod[i][j] = aij[i]*aij[j]
+
+_aij_prod_deriv = np.zeros((_num_element**2,_num_react))
+for k in xrange(_num_element**2):
+    for l in xrange(_num_react):
+        i = k/_num_element
+        j = np.mod(k,_num_element)
+        _aij_prod_deriv[k][l] = _aij_prod[i][j][l]
+
+
+class CEA_tp(Component):
+
+    nguess = Array(iotype="in", desc="molar concentration of the mixtures, \
+                   last element is the total molar concentration")
+    T = Float(iotype="in", desc="Temperature")
+    P = Float(iotype="in", desc="Pressure")
+
+    rhs = Array(np.zeros(_num_element + 1), iotype="out", desc="Right-hand side of resulting system")
+    mu = Array(iotype="out", desc="A chemical thing")
+    ch = Array(np.zeros((_num_element+1, _num_element+1)), iotype="out", 
+               desc="Left-hand side of resulting system")
+
+    def H0(self):
+        ai = a.T
+        T = self.T
+        return (-ai[0]/T**2 + ai[1]/T*np.log(T) + ai[2] + ai[3]*T/2. + ai[4]*T**2/3. + ai[5]*T**3/4. + ai[6]*T**4/5.+ai[7]/T)
+
+    def S0(self):
+        ai = a.T
+        T = self.T
+        return (-ai[0]/(2*T**2) - ai[1]/T + ai[2]*np.log(T) + ai[3]*T + ai[4]*T**2/2. + ai[5]*T**3/3. + ai[6]*T**4/5.+ai[8] )
+
+    def execute(self):
+        """
+        Maps molar concentrations to pi coefficients matrix 
+        and a right-hand-side
+        """ 
+
+        nj = self.nguess[:-1]
+        nmoles = self.nguess[-1]
+        b = np.zeros(_num_element)
+        b0 = np.zeros(_num_element)
+
+        #calculate mu for each reactant
+        self.mu = self.H0() - self.S0() + np.log(nj) + np.log(self.P/nmoles) #pressure in Bars
+
+        #calculate b_i for each element
+        for i in range( 0, _num_element ):
+            b[ i ] =  np.sum(aij[i]*nj) 
+
+        ##determine pi coef for 2.24, 2.56, and 2.64 for each element
+        for i in range( 0, _num_element ):
+            for j in range( 0, _num_element ):
+                tot = np.sum(_aij_prod[i][j]*nj)
+                self.ch[i][j] = tot
+
+
+        #determine the delta n coeff for 2.24, dln/dlnT coeff for 2.56, and dln/dlP coeff 2.64
+        #and pi coef for 2.26,  dpi/dlnT for 2.58, and dpi/dlnP for 2.66
+        #and self.rhs of 2.64
+
+        #determine the delta coeff for 2.24 and pi coef for 2.26\  
+        self.ch[_num_element,:-1]= b
+        self.ch[:-1,_num_element]= b
+
+        #determine delta n coef for eq 2.26
+        sum_nj = np.sum(nj)
+        self.ch[-1,-1] = sum_nj - nmoles
+
+        #determine right side of matrix for eq 2.24
+        for i in range(_num_element ):
+            sum_aij_nj_muj = np.sum(aij[i]*nj*self.mu)
+            self.rhs[i] = sum_aij_nj_muj + b0[i] - b[i]
+
+        #determine the right side of the matrix for eq 2.36
+        sum_nj_muj = np.sum(nj*self.mu)
+        self.rhs[-1] = nmoles - sum_nj + sum_nj_muj
+