GreenPolyProducts/cantera
5
0
Fork 0
mirror of https://github.com/Cantera/cantera.git synced 2025-02-25 18:55:29 -06:00
Code Issues Packages Projects Releases Wiki Activity

Experimental Matlab toolbox now working

Browse Source
This commit is contained in:
Su Sun 2022-01-11 14:36:53 -05:00 committed by Ray Speth
parent 591f3f3539
commit 0189abb705
49 changed files with 2770 additions and 973 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
interfaces/matlab_experimental/1D/AxiStagnFlow.m Normal file
View File

@ -0,0 +1,4 @@
function m = AxiStagnFLow(gas)
% Get an axisymmetric stagnation flow domain.
m = Domain1D(1, gas);
end

11
interfaces/matlab_experimental/1D/AxisymmetricFlow.m Normal file
View File

@ -0,0 +1,11 @@
function m = AxisymmetricFlow(gas, id)
% Create an axisymmetric flow domain.
% :param id:
% String ID of the flow.
m = Domain1D(1, gas);
if nargin == 1
m.setID('flow');
else
m.setID(id);
end
end

531
interfaces/matlab_experimental/1D/Domain1D.m Normal file
View File

@ -0,0 +1,531 @@
classdef Domain1D < handle
properties
dom_id
domain_type
T
end
methods
%% Domain1D class constructor.
function d = Domain1D(a, b, c)
% parameter a:
% Integer type of domain. Possible values are:
% * 1 - Stagnation Flow
% * 2 - Inlet1D
% * 3 - Surf1D
% * 4 - Symm1D
% * 5 - Outlet1D
% * 6 - Reacting Surface
% * 8 - Sim1D
% * -2 - OutletRes
%
% parameter b:
% Instance of class 'Solution' (for a == 1) or 'Interface'
% (for a == 6). Not used for all other valid values of 'a'.
%
% parameter c:
% Integer, either 1 or 2, indicating whether an axisymmetric
% stagnation flow or a free flame should be created. If not
% specified, defaults to 1. Ignored if a!= 1.
checklib;
d.dom_id = 1;
if nargin == 1
if a == 2
calllib(ct, 'inlet_new');
elseif a == 3
calllib(ct, 'surf_new');
elseif a == 4
calllib(ct, 'symm_new');
elseif a == 5
calllib(ct, 'outlet_new');
elseif a == -2
calllib(ct, 'outletres_new');
else
error('Not enough arguments for that job number');
end
elseif nargin == 2
% a stagnation flow
if a == 1
if isa(b, 'Solution')
d.dom_id = calllib(ct, 'stflow_new', ...
b.tp_id, b.kin_id, b.tr_id, 1);
else
error('Wrong argument type. Expecting instance of class Solution.');
end
elseif a == 6
if isa(b, 'Interface')
d.dom_id = calllib(ct, 'reactingsurf_new');
calllib(ct, 'reactingsurf_setkineticsmgr', ...
d.dom_id, b.kin_id);
else
error('Wrong argument type. Expecting instance of class Interface.');
end
else
error('Wrong object type.');
end
elseif nargin == 3
if a == 1
if isa(b, 'Solution')
d.dom_id = calllib(ct, 'stflow_new', ...
b.tp_id, b.kin_id, b.tr_id, c);
else
error('Wrong argument type. Expecting instance of class Solution.');
end
else
error('Unknown domain type.');
end
end
% if d.dom_id < 0
% error(geterr);
% end
d.domain_type = a;
end
%% Utility Methods
function dom_clear(d)
% Delete the Domain1D object
checklib;
calllib(ct, 'domain_del', d.dom_id);
end
%% Domain Methods
function n = componentIndex(d, name)
% Get the index of a component given its name
%
% parameter d:
% Instance of class 'Domain1D'.
% parameter name:
% String name of the component to look up. If a numeric
% value is passed, it will be returned directly.
% return:
% Index of the component.
checklib;
if isa(name, 'double')
n = name;
else
n = calllib(ct, 'domain_componentIndex', ...
d.dom_id, name);
if n >= 0
n = n+1;
end
end
if n <= 0
error('Component not found');
end
end
function n = componentName(d, index)
% Get the name of a component given its index.
%
% parameter d:
% Instance of class 'Domain1D'.
% parameter n:
% Integer or vector of integers of components' names.
% return:
% Cell array of component names.
checklib;
n = length(index);
s = cell(m);
for i = 1:n
id = index(i)-1;
buflen = calllib(ct, 'domain_componentName', ...
d.dom_id, id, 0, 0);
if buflen > 0
aa = char(zeros(1, buflen));
[out_buf, aa] = calllib(ct, 'domain_componentName', ...
d.dom_id, id, buflen, aa);
s{i} = aa;
end
end
end
function d = disableEnergy(d)
% Disable the energy equation.
%
% parameter d:
% Instance of class 'Domain1D'.
checklib;
disp(' ');
disp('Disabling the energy equation...');
calllib(ct, 'stflow_solveEnergyEqn', d.dom_id, 0);
end
function i = domainIndex(d)
% Get the domain index.
% return:
% This function returns an integer flag denoting the
% location of the domain, beginning with 1 at the left.
checklib;
i = calllib(ct, 'domain_index', d.dom_id);
if i >= 0
i = i + 1;
end
if i <= 0
error('Domain not found');
end
end
function i = domainType(d)
% Get the type of domain.
% return:
% This function returns an integer flag denoting the domain
% type.
checklib;
i = calllib(ct, 'domain_type', d.dom_id);
end
function d = enableEnergy(d)
% Enable the energy equation.
% parameter d:
% Instance of class 'Domain1D'.
checklib;
disp(' ');
disp('Enabling the energy equation...');
calllib(ctm 'stflow_solveEnergyEqn', d.dom_id, 1);
end
function zz = gridPoints(d, n)
% Get grid points from a domain.
% parameter d:
% Instance of class 'Domain1D'.
% parameter n:
% Optional, vector of grid points to be retrieved.
% return:
% Vector of grid points.
checklib;
if nargin == 1
np = d.nPoints;
zz = zeros(1, d.np);
for i = 1:np
zz(i) = calllib(ct, 'domain_grid', dom_id, i-1);
end
else
m = length(n);
zz = zeros(1, m);
for i = 1:m
zz(i) = calllib(ct, 'domain_grid', dom_id, n(i)-1);
end
end
end
function a = isFlow(d)
% Determine whether a domain is a flow.
% parameter d:
% Instance of class 'Domain1D'.
% return:
% 1 if the domain is a flow domain, and 0 otherwise.
checklib;
t = d.domainType;
% See Domain1D.h for definitions of constants.
if t < 100
a = 1;
else a = 0;
end
end
function a = isInlet(d)
% Determine whether a domain is an inlet.
% parameter d:
% Instance of class 'Domain1D'.
% return:
% 1 if the domain is an inlet, and 0 otherwise.
checklib;
t = d.domainType;
% See Domain1D.h for definitions of constants.
if t == 104
a = 1;
else a = 0;
end
end
function a = isSurface(d)
% Determine whether a domain is a surface.
% parameter d:
% Instance of class 'Domain1D'.
% return:
% 1 if the domain is a surface, and 0 otherwise.
checklib;
t = d.domainType;
% See Domain1D.h for definitions of constants.
if t == 102
a = 1;
else a = 0;
end
end
function mdot = massFlux(d)
% Get the mass flux.
% return:
% The mass flux in the domain.
checklib;
mdot = calllib(ct, 'bdry_mdot', d.dom_id);
end
function y = massFraction(d, k)
% Get the mass fraction of a species given its integer index.
% This method returns the mass fraction of species 'k', where k
% is the integer index of the species in the flow domain to
% which the boundary domain is attached.
%
% parameter d:
% Instance of class 'Domain1D'
% parameter k:
% Integer species index
% return:
% Mass fraction of species
checklib;
if d.domainIndex == 0
error('No flow domain attached!')
end
if d.isInlet
y = calllib(ct, 'bdry_massFraction', d.dom_id, k-1);
else error('Input domain must be an inlet');
end
end
function n = nComponents(d)
% Get the number of components.
% return:
% Number of variables at each grid point.
checklib;
n = calllib(ct, 'domain_nComponents', d.dom_id);
end
function n = nPoints(d)
% Get the number of grid points.
% return:
% Integer number of grid points.
checklib;
n = calllib(ct, 'domain_nPoints', d.dom_id);
end
function setBounds(d, component, lower, upper)
% Set bounds on the solution components.
% parameter d:
% Instance of class 'Domain1D'.
% parameter component:
% String, component to set the bounds on.
% parameter lower:
% Lower bound.
% parameter upper:
% Upper bound.
checklib;
n = d.componentIndex(component);
calllib(ct, 'domain_setBounds', d.dom_id, n-1, lower, upper);
end
function setCoverageEqs(d, onoff)
% Enable or disable solving the coverage equations.
% parameter d:
% Instance of class 'Domain1D'
% parameter onoff:
% string, one of 'on' or 'yes' to turn solving the coverage
% equations on. One of 'off' or 'no' to turn off the
% coverage equation.
checklib;
if d.domain_type ~= 6
error('Wrong domain type. Expected a reacting surface domain.');
end
ion = -1;
if isa(onoff,'char')
if strcmp(onoff, 'on') || strcmp(onoff, 'yes')
ion = 1;
elseif strcmp(onoff, 'off') || strcmp(onoff, 'no')
ion = 0;
else
error(strcat('unknown option: ', onoff))
end
elseif isa(onoff, 'numeric')
ion = onoff;
end
calllib(ct, 'reactingsurf_enableCoverageEqs', d.dom_id, ion);
end
function setFixedTempProfile(d, profile)
% Set a fixed temperature profile to use when the energy
% equation is not being solved. The profile must be entered as
% an array of positions/temperatures, which may be in rows or
% columns.
%
% parameter d:
% Instance of class 'Domain1D'.
% parameter profile:
% n x 2 or 2 x n array of 'n' points at which the
% temperature is specified.
checklib;
sz = size(profile);
if sz(1) == 2
l = length(1, :);
calllib(ct, '', d.dom_id, l, profile(1, :), l, profile(2, :));
elseif sz(2) == 2
l = length(:, 1);
calllib(ct, '', d.dom_id, l, profile(:, 1); l, profile(:, 2));
else error('Wrong temperature profile array shape.');
end
end
function setID(d, id)
% Set the ID tag for a domain.
% parameter id:
% String ID to assign.
checklib;
calllib(ct, 'domain_setID', d.dom_id, id);
end
function setMdot(d, mdot)
% Set the mass flow rate.
% parameter mdot:
% Mass flow rate.
checklib;
calllib(ct, 'bdry_setTemperature', d.dom_id, mdot);
end
function setMoleFractions(d, x)
% Set the mole fractions.
% parameter x:
% String specifying the species and mole fractions in the
% format "'Spec:X,Spec2:X2'".
checklib;
calllib(ct, 'bdry_setMoleFractions', d.dom_id, x);
end
function setPressure(d, p)
% Set the pressure.
% parameter p:
% Pressure to be set. Unit: Pa.
checklib;
calllib(ct, 'bdry_setPressure', d.dom_id, p);
end
function setProfileD(d, n, p)
% Set the profile of a component.
% Convenience function to allow an instance of 'Domain1D' to
% have a profile of its components set when it is part of a
% 'Stack'.
%
% parameter d:
% Instance of class 'Domain1d'.
% parameter n:
% Integer index of component, vector of component indices,
% string of component name, or cell array of strings of
% component names.
% parameter p:
% n x 2 array, whose columns are the relative (normalized)
% positions and the component values at those points. The
% number of positions 'n' is arbitrary.
checklib;
if d.stack == 0
error('Install domain in stack before calling setProfile.');
end
d.stack.setProfile(d.domainIndex, n, p);
end
function setSteadyTolerances(d, component, rtol, atol)
% Set the steady-state tolerances.
% parameter d:
% Instance of class 'Domain1D'.
% parameter component:
% String or cell array of strings of component values whose
% tolerances should be set. If 'default' is specified, the
% tolerance of all components will be set.
% parameter rtol:
% Relative tolerance.
% parameter atol:
% Absolute tolerance.
checklib;
if strcmp(component, 'default')
nc = d.nComponents;
for ii = 1:nc
calllib(ct, 'domain_setSteadyTolerances', ...
d.dom_id, ii, rtol, atol);
end
elseif iscell(component)
nc = length(component);
for ii = 1:nc
n = d.componentIndex(component{ii});
calllib(ct, 'domain_setSteadyTolerances', ...
d.dom_id, n, rtol, atol);
end
else
n = d.componentIndex(component);
calllib(ct, 'domain_setSteadyTolerances', ...
d.dom_id, ii, rtol, atol);
end
end
function setTransientTolerances(d, component, rtol, atol)
% Set the transient tolerances.
% parameter d:
% Instance of class 'Domain1D'.
% parameter component:
% String or cell array of strings of component values whose
% tolerances should be set. If 'default' is specified, the
% tolerance of all components will be set.
% parameter rtol:
% Relative tolerance.
% parameter atol:
% Absolute tolerance.
checklib;
if strcmp(component, 'default')
nc = d.nComponents;
for ii = 1:nc
calllib(ct, 'domain_setTransientTolerances', ...
d.dom_id, ii, rtol, atol);
end
elseif iscell(component)
nc = length(component);
for ii = 1:nc
n = d.componentIndex(component{ii});
calllib(ct, 'domain_setTransientTolerances', ...
d.dom_id, n, rtol, atol);
end
else
n = d.componentIndex(component);
calllib(ct, 'domain_setTransientTolerances', ...
d.dom_id, ii, rtol, atol);
end
end
function temperature = get.T(d)
% Get the boundary temperature (K).
checklib;
temperature = calllib(ct, 'bdry_temperature', d.dom_id);
end
function set.T(d, t)
% Set the temperature (K).
checklib;
if temperature <= 0
error('The temperature must be positive');
end
calllib(ct, 'bdry_setTemperature', d.dom_id, t);
end
function setupGrid(d, grid)
% Set up the solution grid.
checklib;
calllib(ct, 'domain_setupGrid', d.dom_id, numel(grid), grid);
end
end
end

9
interfaces/matlab_experimental/1D/FreeFlame.m Normal file
View File

@ -0,0 +1,9 @@
function m = FreeFlame(gas, id)
% Create a freely-propagating flat flame.
m = Domain1D(1, gas, 2);
if nargin == 1
m.setID('flame');
else
m.setID(id);
end
end

9
interfaces/matlab_experimental/1D/Inlet.m Normal file
View File

@ -0,0 +1,9 @@
function m = Inlet(id)
% Create an inlet domain.
m = Domain1D(2);
if nargin == 0
m.setID('inlet');
else
m.setID(id);
end
end

9
interfaces/matlab_experimental/1D/Outlet.m Normal file
View File

@ -0,0 +1,9 @@
function m = Outlet(id)
% Create an outlet domain.
m = Domain1D(5);
if nargin == 0
m.setID('outlet');
else
m.setID(id);
end
end

9
interfaces/matlab_experimental/1D/OutletRes.m Normal file
View File

@ -0,0 +1,9 @@
function m = OutletRes(id)
% Create an outlet reservoir domain.
m = Domain1D(-2);
if nargin == 0
m.setID('outletres');
else
m.setID(id);
end
end

458
interfaces/matlab_experimental/1D/Sim1D.m Normal file
View File

@ -0,0 +1,458 @@
classdef Stack < handle
properties
st_id
domains
end
methods
%% Stack class constructor
function s = Stack(domains)
% A stack object is a container for one-dimensional domains,
% which are instances of class Domain1D. The domains are of two
% types - extended domains, and connector domains.
%
% parameter domains:
% Vector of domain instances.
% return:
% Instance of class 'Stack'.
checklib;
s.st_id = 1;
s.domains = domains;
if nargin == 1
nd = length(domains);
ids = zeros(1, nd);
for n=1:nd
ids(n) = domains(n).dom_id;
end
s.st_id = calllib(ct, 'sim1D_new', nd, ids);
else
help(Stack);
error('Wrong number of parameters.');
end
% if s.st_id < 0
% error(geterr);
% end
end
%% Utility Methods
function st_clear(s)
% Delete the Sim1D object
checklib;
calllib(ct, 'sim1D_del', s.st_id);
end
function display(s, fname)
% Show all domains.
%
% parameter s:
% Instance of class 'Stack'.
% parameter fname:
% File to write summary to. If omitted, output is to the
% command window.
checklib;
if nargin == 1
fname = '-';
end
calllib(ct, 'sim1D_showSolution', s.st_id, fname);
end
%% Stack Methods
function n = stackIndex(s, name)
% Get the index of a domain in a stack given its name.
%
% parameter s:
% Instance of class 'Stack'.
% parameter name:
% If double, the value is returned. Otherwise, the name is
% looked up and its index is returned.
% return:
% Index of domain.
checklib;
if isa(name, 'double')
n = name;
else
n = calllib(ct, 'sim1D_domainIndex', s.st_id, name);
if n >= 0
n = n+1;
else
error('Domain not found');
end
end
end
function z = grid(s, name)
% Get the grid in one domain.
%
% parameter s:
% Instance of class 'Stack'.
% parameter name:
% Name of the domain for which the grid should be retrieved.
% return:
% The grid in domain name.
n = s.stackIndex(name);
d = s.domains(n);
z = d.gridPoints;
end
function plotSolution(s, domain, component)
% Plot a specified solution component.
%
% parameter s:
% Instance of class 'Stack'.
% parameter domain:
% Name of domain from which the component should be
% retrieved.
% parameter component:
% Name of the component to be plotted.
n = s.stackIndex(domain);
d = s.domains(n);
z = d.gridPoints;
x = s.solution(domain, component);
plot(z, x);
xlabel('z (m)');
ylabel(component);
end
function r = resid(s, domain, rdt, count)
% Get the residuals.
%
% parameter s:
% Instance of class 'Stack'.
% parameter domain:
% Name of the domain.
% parameter rdt:
% parameter count:
% return:
checklib;
if nargin == 2
rdt = 0.0;
count = 0;
end
idom = s.stackIndex(domain);
d = s.domains(idom);
nc = d.nComponents;
np = d.nPoints;
r = zeros(nc, np);
calllib(ct, 'sim1D_eval', s.st_id, rdt, count);
for m = 1:nc
for n = 1:np
r(m, n) = calllib(ct, 'sim1D_workValue', ...
s.st_id, idom - 1, m - 1, n - 1);
end
end
end
function restore(s, fname, id)
% Restore a previously-saved solution.
% This method can be used ot provide an initial guess for the
% solution.
%
% parameter s:
% Instance of class 'Stack'.
% parameter fname:
% File name of an XML file containing solution info.
% parameter id:
% ID of the element that should be restored.
checklib;
calllib(ct, 'sim1D_restore', s.st_id, fname, id)
end
function saveSoln(s, fname, id, desc)
% Save a solution to a file.
% The output file is in a format that can be used by 'restore'.
%
% parameter s:
% Instance of class 'Stack'.
% parameter fname:
% File name where XML file should be written.
% parameter id:
% ID to be assigned to the XMl element when it is written.
% parameter desc:
% Description to be written to the output file.
checklib;
if nargin == 1
fname = 'soln.xml';
id = 'solution';
desc = '--';
elseif nargin == 2
id = 'solution';
desc = '--';
elseif nargin == 3
desc = '--';
end
calllib(ct, 'sim1D_save', s.st_id, fname, id, desc);
end
function setFlatProfile(s, domain, comp, v)
% Set a component to a value across the entire domain.
%
% parameter s:
% Instance of class 'Stack'.
% parameter domain:
% Integer ID of the domain.
% parameter comp:
% Component to be set.
% parameter v:
% Double value to be set.
checklib;
calllib(ct, 'sim1D_setFlatProfile', s.st_id, ...
domain - 1, comp - 1, v);
end
function setMaxJacAge(s, ss_age, ts_age)
% Set the number of times the Jacobian will be used before it
% is recomputed.
%
% parameter s:
% Instance of class 'Stack'.
% parameter ss_age:
% Maximum age of the Jacobian for steady state analysis.
% parameter ts_age:
% Maximum age of the Jacobian for transient analysis. If not
% specified, deftauls to 'ss_age'.
checklib;
if nargin == 2
ts_age = ss_age;
end
calllib(ct, 'sim1D_setMaxJacAge', s.st_id, ss_age, ts_age);
end
function setProfile(s, name, comp, p)
% Specify a profile for one component,
%
% The solution vector values for this component will be
% linearly interpolated from the discrete function defined by
% p(:, 1) vs p(:, 2).
% Note that "p(1, 1) = 0.0" corresponds to the leftmost grid
% point in the specified domain, and "p(1, n) = 1.0"
% corresponds to the rightmost grid point. This method can be
% called at any time, but is usually used to set the initial
% guess for the solution.
%
% Example (assuming 's' is an instance of class 'Stack'):
% >> zr = [0.0, 0.1, 0.2, 0.4, 0.8, 1.0];
% >> v = [500, 650, 700, 730, 800, 900];
% >> s.setProfile(1, 2, [zr, v]);
%
% parameter s:
% Instance of class 'Stack'.
% parameter name:
% Domain name.
% parameter comp:
% Component number.
% parameter p:
% n x 2 array, whose columns are the relative (normalized)
% positions and the component values at those points. The
% number of positions 'n' is arbitrary.
checklib;
if isa(name, 'double')
n = name;
else
n = s.domainIndex(name);
end
d = s.domains(n);
if isa(comp, 'double') || isa(comp, 'cell')
c = comp;
elseif isa(comp, 'char')
c = {comp};
else
error('Wrong type.');
end
np = length(c);
sz = size(p);
if sz(1) == np + 1;
for j = 1:np
ic = d.componentIndex(c{j});
calllib(ct, 'sim1D_setProfile', s.st_id, ...
n - 1, ic - 1, sz(1), p(1, :), sz(1), p(j+1, :));
end
elseif sz(2) == np + 1;
ic = d.componentIndex(c{j});
calllib(ct, 'sim1D_setProfile', s.st_id, ...
n - 1, ic - 1, sz(2), p(:, 1), sz(2), p(:, j+1));
else
error('Wrong profile shape.');
end
end
function setRefineCriteria(s, n, ratio, slope, curve, prune)
% Set the criteria used to refine the grid.
%
% parameter s:
% Instance of class 'Stack'.
% parameter ratio:
% Maximum size ratio between adjacent cells.
% parameter slope:
% Maximum relative difference in value between adjacent
% points.
% parameter curve:
% Maximum relative difference in slope between adjacent
% cells.
% parameter prune:
% Minimum value for slope or curve for which points will be
% retained or curve value is below prune for all components,
% it will be deleted, unless either neighboring point is
% already marked for deletion.
checklib;
if nargin < 3
ratio = 10.0;
end
if nargin < 4
slope = 0.8;
end
if nargin < 5
curve = 0.8;
end
if nargin < 6
prune = -0.1;
end
calllib(ct, 'sim1D_setRefineCriteria', s.st_id, ...
n - 1, ratio, slope, curve, prune);
end
function setTimeStep(s, stepsize, steps)
% Specify a sequence of time steps.
%
% parameter stepsize:
% Initial step size.
% parameter steps:
% Vector of number of steps to take before re-attempting
% solution of steady-state problem.
% For example, steps = [1, 2, 5, 10] would cause one time
% step to be taken first time the steady-state solution
% attempted. If this failed, two time steps would be taken.
checklib;
calllib(ct, 'sim1D_', s.st_id, ...
stepsize, length(steps), steps);
end
function setValue(s, n, comp, localPoints, v)
% Set the value of a single entry in the solution vector.
%
% Example (assuming 's' is an instance of class 'Stack'):
%
% setValue(s, 3, 5, 1, 5.6);
%
% This sets component 5 at the leftmost point (local point 1)
% in domain 3 to the value 5.6. Note that the local index
% always begins at 1 at the left of each domain, independent of
% the global index of the point, wchih depends on the location
% of this domain in the stack.
%
% parameter s:
% Instance of class 'Stack'.
% parameter n:
% Domain number.
% parameter comp:
% Component number.
% parameter localPoints:
% Local index of the grid point in the domain.
% parameter v:
% Value to be set.
checklib;
calllib(ct, 'sim1D_setValue', s.st_id, ...
n - 1, comp - 1, localPoints - 1, v);
end
function solution(s, domain, component)
% Get a solution component in one domain.
%
% parameter s:
% Instance of class 'Stack'.
% parameter domain:
% String name of the domain from which the solution is
% desired.
% parameter component:
% String component for which the solution is desired. If
% omitted, solution for all of the components will be
% returned in an 'nPoints' x 'nComponnts' array.
% return:
% Either an 'nPoints' x 1 vector, or 'nPoints' x
% 'nCOmponents' array.
checklib;
idom = s.stackIndex(domain);
d = s.domains(idom);
np = d.nPoints;
if nargin == 3
icomp = d.componentIndex(component);
x = zeros(1, np);
for n = 1:np
x(n) = calllib(ct, 'sim1D_Value', s.st_id, ...
idom - 1, icomp - 1, n - 1);
end
else
nc = d.nComponents;
x = zeros(nc, np);
for m = 1:nc
for n = 1:np
x(m, n) = calllib(ct, 'sim1D_Value', s.st_id, ...
idom - 1, m - 1, n - 1);
end
end
end
end
function solve(s, loglevel, refine_grid)
% Solve the problem.
%
% parameter s:
% Instance of class 'Stack'.
% parameter loglevel:
% Integer flag controlling the amount of diagnostic output.
% Zero supresses all output, and 5 produces very verbose
% output.
% parameter refine_grid:
% Integer, 1 to allow grid refinement, 0 to disallow.
checklib;
calllib(ct, 'sim1D_solve', s.st_id, loglevel, refine_grid);
end
function b = subsref(s, index)
% Redefine subscripted references.
switch index.type
case '()'
b = s.domains(index.subs{:});
case '.'
n = s.domainIndex(index.subs);
b = s.domains(n);
otherwise
error('syntax error');
end
end
function writeStats(s)
% Print statistics for the current solution.
% Prints a summary of the number of function and Jacobian
% evaluations for each grid, and the CPU time spent on each
% one.
checklib;
calllib(ct, 'sim1D_writeStats', s.st_id);
end
end
end

17
interfaces/matlab_experimental/1D/Surface.m Normal file
View File

@ -0,0 +1,17 @@
function m = Surface(id, surface_mech)
% Create an surface domain.
% :param surface_mech
% Instance of class 'Interface' defining the surface reaction
% mechanism to be used. Optional.
if nargin < 2
m = Domain1D(3);
if nargin == 0
m.setID('surface');
elseif nargin == 1
m.setID(id);
end
else
m = Domain1D(6, surface_mech);
m.setID(id);
end
end

9
interfaces/matlab_experimental/1D/SymmPlane.m Normal file
View File

@ -0,0 +1,9 @@
function m = SymmPlane(id)
% Create an symmetry plane domain.
m = Domain1D(4);
if nargin == 0
m.setID('symmetry_plane');
else
m.setID(id);
end
end

64
interfaces/matlab_experimental/Base/Interface.m
View File

@ -1,61 +1,55 @@
classdef Interface < handle classdef Interface < handle & ThermoPhase & Kinetics
properties properties
th coverages
kin
end
properties(Constant = true)
lib = 'cantera_shared'
end end
methods methods
%% Interface class constructor %% Interface class constructor
function s = Interface(src, id, p1, p2, p3, p4) function s = Interface(src, id, p1, p2, p3, p4)
% :param src: % parameter src:
% CTI or CTML file containing the interface or edge phase. % CTI or CTML file containing the interface or edge phase.
% :param id: % parameter id:
% Name of the interface or edge phase in the source file. % Name of the interface or edge phase in the source file.
% :param p1/P2/P3/P4: % parameter p1/p2/p3/p4:
% Adjoining phase to the interface; % Adjoining phase to the interface;
% :return: % return:
% Instance of class 'Interface'. % Instance of class 'Interface'.
checklib; checklib;
t = ThermoPhase(src, id); t = ThermoPhase(src, id);
s@ThermoPhase(src, id);
if nargin == 2 if nargin == 2
k = Kinetics(t, src, id); args = {};
elseif nargin == 3 elseif nargin == 3
k = Kinetics(t, src, id, p1); args = {p1};
elseif nargin == 4 elseif nargin == 4
k = Kinetics(t, src, id, p1, p2); args = {p1, p2};
elseif nargin == 5 elseif nargin == 5
k = Kinetics(t, src, id, p1, p2, p3); args = {p1, p2, p3};
elseif nargin == 6 elseif nargin == 6
k = Kinetics(t, src, id, p1, p2, p3, p4); args = {p1, p2, p3, p4};
end end
s@Kinetics(t, src, id, args{:});
s.kin = k;
s.th = t;
end end
%% Interface methods %% Interface methods
function c = coverages(s) function c = get.coverages(s)
% Get the surface coverages of the species on an interface. % Get the surface coverages of the species on an interface.
% %
% :return: % return:
% If no output value is assigned, a bar graph will be % If no output value is assigned, a bar graph will be
% plotted. Otherwise, a vector of length "n_surf_species" % plotted. Otherwise, a vector of length "n_surf_species"
% will be returned. % will be returned.
checklib; checklib;
surf_id = s.th.tr_id; surf_id = s.tr_id;
nsp = s.th.nSpecies; nsp = s.nSpecies;
xx = zeros(1, nsp); xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(s.lib, 'surf_getCoverages', surf_id, xx); calllib(ct, 'surf_getCoverages', surf_id, xx);
c = pt.Value; c = pt.Value;
if nargout == 0 if nargout == 0
@ -63,7 +57,7 @@ classdef Interface < handle
set(gcf, 'Name', 'Coverages') set(gcf, 'Name', 'Coverages')
bar(c); bar(c);
colormap(summer); colormap(summer);
nm = speciesNames(s); nm = s.speciesNames;
set(gca, 'XTickLabel', nm); set(gca, 'XTickLabel', nm);
xlabel('Species Name'); xlabel('Species Name');
ylabel('Coverage'); ylabel('Coverage');
@ -74,17 +68,17 @@ classdef Interface < handle
function c = concentrations(s) function c = concentrations(s)
% Get the concentrations of the species on an interface. % Get the concentrations of the species on an interface.
% %
% :return: % return:
% If no output value is assigned, a bar graph will be % If no output value is assigned, a bar graph will be
% plotted. Otherwise, a vector of length "n_surf_species" % plotted. Otherwise, a vector of length "n_surf_species"
% will be returned. % will be returned.
checklib; checklib;
surf_id = s.th.tr_id; surf_id = s.tr_id;
nsp = s.th.nSpecies; nsp = s.nSpecies;
xx = zeros(1, nsp); xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(s.lib, 'surf_getConcentrations', surf_id, xx); calllib(ct, 'surf_getConcentrations', surf_id, xx);
c = pt.Value; c = pt.Value;
if nargout == 0 if nargout == 0
@ -100,10 +94,10 @@ classdef Interface < handle
end end
end end
function setCoverages(s, cov, norm) function set.coverages(s, cov, norm)
% Set surface coverages of the species on an interface. % Set surface coverages of the species on an interface.
% %
% :param cov: % parameter cov:
% Coverage of the species. "Cov" can be either a vector of % Coverage of the species. "Cov" can be either a vector of
% length "n_surf_species", or a string in the format % length "n_surf_species", or a string in the format
% "Species:Coverage, Species:Coverage". % "Species:Coverage, Species:Coverage".
@ -116,22 +110,22 @@ classdef Interface < handle
norm_flag = 1; norm_flag = 1;
end end
surf_id = s.th.tr_id; surf_id = s.tr_id;
nsp = s.th.nSpecies; nsp = s.nSpecies;
[m, n] = size(cov); [m, n] = size(cov);
if isa(cov, 'double') if isa(cov, 'double')
sz = length(cov); sz = length(cov);
if sz == nsp if sz == nsp
if ((m == nsp && n == 1) || (m == 1 & n == nsp)) if ((m == nsp && n == 1) || (m == 1 & n == nsp))
calllib(s.lib, 'surf_setCoverages', surf_id, cov, norm_flag); calllib(ct, 'surf_setCoverages', surf_id, cov, norm_flag);
else error('wrong size for coverage array'); else error('wrong size for coverage array');
end end
else else
error('wrong size for coverage array'); error('wrong size for coverage array');
end end
elseif isa(cov, 'char') elseif isa(cov, 'char')
calllib(s.lib, 'surf_setCoveragesByName', surf_id, cov); calllib(ct, 'surf_setCoveragesByName', surf_id, cov);
end end
end end

423
interfaces/matlab_experimental/Base/Kinetics.m
View File

@ -1,7 +1,8 @@
classdef Kinetics < handle classdef Kinetics < handle
properties properties
owner owner
id kin_id
Kc % equilibrium constant Kc % equilibrium constant
Kf % forward reaction rate Kf % forward reaction rate
Kr % reverse reaction rate Kr % reverse reaction rate
@ -9,9 +10,7 @@ classdef Kinetics < handle
dS % entropy of reaction dS % entropy of reaction
dG % gibbs free energy of reaction dG % gibbs free energy of reaction
end end
properties(Constant = true)
lib = 'cantera_shared'
end
methods methods
%% Kinetics class constructor %% Kinetics class constructor
@ -41,78 +40,78 @@ classdef Kinetics < handle
end end
end end
end end
kin.id = calllib(kin.lib, 'kin_newFromFile', src, id, ... kin.kin_id = calllib(ct, 'kin_newFromFile', src, id, ...
iph, inb1, inb2, inb3, inb4); iph, inb1, inb2, inb3, inb4);
end end
%% Utility methods %% Utility methods
function clear(kin) function kin_clear(kin)
% Delete the kernel object % Delete the kernel object
checklib; checklib;
calllib(kin.lib, 'kin_del', kin.id); calllib(ct, 'kin_del', kin.kin_id);
end end
%% Get scalar attributes %% Get scalar attributes
function n = nReactions(kin) function n = isReversible(kin, i)
% Get the number of reactions. % Get an array of flags indicating reversibility of a reaction.
% %
% :return: % parameter i:
% Integer number of reactions % Integer reaction number.
% return:
% 1 if reaction number i is reversible. 0 if irreversible.
checklib; checklib;
n = calllib(kin.lib, 'kin_nReactions', kin.id); n = calllib(ct, 'kin_isReversible', kin.kin_id, i);
end end
function n = multiplier(kin, irxn) function n = multiplier(kin, irxn)
% Get the multiplier for reaction rate of progress. % Get the multiplier for reaction rate of progress.
% %
% :param irxn: % parameter irxn:
% Integer reaction number for which the multiplier is % Integer reaction number for which the multiplier is
% desired. % desired.
% :return: % return:
% Multiplier of the rate of progress of reaction irxn. % Multiplier of the rate of progress of reaction irxn.
checklib; checklib;
n = calllib(kin.lib, 'kin_multiplier', kin.id, irxn-1); n = calllib(ct, 'kin_multiplier', kin.kin_id, irxn-1);
end end
function n = nSpecies(kin) function n = nReactions(kin)
% Get the number of reactions.
%
% return:
% Integer number of reactions
checklib;
n = calllib(ct, 'kin_nReactions', kin.kin_id);
end
function n = nSpecies2(kin)
% Get the total number of species. % Get the total number of species.
% %
% :return: % return:
% Integer total number of species. % Integer total number of species.
checklib; checklib;
n = calllib(kin.lib, 'kin_nSpecies', kin.id); n = calllib(ct, 'kin_nSpecies', kin.kin_id);
end
function n = isReversible(kin, i)
% Get an array of flags indicating reversibility of a reaction.
%
% :param i:
% Integer reaction number.
% :return:
% 1 if reaction number i is reversible. 0 if irreversible.
checklib;
n = calllib(kin.lib, 'kin_isReversible', kin.id, i);
end end
function n = stoich_r(kin, species, rxns) function n = stoich_r(kin, species, rxns)
% Get the reactant stoichiometric coefficients. % Get the reactant stoichiometric coefficients.
% %
% :param species: % parameter species:
% Species indices for which reactant stoichiometric % Species indices for which reactant stoichiometric
% coefficients should be retrieved. Optional argument; if % coefficients should be retrieved. Optional argument; if
% specified, "rxns" must be specified as well. % specified, "rxns" must be specified as well.
% :param rxns: % parameter rxns:
% Reaction indicies for which reactant stoichiometric % Reaction indicies for which reactant stoichiometric
% coefficients should be retrieved. Optional argument; if % coefficients should be retrieved. Optional argument; if
% specified, "species" must be specified as well. % specified, "species" must be specified as well.
% :return: % return:
% Returns a sparse matrix of all reactant stoichiometric % Returns a sparse matrix of all reactant stoichiometric
% coefficients. The matrix elements "nu(k, i)" is the % coefficients. The matrix elements "nu(k, i)" is the
% stoichiometric coefficient of species k as a reactant in % stoichiometric coefficient of species k as a reactant in
@ -137,8 +136,8 @@ classdef Kinetics < handle
for k = krange for k = krange
for i = irange for i = irange
t = calllib(kin.lib, 'kin_reactantStoichCoeff', ... t = calllib(ct, 'kin_reactantStoichCoeff', ...
kin.id, k-1, i-1); kin.kin_id, k-1, i-1);
if t ~= 0.0 if t ~= 0.0
temp(k, i) = t; temp(k, i) = t;
end end
@ -151,15 +150,15 @@ classdef Kinetics < handle
function n = stoich_p(kin, species, rxns) function n = stoich_p(kin, species, rxns)
% Get the product stoichiometric coefficients. % Get the product stoichiometric coefficients.
% %
% :param species: % parameter species:
% Species indices for which product stoichiometric % Species indices for which product stoichiometric
% coefficients should be retrieved. Optional argument; if % coefficients should be retrieved. Optional argument; if
% specified, "rxns" must be specified as well. % specified, "rxns" must be specified as well.
% :param rxns: % parameter rxns:
% Reaction indicies for which product stoichiometric % Reaction indicies for which product stoichiometric
% coefficients should be retrieved. Optional argument; if % coefficients should be retrieved. Optional argument; if
% specified, "species" must be specified as well. % specified, "species" must be specified as well.
% :return: % return:
% Returns a sparse matrix of all product stoichiometric % Returns a sparse matrix of all product stoichiometric
% coefficients. The matrix elements "nu(k, i)" is the % coefficients. The matrix elements "nu(k, i)" is the
% stoichiometric coefficient of species k as a product in % stoichiometric coefficient of species k as a product in
@ -184,8 +183,8 @@ classdef Kinetics < handle
for k = krange for k = krange
for i = irange for i = irange
t = calllib(kin.lib, 'kin_productStoichCoeff', ... t = calllib(ct, 'kin_productStoichCoeff', ...
kin.id, k-1, i-1); kin.kin_id, k-1, i-1);
if t ~= 0.0 if t ~= 0.0
temp(k, i) = t; temp(k, i) = t;
end end
@ -198,15 +197,15 @@ classdef Kinetics < handle
function n = stoich_net(kin, species, rxns) function n = stoich_net(kin, species, rxns)
% Get the net stoichiometric coefficients. % Get the net stoichiometric coefficients.
% %
% :param species: % parameter species:
% Species indices for which net stoichiometric coefficients % Species indices for which net stoichiometric coefficients
% should be retrieved. Optional argument; if specified, % should be retrieved. Optional argument; if specified,
% "rxns" must be specified as well. % "rxns" must be specified as well.
% :param rxns: % parameter rxns:
% Reaction indicies for which net stoichiometric % Reaction indicies for which net stoichiometric
% coefficients should be retrieved. Optional argument; if % coefficients should be retrieved. Optional argument; if
% specified, "species" must be specified as well. % specified, "species" must be specified as well.
% :return: % return:
% Returns a sparse matrix of all net stoichiometric % Returns a sparse matrix of all net stoichiometric
% coefficients. The matrix elements "nu(k, i)" is the % coefficients. The matrix elements "nu(k, i)" is the
% stoichiometric coefficient of species k as a net in % stoichiometric coefficient of species k as a net in
@ -227,10 +226,82 @@ classdef Kinetics < handle
%% Get reaction array attributes %% Get reaction array attributes
function cdot = creationRates(kin)
% Get the chemical reaction rates.
%
% return:
% Returns a vector of the creation rates of all species. If
% the output is not assigned to a variable, a bar graph is
% produced. Unit: kmol/m^3-s.
checklib;
nsp = kin.nSpecies;
xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getCreationRates', kin.kin_id, nsp, pt);
cdot = pt.Value;
if nargout == 0
figure
set(gcf, 'Name', 'Creation Rates')
bar(q)
xlabel('Species Number')
ylabel('Creation Rate [kmol/m^3-s]')
title('Species Chemical Reaction Rates')
end
end
function ddot = destructionRates(kin)
% Get the chemical destruction rates.
%
% return:
% Returns a vector of the destruction rates of all species.
% If the output is not assigned to a variable, a bar graph
% is produced. Unit: kmol/m^3-s.
checklib;
nsp = kin.nSpecies;
xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getDestructionRates', kin.kin_id, nsp, pt);
ddot = pt.Value;
if nargout == 0
figure
set(gcf, 'Name', 'Destruction Rates')
bar(q)
xlabel('Species Number')
ylabel('Destruction Rate [kmol/m^3-s]')
title('Species Chemical Reaction Rates')
end
end
function wdot = netProdRates(kin)
% Get the net chemical production rates for all species.
%
% return:
% Returns a vector of the net production (creation-destruction)
% rates of all species. If the output is not assigned to a
% variable, a bar graph is produced. Unit: kmol/m^3-s.
checklib;
nsp = kin.nSpecies;
xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getNetProductionRates', kin.kin_id, nsp, pt);
wdot = pt.Value;
if nargout == 0
figure
set(gcf, 'Name', 'Production Rates')
bar(q)
xlabel('Species Number')
ylabel('Net Production Rate [kmol/m^3-s]')
title('Species Net Chemical Reaction Rates')
end
end
function q = rop_f(kin) function q = rop_f(kin)
% Get the forward rates of progress for all reactions. % Get the forward rates of progress for all reactions.
% %
% :return: % return:
% Returns a column vector of the forward rates of progress % Returns a column vector of the forward rates of progress
% for all reactions. If this function is called without % for all reactions. If this function is called without
% argument, a bar graph is produced. % argument, a bar graph is produced.
@ -239,7 +310,7 @@ classdef Kinetics < handle
nr = kin.nReactions; nr = kin.nReactions;
xx = zeros(1, nr); xx = zeros(1, nr);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getFwdRateOfProgress', kin.id, nr, pt); calllib(ct, 'kin_getFwdRateOfProgress', kin.kin_id, nr, pt);
q = pt.Value; q = pt.Value;
if nargout == 0 if nargout == 0
figure figure
@ -254,7 +325,7 @@ classdef Kinetics < handle
function q = rop_r(kin) function q = rop_r(kin)
% Get the reverse rates of progress for all reactions. % Get the reverse rates of progress for all reactions.
% %
% :return: % return:
% Returns a column vector of the reverse rates of progress % Returns a column vector of the reverse rates of progress
% for all reactions. If this function is called without % for all reactions. If this function is called without
% argument, a bar graph is produced. % argument, a bar graph is produced.
@ -263,7 +334,7 @@ classdef Kinetics < handle
nr = kin.nReactions; nr = kin.nReactions;
xx = zeros(1, nr); xx = zeros(1, nr);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getRevRateOfProgress', kin.id, nr, pt); calllib(ct, 'kin_getRevRateOfProgress', kin.kin_id, nr, pt);
q = pt.Value; q = pt.Value;
if nargout == 0 if nargout == 0
figure figure
@ -278,7 +349,7 @@ classdef Kinetics < handle
function q = rop(kin) function q = rop(kin)
% Get the forward and reverse rates of progress. % Get the forward and reverse rates of progress.
% %
% :return: % return:
% Returns an I x 2 array of reaction rates of progress, % Returns an I x 2 array of reaction rates of progress,
% where I is the number of reactions. The first column % where I is the number of reactions. The first column
% contains the forward rates of progress, and the second % contains the forward rates of progress, and the second
@ -303,7 +374,7 @@ classdef Kinetics < handle
function q = rop_net(kin) function q = rop_net(kin)
% Get the net rates of progress for all reactions. % Get the net rates of progress for all reactions.
% %
% :return: % return:
% Returns a column vector of the net rates of progress % Returns a column vector of the net rates of progress
% for all reactions. If this function is called without % for all reactions. If this function is called without
% argument, a bar graph is produced. % argument, a bar graph is produced.
@ -312,7 +383,7 @@ classdef Kinetics < handle
nr = kin.nReactions; nr = kin.nReactions;
xx = zeros(1, nr); xx = zeros(1, nr);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getNetRateOfProgress', kin.id, nr, pt); calllib(ct, 'kin_getNetRatesOfProgress', kin.kin_id, nr, pt);
q = pt.Value; q = pt.Value;
if nargout == 0 if nargout == 0
figure figure
@ -324,10 +395,91 @@ classdef Kinetics < handle
end end
end end
function rxn = reactionEqn(kin, irxn)
% Get the reaction equation of a reaction
%
% parameter irxn:
% Optional. Integer or vector of reaction numbers.
% return:
% String or cell arrray of strings of the reaction
% equations.
if nargin == 1
m = kin.nReactions;
n = 1
irxn = (n:m)';
elseif nargin == 2
if isa(irxn, 'double')
[m, n] = size(irxn);
else
error('reaction numbers must be numeric');
end
end
rxn = cell(m, n);
for i = 1:m
for j = 1:n
buflen = calllib(ct, 'kin_getReactionString', kin.kin_id, ...
irxn - 1, 0, '');
if buflen > 0
aa = char(zeros(1, buflen));
[~, aa] = calllib(ct, 'kin_getReactionString', ...
kin.kin_id, irxn - 1, buflen, aa);
rxn{i, j} = aa;
end
end
end
end
function enthalpy = get.dH(kin)
% Get the enthalpy of reaction for each reaction.
%
% return:
% Returns a vector of the enthalpy of reaction for each
% reaction. Unit: J/kmol.
checklib;
nr = kin.nReactions;
xx = zeros(1, nr);
pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getDelta', kin.kin_id, 0, nr, pt);
enthalpy = pt.Value;
end
function entropy = get.dS(kin)
% Get the entropy of reaction for each reaction.
%
% return:
% Returns a vector of the entropy of reaction for each
% reaction. Unit: J/kmol-K.
checklib;
nr = kin.nReactions;
xx = zeros(1, nr);
pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getDelta', kin.kin_id, 2, nr, pt);
entropy = pt.Value;
end
function gibbs = get.dG(kin)
% Get the Gibbs free energy of reaction for each reaction.
%
% return:
% Returns a vector of the Gibbs free energy of reaction for
% each reaction. Unit: J/kmol.
checklib;
nr = kin.nReactions;
xx = zeros(1, nr);
pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getDelta', kin.kin_id, 1, nr, pt);
gibbs = pt.Value;
end
function k = get.Kc(kin) function k = get.Kc(kin)
% Get the equilibrium constants for all reactions. % Get the equilibrium constants for all reactions.
% %
% :return: % return:
% Returns a column vector of the equilibrium constants for % Returns a column vector of the equilibrium constants for
% all reactions. The vector has an entry for every reaction, % all reactions. The vector has an entry for every reaction,
% whether reversible or not, but non-zero values occur only % whether reversible or not, but non-zero values occur only
@ -338,12 +490,12 @@ classdef Kinetics < handle
nr = kin.nReactions; nr = kin.nReactions;
xx = zeros(1, nr); xx = zeros(1, nr);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getEquilibriumConstants', kin.id, nr, pt); calllib(ct, 'kin_getEquilibriumConstants', kin.kin_id, nr, pt);
k = pt.Value; k = pt.Value;
if nargout == 0 if nargout == 0
figure figure
set(gcf, 'Name', 'Equilibrium Constants') set(gcf, 'Name', 'Equilibrium Constants')
bar(q) bar(k)
xlabel('Reaction Number') xlabel('Reaction Number')
ylabel('log_{10} Kc [kmol,m, s]') ylabel('log_{10} Kc [kmol,m, s]')
title('Equilibrium Constants') title('Equilibrium Constants')
@ -353,7 +505,7 @@ classdef Kinetics < handle
function k = get.Kf(kin) function k = get.Kf(kin)
% Get the forward reaction rate constants for all reactions. % Get the forward reaction rate constants for all reactions.
% %
% :return: % return:
% Returns a column vector of the forward rates constants of % Returns a column vector of the forward rates constants of
% all of the reactions. % all of the reactions.
@ -361,14 +513,14 @@ classdef Kinetics < handle
nr = kin.nReactions; nr = kin.nReactions;
xx = zeros(1, nr); xx = zeros(1, nr);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getFwdRateConstants', kin.id, nr, pt); calllib(ct, 'kin_getFwdRateConstants', kin.kin_id, nr, pt);
k = pt.Value; k = pt.Value;
end end
function k = get.Kr(kin) function k = get.Kr(kin)
% Get the reverse reaction rate constants for all reactions. % Get the reverse reaction rate constants for all reactions.
% %
% :return: % return:
% Returns a column vector of the reverse rates constants of % Returns a column vector of the reverse rates constants of
% all of the reactions. % all of the reactions.
@ -376,172 +528,33 @@ classdef Kinetics < handle
nr = kin.nReactions; nr = kin.nReactions;
xx = zeros(1, nr); xx = zeros(1, nr);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getRevRateConstants', kin.id, 1, nr, pt); calllib(ct, 'kin_getRevRateConstants', kin.kin_id, 1, nr, pt);
k = pt.Value; k = pt.Value;
end end
function enthalpy = get.dH(kin) function massProdRate = ydot(kin)
% Get the enthalpy of reaction for each reaction.
%
% :return:
% Returns a vector of the enthalpy of reaction for each
% reaction. Unit: J/kmol.
checklib;
nr = kin.nReactions;
xx = zeros(1, nr);
pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getDelta', kin.id, 0, nr, pt);
enthalpy = pt.Value;
end
function entropy = get.dS(kin)
% Get the entropy of reaction for each reaction.
%
% :return:
% Returns a vector of the entropy of reaction for each
% reaction. Unit: J/kmol-K.
checklib;
nr = kin.nReactions;
xx = zeros(1, nr);
pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getDelta', kin.id, 2, nr, pt);
entropy = pt.Value;
end
function gibbs = get.dG(kin)
% Get the Gibbs free energy of reaction for each reaction.
%
% :return:
% Returns a vector of the Gibbs free energy of reaction for
% each reaction. Unit: J/kmol.
checklib;
nr = kin.nReactions;
xx = zeros(1, nr);
pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getDelta', kin.id, 1, nr, pt);
gibbs = pt.Value;
end
function cdot = creationRates(kin)
% Get the chemical reaction rates.
%
% :return:
% Returns a vector of the creation rates of all species. If
% the output is not assigned to a variable, a bar graph is
% produced. Unit: kmol/m^3-s.
checklib;
nsp = kin.nSpecies;
xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getCreationRates', kin.id, nsp, pt);
cdot = pt.Value;
if nargout == 0
figure
set(gcf, 'Name', 'Creation Rates')
bar(q)
xlabel('Species Number')
ylabel('Creation Rate [kmol/m^3-s]')
title('Species Chemical Reaction Rates')
end
end
function ddot = destructionRates(kin)
% Get the chemical destruction rates.
%
% :return:
% Returns a vector of the destruction rates of all species.
% If the output is not assigned to a variable, a bar graph
% is produced. Unit: kmol/m^3-s.
checklib;
nsp = kin.nSpecies;
xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getDestructionRates', kin.id, nsp, pt);
ddot = pt.Value;
if nargout == 0
figure
set(gcf, 'Name', 'Destruction Rates')
bar(q)
xlabel('Species Number')
ylabel('Destruction Rate [kmol/m^3-s]')
title('Species Chemical Reaction Rates')
end
end
function wdot = netProdRates(kin)
% Get the net chemical production rates for all species.
%
% :return:
% Returns a vector of the net production (creation-destruction)
% rates of all species. If the output is not assigned to a
% variable, a bar graph is produced. Unit: kmol/m^3-s.
checklib;
nsp = kin.nSpecies;
xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getNetProductionRates', kin.id, nsp, pt);
wdot = pt.Value;
if nargout == 0
figure
set(gcf, 'Name', 'Production Rates')
bar(q)
xlabel('Species Number')
ylabel('Net Production Rate [kmol/m^3-s]')
title('Species Net Chemical Reaction Rates')
end
end
function ydot = massProdRates(kin)
% Get the mass production rates of the species. % Get the mass production rates of the species.
% %
% :return: % return:
% Returns a vector of the mass production rates. % Returns a vector of the mass production rates.
checklib; checklib;
nsp = kin.nSpecies; nsp = kin.nSpecies;
xx = zeros(1, nsp); xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getSourceTerms', kin.id, nsp, pt); calllib(ct, 'kin_getSourceTerms', kin.kin_id, nsp, pt);
ydot = pt.Value; massProdRate = pt.Value;
end end
% function e = reactionEqn(kin, irxn)
% % Get the reaction equation of a reaction
% %
% % :param irxn:
% % Optional. Integer or vector of reaction numbers.
% % :return:
% % String or cell arrray of strings of the reaction
% % equations.
% Need to resolve string printing issues first!
%
% if nargin == 1
% nr = kin.nReactions;
% irxn = (1:m)';
% elseif nargin == 2
% if isa(irxn, 'double')
% [m, n] = size(irxn);
% else
% error('reaction numbers must be numeric');
% end
% end
% end
%% Set attributes %% Set attributes
function n = setMultiplier(kin, irxn, v) function n = setMultiplier(kin, irxn, v)
% Set the multiplier for the reaction rate of progress. % Set the multiplier for the reaction rate of progress.
% %
% :param irxn: % parameter irxn:
% Integer of vector reaction numbers for which the % Integer of vector reaction numbers for which the
% multiplier should be set. Optional. % multiplier should be set. Optional.
% :param v: % parameter v:
% Value by which the reaction rate of progress should be % Value by which the reaction rate of progress should be
% multiplied. % multiplied.
@ -558,7 +571,7 @@ classdef Kinetics < handle
for i = 1:nr for i = 1:nr
for j = 1:n for j = 1:n
calllib(kin.lib, 'kin_setMultiplier', kin.id, ... calllib(ct, 'kin_setMultiplier', kin.kin_id, ...
irxn(i, j)-1, v); irxn(i, j)-1, v);
end end
end end
@ -567,10 +580,10 @@ classdef Kinetics < handle
function advanceCoeverages(kin, dt) function advanceCoeverages(kin, dt)
% Advance the surface coveages forward in time % Advance the surface coveages forward in time
% %
% :param dt: % parameter dt:
% Time interval by which the coverages should be advanced. % Time interval by which the coverages should be advanced.
calllib(kin.lib, 'kin_advanceCoverages', kin.id, dt); calllib(ct, 'kin_advanceCoverages', kin.kin_id, dt);
end end
end end

84
interfaces/matlab_experimental/Base/Mixture.m
View File

@ -7,10 +7,6 @@ classdef Mixture < handle
P P
end end
properties(Constant = true)
lib = 'cantera_shared'
end
methods methods
%% Mixture class constructor %% Mixture class constructor
@ -38,9 +34,9 @@ classdef Mixture < handle
% store its own state information locally, and syncrhonizes the % store its own state information locally, and syncrhonizes the
% phase objects whenever itrequires phase properties. % phase objects whenever itrequires phase properties.
% %
% :param phases: % parameter phases:
% Cell array of phases and mole numbers. % Cell array of phases and mole numbers.
% :return: % return:
% Instance of class 'Mixture'. % Instance of class 'Mixture'.
checklib; checklib;
@ -50,7 +46,7 @@ classdef Mixture < handle
end end
% Create an empty mixture. % Create an empty mixture.
m.mixindex = calllib(m.lib, 'mix_new'); m.mixindex = calllib(ct, 'mix_new');
m.phases = phases; m.phases = phases;
% If phases are supplied, add them % If phases are supplied, add them
@ -78,7 +74,7 @@ classdef Mixture < handle
function display(m) function display(m)
% Display the state of the mixture on the terminal. % Display the state of the mixture on the terminal.
calllib(m.lib, 'mix_updatePhases', m.mixindex); calllib(ct, 'mix_updatePhases', m.mixindex);
[np, nc] = size(m.phases); [np, nc] = size(m.phases);
for n = 1:np for n = 1:np
s = [sprintf('\n******************* Phase %d', n) ... s = [sprintf('\n******************* Phase %d', n) ...
@ -93,7 +89,7 @@ classdef Mixture < handle
% Delete the MultiPhase object. % Delete the MultiPhase object.
checklib; checklib;
calllib(m.lib, 'mix_del', m.mixindex); calllib(ct, 'mix_del', m.mixindex);
end end
%% Mixture Get methods %% Mixture Get methods
@ -101,11 +97,11 @@ classdef Mixture < handle
function addPhase(m, phase, moles) function addPhase(m, phase, moles)
% Add a phase to the mixture % Add a phase to the mixture
% %
% :param m: % parameter m:
% Instance of class 'Mixture' to which phases is added. % Instance of class 'Mixture' to which phases is added.
% :param phase: % parameter phase:
% Instance of class 'ThermoPhase' which should be added. % Instance of class 'ThermoPhase' which should be added.
% :param moles: % parameter moles:
% Number of moles of the phase to be added. Unit: kmol. % Number of moles of the phase to be added. Unit: kmol.
checklib; checklib;
@ -120,7 +116,7 @@ classdef Mixture < handle
error('Negative moles'); error('Negative moles');
end end
iok = calllib(m.lib, 'mix_addPhase', m.mixindex, phase.tp_id, ... iok = calllib(ct, 'mix_addPhase', m.mixindex, phase.tp_id, ...
moles); moles);
if iok < 0 if iok < 0
error('Error adding phase'); error('Error adding phase');
@ -130,39 +126,39 @@ classdef Mixture < handle
function temperature = get.T(m) function temperature = get.T(m)
% Get the temperature of the mixture. % Get the temperature of the mixture.
% %
% :return: % return:
% Temperature in K. % Temperature in K.
checklib; checklib;
temperature = calllib(m.lib, 'mix_temperature', m.mixindex); temperature = calllib(ct, 'mix_temperature', m.mixindex);
end end
function pressure = get.P(m) function pressure = get.P(m)
% Get the pressure of themixture. % Get the pressure of themixture.
% %
% :return: % return:
% Pressure in Pa. % Pressure in Pa.
checklib; checklib;
pressure = calllib(m.lib, 'mix_pressure', m.mixindex); pressure = calllib(ct, 'mix_pressure', m.mixindex);
end end
function n = nPhases(m) function n = nPhases(m)
% Get the number of phases in the mixture. % Get the number of phases in the mixture.
checklib; checklib;
n = calllib(m.lib, 'mix_nPhases', m.mixindex); n = calllib(ct, 'mix_nPhases', m.mixindex);
end end
function n = nElements(m) function n = nElements(m)
% Get the number of elements in the mixture. % Get the number of elements in the mixture.
checklib; checklib;
n = calllib(m.lib, 'mix_nElements', m.mixindex); n = calllib(ct, 'mix_nElements', m.mixindex);
end end
function n = nSpecies(m) function n = nSpecies(m)
% Get the number of species in the mixture. % Get the number of species in the mixture.
checklib; checklib;
n = calllib(m.lib, 'mix_nSpecies', m.mixindex); n = calllib(ct, 'mix_nSpecies', m.mixindex);
end end
function n = elementIndex(m, name) function n = elementIndex(m, name)
@ -172,7 +168,7 @@ classdef Mixture < handle
% indices start from 1 instead of 0 as in Cantera C++ and % indices start from 1 instead of 0 as in Cantera C++ and
% Python interfaces. % Python interfaces.
checklib; checklib;
n = calllib(m.lib, 'mix_elementIndex', m.mixindex, name) + 1; n = calllib(ct, 'mix_elementIndex', m.mixindex, name) + 1;
end end
function n = speciesIndex(m, k, p) function n = speciesIndex(m, k, p)
@ -182,26 +178,26 @@ classdef Mixture < handle
% indices start from 1 instead of 0 as in Cantera C++ and % indices start from 1 instead of 0 as in Cantera C++ and
% Python interfaces. % Python interfaces.
checklib; checklib;
n = calllib(m.lib, 'mix_speciesIndex', m.mixindex, k-1, p-1) + 1; n = calllib(ct, 'mix_speciesIndex', m.mixindex, k-1, p-1) + 1;
% check back on this one! % check back on this one!
end end
function moles = phaseMoles(m, n) function moles = phaseMoles(m, n)
% Get the number of moles of a phase in a mixture. % Get the number of moles of a phase in a mixture.
% %
% :param n: % parameter n:
% Integer phase number in the input. % Integer phase number in the input.
% :return: % return:
% Moles of phase number 'n'. Unit: kmol. % Moles of phase number 'n'. Unit: kmol.
checklib; checklib;
if nargin == 2 if nargin == 2
moles = calllib(m.lib, 'mix_phaseMoles', m.mixindex, n); moles = calllib(ct, 'mix_phaseMoles', m.mixindex, n);
elseif nargin == 1 elseif nargin == 1
np = m.nPhases; np = m.nPhases;
moles = zeros(1, np); moles = zeros(1, np);
for i = 1:np for i = 1:np
moles(i) = calllib(m.lib, 'mix_phaseMoles', ... moles(i) = calllib(ct, 'mix_phaseMoles', ...
m.mixindex, i); m.mixindex, i);
end end
else error('wrong number of arguments'); else error('wrong number of arguments');
@ -211,14 +207,14 @@ classdef Mixture < handle
function mu = chemPotentials(m) function mu = chemPotentials(m)
% Get the chemical potentials of species in the mixture. % Get the chemical potentials of species in the mixture.
% %
% :return: % return:
% Vector of chemical potentials. Unit: J/kmol. % Vector of chemical potentials. Unit: J/kmol.
checklib; checklib;
nsp = m.nSpecies; nsp = m.nSpecies;
xx = zeros(1, nsp); xx = zeros(1, nsp);
ptr = libpointer('doublePtr', xx); ptr = libpointer('doublePtr', xx);
calllib(m.lib, 'mix_getChemPotential', m.mixindex, nsp, ptr); calllib(ct, 'mix_getChemPotential', m.mixindex, nsp, ptr);
mu = ptr.Value; mu = ptr.Value;
end end
@ -227,30 +223,30 @@ classdef Mixture < handle
function m = set.T(m, temp) function m = set.T(m, temp)
% Set the mixture temperature. % Set the mixture temperature.
% %
% :param temp: % parameter temp:
% Temperature to set. Unit: K. % Temperature to set. Unit: K.
checklib; checklib;
calllib(m.lib, 'mix_setTemperature', m.mixindex, temp); calllib(ct, 'mix_setTemperature', m.mixindex, temp);
end end
function m = set.P(m, pressure) function m = set.P(m, pressure)
% Set the mixture pressure. % Set the mixture pressure.
% %
% :param pressure: % parameter pressure:
% Pressure to set. Unit: Pa. % Pressure to set. Unit: Pa.
checklib; checklib;
calllib(m.lib, 'mix_setPressure', m.mixindex, pressure); calllib(ct, 'mix_setPressure', m.mixindex, pressure);
end end
function setPhaseMoles(m, n, moles) function setPhaseMoles(m, n, moles)
% Set the number of moles of phase n in the mixture. % Set the number of moles of phase n in the mixture.
% %
% :param n: % parameter n:
% Phase number. % Phase number.
% :param moles: % parameter moles:
% Number of moles to set. Unit: kmol. % Number of moles to set. Unit: kmol.
checklib; checklib;
calllib(m.lib, 'mix_setPhaseMoles', m.mixindex, n-1, moles); calllib(ct, 'mix_setPhaseMoles', m.mixindex, n-1, moles);
end end
function setSpeciesMoles(m, moles) function setSpeciesMoles(m, moles)
@ -260,10 +256,10 @@ classdef Mixture < handle
% in the mixture. Note that the species may belong to any % in the mixture. Note that the species may belong to any
% phase, and unspecified species are set to zero. % phase, and unspecified species are set to zero.
% %
% :param moles: % parameter moles:
% Vector or string specifying the moles of species. % Vector or string specifying the moles of species.
checklib; checklib;
calllib(m.lib, 'mix_setMolesByName', m.mixindex, moles); calllib(ct, 'mix_setMolesByName', m.mixindex, moles);
% check back on this one! % check back on this one!
end end
@ -285,25 +281,25 @@ classdef Mixture < handle
% >> equilibrate(mix, 'TP); % >> equilibrate(mix, 'TP);
% >> equilibrate('TP', 1.0e-6, 500); % >> equilibrate('TP', 1.0e-6, 500);
% %
% :param XY: % parameter XY:
% Two-letter string specifying the two properties to hold % Two-letter string specifying the two properties to hold
% fixed. Currently 'TP', 'HP', 'TV', and 'SP' have been % fixed. Currently 'TP', 'HP', 'TV', and 'SP' have been
% implemented. Default: 'TP'. % implemented. Default: 'TP'.
% :param err: % parameter err:
% Error tolerance. Iteration will continue until delta_Mu/RT % Error tolerance. Iteration will continue until delta_Mu/RT
% is less than this value for each reaction. Default: % is less than this value for each reaction. Default:
% 1.0e-9. % 1.0e-9.
% :param maxsteps: % parameter maxsteps:
% Maximum number of steps to take while solving the % Maximum number of steps to take while solving the
% equilibrium problem for specified T and P. Default: 1000. % equilibrium problem for specified T and P. Default: 1000.
% :param maxiter: % parameter maxiter:
% Maximum number of temperature and/or pressure iterations. % Maximum number of temperature and/or pressure iterations.
% This is only relevant if a property pair other than (T, % This is only relevant if a property pair other than (T,
% P)is specified. Default: 200. % P)is specified. Default: 200.
% :param loglevel: % parameter loglevel:
% Set to a value > 0 to write diagnostic output. Larger % Set to a value > 0 to write diagnostic output. Larger
% values generate more detailed information. % values generate more detailed information.
% :return: % return:
% The error in the solution. % The error in the solution.
checklib; checklib;
@ -323,7 +319,7 @@ classdef Mixture < handle
if nargin < 2 if nargin < 2
XY = 'TP' XY = 'TP'
end end
r = calllib(m.lib, 'mix_equilibrate', m.mixindex, XY, err, ... r = calllib(ct, 'mix_equilibrate', m.mixindex, XY, err, ...
maxsteps, maxiter, loglevel); maxsteps, maxiter, loglevel);
end end

27
interfaces/matlab_experimental/Base/Solution.m
View File

@ -1,8 +1,6 @@
classdef Solution < handle classdef Solution < handle & ThermoPhase & Kinetics & Transport
properties properties
thermo tp
kinetics
transport
end end
methods methods
% Solution class constructor % Solution class constructor
@ -11,27 +9,26 @@ classdef Solution < handle
id = '-'; id = '-';
end end
tp = ThermoPhase(src, id); tp = ThermoPhase(src, id);
kin = Kinetics(tp, src, id); s@ThermoPhase(src, id);
s.kinetics = kin; s@Kinetics(tp, src, id);
s.thermo = tp;
if nargin == 3 if nargin == 3
if (strcmp(trans, 'default') || strcmp(trans, 'None')... if ~(strcmp(trans, 'default') || strcmp(trans, 'None')...
|| strcmp(trans, 'Mix') || strcmp(trans, 'Multi')) || strcmp(trans, 'Mix') || strcmp(trans, 'Multi'))
tr = Transport(tp, trans, 0);
else
error('Unknown transport modelling specified.'); error('Unknown transport modelling specified.');
end end
else else
tr = Transport(tp, 'default', 0); trans = 'default';
end end
s.transport = tr; s@Transport(tp, trans, 0);
s.tp_id = tp.tp_id;
end end
% Delete the kernel objects associated with a solution % Delete the kernel objects associated with a solution
function clear(s) function clear(s)
s.thermo.clear; s.tp_clear;
s.kinetics.clear; s.kin_clear;
s.transport.clear; s.tr_clear;
end end
end end
end end

1441
interfaces/matlab_experimental/Base/ThermoPhase.m
View File

File diff suppressed because it is too large Load Diff

56
interfaces/matlab_experimental/Base/Transport.m
View File

@ -1,13 +1,10 @@
classdef Transport < handle classdef Transport < handle
properties properties
th th
tr_id tr_id
end end
properties(Constant = true)
lib = 'cantera_shared'
end
methods methods
%% Transport class constructor %% Transport class constructor
@ -28,22 +25,23 @@ classdef Transport < handle
else else
tr.th = tp; tr.th = tp;
if strcmp(model, 'default') if strcmp(model, 'default')
tr.tr_id = calllib(tr.lib, 'trans_newDefault', ... tr.tr_id = calllib(ct, 'trans_newDefault', ...
tp.tp_id, loglevel); tp.tp_id, loglevel);
else else
tr.tr_id = calllib(tr.lib, 'trans_new', model, ... tr.tr_id = calllib(ct, 'trans_new', model, ...
tp.tp_id, loglevel); tp.tp_id, loglevel);
end end
end end
tr.tp_id = tp.tp_id;
end end
%% Utility methods %% Utility methods
function clear(tr) function tr_clear(tr)
% Delete the kernel object. % Delete the kernel object.
checklib; checklib;
calllib(tr.lib, 'trans_del', tr.tr_id); calllib(ct, 'trans_del', tr.tr_id);
end end
%% Transport Methods %% Transport Methods
@ -51,11 +49,11 @@ classdef Transport < handle
function v = viscosity(tr) function v = viscosity(tr)
% Get the dynamic viscosity. % Get the dynamic viscosity.
% %
% :return: % return:
% Double dynamic viscosity. Unit: Pa*s. % Double dynamic viscosity. Unit: Pa*s.
checklib; checklib;
v = calllib(tr.lib, 'trans_viscosity', tr.tr_id); v = calllib(ct, 'trans_viscosity', tr.tr_id);
if v == -1.0 if v == -1.0
error(geterr); error(geterr);
elseif v < 0.0 elseif v < 0.0
@ -63,14 +61,14 @@ classdef Transport < handle
end end
end end
function v = thermaoConductivity(tr) function v = thermalConductivity(tr)
% Get the thermal conductivity. % Get the thermal conductivity.
% %
% :return: % return:
% Double thermal conductivity. Unit: W/m-K. % Double thermal conductivity. Unit: W/m-K.
checklib; checklib;
v = calllib(tr.lib, 'trans_thermalConductivity', tr.tr_id); v = calllib(ct, 'trans_thermalConductivity', tr.tr_id);
if v == -1.0 if v == -1.0
error(geterr); error(geterr);
elseif v < 0.0 elseif v < 0.0
@ -81,11 +79,11 @@ classdef Transport < handle
function v = electricalConductivity(tr) function v = electricalConductivity(tr)
% Get the electrical conductivity. % Get the electrical conductivity.
% %
% :return: % return:
% Double electrical conductivity. Unit: S/m. % Double electrical conductivity. Unit: S/m.
checklib; checklib;
v = calllib(tr.lib, 'trans_electricalConductivity', tr.tr_id); v = calllib(ct, 'trans_electricalConductivity', tr.tr_id);
if v == -1.0 if v == -1.0
error(geterr); error(geterr);
elseif v < 0.0 elseif v < 0.0
@ -96,7 +94,7 @@ classdef Transport < handle
function v = mixDiffCoeffs(tr) function v = mixDiffCoeffs(tr)
% Get the mixture-averaged diffusion coefficients. % Get the mixture-averaged diffusion coefficients.
% %
% :return: % return:
% Vector of length nSpecies with the mixture-averaged % Vector of length nSpecies with the mixture-averaged
% diffusion coefficients. Unit: m^2/s. % diffusion coefficients. Unit: m^2/s.
@ -104,14 +102,14 @@ classdef Transport < handle
nsp = tr.th.nSpecies; nsp = tr.th.nSpecies;
xx = zeros(1, nsp); xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(tr.lib, 'trans_getMixDiffCoeffs', tr.tr_id, nsp, pt); calllib(ct, 'trans_getMixDiffCoeffs', tr.tr_id, nsp, pt);
v = pt.Value; v = pt.Value;
end end
function v = thermalDiffCoeffs(tr) function v = thermalDiffCoeffs(tr)
% Get the thermal diffusion coefficients. % Get the thermal diffusion coefficients.
% %
% :return: % return:
% Vector of length nSpecies with the thermal diffusion % Vector of length nSpecies with the thermal diffusion
% coefficients. % coefficients.
@ -119,14 +117,14 @@ classdef Transport < handle
nsp = tr.th.nSpecies; nsp = tr.th.nSpecies;
xx = zeros(1, nsp); xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(tr.lib, 'trans_getThermalDiffCoeffs', tr.tr_id, nsp, pt); calllib(ct, 'trans_getThermalDiffCoeffs', tr.tr_id, nsp, pt);
v = pt.Value; v = pt.Value;
end end
function v = binDiffCoeffs(tr) function v = binDiffCoeffs(tr)
% Get the binary diffusion coefficients. % Get the binary diffusion coefficients.
% %
% :return: % return:
% A matrix of binary diffusion coefficients. The matrix is % A matrix of binary diffusion coefficients. The matrix is
% symmetric: d(i, j) = d(j, i). Unit: m^2/s. % symmetric: d(i, j) = d(j, i). Unit: m^2/s.
@ -134,14 +132,14 @@ classdef Transport < handle
nsp = tr.th.nSpecies; nsp = tr.th.nSpecies;
xx = zeros(1, nsp); xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(tr.lib, 'trans_getBinDiffCoeffs', tr.tr_id, nsp, pt); calllib(ct, 'trans_getBinDiffCoeffs', tr.tr_id, nsp, pt);
v = pt.Value; v = pt.Value;
end end
function v = multiDiffCoeffs(tr) function v = multiDiffCoeffs(tr)
% Get the multicomponent diffusion coefficients. % Get the multicomponent diffusion coefficients.
% %
% :return: % return:
% Vector of length nSpecies with the multicomponent % Vector of length nSpecies with the multicomponent
% diffusion coefficients. Unit: m^2/s. % diffusion coefficients. Unit: m^2/s.
@ -149,29 +147,29 @@ classdef Transport < handle
nsp = tr.th.nSpecies; nsp = tr.th.nSpecies;
xx = zeros(1, nsp); xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(tr.lib, 'trans_getMultiDiffCoeffs', tr.tr_id, nsp, pt); calllib(ct, 'trans_getMultiDiffCoeffs', tr.tr_id, nsp, pt);
v = pt.Value; v = pt.Value;
end end
function setParameters(tr, type, k, p) function setParameters(tr, type, k, p)
% Set the parameters. % Set the parameters.
% %
% :param type: % parameter type:
% :param k: % parameter k:
% :param p: % parameter p:
checklib; checklib;
calllib(tr.lib, 'trans_setParameters', tr.tr_id, type, k, p); calllib(ct, 'trans_setParameters', tr.tr_id, type, k, p);
end end
function setThermalConductivity(tr, lam) function setThermalConductivity(tr, lam)
% Set the thermal conductivity. % Set the thermal conductivity.
% %
% :param lam: % parameter lam:
% Thermal conductivity in W/(m-K). % Thermal conductivity in W/(m-K).
checklib; checklib;
setParameters(tr, 1, 0, lam); tr.setParameters(1, 0, lam);
end end
end end

3
interfaces/matlab_experimental/Base/Water.m
View File

@ -1,4 +1,5 @@
function w = Water() function w = Water()
% Return an object representing water. % Return an object representing water.
h = Solution('liquidvapor.yaml', 'water'); w = Solution('liquidvapor.yaml', 'water');
end end

4
interfaces/matlab_experimental/Constants/faradayconstant.m Normal file
View File

@ -0,0 +1,4 @@
function r = faradayconstant
% Get the Faraday constant in C/kmol of electron.
r = 96485336.4595687;
end

4
interfaces/matlab_experimental/Constants/gasconstant.m Normal file
View File

@ -0,0 +1,4 @@
function r = gasconstant
% Universal gas constant in J/kmol-K.
r = 8314.4621;
end

4
interfaces/matlab_experimental/Constants/oneatm.m Normal file
View File

@ -0,0 +1,4 @@
function p = oneatm
% Get 1 atm in Pa.
p = 101325.0;
end

30
interfaces/matlab_experimental/Func/Func.m
View File

@ -8,10 +8,6 @@ classdef Func < handle
typ typ
end end
properties(Constant = true)
lib = 'cantera_shared'
end
methods methods
%% Functor class constructor %% Functor class constructor
@ -44,7 +40,7 @@ classdef Func < handle
% class 'Func1'. See the Cantera C++ documentation for more % class 'Func1'. See the Cantera C++ documentation for more
% details. % details.
% %
% :param typ: % parameter typ:
% String indicating type of functor to create. Possible % String indicating type of functor to create. Possible
% values are: % values are:
% * 'polynomial' % * 'polynomial'
@ -57,11 +53,11 @@ classdef Func < handle
% * 'composite' % * 'composite'
% * 'periodic' % * 'periodic'
% %
% :param n: % parameter n:
% Number of parameters required for the functor % Number of parameters required for the functor
% :param p: % parameter p:
% Vector of parameters % Vector of parameters
% :return: % return:
% Instance of class :mat:func:`Func` % Instance of class :mat:func:`Func`
checklib; checklib;
@ -81,13 +77,13 @@ classdef Func < handle
ptr = libpointer('doublePtr', p); ptr = libpointer('doublePtr', p);
[m, n] = size(p); [m, n] = size(p);
lenp = m * n; lenp = m * n;
nn = calllib(x.lib, 'func_new', type, n, lenp, ptr); nn = calllib(ct, 'func_new', type, n, lenp, ptr);
elseif itype < 45 elseif itype < 45
m = n; m = n;
nn = calllib(x.lib, 'func_new', type, n, m, 0); nn = calllib(ct, 'func_new', type, n, m, 0);
else else
ptr = libpointer('doublePtr', p); ptr = libpointer('doublePtr', p);
nn = calllib(x.lib, 'func_new', type, n, 0, ptr); nn = calllib(ct, 'func_new', type, n, 0, ptr);
end end
end end
@ -124,13 +120,13 @@ classdef Func < handle
function clear(f) function clear(f)
% Clear the functor from memory. % Clear the functor from memory.
checklib; checklib;
calllib(f.lib, 'func_del', f.id); calllib(ct, 'func_del', f.id);
end end
function display(a) function display(a)
% Display the equation of the input function on the terminal. % Display the equation of the input function on the terminal.
% %
% :param a: % parameter a:
% Instance of class 'Func' % Instance of class 'Func'
disp(' '); disp(' ');
@ -162,11 +158,11 @@ classdef Func < handle
function b = subsref(a, s) function b = subsref(a, s)
% Redefine subscripted references for functors. % Redefine subscripted references for functors.
% %
% :param a: % parameter a:
% Instance of class 'Func' % Instance of class 'Func'
% :param s: % parameter s:
% Value at which the function should be evaluated. % Value at which the function should be evaluated.
% :return: % return:
% The value of the function evaluated at 's'. % The value of the function evaluated at 's'.
checklib; checklib;
@ -174,7 +170,7 @@ classdef Func < handle
ind= s.subs{:}; ind= s.subs{:};
b = zeros(1, length(ind)); b = zeros(1, length(ind));
for k = 1:length(ind) for k = 1:length(ind)
b(k) = calllib(a.lib, 'func_value', a.id, ind(k)); b(k) = calllib(ct, 'func_value', a.id, ind(k));
end end
else error('Specify value for x as p(x)'); else error('Specify value for x as p(x)');
end end

12
interfaces/matlab_experimental/Func/gaussian.m Normal file
View File

@ -0,0 +1,12 @@
function g = gaussian(peak, center, width)
% Create a Gaussian functor instance.
% :param peak:
% The peak value.
% :param center:
% Value of x at which the peak is located.
% :param width:
% Full width at half-maximum. The value of the function at center
% +/- (width)/2 is one-half the peak value.
g = Func('gaussian', 0, [peak, center, width]);
end

17
interfaces/matlab_experimental/Func/polynom.m Normal file
View File

@ -0,0 +1,17 @@
function poly = polynom(coeffs)
% Create a polynomial functor instance.
% The polynomial coefficients are specified by a vector [a0, a1, ...
% an]. Examples:
%
% polynom([-2, 6, 3]) % 3x^2+6x-2
% polynom([1.0, -2.5, 0, 0, 2]) % 2x^4-2.5x+1
[n m] = size(coeffs);
if n == 1
poly = Func('polynomial', m - 1, coeffs);
elseif m == 1
poly = Func('polynomial', n - 1, coeffs);
else
error('wrong shape for coefficient array');
end
end

20
interfaces/matlab_experimental/LoadCantera.m
View File

@ -1,8 +1,20 @@
function Load_Cantera function LoadCantera
addpath('Class', 'Utility', 'PresetObjects', 'Examples'); addpath('Class', 'Utility', 'PresetMixtures', 'PresetReactors', ...
if ~libisloaded('cantera_shared') 'PresetObjects', '1D', 'Examples');
if ispc
ctname = 'cantera_shared.dll';
elseif ismac
ctname = 'libcantera_shared.dylib';
elseif isunix
ctname = 'libcantera_shared.so';
else
error('Operating System Not Supported!');
return;
end
if ~libisloaded(ct)
load('Utility/cantera_root.mat'); load('Utility/cantera_root.mat');
[~,warnings] = loadlibrary([cantera_root '/lib/cantera_shared.dll'], ... [~,warnings] = loadlibrary([cantera_root '/lib/' ctname], ...
[cantera_root '/include/cantera/clib/ctmatlab.h'], ... [cantera_root '/include/cantera/clib/ctmatlab.h'], ...
'includepath', [cantera_root '/include'], ... 'includepath', [cantera_root '/include'], ...
'addheader','ct','addheader','ctfunc', ... 'addheader','ct','addheader','ctfunc', ...

14
interfaces/matlab_experimental/Reactor/ConstPressureReactor.m Normal file
View File

@ -0,0 +1,14 @@
function r = ConstPressureReactor(contents)
% Create a constant pressure reactor object.
% Pressure is held constant. The volume is not a state variable, but
% instead takes on whatever value is consistent with holding the
% pressure constant.
%
%: param contents:
% Contents of the reactor of class 'Solution'.
if nargin == 0
contents = 0;
end
r = Reactor(contents, 'ConstPressureReactor');
end

36
interfaces/matlab_experimental/Reactor/FlowDevice.m
View File

@ -7,17 +7,13 @@ classdef FlowDevice < handle
downstream downstream
end end
properties(Constant = true)
lib = 'cantera_shared'
end
methods methods
%% FlowDevice class constructor %% FlowDevice class constructor
function x = FlowDevice(typ) function x = FlowDevice(typ)
% Flow Device class constructor. % Flow Device class constructor.
% %
% :param typ: % parameter typ:
% Type of flow device to be created. Type = % Type of flow device to be created. Type =
% 'MassFlowController', 'PressureController' or 'Valve'. % 'MassFlowController', 'PressureController' or 'Valve'.
@ -36,7 +32,7 @@ classdef FlowDevice < handle
end end
x.type = typ; x.type = typ;
x.id = calllib(x.lib, 'flowdev_new2', typ); x.id = calllib(ct, 'flowdev_new', typ);
% if x.id < 0 % if x.id < 0
% error(geterr); % error(geterr);
% end % end
@ -49,7 +45,7 @@ classdef FlowDevice < handle
function clear(f) function clear(f)
% Clear the specified flow device from memory. % Clear the specified flow device from memory.
checklib; checklib;
calllib(f.lib, 'flowdev_del', f.id); calllib(ct, 'flowdev_del', f.id);
end end
%% FlowDevice methods %% FlowDevice methods
@ -57,9 +53,9 @@ classdef FlowDevice < handle
function install(f, upstream, downstream) function install(f, upstream, downstream)
% Install a flow device between reactors or reservoirs. % Install a flow device between reactors or reservoirs.
% %
% :param upstream: % parameter upstream:
% Upsteram 'Reactor' or 'Reservoir'. % Upsteram 'Reactor' or 'Reservoir'.
% :param downstream: % parameter downstream:
% Downstream 'Reactor' or 'Reservoir'. % Downstream 'Reactor' or 'Reservoir'.
checklib; checklib;
@ -70,7 +66,7 @@ classdef FlowDevice < handle
end end
i = upstream.id; i = upstream.id;
j = downstream.id; j = downstream.id;
ok = calllib(f.lib, 'flowdev_install', f.id, i, j); ok = calllib(ct, 'flowdev_install', f.id, i, j);
% if ok < 0 % if ok < 0
% error(geterr) % error(geterr)
% end % end
@ -81,31 +77,31 @@ classdef FlowDevice < handle
function mdot = massFlowRate(f, time) function mdot = massFlowRate(f, time)
% Get the mass flow rate at a given time. % Get the mass flow rate at a given time.
% %
% :param time: % parameter time:
% Time at which the mass flow rate is desired. % Time at which the mass flow rate is desired.
% :return: % return:
% The mass flow rate through the flow device at the given % The mass flow rate through the flow device at the given
% time. % time.
checklib; checklib;
if nargin == 1 if nargin == 1
mdot = calllib(f.lib, 'flowdev_massFlowRate2', f.id); mdot = calllib(ct, 'flowdev_massFlowRate2', f.id);
else else
warning(['"Time" argument to massFlowRate is deprecated', ... warning(['"Time" argument to massFlowRate is deprecated', ...
'and will be removed after Cantera 2.5.']); 'and will be removed after Cantera 2.5.']);
mdot = calllib(f.lib, 'flowdev_massFlowRate', f.id, time); mdot = calllib(ct, 'flowdev_massFlowRate', f.id, time);
end end
end end
function setFunction(f, mf) function setFunction(f, mf)
% Set the time function with class 'func'. % Set the time function with class 'func'.
% %
% :param mf: % parameter mf:
% Instance of class 'func'. % Instance of class 'func'.
checklib; checklib;
if strcmp(f.type, 'MassFlowController') if strcmp(f.type, 'MassFlowController')
k = calllib(f.lib, 'flowdev_setTimeFunction', f.id, ... k = calllib(ct, 'flowdev_setTimeFunction', f.id, ...
mf.id); mf.id);
% if k < 0 % if k < 0
% error(geterr); % error(geterr);
@ -118,12 +114,12 @@ classdef FlowDevice < handle
function setMassFlowRate(f, mdot) function setMassFlowRate(f, mdot)
% Set the mass flow rate to a constant value. % Set the mass flow rate to a constant value.
% %
% :param mdot: % parameter mdot:
% Mass flow rate % Mass flow rate
checklib; checklib;
if strcmp(f.type, 'MassFlowController') if strcmp(f.type, 'MassFlowController')
k = calllib(f.lib, 'flowdev_setMassFlowCoeff', f.id, mdot); k = calllib(ct, 'flowdev_setMassFlowCoeff', f.id, mdot);
% if k < 0 % if k < 0
% error(geterr); % error(geterr);
% end % end
@ -140,14 +136,14 @@ classdef FlowDevice < handle
% as long as this produces a positive value. If thsi expression % as long as this produces a positive value. If thsi expression
% is negative, zero is returned. % is negative, zero is returned.
% %
% :param k: % parameter k:
% Value fo the valve coefficient. Unit: kg/Pa-s. % Value fo the valve coefficient. Unit: kg/Pa-s.
checklib; checklib;
if ~strcmp(f.type, 'Valve') if ~strcmp(f.type, 'Valve')
error('Valve coefficient can only be set for valves.'); error('Valve coefficient can only be set for valves.');
end end
ok = calllib(f.lib, 'flowdev_setValveCoeff', f.id, k); ok = calllib(ct, 'flowdev_setValveCoeff', f.id, k);
% if k < 0 % if k < 0
% error(geterr); % error(geterr);
% end % end

12
interfaces/matlab_experimental/Reactor/FlowReactor.m Normal file
View File

@ -0,0 +1,12 @@
function r = FlowReactor(contents)
% Create a flow reactor object.
% A reactor representing adiabatic plug flow in a constant-area duct.
%
%: param contents:
% Contents of the reactor of class 'Solution'.
if nargin == 0
contents = 0;
end
r = Reactor(contents, 'FlowReactor');
end

14
interfaces/matlab_experimental/Reactor/IdealGasConstPressureReactor.m Normal file
View File

@ -0,0 +1,14 @@
function r = IdealGasConstPressureReactor(contents)
% Create a constant pressure reactor object with an ideal gas.
% Pressure is held constant. The volume is not a state variable, but
% instead takes on whatever value is consistent with holding the
% pressure constant.
%
%: param contents:
% Contents of the reactor of class 'Solution'.
if nargin == 0
contents = 0;
end
r = Reactor(contents, 'IdealGasConstPressureReactor');
end

13
interfaces/matlab_experimental/Reactor/IdealGasReactor.m Normal file
View File

@ -0,0 +1,13 @@
function r = IdealGasReactor(contents)
% Create a ideal gas reactor object.
% The governing equations are specialized for the ideal gas equation of
% state (and do not work correctly with other thermodynamic models).
%
%: param contents:
% Contents of the reactor of class 'Solution'.
if nargin == 0
contents = 0;
end
r = Reactor(contents, 'IdealGasReactor');
end

21
interfaces/matlab_experimental/Reactor/MassFlowController.m Normal file
View File

@ -0,0 +1,21 @@
function m = MassFlowController(upstream, downstream)
% Create a mass flow controller.
% Creates an instance of class 'FlowDevice' configured to simulate a
% mass flow controller that maintains a constant mass flow rate
% independent of upstream or downstream conditions. If two reactor
% objects are supplied as arguments, the controller is installed
% between the two reactors. Otherwise, the 'install' method should be
% used to install the mass flow controller between reactors.
%
% :param upstream:
% Upstream 'Reactor' or 'Reservoir'.
% :param downstream:
% Downstream 'Reactor' or 'Reservoir.
% :return:
% Instance of class 'FlowDevice'.
m = FlowDevice('MassFlowController');
if nargin == 2
m.install(upstream, downstream)
end
end

86
interfaces/matlab_experimental/Reactor/Reactor.m
View File

@ -17,10 +17,6 @@ classdef Reactor < handle
EnergyFlag EnergyFlag
end end
properties(Constant = true)
lib = 'cantera_shared'
end
methods methods
%% Reactor class constructor %% Reactor class constructor
@ -30,10 +26,10 @@ classdef Reactor < handle
% reactors through flow lines or through walls that may expand % reactors through flow lines or through walls that may expand
% or contract and/orconduct heat. % or contract and/orconduct heat.
% %
% :param contents: % parameter contents:
% Instance of class 'Solution' representing the contents of % Instance of class 'Solution' representing the contents of
% the reactor. % the reactor.
% :param typ: % parameter typ:
% Character array of reactor type. Options are: % Character array of reactor type. Options are:
% 'Reservoir' % 'Reservoir'
% 'Reactor' % 'Reactor'
@ -41,7 +37,7 @@ classdef Reactor < handle
% 'ConstPressureReactor' % 'ConstPressureReactor'
% 'IdealGasReactor' % 'IdealGasReactor'
% 'IdealGasConstPressureReactor' % 'IdealGasConstPressureReactor'
% :return: % return:
% Instance of class 'Reactor'. % Instance of class 'Reactor'.
checklib; checklib;
@ -56,7 +52,7 @@ classdef Reactor < handle
end end
r.type = char(typ); r.type = char(typ);
r.id = calllib(r.lib, 'reactor_new2', typ); r.id = calllib(ct, 'reactor_new', typ);
% if r.id < 0 % if r.id < 0
% error(geterr); % error(geterr);
@ -75,15 +71,15 @@ classdef Reactor < handle
function clear(r) function clear(r)
% Clear the reactor from memory. % Clear the reactor from memory.
checklib; checklib;
calllib(r.lib, 'reactor_del', r.id); calllib(ct, 'reactor_del', r.id);
end end
function insert(r, gas) function insert(r, gas)
% Insert a solution or mixture into a reactor. % Insert a solution or mixture into a reactor.
% %
% :param r: % parameter r:
% Instance of class 'Reactor'. % Instance of class 'Reactor'.
% :param gas: % parameter gas:
% Instance of class 'Solution'. % Instance of class 'Solution'.
r.contents = gas; r.contents = gas;
@ -99,9 +95,9 @@ classdef Reactor < handle
% This method is used internally during Reactor initialization, % This method is used internally during Reactor initialization,
% but is usually not called by users. % but is usually not called by users.
% %
% :param r: % parameter r:
% Instance of class 'Reactor'. % Instance of class 'Reactor'.
% :param t: % parameter t:
% Instance of class 'ThermoPhase' or another object % Instance of class 'ThermoPhase' or another object
% containing an instance of that class. % containing an instance of that class.
checklib; checklib;
@ -110,7 +106,7 @@ classdef Reactor < handle
error('Wrong object type'); error('Wrong object type');
end end
calllib(r.lib, 'reactor_setThermoMgr', r.id, t.tp_id); calllib(ct, 'reactor_setThermoMgr', r.id, t.tp_id);
end end
function setKineticsMgr(r, k) function setKineticsMgr(r, k)
@ -119,9 +115,9 @@ classdef Reactor < handle
% This method is used internally during Reactor initialization, % This method is used internally during Reactor initialization,
% but is usually not called by users. % but is usually not called by users.
% %
% :param r: % parameter r:
% Instance of class 'Reactor'. % Instance of class 'Reactor'.
% :param t: % parameter t:
% Instance of class 'Kinetics' or another object % Instance of class 'Kinetics' or another object
% containing an instance of that class. % containing an instance of that class.
checklib; checklib;
@ -130,7 +126,7 @@ classdef Reactor < handle
error('Wrong object type'); error('Wrong object type');
end end
calllib(r.lib, 'reactor_setKineticsMgr', r.id, k.kin_id); calllib(ct, 'reactor_setKineticsMgr', r.id, k.kin_id);
end end
%% Reactor get methods %% Reactor get methods
@ -138,103 +134,103 @@ classdef Reactor < handle
function temperature = get.T(r) function temperature = get.T(r)
% Get the temperature of the reactor. % Get the temperature of the reactor.
% %
% :return: % return:
% The temperature of the reactor contents at the end of the % The temperature of the reactor contents at the end of the
% last call to 'advance' or 'step'. Unit: K. % last call to 'advance' or 'step'. Unit: K.
checklib; checklib;
temperature = calllib(r.lib, 'reactor_temperature', r.id); temperature = calllib(ct, 'reactor_temperature', r.id);
end end
function pressure = get.P(r) function pressure = get.P(r)
% Get the pressure of the reactor. % Get the pressure of the reactor.
% %
% :return: % return:
% The pressure of the reactor contents at the end of the % The pressure of the reactor contents at the end of the
% last call to 'advance' or 'step'. Unit: Pa. % last call to 'advance' or 'step'. Unit: Pa.
checklib; checklib;
pressure = calllib(r.lib, 'reactor_pressure', r.id); pressure = calllib(ct, 'reactor_pressure', r.id);
end end
function rho = get.D(r) function rho = get.D(r)
% Get the density of the reactor. % Get the density of the reactor.
% %
% :return: % return:
% Density of the phase in the input. Unit: kg/m^3. % Density of the phase in the input. Unit: kg/m^3.
checklib; checklib;
rho = calllib(r.lib, 'reactor_density', r.id); rho = calllib(ct, 'reactor_density', r.id);
end end
function mass = get.M(r) function mass = get.M(r)
% Get the mass of the reactor. % Get the mass of the reactor.
% %
% :return: % return:
% The mass of the reactor contents at the end of the % The mass of the reactor contents at the end of the
% last call to 'advance' or 'step'. The mass is retrieved % last call to 'advance' or 'step'. The mass is retrieved
% from the solution vector. Unit: kg. % from the solution vector. Unit: kg.
checklib; checklib;
mass = calllib(r.lib, 'reactor_mass', r.id); mass = calllib(ct, 'reactor_mass', r.id);
end end
function volume = get.V(r) function volume = get.V(r)
% Get the volume of the reactor. % Get the volume of the reactor.
% %
% :return: % return:
% The volume of the reactor contents at the end of the % The volume of the reactor contents at the end of the
% last call to 'advance' or 'step'. Unit: m^3. % last call to 'advance' or 'step'. Unit: m^3.
checklib; checklib;
volume = calllib(r.lib, 'reactor_volume', r.id); volume = calllib(ct, 'reactor_volume', r.id);
end end
function enthalpy_mass = get.H(r) function enthalpy_mass = get.H(r)
% Get the mass specific enthalpy of the reactor. % Get the mass specific enthalpy of the reactor.
% %
% :return: % return:
% The mass specific enthalpy of the reactor contents at the % The mass specific enthalpy of the reactor contents at the
% end of the last call to 'advance' or 'step'. The enthalpy % end of the last call to 'advance' or 'step'. The enthalpy
% is retrieved from the solution vector. Unit: J/kg. % is retrieved from the solution vector. Unit: J/kg.
checklib; checklib;
enthalpy_mass = calllib(r.lib, 'reactor_enthalpy_mass', r.id); enthalpy_mass = calllib(ct, 'reactor_enthalpy_mass', r.id);
end end
function intEnergy_mass = get.U(r) function intEnergy_mass = get.U(r)
% Get the mass specific internal energy of the reactor. % Get the mass specific internal energy of the reactor.
% %
% :return: % return:
% The mass specific internal energy of the reactor contents % The mass specific internal energy of the reactor contents
% at the end of the last call to 'advance' or 'step'. The % at the end of the last call to 'advance' or 'step'. The
% internal energy is retrieved from the solution vector. % internal energy is retrieved from the solution vector.
% Unit: J/kg. % Unit: J/kg.
checklib; checklib;
intEnergy_mass = calllib(r.lib, 'reactor_intEnergy_mass', r.id); intEnergy_mass = calllib(ct, 'reactor_intEnergy_mass', r.id);
end end
function yi = massFraction(r, species) function yi = massFraction(r, species)
% Get the mass fraction of a species. % Get the mass fraction of a species.
% %
% :param species: % parameter species:
% String or one-based integer id of the species. % String or one-based integer id of the species.
checklib; checklib;
if ischar(species) if ischar(species)
k = r.contents.thermo.speciesIndex(species) - 1; k = r.contents.speciesIndex(species) - 1;
else k = species - 1; else k = species - 1;
end end
yi = calllib(r.lib, 'reactor_massFraction', r.id, k); yi = calllib(ct, 'reactor_massFraction', r.id, k);
end end
function massFractions = get.Y(r) function massFractions = get.Y(r)
% Get the mass fractions of the reactor. % Get the mass fractions of the reactor.
% %
% :return: % return:
% The mass fractions of the reactor contents at the end of % The mass fractions of the reactor contents at the end of
% the last call to 'advance' or 'step'. % the last call to 'advance' or 'step'.
@ -252,21 +248,21 @@ classdef Reactor < handle
function setInitialVolume(r, v0) function setInitialVolume(r, v0)
% Set the initial reactor volume. % Set the initial reactor volume.
% %
% :param v0: % parameter v0:
% Initial volume in m^3. % Initial volume in m^3.
checklib; checklib;
calllib(r.lib, 'reactor_setInitialVolume', r.id, v0); calllib(ct, 'reactor_setInitialVolume', r.id, v0);
end end
function r = set.Mdot(r, MFR) function r = set.Mdot(r, MFR)
% Set the mass flow rate. % Set the mass flow rate.
% %
% :param MFR: % parameter MFR:
% Mass flow rate in kg/s. % Mass flow rate in kg/s.
checklib; checklib;
calllib(r.lib, 'reactor_setMassFlowRate', r.id, MFR); calllib(ct, 'reactor_setMassFlowRate', r.id, MFR);
r.Mdot = MFR; r.Mdot = MFR;
end end
@ -281,9 +277,9 @@ classdef Reactor < handle
% equations enabled if there are reactions present in the % equations enabled if there are reactions present in the
% mechanism file, and disabled otherwise. % mechanism file, and disabled otherwise.
% %
% :param r: % parameter r:
% Instance of class 'Reactor'. % Instance of class 'Reactor'.
% :param flag: % parameter flag:
% String, either "on" or "off" to enable or disable chemical % String, either "on" or "off" to enable or disable chemical
% reactions, respectively. % reactions, respectively.
@ -296,7 +292,7 @@ classdef Reactor < handle
else error('Input must be "on" or "off"'); else error('Input must be "on" or "off"');
end end
calllib(r.lib, 'reactor_setChemistry', r.id, iflag); calllib(ct, 'reactor_setChemistry', r.id, iflag);
r.ChemistryFlag = flag; r.ChemistryFlag = flag;
end end
@ -310,9 +306,9 @@ classdef Reactor < handle
% By default, Reactor objects are created with the energy % By default, Reactor objects are created with the energy
% equation enabled, so usually this method % equation enabled, so usually this method
% %
% :param r: % parameter r:
% Instance of class 'Reactor'. % Instance of class 'Reactor'.
% :param flag: % parameter flag:
% String, either "on" or "off" to enable or disable chemical % String, either "on" or "off" to enable or disable chemical
% reactions, respectively. % reactions, respectively.
@ -326,7 +322,7 @@ classdef Reactor < handle
end end
if iflag >= 0 if iflag >= 0
calllib(r.lib, 'reactor_setEnergy', r.id, iflag); calllib(ct, 'reactor_setEnergy', r.id, iflag);
else error('Input must be "on" or "off".'); else error('Input must be "on" or "off".');
end end

42
interfaces/matlab_experimental/Reactor/ReactorNet.m
View File

@ -8,10 +8,6 @@ classdef ReactorNet < handle
rtol rtol
end end
properties(Constant = true)
lib = 'cantera_shared'
end
methods methods
%% ReactorNet class constructor %% ReactorNet class constructor
@ -21,10 +17,10 @@ classdef ReactorNet < handle
% to simultaneously advance the state of one or more coupled % to simultaneously advance the state of one or more coupled
% reactors. % reactors.
% %
% :param reactors: % parameter reactors:
% Instance of class 'Reactor' or a cell array of instance of % Instance of class 'Reactor' or a cell array of instance of
% 'Reactor'. % 'Reactor'.
% :return: % return:
% Instance of class 'ReactorNet'. % Instance of class 'ReactorNet'.
checklib; checklib;
@ -39,7 +35,7 @@ classdef ReactorNet < handle
reactors = {reactor}; reactors = {reactor};
end end
r.id = calllib(r.lib, 'reactornet_new'); r.id = calllib(ct, 'reactornet_new');
% if r.id < 0 % if r.id < 0
% error(geterr); % error(geterr);
% end % end
@ -56,19 +52,19 @@ classdef ReactorNet < handle
function clear(r) function clear(r)
% Clear the ReactorNet object from the memory. % Clear the ReactorNet object from the memory.
checklib; checklib;
calllib(r.lib, 'reactornet_del', r.id); calllib(ct, 'reactornet_del', r.id);
end end
function addReactor(net, reactor) function addReactor(net, reactor)
% Add a reactor to a network. % Add a reactor to a network.
% %
% :param net: % parameter net:
% Instance of class 'ReactorNet'. % Instance of class 'ReactorNet'.
% :param reactor: % parameter reactor:
% Instance of class 'Solution'. % Instance of class 'Solution'.
checklib; checklib;
calllib(net.lib, 'reactornet_addreactor', net.id, reactor.id); calllib(ct, 'reactornet_addreactor', net.id, reactor.id);
end end
function advance(r, tout) function advance(r, tout)
@ -81,11 +77,11 @@ classdef ReactorNet < handle
% an absolute time, not a time interval.) The integrator may % an absolute time, not a time interval.) The integrator may
% take many internal timesteps before reaching tout. % take many internal timesteps before reaching tout.
% %
% :param tout: % parameter tout:
% End time of the integration. Unit: s. % End time of the integration. Unit: s.
checklib; checklib;
calllib(r.lib, 'reactornet_advance', r.id, tout); calllib(ct, 'reactornet_advance', r.id, tout);
end end
%% ReactorNet set methods %% ReactorNet set methods
@ -93,12 +89,12 @@ classdef ReactorNet < handle
function setInitialTime(r, t) function setInitialTime(r, t)
% Set the initial time of the integration. % Set the initial time of the integration.
% %
% :param t: % parameter t:
% Time at which integration should be restarted, using the % Time at which integration should be restarted, using the
% current state as the initial condition. Unit: s. % current state as the initial condition. Unit: s.
checklib; checklib;
calllib(r.lib, 'reactornet_setInitialTime', r.id, t); calllib(ct, 'reactornet_setInitialTime', r.id, t);
end end
function setMaxTimeStep(r, maxstep) function setMaxTimeStep(r, maxstep)
@ -113,19 +109,19 @@ classdef ReactorNet < handle
% upper bound on the timestep. % upper bound on the timestep.
checklib; checklib;
calllib(r.lib, 'reactornet_setMaxTimeStep', r.id, maxstep); calllib(ct, 'reactornet_setMaxTimeStep', r.id, maxstep);
end end
function setTolerances(r, rerr, aerr) function setTolerances(r, rerr, aerr)
% Set the error tolerance. % Set the error tolerance.
% %
% :param rtol: % parameter rtol:
% Scalar relative error tolerance. % Scalar relative error tolerance.
% :param atol: % parameter atol:
% Scalar absolute error tolerance. % Scalar absolute error tolerance.
checklib; checklib;
calllib(r.lib, 'reactornet_setTolerances', r.id, rerr, aerr); calllib(ct, 'reactornet_setTolerances', r.id, rerr, aerr);
end end
%% ReactorNet get methods %% ReactorNet get methods
@ -139,25 +135,25 @@ classdef ReactorNet < handle
% rapidly changing, the time step becomes smaller to resolve % rapidly changing, the time step becomes smaller to resolve
% the solution. % the solution.
checklib; checklib;
t = calllib(r.lib, 'reactor_step', r.id); t = calllib(ct, 'reactor_step', r.id);
end end
function t = get.time(r) function t = get.time(r)
% Get the current time in s. % Get the current time in s.
checklib; checklib;
t = calllib(r.lib, 'reactornet_time', r.id); t = calllib(ct, 'reactornet_time', r.id);
end end
function t = get.rtol(r) function t = get.rtol(r)
% Get the relative error tolerance % Get the relative error tolerance
checklib; checklib;
t = calllib(r.lib, 'reactornet_rtol', r.id); t = calllib(ct, 'reactornet_rtol', r.id);
end end
function t = get.atol(r) function t = get.atol(r)
% Get the absolute error tolerance % Get the absolute error tolerance
checklib; checklib;
t = calllib(r.lib, 'reactornet_atol', r.id); t = calllib(ct, 'reactornet_atol', r.id);
end end
end end

30
interfaces/matlab_experimental/Reactor/ReactorSurface.m
View File

@ -6,10 +6,6 @@ classdef ReactorSurface < handle
reactor reactor
end end
properties(Constant = true)
lib = 'cantera_shared'
end
methods methods
%% ReactorSurface class constructor %% ReactorSurface class constructor
@ -25,29 +21,29 @@ classdef ReactorSurface < handle
% after initial construction by using the various methods of % after initial construction by using the various methods of
% the 'ReactorSurface' class. % the 'ReactorSurface' class.
% %
% :param kleft: % parameter kleft:
% Surface reaction mechanisms for the left-facing surface. % Surface reaction mechanisms for the left-facing surface.
% This must bean instance of class 'Kinetics', or of a class % This must bean instance of class 'Kinetics', or of a class
% derived from 'Kinetics', such as 'Interface'. % derived from 'Kinetics', such as 'Interface'.
% :param reactor: % parameter reactor:
% Instance of class 'Reactor' to be used as the adjacent % Instance of class 'Reactor' to be used as the adjacent
% bulk phase. % bulk phase.
% :param area: % parameter area:
% The area of the surface in m^2. Defaults to 1.0 m^2 if not % The area of the surface in m^2. Defaults to 1.0 m^2 if not
% specified. % specified.
% :return: % return:
% Instance of class 'ReactorSurface'. % Instance of class 'ReactorSurface'.
checklib; checklib;
s.id = calllib(s.lib, 'reactorsurface_new', 0); s.id = calllib(ct, 'reactorsurface_new', 0);
s.reactor = -1; s.reactor = -1;
% if r.id < 0 % if r.id < 0
% error(geterr); % error(geterr);
% end % end
if nargin >= 1 if nargin >= 1
s.setKinetics(s, kleft); s.setKinetics(kleft);
end end
if nargin >= 2 if nargin >= 2
@ -60,7 +56,7 @@ classdef ReactorSurface < handle
if nargin >= 3 if nargin >= 3
if isnumeric(area) if isnumeric(area)
s.setArea(area); s.area = area;
else else
warning('Area was not a number and was not set'); warning('Area was not a number and was not set');
end end
@ -73,14 +69,14 @@ classdef ReactorSurface < handle
function clear(s) function clear(s)
% Clear the ReactorSurface object from the memory. % Clear the ReactorSurface object from the memory.
checklib; checklib;
calllib(s.lib, 'reactorsurface_del', s.id); calllib(ct, 'reactorsurface_del', s.id);
end end
function install(s, r) function install(s, r)
% Install a ReactorSurface in a Reactor. % Install a ReactorSurface in a Reactor.
checklib; checklib;
s.reactor = r; s.reactor = r;
calllib(s.lib, 'reactorsurface_install', s.id, r.id); calllib(ct, 'reactorsurface_install', s.id, r.id);
end end
%% ReactorSurface get methods %% ReactorSurface get methods
@ -88,7 +84,7 @@ classdef ReactorSurface < handle
function a = get.area(s) function a = get.area(s)
% Get the areaof the reactor surface in m^2. % Get the areaof the reactor surface in m^2.
checklib; checklib;
a = calllib(s.lib, 'reactorsurface_area', s.id); a = calllib(ct, 'reactorsurface_area', s.id);
end end
%% ReactorSurface set methods %% ReactorSurface set methods
@ -96,12 +92,12 @@ classdef ReactorSurface < handle
function s = set.area(s, a) function s = set.area(s, a)
% Set the area of a reactor surface % Set the area of a reactor surface
checklib; checklib;
calllib(s.lib, 'reactorsurface_setArea', s.id, a); calllib(ct, 'reactorsurface_setArea', s.id, a);
end end
function setKinetics(s, kin) function setKinetics(s, kin)
% Setthe surface reaction mechanism on a reactor surface. % Setthe surface reaction mechanism on a reactor surface.
% :param kin: % parameter kin:
% Instance of class 'Kinetics' (or another object derived % Instance of class 'Kinetics' (or another object derived
% from kin) to be used as the kinetic mechanism for this % from kin) to be used as the kinetic mechanism for this
% surface. Typically an instance of class 'Interface'. % surface. Typically an instance of class 'Interface'.
@ -113,7 +109,7 @@ classdef ReactorSurface < handle
ikin = kin.id; ikin = kin.id;
end end
calllib(s.lib, 'reactorsurface_setkinetics', s.id, ikin); calllib(ct, 'reactorsurface_setkinetics', s.id, ikin);
end end
end end

18
interfaces/matlab_experimental/Reactor/Reservoir.m Normal file
View File

@ -0,0 +1,18 @@
function r = Reservoir(contents)
% Create a reservoir object.
% An instance of class 'Reactor' configured so that its intensive state
% is constant in time. A reservoir may be thought of as inifinite in
% extent, perfectly mixed, and non-reacting, so that fluid may be
% extracted or added without changing the composition or thermodynamic
% state. Note that even if the reaction mechanism associated with the
% fluid in the reactor defines reactions, they are disabled within
% reservoirs.
%
%: param contents:
% Contents of the reactor of class 'Solution'.
if nargin == 0
contents = 0;
end
r = Reactor(contents, 'Reservoir');
end

32
interfaces/matlab_experimental/Reactor/Valve.m Normal file
View File

@ -0,0 +1,32 @@
function v = Valve(upstream, downstream)
% Create a valve.
% Creates an instance of class 'FlowDevice' configured to simulate a
% valve that produces a flow rate proportional to the pressure
% difference between the upstream or downstream reactors.
%
% The mass flow rate [kg/s] is computed from the expression:
%
% mdot = K(P_upstream - P_downstream)
%
% as long as this produces a positive value. If this expression is
% negative, zero is returned. Therefore, the 'Valve' object acts as a
% check valve - flow is always from upstream to downstream.
%
% Note: as currently implemented, the valve object does not model real
% valve characteristics - inparticular, it does not model choked flow.
% The mass flow rate is always assumed to be linearly proportional to
% the pressure difference, no matter how large. This MAY change in a
% future release.
%
% :param upstream:
% Upstream 'Reactor' or 'Reservoir'.
% :param downstream:
% Downstream 'Reactor' or 'Reservoir.
% :return:
% Instance of class 'FlowDevice'.
v = FlowDevice('Valve');
if nargin == 2
v.install(upstream, downstream)
end
end

56
interfaces/matlab_experimental/Reactor/Wall.m
View File

@ -8,10 +8,6 @@ classdef Wall < handle
area area
end end
properties(Constant = true)
lib = 'cantera_shared'
end
methods methods
%% Wall class constructor %% Wall class constructor
@ -50,25 +46,25 @@ classdef Wall < handle
% wall can be set by using empty strings or 0.0 for each of the % wall can be set by using empty strings or 0.0 for each of the
% arguments before the velocity with no harm. % arguments before the velocity with no harm.
% %
% :param left: % parameter left:
% Instance of class 'Reactor' to be used as the bulk phase % Instance of class 'Reactor' to be used as the bulk phase
% on the left side of the wall. % on the left side of the wall.
% :param right: % parameter right:
% Instance of class 'Reactor' to be used as the bulk phase % Instance of class 'Reactor' to be used as the bulk phase
% on the right side of the wall. % on the right side of the wall.
% :param area: % parameter area:
% The area of the wall in m^2. Defaults to 1.0 m^2 if not % The area of the wall in m^2. Defaults to 1.0 m^2 if not
% specified. % specified.
% :param k: % parameter k:
% Expansion rate coefficient in m/(s-Pa). Defaults to 0.0 if % Expansion rate coefficient in m/(s-Pa). Defaults to 0.0 if
% not specified. % not specified.
% :param u: % parameter u:
% Heat transfer coefficient in W/(m^2-K). Defaults to 0.0 if % Heat transfer coefficient in W/(m^2-K). Defaults to 0.0 if
% not specified. % not specified.
% :param q: % parameter q:
% Heat flux in W/m^2. Defaults to 0.0 if not specified. Must % Heat flux in W/m^2. Defaults to 0.0 if not specified. Must
% be an instance of 'Func'. % be an instance of 'Func'.
% :param v: % parameter v:
% Velocity of the wall in m/s. Defaults to 0.0 if not % Velocity of the wall in m/s. Defaults to 0.0 if not
% specified. Must be an instance of 'Func'. % specified. Must be an instance of 'Func'.
@ -78,7 +74,7 @@ classdef Wall < handle
typ = 'Wall'; typ = 'Wall';
x.type = char(typ); x.type = char(typ);
x.id = calllib(x.lib, 'wall_new2', x.type); x.id = calllib(ct, 'wall_new', x.type);
% if x.index < 0 % if x.index < 0
% error(geterr); % error(geterr);
% end % end
@ -142,7 +138,7 @@ classdef Wall < handle
function clear(w) function clear(w)
% Clear the Wall object from the memory. % Clear the Wall object from the memory.
checklib; checklib;
calllib(w.lib, 'wall_del', w.id); calllib(ct, 'wall_del', w.id);
end end
function install(w, l, r) function install(w, l, r)
@ -151,51 +147,51 @@ classdef Wall < handle
checklib; checklib;
w.left = l; w.left = l;
w.right = r; w.right = r;
calllib(w.lib, 'wall_install', w.id, l.id, r.id); calllib(ct, 'wall_install', w.id, l.id, r.id);
end end
function ok = ready(w) function ok = ready(w)
% Check whether a wall is ready. % Check whether a wall is ready.
checklib; checklib;
ok = calllib(w.lib, 'wall_ready', w.id); ok = calllib(ct, 'wall_ready', w.id);
end end
%% ReactorNet set methods %% ReactorNet set methods
function w = set.area(w, a) function set.area(w, a)
% Set the area of a wall. % Set the area of a wall.
checklib; checklib;
calllib(w.lib, 'wall_Area', w.id, a); calllib(ct, 'wall_setArea', w.id, a);
end end
function setThermalResistance(w, r) function setThermalResistance(w, r)
% Set the thermal resistance. % Set the thermal resistance.
% %
% :param r: % parameter r:
% Thermal resistance. Unit: K*m^2/W. % Thermal resistance. Unit: K*m^2/W.
checklib; checklib;
calllib(w.lib, 'wall_setThermalResistance', w.id, r); calllib(ct, 'wall_setThermalResistance', w.id, r);
end end
function setHeatTransferCoeff(w, u) function setHeatTransferCoeff(w, u)
% Set the thermal transfer coefficient. % Set the thermal transfer coefficient.
% %
% :param u: % parameter u:
% Heat transfer coefficient. Unit: W/(m^2-K). % Heat transfer coefficient. Unit: W/(m^2-K).
checklib; checklib;
calllib(w.lib, 'wall_setHeatTransferCoeff', w.id, u); calllib(ct, 'wall_setHeatTransferCoeff', w.id, u);
end end
function setExpansionRateCoeff(w, k) function setExpansionRateCoeff(w, k)
% Set the expansion rate coefficient. % Set the expansion rate coefficient.
% %
% :param k: % parameter k:
% Expanstion rate coefficient. Unit: m/(s-Pa). % Expanstion rate coefficient. Unit: m/(s-Pa).
checklib; checklib;
calllib(w.lib, 'wall_setExpansionRateCoeff', w.id, k); calllib(ct, 'wall_setExpansionRateCoeff', w.id, k);
end end
function setHeatFlux(w, f) function setHeatFlux(w, f)
@ -206,11 +202,11 @@ classdef Wall < handle
% to specify a constant heat flux by using the polynomial % to specify a constant heat flux by using the polynomial
% functor with only the first term specified. % functor with only the first term specified.
% %
% :param f: % parameter f:
% Instance of class 'Func'. Unit: W/m^2. % Instance of class 'Func'. Unit: W/m^2.
checklib; checklib;
calllib(w.lib, 'wall_setHeatFlux', w.id, f.id); calllib(ct, 'wall_setHeatFlux', w.id, f.id);
end end
function setVelocity(w, f) function setVelocity(w, f)
@ -221,11 +217,11 @@ classdef Wall < handle
% to specify a constant velocity by using the polynomial % to specify a constant velocity by using the polynomial
% functor with only the first term specified. % functor with only the first term specified.
% %
% :param f: % parameter f:
% Instance of class 'Func'. Unit: m/s. % Instance of class 'Func'. Unit: m/s.
checklib; checklib;
calllib(w.lib, 'wall_setVelocity', w.id, f.id); calllib(ct, 'wall_setVelocity', w.id, f.id);
end end
%% ReactorNet get methods %% ReactorNet get methods
@ -233,19 +229,19 @@ classdef Wall < handle
function a = get.area(w) function a = get.area(w)
% Get the area of the wall in m^2. % Get the area of the wall in m^2.
checklib; checklib;
a = calllib(w.lib, 'wall_area', w.id); a = calllib(ct, 'wall_area', w.id);
end end
function q = qdot(w, t) function q = qdot(w, t)
% Get the total heat transfer through a wall at given time. % Get the total heat transfer through a wall at given time.
checklib; checklib;
q = calllib(w.lib, 'wall_Q', w.id, t); q = calllib(ct, 'wall_Q', w.id, t);
end end
function v = vdot(w, t) function v = vdot(w, t)
% Get the rate of volumetric change at a given time. % Get the rate of volumetric change at a given time.
checklib; checklib;
v = calllib(w.lib, 'wall_vdot', w.id, t); v = calllib(ct, 'wall_vdot', w.id, t);
end end
end end

4
interfaces/matlab_experimental/UnloadCantera.m Normal file
View File

@ -0,0 +1,4 @@
if libisloaded(ct)
unloadlibrary(ct);
end
disp('Cantera has been unloaded');

4
interfaces/matlab_experimental/Utility/UnloadCantera.m
View File

@ -1,4 +0,0 @@
if libisloaded('cantera_shared')
unloadlibrary('cantera_shared');
end
disp('Cantera has been unloaded');

8
interfaces/matlab_experimental/Utility/canteraGitCommit.m Normal file
View File

@ -0,0 +1,8 @@
function v = canteraGitCommit()
% Get Cantera Git commit hash
checklib;
buflen = calllib(ct, 'ct_getGitCommit', 0, '');
aa = char(zeros(1, buflen));
[~, aa] = calllib(ct, 'ct_getGitCommit', buflen, aa);
v = aa;
end

8
interfaces/matlab_experimental/Utility/canteraVersion.m Normal file
View File

@ -0,0 +1,8 @@
function v = canteraVersion()
% Get Cantera version information
checklib;
buflen = calllib(ct, 'ct_getCanteraVersion', 0, '');
aa = char(zeros(1, buflen));
[~, aa] = calllib(ct, 'ct_getCanteraVersion', buflen, aa);
v = aa;
end

2
interfaces/matlab_experimental/Utility/checklib.m
View File

@ -1,3 +1,3 @@
if ~libisloaded('cantera_shared') if ~libisloaded(ct)
error('Cantera is not loaded'); error('Cantera is not loaded');
end end

13
interfaces/matlab_experimental/Utility/cleanup.m
View File

@ -1,12 +1,11 @@
function cleanup() function cleanup()
% Delete all stored Cantera objects and reclaim memory % Delete all stored Cantera objects and reclaim memory
checklib; checklib;
calllib('cantera_shared', 'ct_clearOneDim'); calllib(ct, 'ct_clearOneDim');
calllib('cantera_shared', 'ct_clearMix'); calllib(ct, 'ct_clearMix');
calllib('cantera_shared', 'ct_clearXML'); calllib(ct, 'ct_clearFunc');
calllib('cantera_shared', 'ct_clearFunc'); calllib(ct, 'ct_clearStorage');
calllib('cantera_shared', 'ct_clearStorage'); calllib(ct, 'ct_clearReactors');
calllib('cantera_shared', 'ct_clearReactors'); calllib(ct, 'ct_clearReactionPath');
calllib('cantera_shared', 'ct_clearReactionPath');
clear all clear all
end end

8
interfaces/matlab_experimental/Utility/getDataDirectories.m Normal file
View File

@ -0,0 +1,8 @@
function d = getDataDirectories()
% Get a cell array of the directories searched for data files.
checklib;
buflen = calllib(ct, 'ct_getDataDirectories', 0, '', ';');
aa = char(zeros(1, buflen));
[~, aa, ~] = calllib(ct, 'ct_getDataDirectories', buflen, aa, ';');
d = aa;
end

12
interfaces/matlab_experimental/Utility/geterr.m Normal file
View File

@ -0,0 +1,12 @@
function e = geterr()
% Get the error message from Cantera error.
checklib;
try
buflen = calllib(ct, 'ct_getCanteraError', 0, '');
aa = char(zeros(1, buflen));
[~, aa] = calllib(ct, 'ct_getCanteraError', buflen, aa);
e = aa;
catch ME
e = getReport(ME);
end
end

15
interfaces/matlab_experimental/Utility/importEdge.m Normal file
View File

@ -0,0 +1,15 @@
function s = importEdge(file, name, phase1, phase2, phase3, phase4)
% Import edges between phases.
if nargin == 3
s = Interface(file, name, phase1);
elseif nargin == 4
s = Interface(file, name, phase1, phase2);
elseif nargin == 5
s = Interface(file, name, phase1, phase2, phase3);
elseif nargin == 6
s = Interface(file, name, phase1, phase2, phase3, phase4);
else
error('importEdge only supports 4 neighbor phases.');
end
end

11
interfaces/matlab_experimental/Utility/importInterface.m Normal file
View File

@ -0,0 +1,11 @@
function s = importInterface(file, name, phase1, phase2)
% Import an interface between phases.
if nargin == 3
s = Interface(file, name, phase1);
elseif nargin == 4
s = Interface(file, name, phase1, phase2);
else
error('importInterface only supports 2 bulk phases');
end
end

12
interfaces/matlab_experimental/cleanup.m
View File

@ -1,12 +0,0 @@
function cleanup()
% Delete all stored Cantera objects and reclaim memory
checklib;
calllib('cantera_shared', 'ct_clearOneDim');
calllib('cantera_shared', 'ct_clearMix');
calllib('cantera_shared', 'ct_clearXML');
calllib('cantera_shared', 'ct_clearFunc');
calllib('cantera_shared', 'ct_clearStorage');
calllib('cantera_shared', 'ct_clearReactors');
calllib('cantera_shared', 'ct_clearReactionPath');
clear all
end

8
interfaces/matlab_experimental/ct.m Normal file
View File

@ -0,0 +1,8 @@
function str = ct()
% Return name of 'cantera_shared' library.
if ispc
str = 'cantera_shared';
else
str = 'libcantera_shared';
end
end

9
interfaces/matlab_experimental/readme.txt Normal file
View File

@ -0,0 +1,9 @@
This Matlab Toolbox for Cantera serves as a proof of concept for direct C-library integration into Matlab.
Move the toolbox_C folder into <installation directory to Cantera>\matlab.
For first time users, launch matlab, navigate to the toolbox folder, then run the 'SetCanteraPath(<installation directory to Cantera>)' command.
To start using the toolbox, use the LoadCantera command to load cantera_shared.dll into the memory.
To stop using Cantera and purge all objects from memory, first run cleanup command, then run UnloadCantera command.