isentropic compression process

[Shirin]

Hello,

I am trying to run your code in rcm-example.ipynb to see if I can get the same plot as you did. But I am getting this below error. I have installed Pyked and pyyaml.  

Traceback (most recent call last):

  File "E:/Compression process/Varing reactor volume.py", line 10, in <module>

    ck = ChemKED(dict_input=testfile_rcm)

  File "C:\Program Files\Python37\lib\site-packages\pyked\chemked.py", line 120, in __init__

    self.validate_yaml(self._properties)

  File "C:\Program Files\Python37\lib\site-packages\pyked\chemked.py", line 179, in validate_yaml

    validator = OurValidator(schema)

  File "C:\Program Files\Python37\lib\site-packages\cerberus\validator.py", line 181, in __init__

    self.schema = kwargs.get('schema', None)

  File "C:\Program Files\Python37\lib\site-packages\cerberus\validator.py", line 562, in schema

    self._schema = DefinitionSchema(self, schema)

  File "C:\Program Files\Python37\lib\site-packages\cerberus\schema.py", line 82, in __init__

    self.validate(schema)

  File "C:\Program Files\Python37\lib\site-packages\cerberus\schema.py", line 262, in validate

    self._validate(schema)

  File "C:\Program Files\Python37\lib\site-packages\cerberus\schema.py", line 277, in _validate

    if not self.schema_validator(schema, normalize=False):

  File "C:\Program Files\Python37\lib\site-packages\cerberus\validator.py", line 999, in validate

    self.__validate_unknown_fields(field)

  File "C:\Program Files\Python37\lib\site-packages\cerberus\validator.py", line 1031, in __validate_unknown_fields

    if not validator({field: value}, normalize=False):

  File "C:\Program Files\Python37\lib\site-packages\cerberus\validator.py", line 997, in validate

    self.__validate_definitions(definitions, field)

  File "C:\Program Files\Python37\lib\site-packages\cerberus\validator.py", line 1066, in __validate_definitions

    result = validate_rule(rule)

  File "C:\Program Files\Python37\lib\site-packages\cerberus\validator.py", line 1041, in validate_rule

    return validator(definitions.get(rule, None), field, value)

  File "C:\Program Files\Python37\lib\site-packages\cerberus\validator.py", line 1436, in _validate_schema

    self.__validate_schema_mapping(field, schema, value)

  File "C:\Program Files\Python37\lib\site-packages\cerberus\validator.py", line 1450, in __validate_schema_mapping

    if not validator(value, update=self.update, normalize=False):

  File "C:\Program Files\Python37\lib\site-packages\cerberus\validator.py", line 997, in validate

    self.__validate_definitions(definitions, field)

  File "C:\Program Files\Python37\lib\site-packages\cerberus\validator.py", line 1066, in __validate_definitions

    result = validate_rule(rule)

  File "C:\Program Files\Python37\lib\site-packages\cerberus\validator.py", line 1040, in validate_rule

    validator = self.__get_rule_handler('validate', rule)

  File "C:\Program Files\Python37\lib\site-packages\cerberus\validator.py", line 362, in __get_rule_handler

    "domain.".format(rule, domain)

RuntimeError: There's no handler for 'isvalid_history' in the 'validate' domain.

>>> 

[Bryan W. Weber]

Hi Shirin,

There is currently a problem with PyKED and one of the libraries it uses called Cerberus. You need to install version 1.1 or lower of Cerberus. This should be automatically handled by the installer, but it doesn't work sometimes for some reason. Since this isn't really about Cantera, I'll reply to you off the list to do some more troubleshooting.

Best,

Bryan

[Shirin]

Hi Bryan,

I am trying to write a code for isentropic compression process; I used your example in https://github.com/pr-omethe-us/PyKED/blob/master/docs/rcm-example.ipynb.

But my P-t plot does not make sense, I am sure my time is not write and I don't know how to use time here in this code and in Python.  I am very new in Python. Please help me with this. Thanks,

my code is 

import cantera as ct

import numpy as np

from pyked import ChemKED

temperature = 300 #K

pressure = 1e5  #pa

gas = ct.Solution('gri30.xml')

gas.TP = temperature, pressure

gas.set_equivalence_ratio(1, 'h2:1', 'n2:3.76, o2:1.0')

#gas.set_equivalence_ratio(1, 'ch4:0.09506', 'n2:0.71483, o2:0.19011')

reac = ct.IdealGasReactor(gas)

env = ct.Reservoir(ct.Solution('air.xml'))

#L=2a

#R=l/a

start_crank_angle =180

crank_radius=3.5#0.062

rod_length=5#0.192

rod_radius_ratio= rod_length/crank_radius

rev_per_minute=1500

stroke_length=7#0.124

netw = ct.ReactorNet([reac])

t = 0.0;

dt = 1.0e-5;

#times = np.zeros(100)

tsteps=100

for n in range(1,tsteps):

#for n in range(1,101): #1:100

  t = t + dt;

  netw.advance(t)

#tim[n] = time(netw)

  

# Angular velocity, rad/s

omega = rev_per_minute*np.pi/30

# Start angle, rad

start_crank_rad = start_crank_angle/180.0*np.pi  #deg to rad

#theta = start_crank_rad - omega * time   #rad-rad/s*s

def __call__(self, time):

  theta = self.start_crank_rad - self.omega * time

  velocity= (self.omega*self.stroke_length/2*np.sin(theta) *

                (1 + np.cos(theta)/np.sqrt(self.rod_radius_ratio**2 -

                                           np.sin(theta)**2)))

  return velocity

  

#netw.set_max_time_step(np.min(np.diff(t)))

time = []

temperature = []

pressure = []

#volume = []

mass_fractions = []

while netw.time < 0.05:

    time.append(netw.time)

    temperature.append(reac.T)

    pressure.append(reac.thermo.P)

    #volume.append(reac.volume)

    mass_fractions.append(reac.Y)

    netw.step()

print (time)

print (pressure)

    #%matplotlib notebook

import matplotlib.pyplot as plt

plt.figure()

plt.plot(time, pressure)

plt.ylabel('Pressure [Pa]')

#plt.xlabel('Time [s]')

plt.show()

[Shirin]

Sorry I forgot to add the wall in my code 

w=ct.Wall(reac, env, velocity);

On Friday, November 15, 2019 at 11:27:55 AM UTC-5, Bryan W. Weber wrote:

[Bryan W. Weber]

Hi Shirin,

Thanks for your patience here. There's a few things. First, the example that you copy-pasted into your post here will attempt to simulate an IC engine, including the compression and expansion strokes, so that's not what you want. What you need to do is define a velocity using the volume history that you have.

Assuming that your volume history looks like

Time Volume

0.0. 0.0

1.0. 0.6

2.0. 1.1

...

(Obviously, those numbers are made up), you need to compute the velocity of the piston. Since the piston has constant area, the velocity is the derivative of the volume with respect to time, dV/dt. You can calculate this by finite difference, (V_2 - V_1)/(t_2 - t_1). You first need to load the file that has the data into two arrays, "volume" and "time". Then, you can calculate the velocity as follows:

import numpy as np
velocity = np.diff(volume) / np.diff(time)

Now, Cantera needs a way to determine the velocity of the piston at a particular time during the simulation. My suggestion is to do this with a function:

def velocity_function(sim_t):
    """Compute the velocity of the piston at the simulation time, sim_t"""
    if t > self.time[-1]:
         # Return 0 if the simulation time is greater than a velocity we've computed
        return 0
    else:
        # prev_time_point is the previous value in the time array
        # after the current simulation time
        prev_time_point = self.time[self.time <= t][-1]
        # index is the index of the time array where
        # prev_time_point occurs
        index = np.where(self.time == prev_time_point)[0][0]
        return velocity[index]

 Then, you need to create your simulation:

gas = ct.Solution(input_file, transport_model=None)
gas.TPX = 325, 101325, YOUR_COMPOSITION
env = ct.Reservoir(ct.Solution("air.xml", transport_model=None))
reac = ct.IdealGasReactor(gas)
ct.Wall(reac, env, velocity=velocity_function)

netw = ct.ReactorNet([reac])

while netw.time < end_time:
    netw.step()

Hope that helps!

Bryan

[Govind Bubloo]

Hi Bryan

This post is exactly what i am looking for. I need to compute the velocity from the Time-VOlume history. 

Here is my code.

# 1st i just read the Time_volume history from the text file and stored it in a keywords dictionary.

Exp_time = []

Exp_vol = []

keywords = {}

for line in open("test.txt", "r").readlines():

    line = line.split()

    if len(line)>1 and line[0] == 'VPRO':

        Exp_time.append(line[1])

        Exp_vol.append(line[2])

        keywords = {'vproTime': Exp_time, 'vproVol': Exp_vol}

# Then am computing the velocity for the wall from the volume profile

class VolumeProfile():

       volume = []

       time = []

       def __init__(self, keywords):

              self.time = np.array(keywords['vproTime'])

              self.volume = np.array(keywords['vproVol'])/keywords['vproVol'][0]

              self.velocity = np.diff(self.volume)/np.diff(self.time)

              self.velocity = np.append(self.velocity, 0)

              

       def __call__(self, t):

              if t < self.time[-1]:

                     prev_time_point = self.time[self.time <= t][-1]

                     index = np.where(self.time == prev_time_point)[0][0]

                     return self.velocity[index]

              else:

                     return 0

net = ct.Wall(reac, env, velocity=VolumeProfile(keywords))

But i am getting this error.

((net = ct.Wall(reac, env, velocity=VolumeProfile(keywords))

 self.volume = np.array(keywords['vproVol'])/keywords['vproVol'][0]

TypeError: ufunc 'true_divide' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe'' ))

Could you please explain what mistake i did. Thank you in advance

  

        

[Bryan W. Weber]

Govind,

This is because keywords["vproVol"][0] is a string type, whereas np.array(...) is an array of strings. I think you need to wrap the line[X] in the init function with a float() call to fix this.

Best,

Bryan

[Govindarajan Kumar]

Hi Bryan

Yeah I did it. Thank you so much for your help. You are genius man

Govind

--
You received this message because you are subscribed to the Google Groups "Cantera Users' Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cantera-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cantera-users/a4b5301d-368b-4788-a13d-d341ba4d4592%40googlegroups.com.

[Shirin]

Hi Bryan,

Thank you for your reply. I have a question here, what is the end_time here?

Best Regards,

Shirin

[Bryan W. Weber]

Hi Shirin,

The end_time is whatever time you want to end the simulation, say 100 ms or 400 ms. It is a variable you have to define.

The cases where self was included was a typo, I was copy-pasting from some other code that used a class. In those cases, you just need to replace them with whatever array (or column in an array) is storing the time values. Also, t should become sim_t. Sorry for the typos!

You do not need to add a number for the compression time.

Best,

Bryan

[Shirin Jouzdani]

Hi Bryan,

I have one more question here, I need to define the sim_t right? if so, is the below part correct?

time = 0.0;
dt = 1.0e-6;
times = np.zeros(1000)
tsteps=1000
for n in range(1,tsteps):
 time = time + dt;
 netw.advance(time)
 times[n] = netw.time

I really appreciate your help,

Shirin

--
You received this message because you are subscribed to the Google Groups "Cantera Users' Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cantera-user...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/cantera-users/28d900f5-e041-425a-8dcd-6a0fddf270bc%40googlegroups.com.

[Bryan W. Weber]

Hi Shirin,

I would improve your code like this:

time = 0.0
dt = 1.0e-6

tsteps = 1000
times = np.zeros(tsteps)
for n in range(tsteps):

    time = time + dt
    netw.advance(time)
    times[n] = netw.time

I removed the semicolons and the loop now goes from 0->999 so that you do not end up with an index error when assigning into the times array. Python array indices start at 0, not 1 as they do in MATLAB.

Best,

Bryan

To unsubscribe from this group and stop receiving emails from it, send an email to cantera-users+unsubscribe@googlegroups.com.

[Shirin]

Hi Bryan,

Thank you again for your help and time.

Just to make sure I understand this, the time in your velocity_function is the time values in volume history, correct?

if so, what you just improve is the sim_t?

I think I am confused about the time and sim_t.

Sorry for asking many questions and thank you.

Shirin

On Tuesday, January 7, 2020 at 3:03:47 PM UTC-5, Bryan W. Weber wrote:

Hi Shirin,

To unsubscribe from this group and stop receiving emails from it, send an email to canter...@googlegroups.com.

[Bryan W. Weber]

Hi Shirin,

Yes, the "time" in the function is the array of time values in the volume trace. The sim_t is an argument that is automatically passed to the function by Cantera when the simulation is running. It is equal to the current value of the time during the simulation.

Best,

Bryan

[Shirin]

Hi Bryan,

Thank you again. I ran my simulation but the pressure is constant and does not change. My code is as below:

volume=first_col

time=second_col

velocity = np.diff(volume) / np.diff(time)

velocity = np.append(velocity, 0)

print("Velocity",velocity)

gas = ct.Solution('mech_DME_2000.xml')

T_in = 300 #K

P_in = 1e5  #pa

# to set the composition of the gas

fuel_species = 'ch4'

ifuel = gas.species_index(fuel_species)

io2 = gas.species_index('o2')

in2 = gas.species_index('n2')

#create an array to hold the composition

X = np.zeros(gas.n_species)

air_N2_O2_molar_ratio = 3.76

nu=2

phi=1

X[ifuel] =1

X[io2] =nu/phi

X[in2] = X[io2]*air_N2_O2_molar_ratio

#gas.X={'CH4':1, 'O2':2/phi, 'N2':2*3.76/phi}

gas.TPX = T_in, P_in, X

gas()

env = ct.Reservoir(ct.Solution('air.xml'))

reac = ct.IdealGasReactor(gas)

netw = ct.ReactorNet([reac])

time = 0.0

dt = 1.0e-6

tsteps=2000

times = np.zeros(tsteps)

for n in range(tsteps):

 time = time + dt;

 netw.advance(time)

 times[n] = netw.time

def velocity_function(sim_t):

    """Compute the velocity of the piston at the simulation time, sim_t"""

    if sim_t > time[-1]:

         # Return 0 if the simulation time is greater than a velocity we've computed

        return 0

    else:

        # prev_time_point is the previous value in the time array

        # after the current simulation time

        prev_time_point = time[time <= sim_t][-1]

        # index is the index of the time array where

        # prev_time_point occurs

        index = np.where(time == prev_time_point)[0][0]

        return velocity[index]

ct.Wall(reac, env, velocity=velocity_function)

time = []

temperature = []

pressure = []

volume = []

while netw.time < 0.04:

    time.append(netw.time)

    temperature.append(reac.T)

    pressure.append(reac.thermo.P)

    volume.append(reac.volume)

    netw.step()

print (pressure)

    

#%matplotlib notebook

import matplotlib.pyplot as plt

plt.figure()

plt.plot(time, pressure)

plt.ylabel('Pressure [Pa]')

#plt.xlabel('Time [s]')

plt.show()

Thanks,

Shirin

Traceback (most recent call last):

  File "E:/Compression process/isentropic_compression_ic_new.py", line 92, in <module>

    netw.step()

  File "interfaces\cython\cantera\reactor.pyx", line 862, in cantera._cantera.ReactorNet.step

cantera._cantera.CanteraError: 

***********************************************************************

CanteraError thrown by CVodesIntegrator::step:

CVodes error encountered. Error code: -9

Exceptions caught during RHS evaluation:

***********************************************************************

Exception thrown by Python callback function:

TypeError: '<=' not supported between instances of 'list' and 'float'

***********************************************************************

The right-hand side routine failed at the first call.

Components with largest weighted error estimates:

26: inf

22: 5.95571e+304

28: 2.63428e+304

24: 6.26057e+265

39: 1.6805e+256

36: 1.3336e+256

47: 3.52678e+248

14: 1.92643e+238

44: -7.22948e+236

38: 1.91295e+229

***********************************************************************

[Bryan W. Weber]

Hi Shirin,

You have one too many instances of running the calculation. You need to delete this code:

time = 0.0

dt = 1.0e-6

tsteps=2000

times = np.zeros(tsteps)

for n in range(tsteps):

 time = time + dt;

 netw.advance(time)

 times[n] = netw.time

Best,

Bryan

[Shirin]

Hi Bryan,

Thank you for your reply. When I delete that part, I get the below error:

Traceback (most recent call last):

  File "E:/Compression process/isentropic_compression_ic_new.py", line 101, in <module>

    netw.step()

  File "interfaces\cython\cantera\reactor.pyx", line 862, in cantera._cantera.ReactorNet.step

cantera._cantera.CanteraError: 

***********************************************************************

CanteraError thrown by CVodesIntegrator::step:

CVodes error encountered. Error code: -9

Exceptions caught during RHS evaluation:

***********************************************************************

Exception thrown by Python callback function:

TypeError: '<=' not supported between instances of 'list' and 'float'

***********************************************************************

The right-hand side routine failed at the first call.

Components with largest weighted error estimates:

3: -nan

81: 5.38451e+243

28: -3.43331e-12

79: 5.88753e-56

70: 3.19164e-75

65: 7.43111e-85

59: 4.26401e-99

60: 4.26393e-99

34: 3.53633e-123

38: 9.22873e-128

***********************************************************************

Thanks,

Shirin

[Bryan W. Weber]

Hi Shirin,

You didn't show how you loaded the volume trace. You have to have time be a NumPy array, it can't be a plain Python list.

Best,

Bryan

[Shirin]

Hi Bryan,

Thank you again.I have my volume and time in a excel and I used np.array for both of them but steel I get the same error. 

import cantera as ct

import numpy as np

import xlrd 

  

# Give the location of the file 

loc = ("vol_time_engine.xlsx") 

  

# To open Workbook 

wb = xlrd.open_workbook(loc) 

sheet = wb.sheet_by_index(0) 

  

# row 0 and column 0 

first_row= sheet.row_values(0)

first_col=sheet.col_values(0)

print("First Row ",first_row) 

print("First Column ",first_col )

print("Printing first column elements")

for element in first_col:

    print(element)

second_col=sheet.col_values(1)

print("Second Column ",second_col )

print("Printing second column elements")

for element in second_col:

    print(element)

volume= np.array(first_col)

time = np.array(second_col)

[Shirin]

Hi Bryan,

In the line :

ct.Wall(reac, env, velocity=velocity_function)

don't we need to add input for velocity_function as it is in the definition or something else?

Thanks,

Shirin

[Govindarajan Kumar]

Hi Bryan

Greetings of the day.

I have to make code for finding the Ignition delay time. I was able to make the graph of Volume over time. But dont know how to proceed after this. Can you please help me

from __future__ import division
import sys

import cantera as ct
import numpy as np

import matplotlib.pyplot as plt



#Reading the temperature, Pressure and ratio of species from the RCM input files(VOlume-TIme history file)

for line in open("test.txt","r").readlines():

       data = line.split()
       if len(data)>1 and data[1] == 'ch4':
              a = float(data[2])
       if len(data)>1 and data[1] == 'ar':
              b = float(data[2])
       if len(data)>1 and data[1] == 'o2':
              c = float(data[2])
       if len(data)>1 and data[1] == 'h2o':
              d = float(data[2])
       if len(data)>1 and data[0] == 'TEMP':
           T_initial = float(data[1])
       if len(data)>1 and data[0] == 'PRES':
              P_initial = float(data[1])*101325
gas=ct.Solution("projectMech.cti")
gas.TPX = T_initial, P_initial, {'CH4':a, 'Ar':b, 'O2':c, 'H2O':d}
#print(gas())


reac = ct.IdealGasReactor(gas)

env = ct.Reservoir(ct.Solution('air.xml'))

#Reading the Experimental Volume and TIme histories from the same input file and assigning it to a Dictionay Keywords

Exp_time = []
Exp_vol = []
keywords = {}
for line in open("test.txt", "r").readlines():
    line = line.split()
    if len(line)>1 and line[0] == 'VPRO':

        Exp_time.append(float(line[1]))
        Exp_vol.append(float(line[2]))
        keywords = {'Time': Exp_time, 'Vol': Exp_vol}
#print(Exp_vol)

class VolumeProfile():
       def __init__(self, keywords):
              self.time = np.array(keywords['Time'])
              self.volume = np.array(keywords['Vol'])/keywords['Vol'][0]
              self.velocity = np.diff(Exp_vol)/np.diff(Exp_time)
              self.velocity = np.append(self.velocity, 0)
       def __call__(self, t): #Return the velocity when called during a time step.
              if t < self.time[-1]:
                     prev_time_point = self.time[self.time <= t][-1] #prev_time_point is the previous value in the time array after the current simulation time
                     index = np.where(self.time == prev_time_point)[0][0] #index is the index of the time array where prev_time_point occurs

                     return self.velocity[index]
              else:
                     return 0
net = ct.Wall(reac, env, velocity = VolumeProfile(keywords))

netw = ct.ReactorNet([reac])

time = []
temperature = []
pressure = []
volume = []

mass_fractions = []
while netw.time < 0.6:
       time.append(netw.time)
       temperature.append(reac.T)
       pressure.append(reac.thermo.P)
       volume.append(reac.volume)
       mass_fractions.append(reac.Y)
       netw.step()

  
       
#plt.plot(time, pressure)
#plt.ylabel('Pressure [Pa]')
#plt.xlabel('Time [s]')
#plt.plot(time, temperature)
#plt.ylabel('temperature [K]')
#plt.xlabel('Time [s]')
#plt.plot(Exp_time, Exp_vol, label='Experimental Volume', color='k', linestyle=':')
#plt.plot(time, volume, label='Simulated Volume',linestyle='-', markersize=15)
#plt.legend()
#plt.ylabel('Volume [cm^3]')
#plt.xlabel('Time [s]')
#plt.xlim(0,1)
#plt.ylim(0,1.1)
#plt.show()

--

You received this message because you are subscribed to the Google Groups "Cantera Users' Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email to cantera-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cantera-users/169d4265-bab4-4eb3-82c0-fc8e3915e397%40googlegroups.com.

[Shirin]

Hi Bryan,

It seems the error that I get it is related to the line :ct.Wall(reac, env, velocity=velocity_function) .

Because when I change my code to the class VolumeProfile() version it gives me different error as below: (I also attached my code following the error)

File "interfaces\cython\cantera\reactor.pyx", line 862, in cantera._cantera.ReactorNet.step

cantera._cantera.CanteraError: 

***********************************************************************

CanteraError thrown by CVodesIntegrator::step:

CVodes error encountered. Error code: -3

At t = 0.00039434 and h = 3.2007e-26, the error test failed repeatedly or with |h| = hmin.

Components with largest weighted error estimates:

2: -104.391

34: 2.3888e-07

35: 1.28338e-08

28: -3.77103e-09

26: -2.32004e-09

25: 9.15911e-10

46: -4.49081e-10

20: 3.79094e-10

17: 1.95713e-10

27: 1.52556e-10

***********************************************************************

My code now is as below:

import cantera as ct

import numpy as np

import xlrd 

  

# Give the location of the file 

loc = ("vol_time_engine.xlsx") 

  

# To open Workbook 

wb = xlrd.open_workbook(loc) 

sheet = wb.sheet_by_index(0) 

  

# row 0 and column 0 

first_row= sheet.row_values(0)

first_col=sheet.col_values(0)

print("First Row ",first_row) 

print("First Column ",first_col )

print("Printing first column elements")

for element in first_col:

    print(element)

second_col=sheet.col_values(1)

print("Second Column ",second_col )

print("Printing second column elements")

for element in second_col:

    print(element)

exp_volume= np.array(first_col)

exp_time = np.array(second_col)

#print (time)

#velocity = np.diff(volume) / np.diff(time)

#velocity = np.append(velocity, 0)

#print("Velocity",velocity)

gas = ct.Solution('mech_DME_2000.xml')

T_in = 300 #K

P_in = 1e5  #pa

# to set the composition of the gas

fuel_species = 'ch4'

ifuel = gas.species_index(fuel_species)

io2 = gas.species_index('o2')

in2 = gas.species_index('n2')

#create an array to hold the composition

X = np.zeros(gas.n_species)

air_N2_O2_molar_ratio = 3.76

nu=2

phi=1

X[ifuel] =1

X[io2] =nu/phi

X[in2] = X[io2]*air_N2_O2_molar_ratio

#gas.X={'CH4':1, 'O2':2/phi, 'N2':2*3.76/phi}

gas.TPX = T_in, P_in, {'CH4':1, 'O2':2/phi, 'N2':2*3.76/phi}

gas()

env = ct.Reservoir(ct.Solution('air.xml'))

reac = ct.IdealGasReactor(gas)

netw = ct.ReactorNet([reac])

keywords = {'Time': exp_time, 'Vol': exp_volume}

class VolumeProfile():

       def __init__(self, keywords):

              self.time = np.array(keywords['Time'])

              self.volume = np.array(keywords['Vol'])/keywords['Vol'][0]

              self.velocity = np.diff(exp_volume)/np.diff(exp_time)

              self.velocity = np.append(self.velocity, 0)

       def __call__(self, t): #Return the velocity when called during a time step.

              if t < self.time[-1]:

                     prev_time_point = self.time[self.time <= t][-1] #prev_time_point is the previous value in the time array after the current simulation time

                     index = np.where(self.time == prev_time_point)[0][0] #index is the index of the time array where prev_time_point occurs

                     return self.velocity[index]

              else:

                     return 0

ct.Wall(reac, env, velocity=VolumeProfile(keywords))

time = []

temperature = []

pressure = []

volume = []

while netw.time < 0.02:

    time.append(netw.time)

    temperature.append(reac.T)

    pressure.append(reac.thermo.P)

    volume.append(reac.volume)

    netw.step()

print (pressure)

    

#%matplotlib notebook

import matplotlib.pyplot as plt

plt.figure()

plt.plot(time, pressure)

plt.ylabel('Pressure [Pa]')

#plt.xlabel('Time [s]')

plt.show()

    

Could you please help me with the new error?

Thanks,

Shirin

[Shirin]

Hello Bryan,

I have posted the error that I get in my previous post. I just wanted to let you know that I am storing my volume and time in an excel with two columns.

Thanks,

Shirin

[Shirin]

Hello,

I have a question here, what is VPRO?

Thanks,

Shirin

On Sat, Jan 11, 2020 at 8:49 PM Bryan W. Weber <bryan...@gmail.com> wrote:

Hi Shirin,

You didn't show how you loaded the volume trace. You have to have time be a NumPy array, it can't be a plain Python list.

Best,

Bryan

--
You received this message because you are subscribed to the Google Groups "Cantera Users' Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email to canter...@googlegroups.com.

[Shirin]

Hi Bryan,

I have been trying to find the error and changed some parts in my code, and now I am getting the new error. I attached the error and then my code. I think the error is now from reading my data. Please help me to find the error.

Thanks,

Shirin

Traceback (most recent call last):

  File "E:/Compression process/new_test.py", line 102, in <module>

    netw.step()

  File "interfaces\cython\cantera\reactor.pyx", line 862, in cantera._cantera.ReactorNet.step

cantera._cantera.CanteraError: 

***********************************************************************

CanteraError thrown by CVodesIntegrator::step:

CVodes error encountered. Error code: -9

Exceptions caught during RHS evaluation:

***********************************************************************

Exception thrown by Python callback function:

IndexError: index -1 is out of bounds for axis 0 with size 0

***********************************************************************

The right-hand side routine failed at the first call.

Components with largest weighted error estimates:

0: 0

1: 0

2: 0

3: 0

4: 0

5: 0

6: 0

7: 0

__________________________________________________________________________________________

import cantera as ct

import numpy as np

import sys

import xlrd 

  

# Give the location of the file 

loc = ("vol_time_engine.xlsx") 

  

# To open Workbook 

wb = xlrd.open_workbook(loc) 

sheet = wb.sheet_by_index(0) 

  

gas = ct.Solution('mech_DME_2000.xml')

T_in = 300 #K

P_in = 1e5  #pa

# to set the composition of the gas

fuel_species = 'ch4'

ifuel = gas.species_index(fuel_species)

io2 = gas.species_index('o2')

in2 = gas.species_index('n2')

#create an array to hold the composition

X = np.zeros(gas.n_species)

air_N2_O2_molar_ratio = 3.76

nu=2

phi=1

X[ifuel] =1

X[io2] =nu/phi

X[in2] = X[io2]*air_N2_O2_molar_ratio

#gas.X={'CH4':1, 'O2':2/phi, 'N2':2*3.76/phi}

gas.TPX = T_in, P_in, {'CH4':1, 'O2':2/phi, 'N2':2*3.76/phi}

gas()

env = ct.Reservoir(ct.Solution('air.xml'))

reac = ct.IdealGasReactor(gas)

Exp_time = []

Exp_vol = []

keywords={}

for i in range(0,sheet.nrows,1):

    Exp_vol.append(sheet.cell_value(i, 0))

    Exp_time.append(sheet.cell_value(i, 1))

    keywords = {'Time': Exp_time, 'Vol': Exp_vol}

   

print("Volume: ", len(Exp_vol))

print("Time: ", len(Exp_time))

print("keywords: ", keywords)

class VolumeProfile():

       def __init__(self, keywords):

              self.time = np.array(keywords['Time'])

              self.volume = np.array(keywords['Vol'])/keywords['Vol'][0]

              self.velocity = np.diff(Exp_vol)/np.diff(Exp_time)

              self.velocity = np.append(self.velocity, 0)

       def __call__(self, t): #Return the velocity when called during a time step.

              if t < self.time[-1]:

                     prev_time_point = self.time[self.time <= t][-1] #prev_time_point is the previous value in the time array after the current simulation time

                     index = np.where(self.time == prev_time_point)[0][0] #index is the index of the time array where prev_time_point occurs

                     return self.velocity[index]

              else:

                     return 0

net=ct.Wall(reac, env, velocity=VolumeProfile(keywords))

netw = ct.ReactorNet([reac])

time = []

temperature = []

pressure = []

volume = []

while netw.time < 0.02:

    time.append(netw.time)

    temperature.append(reac.T)

    pressure.append(reac.thermo.P)

    volume.append(reac.volume)

    netw.step()

print (pressure)

    

#%matplotlib notebook

import matplotlib.pyplot as plt

plt.figure()

plt.plot(time, pressure)

plt.ylabel('Pressure [Pa]')

#plt.xlabel('Time [s]')

plt.show()

    

[Shirin]

Hello,

I have been trying to find the isentropic compression process in ic engine. I have been getting an error.  I have tried everything to solve my problem but I can not find the error.

It seems it does not go to the next step and when I asked to print time and pressure and volume in while loop it only shows the initial value.

I attached my code and error, please help me to solve this problem. Thank you.

import cantera as ct
import numpy as np
import sys
import xlrd
 
# Give the location of the file
loc = ("vol_time_engine.xlsx")
 
# To open Workbook
wb = xlrd.open_workbook(loc)
sheet = wb.sheet_by_index(0)

 

gas = ct.Solution('gri30.xml')
print(gas.species_names)

T_in = 300 #K
P_in = 1e5  #pa
# to set the composition of the gas
fuel_species = 'ch4'
ifuel = gas.species_index(fuel_species)

print(ifuel)

io2 = gas.species_index('o2')

print(io2 )

in2 = gas.species_index('n2')

print(in2)
phi=1

gas.TPX = T_in, P_in, {'CH4':0.5, 'O2':2/phi, 'N2':2*3.76/phi}
gas()
print("GasTPX: ", gas.TPX)

reac = ct.IdealGasReactor(gas)
env = ct.Reservoir(ct.Solution('air.xml'))

Exp_time = []
Exp_vol = []
keywords={}

for i in range(1,sheet.nrows,1):

       Exp_vol.append(sheet.cell_value(i, 0))
       Exp_time.append(sheet.cell_value(i, 1))
       keywords = {'Time': Exp_time, 'Vol': Exp_vol}
   

print("Volume: ", Exp_vol)

print("Time: ", len(Exp_time))
print("keywords: ", keywords)




class VolumeProfile():
       def __init__(self, keywords):
              self.time = np.array(keywords['Time'])              
              self.volume = np.array(keywords['Vol'])/keywords['Vol'][0]              
              self.velocity = np.diff(Exp_vol)/np.diff(Exp_time)              
              self.velocity = np.append(self.velocity, 0)
       def __call__(self, t): #Return the velocity when called during a time step.
              if t < self.time[-1]:
                     prev_time_point = self.time[self.time <= t][-1] #prev_time_point is the previous value in the time array after the current simulation time                    
                     index = np.where(self.time == prev_time_point)[0][0] #index is the index of the time array where prev_time_point occurs
                     return self.velocity[index]
              else:
                     return 0

net=ct.Wall(reac, env, velocity=VolumeProfile(keywords))
netw = ct.ReactorNet([reac])

time = []
temperature = []
pressure = []
volume = []

while netw.time < 0.05:
    time.append(netw.time)
    print(time)
    temperature.append(reac.T)
    print(temperature)
    pressure.append(reac.thermo.P)
    print(pressure)
    volume.append(reac.volume)
    print(volume)

    netw.step()


print (pressure)
   
#%matplotlib notebook
import matplotlib.pyplot as plt

plt.figure()
plt.plot(time, pressure)
plt.ylabel('Pressure [Pa]')
#plt.xlabel('Time [s]')
plt.show()

 _________________________________________________________________________________

my error is:   

['H2', 'H', 'O', 'O2', 'OH', 'H2O', 'HO2', 'H2O2', 'C', 'CH', 'CH2', 'CH2(S)', 'CH3', 'CH4', 'CO', 'CO2', 'HCO', 'CH2O', 'CH2OH', 'CH3O', 'CH3OH', 'C2H', 'C2H2', 'C2H3', 'C2H4', 'C2H5', 'C2H6', 'HCCO', 'CH2CO', 'HCCOH', 'N', 'NH', 'NH2', 'NH3', 'NNH', 'NO', 'NO2', 'N2O', 'HNO', 'CN', 'HCN', 'H2CN', 'HCNN', 'HCNO', 'HOCN', 'HNCO', 'NCO', 'N2', 'AR', 'C3H7', 'C3H8', 'CH2CHO', 'CH3CHO']
13
3
47

  gri30:

       temperature             300  K
          pressure          100000  Pa
           density         1.13103  kg/m^3
  mean mol. weight         28.2116  amu

                          1 kg            1 kmol
                       -----------      ------------
          enthalpy     -1.2998e+05       -3.667e+06     J
   internal energy      -2.184e+05       -6.161e+06     J
           entropy          7088.5            2e+05     J/K
    Gibbs function     -2.2565e+06       -6.366e+07     J
 heat capacity c_p          1044.7        2.947e+04     J/K
 heat capacity c_v          749.94        2.116e+04     J/K

                           X                 Y          Chem. Pot. / RT
                     -------------     ------------     ------------
                O2       0.199601         0.226396         -26.2983
               CH4      0.0499002        0.0283761         -55.3337
                N2       0.750499         0.745228         -23.3333
     [  +50 minor]              0                0

GasTPX:  (300.0, 100000.0, array([0.       , 0.       , 0.       , 0.1996008, 0.       , 0.       ,
       0.       , 0.       , 0.       , 0.       , 0.       , 0.       ,
       0.       , 0.0499002, 0.       , 0.       , 0.       , 0.       ,
       0.       , 0.       , 0.       , 0.       , 0.       , 0.       ,
       0.       , 0.       , 0.       , 0.       , 0.       , 0.       ,
       0.       , 0.       , 0.       , 0.       , 0.       , 0.       ,
       0.       , 0.       , 0.       , 0.       , 0.       , 0.       ,
       0.       , 0.       , 0.       , 0.       , 0.       , 0.750499 ,
       0.       , 0.       , 0.       , 0.       , 0.       ]))
Volume:  
Time:  720
keywords:  
[0.0]
[300.0]
[99999.99999999997]
[1.0]

Traceback (most recent call last):

  File "E:/Compression process/compression.py", line 86, in <module>

    netw.step()
  File "interfaces\cython\cantera\reactor.pyx", line 862, in cantera._cantera.ReactorNet.step
cantera._cantera.CanteraError:
***********************************************************************
CanteraError thrown by CVodesIntegrator::step:
CVodes error encountered. Error code: -9
Exceptions caught during RHS evaluation:

***********************************************************************
Exception thrown by Python callback function:
IndexError: index -1 is out of bounds for axis 0 with size 0
***********************************************************************


The right-hand side routine failed at the first call.
Components with largest weighted error estimates:

55: 2.44141e+17
47: 9.76562e+16
54: 6.4e+16
39: 3.90625e+16
46: 3.2e+16
38: 1.6e+16
31: 1.5625e+16
53: 1.13906e+16
30: 8e+15
45: 7.59375e+15
***********************************************************************

[Shirin]

Hello,

I am trying to solve to isentropic compression process in IC engine. I have had different problems so far and I could solve them but now I am getting a new error. Please help me to solve this one. I think the problem is from my data which I attached here.

My code is:

import cantera as ct

import numpy as np

import sys

import xlrd 

  

# Give the location of the file 

loc = ("vol_time_engine.xlsx") 

  

# To open Workbook 

wb = xlrd.open_workbook(loc) 

sheet = wb.sheet_by_index(0)

  

gas = ct.Solution('gri30.xml')

print(gas.species_names)

T_in = 300 #K

P_in = 1e5  #pa

# to set the composition of the gas

fuel_species = 'ch4'

ifuel = gas.species_index(fuel_species)

print(ifuel)

io2 = gas.species_index('o2')

print(io2 )

in2 = gas.species_index('n2')

print(in2)

phi=1

gas.TPX = T_in, P_in, {'CH4':1, 'O2':2/phi, 'N2':2*3.76/phi}

gas()

print("GasTPX: ", gas.TPX)

reac = ct.IdealGasReactor(gas)

env = ct.Reservoir(ct.Solution('air.xml'))

Exp_time = []

Exp_vol = []

keywords={}

for i in range(0,sheet.nrows,1):

       Exp_vol.append(sheet.cell_value(i, 0))

       Exp_time.append(sheet.cell_value(i, 1))

       keywords = {'Time': Exp_time, 'Vol': Exp_vol}

   

print("Volume: ", Exp_vol)

print("Time: ", len(Exp_time))

print("keywords: ", keywords)

class VolumeProfile():

       def __init__(self, keywords):

              self.time = np.array(keywords['Time'])              

              self.volume = np.array(keywords['Vol'])/keywords['Vol'][0]              

              self.velocity = np.diff(Exp_vol)/np.diff(Exp_time)              

              self.velocity = np.append(self.velocity, 0)

       def __call__(self, t): #Return the velocity when called during a time step.

              if t < self.time[-1]:

                     prev_time_point = self.time[self.time <= t][-1] #prev_time_point is the previous value in the time array after the current simulation time                     

                     index = np.where(self.time == prev_time_point)[0][0] #index is the index of the time array where prev_time_point occurs

                     return self.velocity[index]

              else:

                     return 0

net=ct.Wall(reac, env, velocity=VolumeProfile(keywords))

netw = ct.ReactorNet([reac])

time = []

temperature = []

pressure = []

volume = []

while netw.time < 0.001:

    time.append(netw.time)

    print(time)

    temperature.append(reac.T)

    print(temperature)

    pressure.append(reac.thermo.P)

    print(pressure)

    volume.append(reac.volume)

    print(volume)

    netw.step()

print (pressure)

    

#%matplotlib notebook

import matplotlib.pyplot as plt

plt.figure()

plt.plot(time, pressure)

plt.ylabel('Pressure [Pa]')

#plt.xlabel('Time [s]')

plt.show()

    

___________________________________________________

error is :

Traceback (most recent call last):

  File "E:\Compression process\isentropic_compression_ic_new_class.py", line 125, in <module>

    netw.step()

  File "interfaces\cython\cantera\reactor.pyx", line 862, in cantera._cantera.ReactorNet.step

cantera._cantera.CanteraError: 

***********************************************************************

CanteraError thrown by CVodesIntegrator::step:

CVodes error encountered. Error code: -3

Exceptions caught during RHS evaluation:

***********************************************************************

CanteraError thrown by checkFinite:

ydot contains non-finite elements:

ydot[2] = -nan

ydot[3] = -nan

ydot[4] = -nan

ydot[5] = -nan

ydot[6] = -nan

ydot[7] = -nan

ydot[8] = -nan

ydot[9] = -nan

ydot[10] = -nan

ydot[11] = -nan

ydot[12] = -nan

ydot[13] = -nan

ydot[14] = -inf

ydot[15] = -nan

ydot[16] = -nan

ydot[17] = -nan

ydot[18] = -nan

ydot[19] = -nan

ydot[20] = -nan

ydot[21] = -nan

ydot[22] = -nan

ydot[23] = -inf

ydot[24] = -nan

ydot[25] = -nan

ydot[26] = -nan

ydot[27] = -nan

ydot[28] = -nan

ydot[29] = inf

ydot[30] = -nan

ydot[31] = -nan

ydot[32] = -nan

ydot[33] = -nan

ydot[34] = -inf

ydot[35] = -inf

ydot[37] = inf

ydot[38] = -nan

ydot[39] = -nan

ydot[40] = -nan

ydot[41] = -nan

ydot[42] = -nan

ydot[43] = -nan

ydot[44] = inf

ydot[45] = -nan

ydot[46] = -nan

ydot[47] = inf

ydot[48] = inf

ydot[49] = -nan

ydot[50] = -nan

ydot[52] = -nan

ydot[53] = -inf

ydot[54] = -nan

***********************************************************************

***********************************************************************

CanteraError thrown by checkFinite:

ydot contains non-finite elements:

ydot[2] = -nan

ydot[3] = -nan

ydot[4] = -nan

ydot[5] = -nan

ydot[6] = -nan

ydot[7] = -nan

ydot[8] = -nan

ydot[9] = -nan

ydot[10] = -nan

ydot[11] = -nan

ydot[12] = -nan

ydot[13] = -nan

ydot[14] = -nan

ydot[15] = -nan

ydot[16] = -nan

ydot[17] = -nan

ydot[18] = -nan

ydot[19] = inf

ydot[20] = -nan

ydot[21] = -nan

ydot[22] = -nan

ydot[23] = inf

ydot[24] = -inf

ydot[25] = -nan

ydot[26] = -nan

ydot[27] = -nan

ydot[28] = -nan

ydot[29] = inf

ydot[30] = -nan

ydot[31] = inf

ydot[32] = inf

ydot[33] = -nan

ydot[34] = -inf

ydot[35] = -inf

ydot[37] = inf

ydot[38] = -nan

ydot[39] = -nan

ydot[40] = -nan

ydot[41] = -nan

ydot[42] = -nan

ydot[43] = -nan

ydot[44] = inf

ydot[45] = -nan

ydot[46] = -nan

ydot[47] = inf

ydot[48] = inf

ydot[49] = -nan

ydot[50] = -nan

ydot[52] = inf

ydot[53] = -inf

ydot[54] = -inf

***********************************************************************

***********************************************************************

CanteraError thrown by checkFinite:

ydot contains non-finite elements:

ydot[2] = -nan

ydot[3] = -nan

ydot[4] = -nan

ydot[5] = -nan

ydot[6] = -nan

ydot[7] = -nan

ydot[8] = -nan

ydot[9] = -nan

ydot[10] = -inf

ydot[11] = -nan

ydot[12] = -nan

ydot[13] = -nan

ydot[14] = -nan

ydot[15] = -nan

ydot[16] = -nan

ydot[17] = -nan

ydot[18] = -nan

ydot[19] = -nan

ydot[20] = -nan

ydot[21] = -nan

ydot[22] = -nan

ydot[24] = -nan

ydot[25] = -nan

ydot[26] = inf

ydot[27] = inf

ydot[28] = -nan

ydot[29] = -inf

ydot[30] = -nan

ydot[31] = -nan

ydot[32] = -nan

ydot[33] = -nan

ydot[34] = -inf

ydot[35] = -inf

ydot[37] = -nan

ydot[38] = -nan

ydot[39] = -nan

ydot[40] = -nan

ydot[41] = -nan

ydot[42] = -nan

ydot[43] = -nan

ydot[44] = inf

ydot[45] = -nan

ydot[46] = -nan

ydot[47] = inf

ydot[48] = inf

ydot[49] = -nan

ydot[50] = -nan

ydot[52] = -nan

ydot[53] = -inf

ydot[54] = -nan

***********************************************************************

At t = 1.21e-07 and h = 1.01169e-38, the error test failed repeatedly or with |h| = hmin.

Components with largest weighted error estimates:

15: 128.324

52: -54.9843

4: -7.3937

28: -0.104056

2: -0.0750112

5: -0.019342

19: 0.000384387

20: 0.000240976

12: 7.90863e-05

8: 6.67031e-06

***********************************************************************

>>> 

>>> 

[Bryan W. Weber]

Hi Shirin,

Use the volume_history.xlsx file you attached to your post, I could not reproduce the error that you get, and the script completed successfully. Can you please post which version of Cantera you're using? It also seems the script you're running is not the same one in your post, since the one in your post has less than 100 lines but the error is on line 125. Please post (attach) the exact script and input file you're using.

Best,

Bryan

[Shirin]

Hi Bryan,

Thank you for your reply. When I use the below code with the attached data I get the below error. However, when I change the wall position and put the reac first and then env in line :net=ct.Wall(reac, env, velocity=VolumeProfile(keywords)) the code runs but the results are not correct.

I am using Cantera 2.4.

In addition, when I use different mechanism I get different error. I dont want to confuse you so I go with this error and after solving this I will post the next error. 

Can I ask how to see the line numbers in the code? I am using Python 3.7.5 Shell to run the script.

Thanks,

Shirin

My code:

import cantera as ct

import numpy as np

import sys

import xlrd 

  

# Give the location of the file 

loc = ("vol_time_engine_.xlsx") 

  

# To open Workbook 

wb = xlrd.open_workbook(loc) 

sheet = wb.sheet_by_index(0)

  

#gas = ct.Solution('mech_DME_2000.xml')

gas = ct.Solution('gri30.xml')

print(gas.species_names)

T_in = 300 #K

P_in = 1e5  #pa

# to set the composition of the gas

fuel_species = 'ch4'

ifuel = gas.species_index(fuel_species)

print(ifuel)

io2 = gas.species_index('o2')

print(io2 )

in2 = gas.species_index('n2')

print(in2)

#air_N2_O2_molar_ratio = 3.76

nu=2

phi=1

#X[ifuel] =1

#X[io2] =nu/phi

#X[in2] = X[io2]*air_N2_O2_molar_ratio

#gas.X={'CH4':1, 'O2':2/phi, 'N2':2*3.76/phi}

#gas.TPX = T_in, P_in, {'ch4':1, 'o2':2/phi, 'n2':2*3.76/phi}

gas.TPX = T_in, P_in, {'CH4':0.5, 'O2':nu/phi, 'N2':nu*3.76/phi}

#gas.TPX = T_in, P_in, {'ch4':0.5, 'o2':nu/phi, 'n2':nu*3.76/phi}

gas()

print("GasTPX: ", gas.TPX)

print("gas.n_species: ",gas.n_species)

reac = ct.IdealGasReactor(gas)

env = ct.Reservoir(ct.Solution('air.xml'))

Exp_time = []

Exp_vol = []

keywords={}

for i in range(0,sheet.nrows,1):

       Exp_vol.append(sheet.cell_value(i, 0))

       Exp_time.append(sheet.cell_value(i, 1))

       keywords = {'Time': Exp_time, 'Vol': Exp_vol}

   

print("Volume: ", Exp_vol)

print("Time: ", len(Exp_time))

print("keywords: ", keywords)

class VolumeProfile():

       def __init__(self, keywords):

              self.time = np.array(keywords['Time'])

              #print("self.time :",self.time)

              self.volume = np.array(keywords['Vol'])/keywords['Vol'][0]

              #print("self.volume: ", self.volume)

              self.velocity = np.diff(Exp_vol)/np.diff(Exp_time)

              print("self.velocity: ",self.velocity )

              self.velocity = np.append(self.velocity, 0)

       def __call__(self, t): #Return the velocity when called during a time step.

              if t < self.time[-1]:

                     prev_time_point = self.time[self.time <= t][-1] #prev_time_point is the previous value in the time array after the current simulation time                     

                     index = np.where(self.time == prev_time_point)[0][0] #index is the index of the time array where prev_time_point occurs

                     return self.velocity[index]

              else:

                     return 0

net=ct.Wall(env, reac, velocity=VolumeProfile(keywords))

netw = ct.ReactorNet([reac])

time = []

temperature = []

pressure = []

volume = []

while netw.time < 0.08:

      time.append(netw.time)

      temperature.append(reac.T)    

      pressure.append(reac.thermo.P)    

      volume.append(reac.volume)    

      netw.step()

print (time)

print (pressure)

    

#%matplotlib notebook

import matplotlib.pyplot as plt

plt.figure()

plt.plot(time, pressure)

plt.ylabel('Pressure [Pa]')

plt.xlabel('Time [s]')

plt.figure()

plt.plot(Exp_time, Exp_vol, label='Experimental Volume', color='k', linestyle=':')

plt.plot(time, volume, label='Simulated Volume',linestyle='-', markersize=15)

plt.legend()

plt.ylabel('Volume [cm^3]')

plt.xlabel('Time [s]')

plt.figure()

plt.plot(time, temperature)

plt.ylabel('temperature [K]')

plt.show()

    

______________________________________________

Error:

Traceback (most recent call last):

  File "E:\Compression process\isentropic_compression_ic_new_class.py", line 123, in <module>

    netw.step()

  File "interfaces\cython\cantera\reactor.pyx", line 862, in cantera._cantera.ReactorNet.step

cantera._cantera.CanteraError: 

***********************************************************************

CanteraError thrown by CVodesIntegrator::step:

CVodes error encountered. Error code: -3

At t = 0.000661096 and h = 2.60813e-34, the error test failed repeatedly or with |h| = hmin.

Components with largest weighted error estimates:

52: -96.0805

20: 82.0395

28: 68.9422

5: -0.0106475

2: 0.00604313

12: 0.00125758

4: -0.000951175

8: 1.99692e-06

17: 9.51322e-07

22: 3.53662e-08

***********************************************************************

[Shirin]

Hello,

I am trying to solve  isentropic compression process in IC engine. When I use the below code with the attached data I get the below error. However, when I change the wall position and put the reac first and then env in line :net=ct.Wall(reac, env, velocity=VolumeProfile(keywords)) the code runs but the results are not correct.

I am using Cantera 2.4.

I think this code is applicable only to RCM , any thought? 

I would appreciate your help,

[Bryan W. Weber]

Hi Shirin,

I think the change you need to make to your code is to limit the step size of the integrator. I think what's happening is that it is trying to step too far in time and then it runs into errors that are too big and it can't recover. Please see the code below which works for me with Cantera 2.4.0 and Python 3.7.0:

import cantera as ct
import numpy as np
import sys
import xlrd
 
# Give the location of the file
loc = ("vol_time_engine_.xlsx")
 
# To open Workbook
wb = xlrd.open_workbook(loc)
sheet = wb.sheet_by_index(0)

#gas = ct.Solution('mech_DME_2000.xml')
gas = ct.Solution('gri30.xml')

# print(gas.species_names)

T_in = 300 #K
P_in = 1e5  #pa
# to set the composition of the gas
fuel_species = 'ch4'
ifuel = gas.species_index(fuel_species)

# print(ifuel)

io2 = gas.species_index('o2')

#print(io2 )

in2 = gas.species_index('n2')

#print(in2)

#air_N2_O2_molar_ratio = 3.76
nu=2
phi=1
#X[ifuel] =1
#X[io2] =nu/phi
#X[in2] = X[io2]*air_N2_O2_molar_ratio
#gas.X={'CH4':1, 'O2':2/phi, 'N2':2*3.76/phi}


#gas.TPX = T_in, P_in, {'ch4':1, 'o2':2/phi, 'n2':2*3.76/phi}
gas.TPX = T_in, P_in, {'CH4':0.5, 'O2':nu/phi, 'N2':nu*3.76/phi}
#gas.TPX = T_in, P_in, {'ch4':0.5, 'o2':nu/phi, 'n2':nu*3.76/phi}

#gas()
print("GasTPX: ", gas.TPX)
#print("gas.n_species: ",gas.n_species)

reac = ct.IdealGasReactor(gas)
env = ct.Reservoir(ct.Solution('air.xml'))

Exp_time = []
Exp_vol = []

# keywords={}

for i in range(0,sheet.nrows,1):
       Exp_vol.append(sheet.cell_value(i, 0))
       Exp_time.append(sheet.cell_value(i, 1))
       keywords = {'Time': Exp_time, 'Vol': Exp_vol}
   

#print("Volume: ", Exp_vol)
#print("Time: ", len(Exp_time))
#print("keywords: ", keywords)

class VolumeProfile():
       def __init__(self, keywords):
              self.time = np.array(keywords['Time'])
              #print("self.time :",self.time)
              self.volume = np.array(keywords['Vol'])/keywords['Vol'][0]
              #print("self.volume: ", self.volume)

              self.velocity = np.diff(self.volume)/np.diff(self.time)
              # print("self.velocity: ",self.velocity )

              self.velocity = np.append(self.velocity, 0)
       def __call__(self, t): #Return the velocity when called during a time step.
              if t < self.time[-1]:
                     prev_time_point = self.time[self.time <= t][-1] #prev_time_point is the previous value in the time array after the current simulation time                    
                     index = np.where(self.time == prev_time_point)[0][0] #index is the index of the time array where prev_time_point occurs
                     return self.velocity[index]
              else:
                     return 0

vp = VolumeProfile(keywords)
net=ct.Wall(reac, env, velocity=vp)
print(np.min(np.diff(vp.time)))

netw = ct.ReactorNet([reac])

netw.set_max_time_step(np.min(np.diff(vp.time)))

time = []
temperature = []
pressure = []
volume = []

while netw.time < 0.08:
      time.append(netw.time)
      temperature.append(reac.T)    
      pressure.append(reac.thermo.P)    
      volume.append(reac.volume)    
      netw.step()

#print (time)
#print (pressure)

import matplotlib.pyplot as plt

plt.figure()
plt.plot(time, pressure)
plt.ylabel('Pressure [Pa]')
plt.xlabel('Time [s]')
plt.figure()
plt.plot(Exp_time, Exp_vol, label='Experimental Volume', color='k', linestyle=':')
plt.plot(time, volume, label='Simulated Volume',linestyle='-', markersize=15)
plt.legend()
plt.ylabel('Volume [cm^3]')
plt.xlabel('Time [s]')
plt.figure()
plt.plot(time, temperature)
plt.ylabel('temperature [K]')

plt.savefig('temperature.png')

The main section that changed was this one:

vp = VolumeProfile(keywords)
net=ct.Wall(reac, env, velocity=vp)
print(np.min(np.diff(vp.time)))

netw = ct.ReactorNet([reac])

netw.set_max_time_step(np.min(np.diff(vp.time)))

With this change, I get the attached for the temperature.

The reason everything was backwards before was because the order of the arguments to ct.Wall affects which direction the wall is moving.

Hope it helps,

Bryan

[Shirin]

Hi Bryan,

Thank you so much for revising the code. Actually it is running now but the simulation of volume is not. I'v attached the figure for volume.

Best Regards,

Shirin

On Friday, February 7, 2020 at 11:33:21 AM UTC-5, Bryan W. Weber wrote:

Hi Shirin,

plt.plot(Exp_time,<span

[Bryan W. Weber]

Shirin,

The volume is normalized to be between 0 and 1. That's why it appears as a flat line. I bet if you multiply the output volume by 600, you will find the same trace.

Best,

Bryan

[Shirin]

Bryan,

That is correct, sorry for that question. I appreciate your help.

Best Regards,

Shirin

[Shirin]

Hello,

I need to see the mole fraction of some species but I don't know how to call it in python. I am simulating the compression process of a mixture of two fuels and compare the result with what I get in GT_Power. But my temperature and pressure profiles dont show their maximum and I dont know why. I have attached the temp profile. 

Thanks,

Shirin

[Anant Girdhar]

Hi Shirin, 

You can get a list of the mole fractions using: 

    gas.X

The order of the species in this list is given by: 

    gas.species_names

So, one way to get the mole fractions for specific species (say, H2), would be:

    gas.X[gas.species_names.index("H2")]

I hope this helps. 

anant

--

You received this message because you are subscribed to the Google Groups "Cantera Users' Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email to cantera-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cantera-users/014280c1-bd2e-4ba1-a34b-b63533649701%40googlegroups.com.

[Lauge Thorsen]

Hello Cantera people

Is this requirement of Cerberus version 1.1 still a requirement? 

I want to run the rcm-example (https://github.com/pr-omethe-us/PyKED/blob/master/docs/rcm-example.ipynb) in my jupyter notebook setup. I follow the example and get a similar error to what Shirin is showing (see below).

However, when I try to reinstall Cerberus with a 1.1 version it fails (pip install --force-reinstall -v "cerberus==1.1") - see attached picture.Is there some problems with my installation or how do I move forward with this issue? Any help would be appreciated!

Best

Lauge

Input:

import cantera as ct
import numpy as np

from pyked import ChemKED

from urllib.request import urlopen
import yaml
rcm_link = 'https://raw.githubusercontent.com/pr-omethe-us/PyKED/master/pyked/tests/testfile_rcm.yaml'
with urlopen(rcm_link) as response:
    testfile_rcm = yaml.safe_load(response.read())
ck = ChemKED(dict_input=testfile_rcm)
dp = ck.datapoints[0]

Output:

--------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) Cell In[19], line 6 4 with urlopen(rcm_link) as response: 5 testfile_rcm = yaml.safe_load(response.read()) ----> 6 ck = ChemKED(dict_input=testfile_rcm) 7 dp = ck.datapoints[0] File ~\ct-env\lib\site-packages\pyked\chemked.py:120, in ChemKED.__init__(self, yaml_file, dict_input, skip_validation) 117 raise NameError("ChemKED needs either a YAML filename or dictionary as input.") 119 if not skip_validation: --> 120 self.validate_yaml(self._properties) 122 self.datapoints = [] 123 for point in self._properties['datapoints']: File ~\ct-env\lib\site-packages\pyked\chemked.py:179, in ChemKED.validate_yaml(self, properties) 169 def validate_yaml(self, properties): 170 """Validate the parsed YAML file for adherance to the ChemKED format. 171 172 Arguments: (...) 177 string contains the errors that are present. 178 """ --> 179 validator = OurValidator(schema) 180 if not validator.validate(properties): 181 for key, value in validator.errors.items(): File ~\ct-env\lib\site-packages\cerberus\validator.py:183, in BareValidator.__init__(self, *args, **kwargs) 178 """ The error handler used to format :attr:`~cerberus.Validator.errors` 179 and process submitted errors with 180 :meth:`~cerberus.Validator._error`. 181 Type: :class:`~cerberus.errors.BaseErrorHandler` """ 182 self.__store_config(args, kwargs) --> 183 self.schema = kwargs.get('schema', None) 184 self.allow_unknown = kwargs.get('allow_unknown', False) 185 self.require_all = kwargs.get('require_all', False) File ~\ct-env\lib\site-packages\cerberus\validator.py:594, in BareValidator.schema(self, schema) 592 self._schema = schema 593 else: --> 594 self._schema = DefinitionSchema(self, schema) File ~\ct-env\lib\site-packages\cerberus\schema.py:82, in DefinitionSchema.__init__(self, validator, schema) 73 self.schema_validator = SchemaValidator( 74 None, 75 allow_unknown=self.validation_schema, (...) 78 target_validator=validator, 79 ) 81 schema = self.expand(schema) ---> 82 self.validate(schema) 83 self.schema = schema File ~\ct-env\lib\site-packages\cerberus\schema.py:265, in DefinitionSchema.validate(self, schema) 263 _hash = (mapping_hash(schema), mapping_hash(self.validator.types_mapping)) 264 if _hash not in self.validator._valid_schemas: --> 265 self._validate(schema) 266 self.validator._valid_schemas.add(_hash) File ~\ct-env\lib\site-packages\cerberus\schema.py:282, in DefinitionSchema._validate(self, schema) 279 test_rules[rule.replace(" ", "_")] = constraint 280 test_schema[field] = test_rules --> 282 if not self.schema_validator(test_schema, normalize=False): 283 raise SchemaError(self.schema_validator.errors) File ~\ct-env\lib\site-packages\cerberus\validator.py:1042, in BareValidator.validate(self, document, schema, update, normalize) 1040 self.__validate_definitions(definitions, field) 1041 else: -> 1042 self.__validate_unknown_fields(field) 1044 if not self.update: 1045 self.__validate_required_fields(self.document) File ~\ct-env\lib\site-packages\cerberus\validator.py:1075, in BareValidator.__validate_unknown_fields(self, field) 1071 schema_crumb = 'allow_unknown' if self.is_child else '__allow_unknown__' 1072 validator = self._get_child_validator( 1073 schema_crumb=schema_crumb, schema={field: self.allow_unknown} 1074 ) -> 1075 if not validator({field: value}, normalize=False): 1076 self._error(validator._errors) 1077 else: File ~\ct-env\lib\site-packages\cerberus\validator.py:1040, in BareValidator.validate(self, document, schema, update, normalize) 1038 definitions = self.schema.get(field) 1039 if definitions is not None: -> 1040 self.__validate_definitions(definitions, field) 1041 else: 1042 self.__validate_unknown_fields(field) File ~\ct-env\lib\site-packages\cerberus\validator.py:1110, in BareValidator.__validate_definitions(self, definitions, field) 1108 rule = self._remaining_rules.pop(0) 1109 try: -> 1110 result = validate_rule(rule) 1111 # TODO remove on next breaking release 1112 if result: File ~\ct-env\lib\site-packages\cerberus\validator.py:1085, in BareValidator.__validate_definitions.<locals>.validate_rule(rule) 1083 def validate_rule(rule): 1084 validator = self.__get_rule_handler('validate', rule) -> 1085 return validator(definitions.get(rule, None), field, value) File ~\ct-env\lib\site-packages\cerberus\validator.py:1491, in BareValidator._validate_schema(self, schema, field, value) 1489 self.__validate_schema_sequence(field, schema, value) 1490 elif isinstance(value, Mapping): -> 1491 self.__validate_schema_mapping(field, schema, value) File ~\ct-env\lib\site-packages\cerberus\validator.py:1504, in BareValidator.__validate_schema_mapping(self, field, schema, value) 1496 validator = self._get_child_validator( 1497 document_crumb=field, 1498 schema_crumb=(field, 'schema'), (...) 1501 require_all=field_rules.get('require_all', self.require_all), 1502 ) 1503 try: -> 1504 if not validator(value, update=self.update, normalize=False): 1505 self._error(field, errors.MAPPING_SCHEMA, validator._errors) 1506 except _SchemaRuleTypeError: File ~\ct-env\lib\site-packages\cerberus\validator.py:1040, in BareValidator.validate(self, document, schema, update, normalize) 1038 definitions = self.schema.get(field) 1039 if definitions is not None: -> 1040 self.__validate_definitions(definitions, field) 1041 else: 1042 self.__validate_unknown_fields(field) File ~\ct-env\lib\site-packages\cerberus\validator.py:1110, in BareValidator.__validate_definitions(self, definitions, field) 1108 rule = self._remaining_rules.pop(0) 1109 try: -> 1110 result = validate_rule(rule) 1111 # TODO remove on next breaking release 1112 if result: File ~\ct-env\lib\site-packages\cerberus\validator.py:1084, in BareValidator.__validate_definitions.<locals>.validate_rule(rule) 1083 def validate_rule(rule): -> 1084 validator = self.__get_rule_handler('validate', rule) 1085 return validator(definitions.get(rule, None), field, value) File ~\ct-env\lib\site-packages\cerberus\validator.py:366, in BareValidator.__get_rule_handler(self, domain, rule) 364 result = getattr(self, methodname, None) 365 if result is None: --> 366 raise RuntimeError( 367 "There's no handler for '{}' in the '{}' " 368 "domain.".format(rule, domain) 369 ) 370 return result RuntimeError: There's no handler for 'isvalid_history' in the 'validate' domain.

On Friday 15 November 2019 at 17:27:55 UTC+1 bryan....@gmail.com wrote:

Hi Shirin,

There is currently a problem with PyKED and one of the libraries it uses called Cerberus. You need to install version 1.1 or lower of Cerberus. This should be automatically handled by the installer, but it doesn't work sometimes for some reason. Since this isn't really about Cantera, I'll reply to you off the list to do some more troubleshooting.

Best,

Bryan

On Friday, November 15, 2019 at 2:53:27 AM UTC-5, Shirin wrote:

Hello,

I am trying to run your code in rcm-example.ipynb to see if I can get the same plot as you did. But I am getting this below error. I have installed Pyked and pyyaml.  

Traceback (most recent call last):

[Bryan Weber]

Hi Lauge,

Sorry for the delay! Unfortunately, PyKED is essentially unmaintained as far as I know. The reason for the `pip` install error is because of some underlying changes in Python itself, so I suppose Cerberus 1.1 is unsupported on recent Python as well. This was a change made in Python 3.10, so installing Python 3.9 may work. If you encounter further problems, could you post them on the PyKED GitHub page: https://github.com/pr-omethe-us/PyKED/issues

Thank you!

Bryan