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

Made changes based on Feedback.

Browse Source
This commit is contained in:
ssun30 2022-03-22 13:51:03 -04:00 committed by Ray Speth
parent 0294939a94
commit 21980e6e6f
18 changed files with 732 additions and 1028 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
interfaces/matlab_experimental/1D/AxiStagnFlow.m
View File

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

302
interfaces/matlab_experimental/1D/Domain1D.m
View File

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

202
interfaces/matlab_experimental/1D/Sim1D.m
View File

@ -1,7 +1,7 @@
classdef Stack < handle classdef Stack < handle
properties properties
st_id stID
domains domains
end end
@ -13,27 +13,27 @@ classdef Stack < handle
% which are instances of class Domain1D. The domains are of two % which are instances of class Domain1D. The domains are of two
% types - extended domains, and connector domains. % types - extended domains, and connector domains.
% %
% parameter domains: % :parameter domains:
% Vector of domain instances. % Vector of domain instances.
% return: % :return:
% Instance of class 'Stack'. % Instance of class 'Stack'.
checklib; checklib;
s.st_id = -1; s.stID = -1;
s.domains = domains; s.domains = domains;
if nargin == 1 if nargin == 1
nd = length(domains); nd = length(domains);
ids = zeros(1, nd); ids = zeros(1, nd);
for n=1:nd for n=1:nd
ids(n) = domains(n).dom_id; ids(n) = domains(n).domID;
end end
s.st_id = calllib(ct, 'sim1D_new', nd, ids); s.stID = calllib(ct, 'sim1D_new', nd, ids);
else else
help(Stack); help(Stack);
error('Wrong number of parameters.'); error('Wrong number of :parameters.');
end end
% if s.st_id < 0 % if s.stID < 0
% error(geterr); % error(geterr);
% end % end
end end
@ -42,24 +42,23 @@ classdef Stack < handle
function st_clear(s) function st_clear(s)
% Delete the Sim1D object % Delete the Sim1D object
checklib;
calllib(ct, 'sim1D_del', s.st_id); calllib(ct, 'sim1D_del', s.stID);
end end
function display(s, fname) function display(s, fname)
% Show all domains. % Show all domains.
% %
% parameter s: % :parameter s:
% Instance of class 'Stack'. % Instance of class 'Stack'.
% parameter fname: % :parameter fname:
% File to write summary to. If omitted, output is to the % File to write summary to. If omitted, output is to the
% command window. % command window.
checklib;
if nargin == 1 if nargin == 1
fname = '-'; fname = '-';
end end
calllib(ct, 'sim1D_showSolution', s.st_id, fname); calllib(ct, 'sim1D_showSolution', s.stID, fname);
end end
%% Stack Methods %% Stack Methods
@ -67,19 +66,18 @@ classdef Stack < handle
function n = stackIndex(s, name) function n = stackIndex(s, name)
% Get the index of a domain in a stack given its name. % Get the index of a domain in a stack given its name.
% %
% parameter s: % :parameter s:
% Instance of class 'Stack'. % Instance of class 'Stack'.
% parameter name: % :parameter name:
% If double, the value is returned. Otherwise, the name is % If double, the value is :returned. Otherwise, the name is
% looked up and its index is returned. % looked up and its index is :returned.
% return: % :return:
% Index of domain. % Index of domain.
checklib;
if isa(name, 'double') if isa(name, 'double')
n = name; n = name;
else else
n = calllib(ct, 'sim1D_domainIndex', s.st_id, name); n = calllib(ct, 'sim1D_domainIndex', s.stID, name);
if n >= 0 if n >= 0
n = n+1; n = n+1;
else else
@ -91,11 +89,11 @@ classdef Stack < handle
function z = grid(s, name) function z = grid(s, name)
% Get the grid in one domain. % Get the grid in one domain.
% %
% parameter s: % :parameter s:
% Instance of class 'Stack'. % Instance of class 'Stack'.
% parameter name: % :parameter name:
% Name of the domain for which the grid should be retrieved. % Name of the domain for which the grid should be retrieved.
% return: % :return:
% The grid in domain name. % The grid in domain name.
n = s.stackIndex(name); n = s.stackIndex(name);
@ -106,12 +104,12 @@ classdef Stack < handle
function plotSolution(s, domain, component) function plotSolution(s, domain, component)
% Plot a specified solution component. % Plot a specified solution component.
% %
% parameter s: % :parameter s:
% Instance of class 'Stack'. % Instance of class 'Stack'.
% parameter domain: % :parameter domain:
% Name of domain from which the component should be % Name of domain from which the component should be
% retrieved. % retrieved.
% parameter component: % :parameter component:
% Name of the component to be plotted. % Name of the component to be plotted.
n = s.stackIndex(domain); n = s.stackIndex(domain);
@ -126,15 +124,13 @@ classdef Stack < handle
function r = resid(s, domain, rdt, count) function r = resid(s, domain, rdt, count)
% Get the residuals. % Get the residuals.
% %
% parameter s: % :parameter s:
% Instance of class 'Stack'. % Instance of class 'Stack'.
% parameter domain: % :parameter domain:
% Name of the domain. % Name of the domain.
% parameter rdt: % :parameter rdt:
% parameter count: % :parameter count:
% return: % :return:
checklib;
if nargin == 2 if nargin == 2
rdt = 0.0; rdt = 0.0;
@ -148,11 +144,11 @@ classdef Stack < handle
np = d.nPoints; np = d.nPoints;
r = zeros(nc, np); r = zeros(nc, np);
calllib(ct, 'sim1D_eval', s.st_id, rdt, count); calllib(ct, 'sim1D_eval', s.stID, rdt, count);
for m = 1:nc for m = 1:nc
for n = 1:np for n = 1:np
r(m, n) = calllib(ct, 'sim1D_workValue', ... r(m, n) = calllib(ct, 'sim1D_workValue', ...
s.st_id, idom - 1, m - 1, n - 1); s.stID, idom - 1, m - 1, n - 1);
end end
end end
end end
@ -162,29 +158,28 @@ classdef Stack < handle
% This method can be used ot provide an initial guess for the % This method can be used ot provide an initial guess for the
% solution. % solution.
% %
% parameter s: % :parameter s:
% Instance of class 'Stack'. % Instance of class 'Stack'.
% parameter fname: % :parameter fname:
% File name of an XML file containing solution info. % File name of an XML file containing solution info.
% parameter id: % :parameter id:
% ID of the element that should be restored. % ID of the element that should be restored.
checklib;
calllib(ct, 'sim1D_restore', s.st_id, fname, id) calllib(ct, 'sim1D_restore', s.stID, fname, id)
end end
function saveSoln(s, fname, id, desc) function saveSoln(s, fname, id, desc)
% Save a solution to a file. % Save a solution to a file.
% The output file is in a format that can be used by 'restore'. % The output file is in a format that can be used by 'restore'.
% %
% parameter s: % :parameter s:
% Instance of class 'Stack'. % Instance of class 'Stack'.
% parameter fname: % :parameter fname:
% File name where XML file should be written. % File name where XML file should be written.
% parameter id: % :parameter id:
% ID to be assigned to the XMl element when it is written. % ID to be assigned to the XMl element when it is written.
% parameter desc: % :parameter desc:
% Description to be written to the output file. % Description to be written to the output file.
checklib;
if nargin == 1 if nargin == 1
fname = 'soln.xml'; fname = 'soln.xml';
@ -196,22 +191,22 @@ classdef Stack < handle
elseif nargin == 3 elseif nargin == 3
desc = '--'; desc = '--';
end end
calllib(ct, 'sim1D_save', s.st_id, fname, id, desc); calllib(ct, 'sim1D_save', s.stID, fname, id, desc);
end end
function setFlatProfile(s, domain, comp, v) function setFlatProfile(s, domain, comp, v)
% Set a component to a value across the entire domain. % Set a component to a value across the entire domain.
% %
% parameter s: % :parameter s:
% Instance of class 'Stack'. % Instance of class 'Stack'.
% parameter domain: % :parameter domain:
% Integer ID of the domain. % Integer ID of the domain.
% parameter comp: % :parameter comp:
% Component to be set. % Component to be set.
% parameter v: % :parameter v:
% Double value to be set. % Double value to be set.
checklib;
calllib(ct, 'sim1D_setFlatProfile', s.st_id, ... calllib(ct, 'sim1D_setFlatProfile', s.stID, ...
domain - 1, comp - 1, v); domain - 1, comp - 1, v);
end end
@ -219,19 +214,18 @@ classdef Stack < handle
% Set the number of times the Jacobian will be used before it % Set the number of times the Jacobian will be used before it
% is recomputed. % is recomputed.
% %
% parameter s: % :parameter s:
% Instance of class 'Stack'. % Instance of class 'Stack'.
% parameter ss_age: % :parameter ss_age:
% Maximum age of the Jacobian for steady state analysis. % Maximum age of the Jacobian for steady state analysis.
% parameter ts_age: % :parameter ts_age:
% Maximum age of the Jacobian for transient analysis. If not % Maximum age of the Jacobian for transient analysis. If not
% specified, deftauls to 'ss_age'. % specified, deftauls to 'ss_age'.
checklib;
if nargin == 2 if nargin == 2
ts_age = ss_age; ts_age = ss_age;
end end
calllib(ct, 'sim1D_setMaxJacAge', s.st_id, ss_age, ts_age); calllib(ct, 'sim1D_setMaxJacAge', s.stID, ss_age, ts_age);
end end
function setProfile(s, name, comp, p) function setProfile(s, name, comp, p)
@ -251,19 +245,17 @@ classdef Stack < handle
% >> v = [500, 650, 700, 730, 800, 900]; % >> v = [500, 650, 700, 730, 800, 900];
% >> s.setProfile(1, 2, [zr, v]); % >> s.setProfile(1, 2, [zr, v]);
% %
% parameter s: % :parameter s:
% Instance of class 'Stack'. % Instance of class 'Stack'.
% parameter name: % :parameter name:
% Domain name. % Domain name.
% parameter comp: % :parameter comp:
% Component number. % Component number.
% parameter p: % :parameter p:
% n x 2 array, whose columns are the relative (normalized) % n x 2 array, whose columns are the relative (normalized)
% positions and the component values at those points. The % positions and the component values at those points. The
% number of positions 'n' is arbitrary. % number of positions 'n' is arbitrary.
checklib;
if isa(name, 'double') if isa(name, 'double')
n = name; n = name;
else else
@ -285,12 +277,12 @@ classdef Stack < handle
if sz(1) == np + 1; if sz(1) == np + 1;
for j = 1:np for j = 1:np
ic = d.componentIndex(c{j}); ic = d.componentIndex(c{j});
calllib(ct, 'sim1D_setProfile', s.st_id, ... calllib(ct, 'sim1D_setProfile', s.stID, ...
n - 1, ic - 1, sz(1), p(1, :), sz(1), p(j+1, :)); n - 1, ic - 1, sz(1), p(1, :), sz(1), p(j+1, :));
end end
elseif sz(2) == np + 1; elseif sz(2) == np + 1;
ic = d.componentIndex(c{j}); ic = d.componentIndex(c{j});
calllib(ct, 'sim1D_setProfile', s.st_id, ... calllib(ct, 'sim1D_setProfile', s.stID, ...
n - 1, ic - 1, sz(2), p(:, 1), sz(2), p(:, j+1)); n - 1, ic - 1, sz(2), p(:, 1), sz(2), p(:, j+1));
else else
error('Wrong profile shape.'); error('Wrong profile shape.');
@ -300,24 +292,22 @@ classdef Stack < handle
function setRefineCriteria(s, n, ratio, slope, curve, prune) function setRefineCriteria(s, n, ratio, slope, curve, prune)
% Set the criteria used to refine the grid. % Set the criteria used to refine the grid.
% %
% parameter s: % :parameter s:
% Instance of class 'Stack'. % Instance of class 'Stack'.
% parameter ratio: % :parameter ratio:
% Maximum size ratio between adjacent cells. % Maximum size ratio between adjacent cells.
% parameter slope: % :parameter slope:
% Maximum relative difference in value between adjacent % Maximum relative difference in value between adjacent
% points. % points.
% parameter curve: % :parameter curve:
% Maximum relative difference in slope between adjacent % Maximum relative difference in slope between adjacent
% cells. % cells.
% parameter prune: % :parameter prune:
% Minimum value for slope or curve for which points will be % Minimum value for slope or curve for which points will be
% retained or curve value is below prune for all components, % retained or curve value is below prune for all components,
% it will be deleted, unless either neighboring point is % it will be deleted, unless either neighboring point is
% already marked for deletion. % already marked for deletion.
checklib;
if nargin < 3 if nargin < 3
ratio = 10.0; ratio = 10.0;
end end
@ -330,23 +320,23 @@ classdef Stack < handle
if nargin < 6 if nargin < 6
prune = -0.1; prune = -0.1;
end end
calllib(ct, 'sim1D_setRefineCriteria', s.st_id, ... calllib(ct, 'sim1D_setRefineCriteria', s.stID, ...
n - 1, ratio, slope, curve, prune); n - 1, ratio, slope, curve, prune);
end end
function setTimeStep(s, stepsize, steps) function setTimeStep(s, stepsize, steps)
% Specify a sequence of time steps. % Specify a sequence of time steps.
% %
% parameter stepsize: % :parameter stepsize:
% Initial step size. % Initial step size.
% parameter steps: % :parameter steps:
% Vector of number of steps to take before re-attempting % Vector of number of steps to take before re-attempting
% solution of steady-state problem. % solution of steady-state problem.
% For example, steps = [1, 2, 5, 10] would cause one time % For example, steps = [1, 2, 5, 10] would cause one time
% step to be taken first time the steady-state solution % step to be taken first time the steady-state solution
% attempted. If this failed, two time steps would be taken. % attempted. If this failed, two time steps would be taken.
checklib;
calllib(ct, 'sim1D_TimeStep', s.st_id, ... calllib(ct, 'sim1D_TimeStep', s.stID, ...
stepsize, length(steps), steps); stepsize, length(steps), steps);
end end
@ -363,37 +353,36 @@ classdef Stack < handle
% the global index of the point, wchih depends on the location % the global index of the point, wchih depends on the location
% of this domain in the stack. % of this domain in the stack.
% %
% parameter s: % :parameter s:
% Instance of class 'Stack'. % Instance of class 'Stack'.
% parameter n: % :parameter n:
% Domain number. % Domain number.
% parameter comp: % :parameter comp:
% Component number. % Component number.
% parameter localPoints: % :parameter localPoints:
% Local index of the grid point in the domain. % Local index of the grid point in the domain.
% parameter v: % :parameter v:
% Value to be set. % Value to be set.
checklib;
calllib(ct, 'sim1D_setValue', s.st_id, ... calllib(ct, 'sim1D_setValue', s.stID, ...
n - 1, comp - 1, localPoints - 1, v); n - 1, comp - 1, localPoints - 1, v);
end end
function x = solution(s, domain, component) function x = solution(s, domain, component)
% Get a solution component in one domain. % Get a solution component in one domain.
% %
% parameter s: % :parameter s:
% Instance of class 'Stack'. % Instance of class 'Stack'.
% parameter domain: % :parameter domain:
% String name of the domain from which the solution is % String name of the domain from which the solution is
% desired. % desired.
% parameter component: % :parameter component:
% String component for which the solution is desired. If % String component for which the solution is desired. If
% omitted, solution for all of the components will be % omitted, solution for all of the components will be
% returned in an 'nPoints' x 'nComponnts' array. % :returned in an 'nPoints' x 'nComponnts' array.
% return: % :return:
% Either an 'nPoints' x 1 vector, or 'nPoints' x % Either an 'nPoints' x 1 vector, or 'nPoints' x
% 'nCOmponents' array. % 'nCOmponents' array.
checklib;
idom = s.stackIndex(domain); idom = s.stackIndex(domain);
d = s.domains(idom); d = s.domains(idom);
@ -402,7 +391,7 @@ classdef Stack < handle
icomp = d.componentIndex(component); icomp = d.componentIndex(component);
x = zeros(1, np); x = zeros(1, np);
for n = 1:np for n = 1:np
x(n) = calllib(ct, 'sim1D_value', s.st_id, ... x(n) = calllib(ct, 'sim1D_value', s.stID, ...
idom - 1, icomp - 1, n - 1); idom - 1, icomp - 1, n - 1);
end end
else else
@ -410,7 +399,7 @@ classdef Stack < handle
x = zeros(nc, np); x = zeros(nc, np);
for m = 1:nc for m = 1:nc
for n = 1:np for n = 1:np
x(m, n) = calllib(ct, 'sim1D_value', s.st_id, ... x(m, n) = calllib(ct, 'sim1D_value', s.stID, ...
idom - 1, m - 1, n - 1); idom - 1, m - 1, n - 1);
end end
end end
@ -420,38 +409,25 @@ classdef Stack < handle
function solve(s, loglevel, refine_grid) function solve(s, loglevel, refine_grid)
% Solve the problem. % Solve the problem.
% %
% parameter s: % :parameter s:
% Instance of class 'Stack'. % Instance of class 'Stack'.
% parameter loglevel: % :parameter loglevel:
% Integer flag controlling the amount of diagnostic output. % Integer flag controlling the amount of diagnostic output.
% Zero supresses all output, and 5 produces very verbose % Zero supresses all output, and 5 produces very verbose
% output. % output.
% parameter refine_grid: % :parameter refine_grid:
% Integer, 1 to allow grid refinement, 0 to disallow. % 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) calllib(ct, 'sim1D_solve', s.stID, loglevel, refine_grid);
% % Redefine subscripted references. end
% 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) function writeStats(s)
% Print statistics for the current solution. % Print statistics for the current solution.
% Prints a summary of the number of function and Jacobian % Prints a summary of the number of function and Jacobian
% evaluations for each grid, and the CPU time spent on each % evaluations for each grid, and the CPU time spent on each
% one. % one.
checklib;
calllib(ct, 'sim1D_writeStats', s.st_id, 1); calllib(ct, 'sim1D_writeStats', s.stID, 1);
end end
end end

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

@ -8,13 +8,13 @@ classdef Interface < handle & ThermoPhase & Kinetics
%% Interface class constructor %% Interface class constructor
function s = Interface(src, id, p1, p2, p3, p4) function s = Interface(src, id, p1, p2, p3, p4)
% parameter src: % :parameter src:
% CTI or CTML file containing the interface or edge phase. % CTI or CTML file containing the interface or edge phase.
% parameter 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.
% parameter 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;
@ -32,7 +32,7 @@ classdef Interface < handle & ThermoPhase & Kinetics
args = {p1, p2, p3, p4}; args = {p1, p2, p3, p4};
end end
s@Kinetics(t, src, id, args{:}); s@Kinetics(t, src, id, args{:});
s.tp_id = t.tp_id; s.tpID = t.tpID;
end end
%% Interface methods %% Interface methods
@ -40,78 +40,46 @@ classdef Interface < handle & ThermoPhase & Kinetics
function c = get.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 % Vector of length "n_surf_species" for coverage.
% plotted. Otherwise, a vector of length "n_surf_species"
% will be returned.
checklib; surfID = s.tpID;
surf_id = s.tp_id;
nsp = s.nSpecies; nsp = s.nSpecies;
xx = zeros(1, nsp); xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(ct, 'surf_getCoverages', surf_id, xx); calllib(ct, 'surf_getCoverages', surfID, xx);
c = pt.Value; c = pt.Value;
% if nargout == 0
% figure
% set(gcf, 'Name', 'Coverages')
% bar(c);
% colormap(summer);
% nm = s.speciesNames;
% set(gca, 'XTickLabel', nm);
% xlabel('Species Name');
% ylabel('Coverage');
% title('Surface Species Coverages');
% end
end end
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 % Vector of length "n_surf_species" for concentration.
% plotted. Otherwise, a vector of length "n_surf_species"
% will be returned.
checklib; surfID = s.tr_id;
surf_id = s.tr_id;
nsp = s.nSpecies; nsp = s.nSpecies;
xx = zeros(1, nsp); xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(ct, 'surf_getConcentrations', surf_id, xx); calllib(ct, 'surf_getConcentrations', surfID, xx);
c = pt.Value; c = pt.Value;
% if nargout == 0
% figure
% set(gcf, 'Name', 'Concentrations')
% bar(c);
% colormap(summer);
% nm = speciesNames(s);
% set(gca, 'XTickLabel', nm);
% xlabel('Species Name');
% ylabel('Concentration [kmol/m^2]');
% title('Surface Species Concentrations');
% end
end end
function set.coverages(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.
% %
% parameter 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".
checklib;
if nargin == 3 && strcmp(norm, 'nonorm') if nargin == 3 && strcmp(norm, 'nonorm')
norm_flag = 0; norm_flag = 0;
else else
norm_flag = 1; norm_flag = 1;
end end
surf_id = s.tr_id; surfID = s.tr_id;
nsp = s.nSpecies; nsp = s.nSpecies;
[m, n] = size(cov); [m, n] = size(cov);
@ -119,14 +87,14 @@ classdef Interface < handle & ThermoPhase & Kinetics
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(ct, 'surf_setCoverages', surf_id, cov, norm_flag); calllib(ct, 'surf_setCoverages', surfID, 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(ct, 'surf_setCoveragesByName', surf_id, cov); calllib(ct, 'surf_setCoveragesByName', surfID, cov);
end end
end end

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

@ -1,8 +1,7 @@
classdef Kinetics < handle classdef Kinetics < handle
properties properties
owner kinID
kin_id
Kc % equilibrium constant Kc % equilibrium constant
Kf % forward reaction rate Kf % forward reaction rate
Kr % reverse reaction rate Kr % reverse reaction rate
@ -15,6 +14,7 @@ classdef Kinetics < handle
%% Kinetics class constructor %% Kinetics class constructor
function kin = Kinetics(ph, src, id, n1, n2, n3, n4) function kin = Kinetics(ph, src, id, n1, n2, n3, n4)
checklib; checklib;
% indices for bulk phases in a heterogeneous mechanism % indices for bulk phases in a heterogeneous mechanism
inb1 = -1; inb1 = -1;
@ -24,105 +24,88 @@ classdef Kinetics < handle
if nargin == 2 if nargin == 2
id = '-'; id = '-';
end end
kin.owner = 1;
% get the integer indices used to find the stored objects % get the integer indices used to find the stored objects
% representing the phases participating in the mechanism % representing the phases participating in the mechanism
iph = ph.tp_id; iph = ph.tpID;
if nargin > 3 if nargin > 3
inb1 = n1.tp_id; inb1 = n1.tpID;
if nargin > 4 if nargin > 4
inb2 = n2.tp_id; inb2 = n2.tpID;
if nargin > 5 if nargin > 5
inb3 = n3.tp_id; inb3 = n3.tpID;
if nargin > 6 if nargin > 6
inb4 = n4.tp_id; inb4 = n4.tpID;
end end
end end
end end
end end
kin.kin_id = calllib(ct, 'kin_newFromFile', src, id, ... kin.kinID = calllib(ct, 'kin_newFromFile', src, id, ...
iph, inb1, inb2, inb3, inb4); iph, inb1, inb2, inb3, inb4);
end end
%% Utility methods %% Utility methods
function kin_clear(kin) function kinClear(kin)
% Delete the kernel object % Delete the kernel object
checklib; calllib(ct, 'kin_del', kin.kinID);
calllib(ct, 'kin_del', kin.kin_id);
end end
%% Get scalar attributes %% Get scalar attributes
function n = isReversible(kin, i)
% Get an array of flags indicating reversibility of a reaction.
%
% parameter i:
% Integer reaction number.
% return:
% 1 if reaction number i is reversible. 0 if irreversible.
checklib;
n = calllib(ct, 'kin_isReversible', kin.kin_id, i);
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.
% %
% parameter 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; n = calllib(ct, 'kin_multiplier', kin.kinID, irxn-1);
n = calllib(ct, 'kin_multiplier', kin.kin_id, irxn-1);
end end
function n = nReactions(kin) function n = nReactions(kin)
% Get the number of reactions. % Get the number of reactions.
% %
% return: % :return:
% Integer number of reactions % Integer number of reactions
checklib; n = calllib(ct, 'kin_nReactions', kin.kinID);
n = calllib(ct, 'kin_nReactions', kin.kin_id);
end end
function n = nSpecies2(kin) function n = nTotalSpecies(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; n = calllib(ct, 'kin_nSpecies', kin.kinID);
n = calllib(ct, 'kin_nSpecies', kin.kin_id);
end end
function n = stoich_r(kin, species, rxns) function n = stoichReactant(kin, species, rxns)
% Get the reactant stoichiometric coefficients. % Get the reactant stoichiometric coefficients.
% %
% parameter 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.
% parameter 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
% reaction i. If "species" and "rxns" are specified, the % reaction i. If "species" and "rxns" are specified, the
% matrix will contain only entries for the specified species % matrix will contain only entries for the specified species
% and reactions. For example, "stoich_p(a, 3, [1, 3, 5, % and reactions. For example, "stoich_p(a, 3, [1, 3, 5,
% 7])" returns a sparse matrix containing only the % 7])" :returns a sparse matrix containing only the
% coefficients for specis 3 in reactions 1, 3, 5, and 7. % coefficients for specis 3 in reactions 1, 3, 5, and 7.
checklib; nsp = kin.nTotalSpecies;
nsp = kin.nSpecies;
nr = kin.nReactions; nr = kin.nReactions;
temp = sparse(nsp, nr); temp = sparse(nsp, nr);
if nargin == 1 if nargin == 1
@ -137,7 +120,7 @@ classdef Kinetics < handle
for k = krange for k = krange
for i = irange for i = irange
t = calllib(ct, 'kin_reactantStoichCoeff', ... t = calllib(ct, 'kin_reactantStoichCoeff', ...
kin.kin_id, k-1, i-1); kin.kinID, k-1, i-1);
if t ~= 0.0 if t ~= 0.0
temp(k, i) = t; temp(k, i) = t;
end end
@ -147,29 +130,28 @@ classdef Kinetics < handle
n = temp; n = temp;
end end
function n = stoich_p(kin, species, rxns) function n = stoichProduct(kin, species, rxns)
% Get the product stoichiometric coefficients. % Get the product stoichiometric coefficients.
% %
% parameter 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.
% parameter 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
% reaction i. If "species" and "rxns" are specified, the % reaction i. If "species" and "rxns" are specified, the
% matrix will contain only entries for the specified species % matrix will contain only entries for the specified species
% and reactions. For example, "stoich_p(a, 3, [1, 3, 5, % and reactions. For example, "stoich_p(a, 3, [1, 3, 5,
% 7])" returns a sparse matrix containing only the % 7])" :returns a sparse matrix containing only the
% coefficients for specis 3 in reactions 1, 3, 5, and 7. % coefficients for specis 3 in reactions 1, 3, 5, and 7.
checklib; nsp = kin.nTotalSpecies;
nsp = kin.nSpecies;
nr = kin.nReactions; nr = kin.nReactions;
temp = sparse(nsp, nr); temp = sparse(nsp, nr);
if nargin == 1 if nargin == 1
@ -184,38 +166,37 @@ classdef Kinetics < handle
for k = krange for k = krange
for i = irange for i = irange
t = calllib(ct, 'kin_productStoichCoeff', ... t = calllib(ct, 'kin_productStoichCoeff', ...
kin.kin_id, k-1, i-1); kin.kinID, k-1, i-1);
if t ~= 0.0 if t ~= 0.0
temp(k, i) = t; temp(k, i) = t;
end end
end end
end end
n = temp; n = temp;
end end
function n = stoich_net(kin, species, rxns) function n = stoichNet(kin, species, rxns)
% Get the net stoichiometric coefficients. % Get the net stoichiometric coefficients.
% %
% parameter 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.
% parameter 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 % A sparse matrix of all net stoichiometric coefficients.
% coefficients. The matrix elements "nu(k, i)" is the % The matrix elements "nu(k, i)" is the stoichiometric
% stoichiometric coefficient of species k as a net in % coefficient of species k as a net in reaction.
% reaction i. If "species" and "rxns" are specified, the % If "species" and "rxns" are specified, the matrix will
% matrix will contain only entries for the specified species % contain only entries for the specified species and reactions.
% and reactions. For example, "stoich_net(a, 3, [1, 3, 5, % For example, "stoich_net(a, 3, [1, 3, 5, 7])" returns a
% 7])" returns a sparse matrix containing only the % sparse matrix containing only the coefficients for
% coefficients for specis 3 in reactions 1, 3, 5, and 7. % specis 3 in reactions 1, 3, 5, and 7.
checklib;
if nargin == 1 if nargin == 1
n = stoich_p(kin)-stoich_r(kin); n = stoich_p(kin)-stoich_r(kin);
elseif nargin == 3 elseif nargin == 3
@ -229,178 +210,116 @@ classdef Kinetics < handle
function cdot = creationRates(kin) function cdot = creationRates(kin)
% Get the chemical reaction rates. % Get the chemical reaction rates.
% %
% return: % :return:
% Returns a vector of the creation rates of all species. If % A vector of the creation rates of all species. Unit: kmol/m^3-s.
% the output is not assigned to a variable, a bar graph is
% produced. Unit: kmol/m^3-s.
checklib; nsp = kin.nTotalSpecies;
nsp = kin.nSpecies;
xx = zeros(1, nsp); xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getCreationRates', kin.kin_id, nsp, pt); calllib(ct, 'kin_getCreationRates', kin.kinID, nsp, pt);
cdot = pt.Value; 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 end
function ddot = destructionRates(kin) function ddot = destructionRates(kin)
% Get the chemical destruction rates. % Get the chemical destruction rates.
% %
% return: % :return:
% Returns a vector of the destruction rates of all species. % A vector of the destruction rates of all species. Unit: kmol/m^3-s.
% If the output is not assigned to a variable, a bar graph
% is produced. Unit: kmol/m^3-s.
checklib; nsp = kin.nTotalSpecies;
nsp = kin.nSpecies;
xx = zeros(1, nsp); xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getDestructionRates', kin.kin_id, nsp, pt); calllib(ct, 'kin_getDestructionRates', kin.kinID, nsp, pt);
ddot = pt.Value; ddot = pt.Value;
% if nargout == 0 end
% figure
% set(gcf, 'Name', 'Destruction Rates') function n = isReversible(kin, i)
% bar(q) % Get an array of flags indicating reversibility of a reaction.
% xlabel('Species Number') %
% ylabel('Destruction Rate [kmol/m^3-s]') % :parameter i:
% title('Species Chemical Reaction Rates') % Integer reaction number.
% end % :return:
% 1 if reaction number i is reversible. 0 if irreversible.
n = calllib(ct, 'kin_isReversible', kin.kinID, i);
end end
function wdot = netProdRates(kin) function wdot = netProdRates(kin)
% Get the net chemical production rates for all species. % Get the net chemical production rates for all species.
% %
% return: % :return:
% Returns a vector of the net production (creation-destruction) % A vector of the net production (creation-destruction)
% rates of all species. If the output is not assigned to a % rates of all species. Unit: kmol/m^3-s.
% variable, a bar graph is produced. Unit: kmol/m^3-s.
checklib; nsp = kin.nTotalSpecies;
nsp = kin.nSpecies;
xx = zeros(1, nsp); xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getNetProductionRates', kin.kin_id, nsp, pt); calllib(ct, 'kin_getNetProductionRates', kin.kinID, nsp, pt);
wdot = pt.Value; 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 end
function q = rop_f(kin) function q = ropForward(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 % A column vector of the forward rates of progress for all
% for all reactions. If this function is called without % reactions.
% argument, a bar graph is produced.
checklib;
nr = kin.nReactions; nr = kin.nReactions;
xx = zeros(1, nr); xx = zeros(1, nr);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getFwdRateOfProgress', kin.kin_id, nr, pt); calllib(ct, 'kin_getFwdRateOfProgress', kin.kinID, nr, pt);
q = pt.Value; q = pt.Value;
% if nargout == 0
% figure
% set(gcf, 'Name', 'Rates of Progress')
% bar(q)
% xlabel('Reaction Number')
% ylabel('Forward Rate of Progress [kmol/m^3]')
% title('Forward Rates of Progress')
% end
end end
function q = rop_r(kin) function q = ropReverse(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 % A column vector of the reverse rates of progress for all
% for all reactions. If this function is called without % reactions.
% argument, a bar graph is produced.
checklib;
nr = kin.nReactions; nr = kin.nReactions;
xx = zeros(1, nr); xx = zeros(1, nr);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getRevRateOfProgress', kin.kin_id, nr, pt); calllib(ct, 'kin_getRevRateOfProgress', kin.kinID, nr, pt);
q = pt.Value; q = pt.Value;
% if nargout == 0
% figure
% set(gcf, 'Name', 'Rates of Progress')
% bar(q)
% xlabel('Reaction Number')
% ylabel('Reverse Rate of Progress [kmol/m^3]')
% title('Reverse Rates of Progress')
% end
end end
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, % An I x 2 array of reaction rates of progress, where I is
% where I is the number of reactions. The first column % the number of reactions. The first column contains the
% contains the forward rates of progress, and the second % forward rates of progress, and the second column the
% column the reverse rates. If this function is called % reverse rates.
% without arguments, a bar graph is produced.
checklib;
f = rop_f(kin); f = rop_f(kin);
r = rop_r(kin); r = rop_r(kin);
q = [f, r]; q = [f, r];
% if nargout == 0
% figure
% set(gcf, 'Name', 'Rates of Progress')
% bar(q)
% xlabel('Reaction Number')
% ylabel('Rate of Progress [kmol/m^3]')
% title('Rates of Progress')
% legend('Forward', 'Reverse')
% end
end end
function q = rop_net(kin) function q = ropNet(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 % A column vector of the net rates of progress for all
% for all reactions. If this function is called without % reactions.
% argument, a bar graph is produced.
checklib;
nr = kin.nReactions; nr = kin.nReactions;
xx = zeros(1, nr); xx = zeros(1, nr);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getNetRatesOfProgress', kin.kin_id, nr, pt); calllib(ct, 'kin_getNetRatesOfProgress', kin.kinID, nr, pt);
q = pt.Value; q = pt.Value;
% if nargout == 0
% figure
% set(gcf, 'Name', 'Rates of Progress')
% bar(q)
% xlabel('Reaction Number')
% ylabel('Net Rate of Progress [kmol/m^3]')
% title('Net Rates of Progress')
% end
end end
function rxn = reactionEqn(kin, irxn) function rxn = reactionEqn(kin, irxn)
% Get the reaction equation of a reaction % Get the reaction equation of a reaction
% %
% parameter irxn: % :parameter irxn:
% Optional. Integer or vector of reaction numbers. % Optional. Integer or vector of reaction numbers.
% return: % :return:
% String or cell arrray of strings of the reaction % String or cell arrray of strings of the reaction
% equations. % equations.
@ -419,12 +338,12 @@ classdef Kinetics < handle
rxn = cell(m, n); rxn = cell(m, n);
for i = 1:m for i = 1:m
for j = 1:n for j = 1:n
buflen = calllib(ct, 'kin_getReactionString', kin.kin_id, ... buflen = calllib(ct, 'kin_getReactionString', kin.kinID, ...
irxn - 1, 0, ''); irxn - 1, 0, '');
if buflen > 0 if buflen > 0
aa = char(zeros(1, buflen)); aa = char(zeros(1, buflen));
[~, aa] = calllib(ct, 'kin_getReactionString', ... [~, aa] = calllib(ct, 'kin_getReactionString', ...
kin.kin_id, irxn - 1, buflen, aa); kin.kinID, irxn - 1, buflen, aa);
rxn{i, j} = aa; rxn{i, j} = aa;
end end
end end
@ -434,115 +353,99 @@ classdef Kinetics < handle
function enthalpy = get.dH(kin) function enthalpy = get.dH(kin)
% Get the enthalpy of reaction for each reaction. % Get the enthalpy of reaction for each reaction.
% %
% return: % :return:
% Returns a vector of the enthalpy of reaction for each % A vector of the enthalpy of reaction for each reaction.
% reaction. Unit: J/kmol. % Unit: J/kmol.
checklib;
nr = kin.nReactions; nr = kin.nReactions;
xx = zeros(1, nr); xx = zeros(1, nr);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getDelta', kin.kin_id, 0, nr, pt); calllib(ct, 'kin_getDelta', kin.kinID, 0, nr, pt);
enthalpy = pt.Value; enthalpy = pt.Value;
end end
function entropy = get.dS(kin) function entropy = get.dS(kin)
% Get the entropy of reaction for each reaction. % Get the entropy of reaction for each reaction.
% %
% return: % :return:
% Returns a vector of the entropy of reaction for each % A vector of the entropy of reaction for each reaction.
% reaction. Unit: J/kmol-K. % Unit: J/kmol-K.
checklib;
nr = kin.nReactions; nr = kin.nReactions;
xx = zeros(1, nr); xx = zeros(1, nr);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getDelta', kin.kin_id, 2, nr, pt); calllib(ct, 'kin_getDelta', kin.kinID, 2, nr, pt);
entropy = pt.Value; entropy = pt.Value;
end end
function gibbs = get.dG(kin) function gibbs = get.dG(kin)
% Get the Gibbs free energy of reaction for each reaction. % Get the Gibbs free energy of reaction for each reaction.
% %
% return: % :return:
% Returns a vector of the Gibbs free energy of reaction for % A vector of the Gibbs free energy of reaction for each
% each reaction. Unit: J/kmol. % reaction. Unit: J/kmol.
checklib;
nr = kin.nReactions; nr = kin.nReactions;
xx = zeros(1, nr); xx = zeros(1, nr);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getDelta', kin.kin_id, 1, nr, pt); calllib(ct, 'kin_getDelta', kin.kinID, 1, nr, pt);
gibbs = pt.Value; gibbs = pt.Value;
end 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 % A column vector of the equilibrium constants for all
% all reactions. The vector has an entry for every reaction, % 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
% for the reversible reactions. If the output is not % for the reversible reactions.
% assigned to a variable, a bar graph is produced.
checklib;
nr = kin.nReactions; nr = kin.nReactions;
xx = zeros(1, nr); xx = zeros(1, nr);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getEquilibriumConstants', kin.kin_id, nr, pt); calllib(ct, 'kin_getEquilibriumConstants', kin.kinID, nr, pt);
k = pt.Value; k = pt.Value;
% if nargout == 0
% figure
% set(gcf, 'Name', 'Equilibrium Constants')
% bar(k)
% xlabel('Reaction Number')
% ylabel('log_{10} Kc [kmol,m, s]')
% title('Equilibrium Constants')
% end
end end
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 % A column vector of the forward rates constants of all
% all of the reactions. % reactions.
checklib;
nr = kin.nReactions; nr = kin.nReactions;
xx = zeros(1, nr); xx = zeros(1, nr);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getFwdRateConstants', kin.kin_id, nr, pt); calllib(ct, 'kin_getFwdRateConstants', kin.kinID, 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 % A column vector of the reverse rates constants of all
% all of the reactions. % reactions.
checklib;
nr = kin.nReactions; nr = kin.nReactions;
xx = zeros(1, nr); xx = zeros(1, nr);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getRevRateConstants', kin.kin_id, 1, nr, pt); calllib(ct, 'kin_getRevRateConstants', kin.kinID, 1, nr, pt);
k = pt.Value; k = pt.Value;
end end
function massProdRate = ydot(kin) function massProdRate = ydot(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. % A vector of the mass production rates.
checklib; nsp = kin.nTotalSpecies;
nsp = kin.nSpecies;
xx = zeros(1, nsp); xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getSourceTerms', kin.kin_id, nsp, pt); calllib(ct, 'kin_getSourceTerms', kin.kinID, nsp, pt);
massProdRate = pt.Value; massProdRate = pt.Value;
end end
@ -551,10 +454,10 @@ classdef Kinetics < handle
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.
% %
% parameter 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.
% parameter 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.
@ -571,7 +474,7 @@ classdef Kinetics < handle
for i = 1:nr for i = 1:nr
for j = 1:n for j = 1:n
calllib(ct, 'kin_setMultiplier', kin.kin_id, ... calllib(ct, 'kin_setMultiplier', kin.kinID, ...
irxn(i, j)-1, v); irxn(i, j)-1, v);
end end
end end
@ -580,10 +483,10 @@ classdef Kinetics < handle
function advanceCoverages(kin, dt) function advanceCoverages(kin, dt)
% Advance the surface coveages forward in time % Advance the surface coveages forward in time
% %
% parameter dt: % :parameter dt:
% Time interval by which the coverages should be advanced. % Time interval by which the coverages should be advanced.
calllib(ct, 'kin_advanceCoverages', kin.kin_id, dt); calllib(ct, 'kin_advanceCoverages', kin.kinID, dt);
end end
end end

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

@ -1,7 +1,7 @@
classdef Mixture < handle classdef Mixture < handle
properties properties
mixindex mixID
phases phases
T T
P P
@ -13,30 +13,30 @@ classdef Mixture < handle
function m = Mixture(phases) function m = Mixture(phases)
% To construct a mixture, supply a cell array of phases and mole % To construct a mixture, supply a cell array of phases and mole
% numbers: % numbers:
% >> gas = Solution('gas.cti'); % >> air = Solution('air.yaml');
% >> graphite = Solution('graphite.cti'); % >> graphite = Solution('graphite.yaml');
% >> mix = Mixture({gas, 1.0; graphite, 0.1}); % >> mix = Mixture({air, 1.0; graphite, 0.1});
% %
% Phases may also be added later using the addPhase method: % Phases may also be added later using the addPhase method:
% >> water = Solution('water.cti'); % >> water = Solution('water.cti');
% >> addPhase(mix, water, 3.0); % >> mix.addPhase(water, 3.0);
% %
% Note that the objects representing each phase compute only % Note that the objects representing each phase compute only
% the intensive state of the phase - they do not store any % the intensive state of the phase - they do not store any
% information on the amount of this phase. Mixture objects, on % information on the amount of this phase. Mixture objects, on
% the other hand, represent full extensive state. % the other hand, represent full extensive state.
% %
% Mixture objects are 'lightweight'in the sense that they do % Mixture objects are lightweight in the sense that they do
% not store parameters needed to compute thermodynamic or % not store :parameters needed to compute thermodynamic or
% kinetic properties of the phases. These are contained in the % kinetic properties of the phases. These are contained in the
% ('heavyweight') phase objects. Multiple mixture objects are % ('heavyweight') phase objects. Multiple mixture objects are
% constructed using the same set of phase objects. Each one % constructed using the same set of phase objects. Each one
% store its own state information locally, and syncrhonizes the % stores its own state information locally, and synchronizes the
% phase objects whenever itrequires phase properties. % phase objects whenever it requires phase properties.
% %
% parameter 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;
@ -46,13 +46,13 @@ classdef Mixture < handle
end end
% Create an empty mixture. % Create an empty mixture.
m.mixindex = calllib(ct, 'mix_new'); m.mixID = calllib(ct, 'mix_new');
m.phases = phases; m.phases = phases;
% If phases are supplied, add them % If phases are supplied, add them
if nargin == 1 if nargin == 1
if ~isa(phases, 'cell') if ~isa(phases, 'cell')
error('Enter phasesas a cell array.'); error('Enter phases as a cell array.');
end end
% First column contains the phase objects, and the second % First column contains the phase objects, and the second
@ -74,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(ct, 'mix_updatePhases', m.mixindex); calllib(ct, 'mix_updatePhases', m.mixID);
[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) ...
@ -88,24 +88,19 @@ classdef Mixture < handle
function clear(m) function clear(m)
% Delete the MultiPhase object. % Delete the MultiPhase object.
checklib; calllib(ct, 'mix_del', m.mixID);
calllib(ct, 'mix_del', m.mixindex);
end end
%% Mixture Get methods
function addPhase(m, phase, moles) function addPhase(m, phase, moles)
% Add a phase to the mixture % Add a phase to the mixture
% %
% parameter m: % :parameter m:
% Instance of class 'Mixture' to which phases is added. % Instance of class 'Mixture' to which phases is added.
% parameter phase: % :parameter phase:
% Instance of class 'ThermoPhase' which should be added. % Instance of class 'ThermoPhase' which should be added.
% parameter 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;
if ~isa(phase, 'ThermoPhase') if ~isa(phase, 'ThermoPhase')
error('Phase object of wrong type.'); error('Phase object of wrong type.');
end end
@ -116,49 +111,49 @@ classdef Mixture < handle
error('Negative moles'); error('Negative moles');
end end
iok = calllib(ct, 'mix_addPhase', m.mixindex, phase.tp_id, ... iok = calllib(ct, 'mix_addPhase', m.mixID, phase.tp_id, ...
moles); moles);
if iok < 0 if iok < 0
error('Error adding phase'); error('Error adding phase');
end end
end end
%% Mixture Get methods
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; temperature = calllib(ct, 'mix_temperature', m.mixID);
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; pressure = calllib(ct, 'mix_pressure', m.mixID);
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;
n = calllib(ct, 'mix_nPhases', m.mixindex); n = calllib(ct, 'mix_nPhases', m.mixID);
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;
n = calllib(ct, 'mix_nElements', m.mixindex); n = calllib(ct, 'mix_nElements', m.mixID);
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;
n = calllib(ct, 'mix_nSpecies', m.mixindex); n = calllib(ct, 'mix_nSpecies', m.mixID);
end end
function n = elementIndex(m, name) function n = elementIndex(m, name)
@ -167,8 +162,8 @@ classdef Mixture < handle
% Note: In keeping with the conventions used by Matlab, the % Note: In keeping with the conventions used by Matlab, the
% 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;
n = calllib(ct, 'mix_elementIndex', m.mixindex, name) + 1; n = calllib(ct, 'mix_elementIndex', m.mixID, name) + 1;
end end
function n = speciesIndex(m, k, p) function n = speciesIndex(m, k, p)
@ -177,28 +172,27 @@ classdef Mixture < handle
% Note: In keeping with the conventions used by Matlab, the % Note: In keeping with the conventions used by Matlab, the
% 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;
n = calllib(ct, 'mix_speciesIndex', m.mixindex, k-1, p-1) + 1; n = calllib(ct, 'mix_speciesIndex', m.mixID, 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.
% %
% parameter 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;
if nargin == 2 if nargin == 2
moles = calllib(ct, 'mix_phaseMoles', m.mixindex, n); moles = calllib(ct, 'mix_phaseMoles', m.mixID, 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(ct, 'mix_phaseMoles', ... moles(i) = calllib(ct, 'mix_phaseMoles', ...
m.mixindex, i); m.mixID, i);
end end
else error('wrong number of arguments'); else error('wrong number of arguments');
end end
@ -207,14 +201,13 @@ 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;
nsp = m.nSpecies; nsp = m.nSpecies;
xx = zeros(1, nsp); xx = zeros(1, nsp);
ptr = libpointer('doublePtr', xx); ptr = libpointer('doublePtr', xx);
calllib(ct, 'mix_getChemPotential', m.mixindex, nsp, ptr); calllib(ct, 'mix_getChemPotential', m.mixID, nsp, ptr);
mu = ptr.Value; mu = ptr.Value;
end end
@ -223,30 +216,30 @@ classdef Mixture < handle
function m = set.T(m, temp) function m = set.T(m, temp)
% Set the mixture temperature. % Set the mixture temperature.
% %
% parameter temp: % :parameter temp:
% Temperature to set. Unit: K. % Temperature to set. Unit: K.
checklib;
calllib(ct, 'mix_setTemperature', m.mixindex, temp); calllib(ct, 'mix_setTemperature', m.mixID, temp);
end end
function m = set.P(m, pressure) function m = set.P(m, pressure)
% Set the mixture pressure. % Set the mixture pressure.
% %
% parameter pressure: % :parameter pressure:
% Pressure to set. Unit: Pa. % Pressure to set. Unit: Pa.
checklib;
calllib(ct, 'mix_setPressure', m.mixindex, pressure); calllib(ct, 'mix_setPressure', m.mixID, 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.
% %
% parameter n: % :parameter n:
% Phase number. % Phase number.
% parameter moles: % :parameter moles:
% Number of moles to set. Unit: kmol. % Number of moles to set. Unit: kmol.
checklib;
calllib(ct, 'mix_setPhaseMoles', m.mixindex, n-1, moles); calllib(ct, 'mix_setPhaseMoles', m.mixID, n-1, moles);
end end
function setSpeciesMoles(m, moles) function setSpeciesMoles(m, moles)
@ -256,10 +249,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.
% %
% parameter moles: % :parameter moles:
% Vector or string specifying the moles of species. % Vector or string specifying the moles of species.
checklib;
calllib(ct, 'mix_setMolesByName', m.mixindex, moles); calllib(ct, 'mix_setMolesByName', m.mixID, moles);
% check back on this one! % check back on this one!
end end
@ -281,29 +274,27 @@ classdef Mixture < handle
% >> equilibrate(mix, 'TP); % >> equilibrate(mix, 'TP);
% >> equilibrate('TP', 1.0e-6, 500); % >> equilibrate('TP', 1.0e-6, 500);
% %
% parameter 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'.
% parameter 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.
% parameter 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.
% parameter 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.
% parameter 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;
if nargin < 6 if nargin < 6
loglevel = 0; loglevel = 0;
end end
@ -319,7 +310,7 @@ classdef Mixture < handle
if nargin < 2 if nargin < 2
XY = 'TP' XY = 'TP'
end end
r = calllib(ct, 'mix_equilibrate', m.mixindex, XY, err, ... r = calllib(ct, 'mix_equilibrate', m.mixID, XY, err, ...
maxsteps, maxiter, loglevel); maxsteps, maxiter, loglevel);
end end

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

@ -20,14 +20,14 @@ classdef Solution < handle & ThermoPhase & Kinetics & Transport
trans = 'default'; trans = 'default';
end end
s@Transport(tp, trans, 0); s@Transport(tp, trans, 0);
s.tp_id = tp.tp_id; s.tpID = tp.tpID;
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.tp_clear; s.tpClear;
s.kin_clear; s.kinClear;
s.tr_clear; s.trClear;
end end
end end

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

@ -1,8 +1,7 @@
classdef ThermoPhase < handle classdef ThermoPhase < handle
properties properties
tp_owner tpID
tp_id
T % temperature T % temperature
P % pressure P % pressure
D % density D % density
@ -72,8 +71,7 @@ classdef ThermoPhase < handle
if nargin == 1 if nargin == 1
id = '-'; id = '-';
end end
tp.tp_owner = 1; tp.tpID = calllib(ct, 'thermo_newFromFile', src, id);
tp.tp_id = calllib(ct, 'thermo_newFromFile', src, id);
tp.basis = 'molar'; tp.basis = 'molar';
end end
@ -85,14 +83,13 @@ classdef ThermoPhase < handle
if nargin < 2 || ~isnumeric(threshold) if nargin < 2 || ~isnumeric(threshold)
threshold = 1e-14; threshold = 1e-14;
end end
calllib(ct, 'thermo_print', tp.tp_id, 1, threshold); calllib(ct, 'thermo_print', tp.tpID, 1, threshold);
end end
function tp_clear(tp) function tpClear(tp)
% Delete the kernel object. % Delete the kernel object.
checklib; calllib(ct, 'thermo_del', tp.tpID);
calllib(ct, 'thermo_del', tp.tp_id);
end end
function tp = set.basis(tp, b) function tp = set.basis(tp, b)
@ -110,33 +107,31 @@ classdef ThermoPhase < handle
%% PhaseGet single methods %% PhaseGet single methods
function amu = atomicMasses(tp) function amu = atomicWeights(tp)
% Get the atomic masses of the elements. % Get the atomic masses of the elements.
% %
% return: % :return:
% Vector of element atomic masses. Unit: kg/kmol % Vector of element atomic masses. Unit: kg/kmol
checklib;
nel = tp.nElements; nel = tp.nElements;
aa = zeros(1, nel); aa = zeros(1, nel);
pt = libpointer('doublePtr', aa); pt = libpointer('doublePtr', aa);
calllib(ct, 'thermo_getAtomicWeights', ... calllib(ct, 'thermo_getAtomicWeights', ...
tp.tp_id, nel, pt); tp.tpID, nel, pt);
amu = pt.Value; amu = pt.Value;
end end
function e = charges(tp) function e = charges(tp)
% Get the array of species charges. % Get the array of species charges.
% %
% return: % :return:
% Vector of species charges. Unit: elem. charge % Vector of species charges. Unit: elem. charge
checklib;
nsp = tp.nSpecies; nsp = tp.nSpecies;
yy = zeros(1, nsp); yy = zeros(1, nsp);
pt = libpointer('doublePtr', yy); pt = libpointer('doublePtr', yy);
calllib(ct, 'thermo_getCharges', ... calllib(ct, 'thermo_getCharges', ...
tp.tp_id, nsp, pt); tp.tpID, nsp, pt);
e = pt.Value; e = pt.Value;
end end
@ -149,19 +144,18 @@ classdef ThermoPhase < 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.
% %
% parameter name: % :parameter name:
% String or cell array of elements whose index is requested % String or cell array of elements whose index is requested
% return: % :return:
% Integer number of elements in the phase. % Integer number of elements in the phase.
checklib;
if iscell(name) if iscell(name)
[m, n] = size(name); [m, n] = size(name);
k = zeros(m, n); k = zeros(m, n);
for i = 1:m for i = 1:m
for j = 1:n for j = 1:n
k(i, j) = calllib(ct, 'thermo_elementIndex', ... k(i, j) = calllib(ct, 'thermo_elementIndex', ...
tp.tp_id, name{i, j}) + 1; tp.tpID, name{i, j}) + 1;
if k(i, j) > 1e3 if k(i, j) > 1e3
warning(['Element ', name{i, j}, ... warning(['Element ', name{i, j}, ...
' does not exist in the phase']); ' does not exist in the phase']);
@ -171,7 +165,7 @@ classdef ThermoPhase < handle
end end
elseif ischar(name) elseif ischar(name)
k = calllib(ct, 'thermo_elementIndex', ... k = calllib(ct, 'thermo_elementIndex', ...
tp.tp_id, name) + 1; tp.tpID, name) + 1;
if k > 1e3 if k > 1e3
warning(['Element ', name, ... warning(['Element ', name, ...
' does not exist in the phase']); ' does not exist in the phase']);
@ -184,9 +178,7 @@ classdef ThermoPhase < handle
function elMassFrac = elementalMassFraction(tp, element) function elMassFrac = elementalMassFraction(tp, element)
% Determine the elemental mass fraction in gas object. % Determine the elemental mass fraction in gas object.
checklib; % Check input :parameters.
% Check input parameters.
if nargin ~= 2 if nargin ~= 2
error('elementalMassFraction expects two input arguments.'); error('elementalMassFraction expects two input arguments.');
end end
@ -226,43 +218,42 @@ classdef ThermoPhase < handle
function mmw = meanMolecularWeight(tp) function mmw = meanMolecularWeight(tp)
% Get the mean molecular weight. % Get the mean molecular weight.
% %
% return: % :return:
% Double mean molecular weight. Unit: kg/kmol % Double mean molecular weight. Unit: kg/kmol
checklib; mmw = calllib(ct, 'thermo_meanMolecularWeight', tp.tpID);
mmw = calllib(ct, 'thermo_meanMolecularWeight', tp.tp_id);
end end
function density = molarDensity(tp) function density = molarDensity(tp)
% Get the molar basis density in kmol/m^3. % Get the molar basis density in kmol/m^3.
checklib;
density = calllib(ct, 'thermo_molarDensity', tp.tp_id); density = calllib(ct, 'thermo_molarDensity', tp.tpID);
end end
function mw = MolecularWeights(tp) function mw = MolecularWeights(tp)
% Get the array of molecular weights of all species. % Get the array of molecular weights of all species.
% %
% return: % :return:
% Vector of species molecular weights. Unit: kg/kmol % Vector of species molecular weights. Unit: kg/kmol
checklib;
nsp = tp.nSpecies; nsp = tp.nSpecies;
yy = zeros(1, nsp); yy = zeros(1, nsp);
pt = libpointer('doublePtr', yy); pt = libpointer('doublePtr', yy);
calllib(ct, 'thermo_getMolecularWeights', ... calllib(ct, 'thermo_getMolecularWeights', ...
tp.tp_id, nsp, pt); tp.tpID, nsp, pt);
mw = pt.Value; mw = pt.Value;
end end
function n = nAtoms(tp, species, element) function n = nAtoms(tp, species, element)
% Get the number of atoms of an element in a species. % Get the number of atoms of an element in a species.
% %
% parameter k: % :parameter k:
% String species name or integer species number. % String species name or integer species number.
% parameter m: % :parameter m:
% String element name or integer element number. % String element name or integer element number.
% return: % :return:
% Integer number of atoms of the element in the species. % Integer number of atoms of the element in the species.
if nargin == 3 if nargin == 3
if ischar(species) if ischar(species)
k = tp.speciesIndex(species); k = tp.speciesIndex(species);
@ -280,7 +271,7 @@ classdef ThermoPhase < handle
n = -1; n = -1;
return return
end end
n = calllib(ct, 'thermo_nAtoms', tp.tp_id, k-1, m-1); n = calllib(ct, 'thermo_nAtoms', tp.tpID, k-1, m-1);
else else
error('Two input arguments required.') error('Two input arguments required.')
end end
@ -289,21 +280,19 @@ classdef ThermoPhase < handle
function nel = nElements(tp) function nel = nElements(tp)
% Get the number of elements. % Get the number of elements.
% %
% return: % :return:
% Integer number of elements in the phase. % Integer number of elements in the phase.
checklib; nel = calllib(ct, 'thermo_nElements', tp.tpID);
nel = calllib(ct, 'thermo_nElements', tp.tp_id);
end end
function nsp = nSpecies(tp) function nsp = nSpecies(tp)
% Get the number of species. % Get the number of species.
% %
% return: % :return:
% Integer number of species in the phase. % Integer number of species in the phase.
checklib; nsp = calllib(ct, 'thermo_nSpecies', tp.tpID);
nsp = calllib(ct, 'thermo_nSpecies', tp.tp_id);
end end
function k = speciesIndex(tp, name) function k = speciesIndex(tp, name)
@ -315,19 +304,18 @@ classdef ThermoPhase < 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.
% %
% parameter name: % :parameter name:
% String or cell array of species whose index is requested. % String or cell array of species whose index is requested.
% return: % :return:
% Integer number of species in the phase. % Integer number of species in the phase.
checklib;
if iscell(name) if iscell(name)
[m, n] = size(name); [m, n] = size(name);
k = zeros(m, n); k = zeros(m, n);
for i = 1:m for i = 1:m
for j = 1:n for j = 1:n
k(i, j) = calllib(ct, 'thermo_speciesIndex', ... k(i, j) = calllib(ct, 'thermo_speciesIndex', ...
tp.tp_id, name{i, j}) + 1; tp.tpID, name{i, j}) + 1;
if k(i, j) > 1e3 if k(i, j) > 1e3
warning(['Species ', name{i, j}, ... warning(['Species ', name{i, j}, ...
' does not exist in the phase']); ' does not exist in the phase']);
@ -337,7 +325,7 @@ classdef ThermoPhase < handle
end end
elseif ischar(name) elseif ischar(name)
k = calllib(ct, 'thermo_speciesIndex', ... k = calllib(ct, 'thermo_speciesIndex', ...
tp.tp_id, name) + 1; tp.tpID, name) + 1;
if k > 1e3 if k > 1e3
warning(['Species ', name, ... warning(['Species ', name, ...
' does not exist in the phase.']); ' does not exist in the phase.']);
@ -351,21 +339,22 @@ classdef ThermoPhase < handle
function nm = speciesName(tp, k) function nm = speciesName(tp, k)
% Get the name of a species given the index. % Get the name of a species given the index.
% %
% parameter k: % :parameter k:
% Scalar or array of integer species index. % Scalar or array of integer species index.
% return: % :return:
% Cell array of strings species name. % Cell array of strings species name.
[m, n] = size(k); [m, n] = size(k);
nm = cell(m, n); nm = cell(m, n);
for i = 1:m for i = 1:m
for j = 1:n for j = 1:n
ksp = k(i, j) - 1; ksp = k(i, j) - 1;
buflen = calllib(ct, 'thermo_getSpeciesName', ... buflen = calllib(ct, 'thermo_getSpeciesName', ...
tp.tp_id, ksp, 0, ''); tp.tpID, ksp, 0, '');
if buflen > 0 if buflen > 0
aa = char(zeros(1, buflen)); aa = char(zeros(1, buflen));
[~, aa] = calllib(ct, 'thermo_getSpeciesName', ... [~, aa] = calllib(ct, 'thermo_getSpeciesName', ...
tp.tp_id, ksp, buflen, aa); tp.tpID, ksp, buflen, aa);
nm{i, j} = aa; nm{i, j} = aa;
end end
end end
@ -374,56 +363,55 @@ classdef ThermoPhase < handle
function n = speciesNames(tp) function n = speciesNames(tp)
% Get all species names. % Get all species names.
n = tp.speciesName(1:tp.nSpecies); n = tp.speciesName(1:tp.nSpecies);
end end
function temperature = get.T(tp) function temperature = get.T(tp)
% Get the temperature. % Get the temperature.
% %
% return: % :return:
% Double temperature. Unit: K % Double temperature. Unit: K
checklib; temperature = calllib(ct, 'thermo_temperature', tp.tpID);
temperature = calllib(ct, 'thermo_temperature', tp.tp_id);
end end
function pressure = get.P(tp) function pressure = get.P(tp)
% Get the pressure. % Get the pressure.
% %
% return: % :return:
% Double pressure. Unit: Pa % Double pressure. Unit: Pa
checklib; pressure = calllib(ct, 'thermo_pressure', tp.tpID);
pressure = calllib(ct, 'thermo_pressure', tp.tp_id);
end end
function density = get.D(tp) function density = get.D(tp)
% Get the mass basis density in kg/m^3. % Get the mass basis density in kg/m^3.
checklib;
density = calllib(ct, 'thermo_density', tp.tp_id); density = calllib(ct, 'thermo_density', tp.tpID);
end end
function volume = get.V(tp) function volume = get.V(tp)
% Get the specific volume depending on the basis. % Get the specific volume depending on the basis.
% %
% return: % :return:
% Density depending on the basis. Units: % Density depending on the basis. Units:
% m^3/kmol (molar) m^3/kg (mass). % m^3/kmol (molar) m^3/kg (mass).
volume = 1/tp.D; volume = 1/tp.D;
end end
function moleFractions = get.X(tp) function moleFractions = get.X(tp)
% Get the mole fractions of all species. % Get the mole fractions of all species.
% %
% return: % :return:
% Vector of species mole fractions. % Vector of species mole fractions.
checklib;
nsp = tp.nSpecies; nsp = tp.nSpecies;
xx = zeros(1, nsp); xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(ct, 'thermo_getMoleFractions', ... calllib(ct, 'thermo_getMoleFractions', ...
tp.tp_id, nsp, pt); tp.tpID, nsp, pt);
moleFractions = pt.Value; moleFractions = pt.Value;
% if no output argument is specified, a bar plot is produced. % if no output argument is specified, a bar plot is produced.
@ -440,10 +428,10 @@ classdef ThermoPhase < handle
function x = moleFraction(tp, species) function x = moleFraction(tp, species)
% Get the mole fraction of one or a list of species. % Get the mole fraction of one or a list of species.
% %
% parameter species: % :parameter species:
% String or cell array of species whose mole fraction is % String or cell array of species whose mole fraction is
% requested. % requested.
% return: % :return:
% Scalar or vector of species mole fractions. % Scalar or vector of species mole fractions.
xarray = tp.X; xarray = tp.X;
@ -469,15 +457,14 @@ classdef ThermoPhase < handle
function massFractions = get.Y(tp) function massFractions = get.Y(tp)
% Get the mass fractions of all species. % Get the mass fractions of all species.
% %
% return: % :return:
% Vector of species mass fractions. % Vector of species mass fractions.
checklib;
nsp = tp.nSpecies; nsp = tp.nSpecies;
yy = zeros(1, nsp); yy = zeros(1, nsp);
pt = libpointer('doublePtr', yy); pt = libpointer('doublePtr', yy);
calllib(ct, 'thermo_getMassFractions', ... calllib(ct, 'thermo_getMassFractions', ...
tp.tp_id, nsp, pt); tp.tpID, nsp, pt);
massFractions = pt.Value; massFractions = pt.Value;
% If no output argument is specified, a bar plot is produced. % If no output argument is specified, a bar plot is produced.
@ -494,10 +481,10 @@ classdef ThermoPhase < handle
function y = massFraction(tp, species) function y = massFraction(tp, species)
% Get the mass fraction of one or a list of species. % Get the mass fraction of one or a list of species.
% %
% parameter species: % :parameter species:
% String or cell array of species whose mass fraction is % String or cell array of species whose mass fraction is
% requested. % requested.
% return: % :return:
% Scalar or vector of species mass fractions. % Scalar or vector of species mass fractions.
yy = tp.Y; yy = tp.Y;
@ -525,30 +512,28 @@ classdef ThermoPhase < handle
function mu = chemical_potentials(tp) function mu = chemical_potentials(tp)
% Get the chemical potentials of the species. % Get the chemical potentials of the species.
% %
% return: % :return:
% Vector of species chemical potentials. Unit: J/kmol. % Vector of species chemical potentials. Unit: J/kmol.
checklib;
nsp = tp.nSpecies; nsp = tp.nSpecies;
xx = zeros(1, nsp); xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(ct, 'thermo_chemPotentials', ... calllib(ct, 'thermo_chemPotentials', ...
tp.tp_id, nsp, pt); tp.tpID, nsp, pt);
mu = pt.Value; mu = pt.Value;
end end
function c = cv(tp) function c = cv(tp)
% Get the specific heat at constant volume. % Get the specific heat at constant volume.
% %
% return: % :return:
% Specific heat of the mixture at constant volume depending % Specific heat of the mixture at constant volume depending
% on the basis. Units: J/kmol-K (molar) J/kg-K (mass). % on the basis. Units: J/kmol-K (molar) J/kg-K (mass).
checklib;
if strcmp(tp.basis, 'molar') if strcmp(tp.basis, 'molar')
c = calllib(ct, 'thermo_cv_mole', tp.tp_id); c = calllib(ct, 'thermo_cv_mole', tp.tpID);
elseif strcmp(tp.basis, 'mass') elseif strcmp(tp.basis, 'mass')
c = calllib(ct, 'thermo_cv_mass', tp.tp_id); c = calllib(ct, 'thermo_cv_mass', tp.tpID);
else error("basis not specified"); else error("basis not specified");
end end
end end
@ -556,15 +541,14 @@ classdef ThermoPhase < handle
function c = cp(tp) function c = cp(tp)
% Get the specific heat at constant pressure. % Get the specific heat at constant pressure.
% %
% return: % :return:
% Specific heat of the mixture at constant pressure depending % Specific heat of the mixture at constant pressure depending
% on the basis. Units: J/kmol-K (molar) J/kg-K (mass). % on the basis. Units: J/kmol-K (molar) J/kg-K (mass).
checklib;
if strcmp(tp.basis, 'molar') if strcmp(tp.basis, 'molar')
c = calllib(ct, 'thermo_cp_mole', tp.tp_id); c = calllib(ct, 'thermo_cp_mole', tp.tpID);
elseif strcmp(tp.basis, 'mass') elseif strcmp(tp.basis, 'mass')
c = calllib(ct, 'thermo_cp_mass', tp.tp_id); c = calllib(ct, 'thermo_cp_mass', tp.tpID);
else error("basis not specified"); else error("basis not specified");
end end
end end
@ -572,53 +556,49 @@ classdef ThermoPhase < handle
function d = critDensity(tp) function d = critDensity(tp)
% Get the critical density. % Get the critical density.
% %
% return: % :return:
% Critical density. Unit: K. % Critical density. Unit: K.
checklib; d = calllib(ct, 'thermo_critDensity', tp.tpID);
d = calllib(ct, 'thermo_critDensity', tp.tp_id);
end end
function p = critPressure(tp) function p = critPressure(tp)
% Get the critical pressure. % Get the critical pressure.
% %
% return: % :return:
% Critical temperature. Unit: Pa. % Critical temperature. Unit: Pa.
checklib; p = calllib(ct, 'thermo_critPressure', tp.tpID);
p = calllib(ct, 'thermo_critPressure', tp.tp_id);
end end
function t = critTemperature(tp) function t = critTemperature(tp)
% Get the critical temperature. % Get the critical temperature.
% %
% return: % :return:
% Critical temperature. Unit: K. % Critical temperature. Unit: K.
checklib;
t = calllib(ct, 'thermo_critTemperature', tp.tp_id); t = calllib(ct, 'thermo_critTemperature', tp.tpID);
end end
function v = electricPotential(tp) function v = electricPotential(tp)
% Get the electric potential % Get the electric potential
% %
% return: % :return:
% Electric potential of the phase. Unit: V. % Electric potential of the phase. Unit: V.
checklib; v = calllib(ct, 'thermo_electricPotential', tp.tpID);
v = calllib(ct, 'thermo_electricPotential', tp.tp_id);
end end
function e = eosType(tp) function e = eosType(tp)
% Get the type of the equation of state % Get the type of the equation of state
% checklib;
% buflen = calllib(ct, 'thermo_getEosType', tp.tp_id, 0, ''); buflen = calllib(ct, 'thermo_getEosType', tp.tpID, 0, '');
% if buflen > 0 if buflen > 0
% aa = char(zeros(1, buflen)); aa = char(zeros(1, buflen));
% [~, aa] = calllib(ct, 'thermo_getEosType', ... [~, aa] = calllib(ct, 'thermo_getEosType', ...
% tp.tp_id, buflen, aa); tp.tpID, buflen, aa);
% end end
% e = aa; e = aa;
e = 'IdealGas';
end end
function v = isIdealGas(tp) function v = isIdealGas(tp)
@ -635,53 +615,48 @@ classdef ThermoPhase < handle
function b = isothermalCompressibility(tp) function b = isothermalCompressibility(tp)
% Get the isothermal compressibility % Get the isothermal compressibility
% %
% return: % :return:
% Isothermal compressibility. Unit: 1/Pa. % Isothermal compressibility. Unit: 1/Pa.
checklib; b = calllib(ct, 'thermo_isothermalCompressibility', tp.tpID);
b = calllib(ct, 'thermo_isothermalCompressibility', tp.tp_id);
end end
function t = maxTemp(tp) function t = maxTemp(tp)
% Get the maximum temperature of the parameter fits. % Get the maximum temperature of the :parameter fits.
% %
% return: % :return:
% Vector of maximum temperatures of all species. % Vector of maximum temperatures of all species.
checklib; t = calllib(ct, 'thermo_maxTemp', tp.tpID, -1);
t = calllib(ct, 'thermo_maxTemp', tp.tp_id, -1);
end end
function t = minTemp(tp) function t = minTemp(tp)
% Get the minimum temperature of the parameter fits. % Get the minimum temperature of the :parameter fits.
% %
% return: % :return:
% Vector of minimum temperatures of all species. % Vector of minimum temperatures of all species.
checklib; t = calllib(ct, 'thermo_minTemp', tp.tpID, -1);
t = calllib(ct, 'thermo_minTemp', tp.tp_id, -1);
end end
function p = P_sat(tp, t) function p = P_sat(tp, t)
% Get the saturation pressure for a given temperature. % Get the saturation pressure for a given temperature.
% %
% parameter t: % :parameter t:
% Temperature. Unit: K. % Temperature. Unit: K.
% return: % :return:
% Saturation pressure for temperature t. Unit: Pa. % Saturation pressure for temperature t. Unit: Pa.
checklib; p = calllib(ct, 'thermo_satPressure', tp.tpID, t);
p = calllib(ct, 'thermo_satPressure', tp.tp_id, t);
end end
function p = refPressure(tp) function p = refPressure(tp)
% Get the reference pressure. % Get the reference pressure.
% %
% return: % :return:
% Reference pressure. Unit: Pa. % Reference pressure. Unit: Pa.
checklib; p = calllib(ct, 'thermo_refPressure', tp.tpID, -1);
p = calllib(ct, 'thermo_refPressure', tp.tp_id, -1);
end end
function c = soundspeed(tp) function c = soundspeed(tp)
@ -697,10 +672,9 @@ classdef ThermoPhase < handle
% %
% c = sqrt( % c = sqrt(
% %
% return: % :return:
% The speed of sound. Unit: m/s % The speed of sound. Unit: m/s
checklib;
if tp.isIdealGas if tp.isIdealGas
tp.basis = 'mass'; tp.basis = 'mass';
gamma = tp.cp/tp.cv; gamma = tp.cp/tp.cv;
@ -722,47 +696,43 @@ classdef ThermoPhase < handle
function a = thermalExpansionCoeff(tp) function a = thermalExpansionCoeff(tp)
% Get the thermal expansion coefficient. % Get the thermal expansion coefficient.
% %
% return: % :return:
% Thermal expansion coefficient. Unit: 1/K. % Thermal expansion coefficient. Unit: 1/K.
checklib; a = calllib(ct, 'thermo_thermalExpansionCoeff', tp.tpID);
a = calllib(ct, 'thermo_thermalExpansionCoeff', tp.tp_id);
end end
function t = T_sat(tp, p) function t = T_sat(tp, p)
% Get the saturation temperature for a given pressure. % Get the saturation temperature for a given pressure.
% %
% parameter p: % :parameter p:
% Pressure. Unit: Pa. % Pressure. Unit: Pa.
% return: % :return:
% Saturation temperature for pressure p. Unit: K. % Saturation temperature for pressure p. Unit: K.
checklib; t = calllib(ct, 'thermo_satTemperature', tp.tpID, p);
t = calllib(ct, 'thermo_satTemperature', tp.tp_id, p);
end end
function v = vaporFraction(tp) function v = vaporFraction(tp)
% Get the vapor fractions. % Get the vapor fractions.
% %
% return: % :return:
% Vapor fraction. % Vapor fraction.
checklib; v = calllib(ct, 'thermo_vaporFraction', tp.tpID);
v = calllib(ct, 'thermo_vaporFraction', tp.tp_id);
end end
function enthalpy = get.H(tp) function enthalpy = get.H(tp)
% Get the enthalpy. % Get the enthalpy.
% %
% return: % :return:
% Enthalpy of the mixture depending on the basis. % Enthalpy of the mixture depending on the basis.
% Units: J/kmol (molar) J/kg (mass). % Units: J/kmol (molar) J/kg (mass).
checklib;
if strcmp(tp.basis, 'molar') if strcmp(tp.basis, 'molar')
enthalpy = calllib(ct, 'thermo_enthalpy_mole', tp.tp_id); enthalpy = calllib(ct, 'thermo_enthalpy_mole', tp.tpID);
elseif strcmp(tp.basis, 'mass') elseif strcmp(tp.basis, 'mass')
enthalpy = calllib(ct, 'thermo_enthalpy_mass', tp.tp_id); enthalpy = calllib(ct, 'thermo_enthalpy_mass', tp.tpID);
else error("basis not specified"); else error("basis not specified");
end end
end end
@ -770,29 +740,27 @@ classdef ThermoPhase < handle
function enthalpy = enthalpies_RT(tp) function enthalpy = enthalpies_RT(tp)
% Get the non-dimensional enthalpy. % Get the non-dimensional enthalpy.
% %
% return: % :return:
% Vector of standard-state species enthalpies divided by RT. % Vector of standard-state species enthalpies divided by RT.
checklib;
nsp = tp.nSpecies; nsp = tp.nSpecies;
xx = zeros(1, nsp); xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx); pt = libpointer('doublePtr', xx);
calllib(ct, 'thermo_getEnthalpies_RT', tp.tp_id, nsp, pt); calllib(ct, 'thermo_getEnthalpies_RT', tp.tpID, nsp, pt);
enthalpy = pt.Value; enthalpy = pt.Value;
end end
function entropy = get.S(tp) function entropy = get.S(tp)
% Get the entropy. % Get the entropy.
% %
% return: % :return:
% Entropy of the mixture depending on the basis. % Entropy of the mixture depending on the basis.
% Units: J/kmol-K (molar) J/kg-K (mass). % Units: J/kmol-K (molar) J/kg-K (mass).
checklib;
if strcmp(tp.basis, 'molar') if strcmp(tp.basis, 'molar')
entropy = calllib(ct, 'thermo_entropy_mole', tp.tp_id); entropy = calllib(ct, 'thermo_entropy_mole', tp.tpID);
elseif strcmp(tp.basis, 'mass') elseif strcmp(tp.basis, 'mass')
entropy = calllib(ct, 'thermo_entropy_mass', tp.tp_id); entropy = calllib(ct, 'thermo_entropy_mass', tp.tpID);
else error("basis not specified"); else error("basis not specified");
end end
end end
@ -800,15 +768,14 @@ classdef ThermoPhase < handle
function intEnergy = get.U(tp) function intEnergy = get.U(tp)
% Get the internal energy. % Get the internal energy.
% %
% return: % :return:
% Internal energy of the mixture depending on the basis. % Internal energy of the mixture depending on the basis.
% Units: J/kmol (molar) J/kg (mass). % Units: J/kmol (molar) J/kg (mass).
checklib;
if strcmp(tp.basis, 'molar') if strcmp(tp.basis, 'molar')
intEnergy = calllib(ct, 'thermo_intEnergy_mole', tp.tp_id); intEnergy = calllib(ct, 'thermo_intEnergy_mole', tp.tpID);
elseif strcmp(tp.basis, 'mass') elseif strcmp(tp.basis, 'mass')
intEnergy = calllib(ct, 'thermo_intEnergy_mass', tp.tp_id); intEnergy = calllib(ct, 'thermo_intEnergy_mass', tp.tpID);
else error("basis not specified"); else error("basis not specified");
end end
end end
@ -816,15 +783,14 @@ classdef ThermoPhase < handle
function gibbs = get.G(tp) function gibbs = get.G(tp)
% Get the Gibss free energy. % Get the Gibss free energy.
% %
% return: % :return:
% Gibbs free energy of the mixture depending on the basis. % Gibbs free energy of the mixture depending on the basis.
% Units: J/kmol (molar) J/kg (mass). % Units: J/kmol (molar) J/kg (mass).
checklib;
if strcmp(tp.basis, 'molar') if strcmp(tp.basis, 'molar')
gibbs = calllib(ct, 'thermo_gibbs_mole', tp.tp_id); gibbs = calllib(ct, 'thermo_gibbs_mole', tp.tpID);
elseif strcmp(tp.basis, 'mass') elseif strcmp(tp.basis, 'mass')
gibbs = calllib(ct, 'thermo_gibbs_mass', tp.tp_id); gibbs = calllib(ct, 'thermo_gibbs_mass', tp.tpID);
else error("basis not specified"); else error("basis not specified");
end end
end end
@ -836,12 +802,6 @@ classdef ThermoPhase < handle
end end
function output = get.DPX(tp) function output = get.DPX(tp)
% Get density, pressure, and mole fractions depending on the basis.
% return:
% Density. Unit: kmol/m^3 (molar) kg/m^3 (mass).
% Pressure. Unit: Pa.
% Mole fractions of all species.
output = {tp.D, tp.P, tp.X}; output = {tp.D, tp.P, tp.X};
end end
@ -854,12 +814,6 @@ classdef ThermoPhase < handle
end end
function output = get.HPX(tp) function output = get.HPX(tp)
% Get enthalpy, pressure, and mole fractions depending on the basis.
% return:
% Enthalpy. Unit: J/kmol (molar) J/kg (mass).
% Pressure. Unit: Pa.
% Mole fractions of all species.
output = {tp.H, tp.P, tp.X}; output = {tp.H, tp.P, tp.X};
end end
@ -1015,47 +969,44 @@ classdef ThermoPhase < handle
function tp = setElectricPotential(tp, phi) function tp = setElectricPotential(tp, phi)
% Set the electric potential in V. % Set the electric potential in V.
calllib(ct, 'thermo_setElectricPotential', tp.tp_id, phi);
calllib(ct, 'thermo_setElectricPotential', tp.tpID, phi);
end end
function tp = setState_Psat(tp, p, q) function tp = setState_Psat(tp, p, q)
% Set saturated vapor % Set saturated vapor
checklib;
calllib(ct, 'thermo_setState_Psat', tp.tp_id, p, q); calllib(ct, 'thermo_setState_Psat', tp.tpID, p, q);
end end
function tp = setState_Tsat(tp, t, q) function tp = setState_Tsat(tp, t, q)
% Set saturated liquid % Set saturated liquid
checklib;
calllib(ct, 'thermo_setState_Tsat', tp.tp_id, t, 1 - q); calllib(ct, 'thermo_setState_Tsat', tp.tpID, t, 1 - q);
end end
function set.T(tp, temperature) function set.T(tp, temperature)
checklib;
if temperature <= 0 if temperature <= 0
error('The temperature must be positive'); error('The temperature must be positive');
end end
calllib(ct, 'thermo_setTemperature', tp.tp_id, temperature); calllib(ct, 'thermo_setTemperature', tp.tpID, temperature);
end end
function set.P(tp, pressure) function set.P(tp, pressure)
checklib;
if pressure <= 0 if pressure <= 0
error('The pressure must be positive'); error('The pressure must be positive');
end end
calllib(ct, 'thermo_setPressure', tp.tp_id, pressure); calllib(ct, 'thermo_setPressure', tp.tpID, pressure);
end end
function set.D(tp, density) function set.D(tp, density)
checklib;
if density <= 0 if density <= 0
error('The density must be positive'); error('The density must be positive');
end end
calllib(ct, 'thermo_setDensity', tp.tp_id, density); calllib(ct, 'thermo_setDensity', tp.tpID, density);
end end
function set.X(tp, xx) function set.X(tp, xx)
checklib;
lim = 1e-9; lim = 1e-9;
if isa(xx, 'double') if isa(xx, 'double')
nsp = tp.nSpecies; nsp = tp.nSpecies;
@ -1063,32 +1014,30 @@ classdef ThermoPhase < handle
norm = 0; norm = 0;
else norm = 1; else norm = 1;
end end
calllib(ct, 'thermo_setMoleFractions', tp.tp_id, ... calllib(ct, 'thermo_setMoleFractions', tp.tpID, ...
nsp, xx, norm); nsp, xx, norm);
elseif isa(xx, 'char') elseif isa(xx, 'char')
calllib(ct, 'thermo_setMoleFractionsByName', tp.tp_id, xx); calllib(ct, 'thermo_setMoleFractionsByName', tp.tpID, xx);
end end
end end
function set.Y(tp, yy) function set.Y(tp, yy)
checklib;
if isa(yy, 'double') if isa(yy, 'double')
nsp = tp.nSpecies; nsp = tp.nSpecies;
if sum(yy) -1 <= 1e-9 if sum(yy) -1 <= 1e-9
norm = 0; norm = 0;
else norm = 1; else norm = 1;
end end
calllib(ct, 'thermo_setMassFractions', tp.tp_id, ... calllib(ct, 'thermo_setMassFractions', tp.tpID, ...
nsp, yy, norm); nsp, yy, norm);
elseif isa(yy, 'char') elseif isa(yy, 'char')
calllib(ct, 'thermo_setMassFractionsByName', tp.tp_id, yy); calllib(ct, 'thermo_setMassFractionsByName', tp.tpID, yy);
end end
end end
%% PhaseSet multi methods %% PhaseSet multi methods
function set.DP(tp, input) function set.DP(tp, input)
checklib;
d = input{1}; d = input{1};
p = input{2}; p = input{2};
if d <= 0 if d <= 0
@ -1097,7 +1046,7 @@ classdef ThermoPhase < handle
if p <= 0 if p <= 0
error('The pressure must be positive'); error('The pressure must be positive');
end end
calllib(ct, 'thermo_set_RP', tp.tp_id, [d, p]); calllib(ct, 'thermo_set_RP', tp.tpID, [d, p]);
end end
function set.DPX(tp, input) function set.DPX(tp, input)
@ -1111,13 +1060,12 @@ classdef ThermoPhase < handle
end end
function set.HP(tp, input) function set.HP(tp, input)
checklib;
h = input{1}; h = input{1};
p = input{2}; p = input{2};
if p <= 0 if p <= 0
error('The pressure must be positive'); error('The pressure must be positive');
end end
calllib(ct, 'thermo_set_HP', tp.tp_id, [h, p]); calllib(ct, 'thermo_set_HP', tp.tpID, [h, p]);
end end
function set.HPX(tp, input) function set.HPX(tp, input)
@ -1131,7 +1079,6 @@ classdef ThermoPhase < handle
end end
function set.PV(tp, input) function set.PV(tp, input)
checklib;
p = input{1}; p = input{1};
v = input{2}; v = input{2};
if p <= 0 if p <= 0
@ -1140,7 +1087,7 @@ classdef ThermoPhase < handle
if v <= 0 if v <= 0
error('The specific volume must be positive'); error('The specific volume must be positive');
end end
calllib(ct, 'thermo_set_PV', tp.tp_id, [p, v]); calllib(ct, 'thermo_set_PV', tp.tpID, [p, v]);
end end
function set.PVX(tp, input) function set.PVX(tp, input)
@ -1154,10 +1101,9 @@ classdef ThermoPhase < handle
end end
function set.SH(tp, input) function set.SH(tp, input)
checklib;
s = input{1}; s = input{1};
h = input{2}; h = input{2};
calllib(ct, 'thermo_set_SH', tp.tp_id, [s, h]); calllib(ct, 'thermo_set_SH', tp.tpID, [s, h]);
end end
function set.SHX(tp, input) function set.SHX(tp, input)
@ -1171,13 +1117,12 @@ classdef ThermoPhase < handle
end end
function set.SP(tp, input) function set.SP(tp, input)
checklib;
s = input{1}; s = input{1};
p = input{2}; p = input{2};
if p <= 0 if p <= 0
error('The pressure must be positive'); error('The pressure must be positive');
end end
calllib(ct, 'thermo_set_SP', tp.tp_id, [s, p]); calllib(ct, 'thermo_set_SP', tp.tpID, [s, p]);
end end
function set.SPX(tp, input) function set.SPX(tp, input)
@ -1191,13 +1136,12 @@ classdef ThermoPhase < handle
end end
function set.ST(tp, input) function set.ST(tp, input)
checklib;
s = input{1}; s = input{1};
t = input{2}; t = input{2};
if t <= 0 if t <= 0
error('The temperature must be positive'); error('The temperature must be positive');
end end
calllib(ct, 'thermo_set_ST', tp.tp_id, [s, t]); calllib(ct, 'thermo_set_ST', tp.tpID, [s, t]);
end end
function set.STX(tp, input) function set.STX(tp, input)
@ -1211,13 +1155,12 @@ classdef ThermoPhase < handle
end end
function set.SV(tp, input) function set.SV(tp, input)
checklib;
s = input{1}; s = input{1};
v = input{2}; v = input{2};
if v <= 0 if v <= 0
error('The specific volume must be positive'); error('The specific volume must be positive');
end end
calllib(ct, 'thermo_set_SV', tp.tp_id, [s, v]); calllib(ct, 'thermo_set_SV', tp.tpID, [s, v]);
end end
function set.SVX(tp, input) function set.SVX(tp, input)
@ -1254,13 +1197,12 @@ classdef ThermoPhase < handle
end end
function set.TH(tp, input) function set.TH(tp, input)
checklib;
t = input{1}; t = input{1};
if t <= 0 if t <= 0
error('The temperature must be positive'); error('The temperature must be positive');
end end
h = input{2}; h = input{2};
calllib(ct, 'thermo_set_TH', tp.tp_id, [t, h]); calllib(ct, 'thermo_set_TH', tp.tpID, [t, h]);
end end
function set.THX(tp, input) function set.THX(tp, input)
@ -1297,7 +1239,6 @@ classdef ThermoPhase < handle
end end
function set.TV(tp, input) function set.TV(tp, input)
checklib;
t = input{1}; t = input{1};
v = input{2}; v = input{2};
if t <= 0 if t <= 0
@ -1306,7 +1247,7 @@ classdef ThermoPhase < handle
if v <= 0 if v <= 0
error('The specific volume must be positive'); error('The specific volume must be positive');
end end
calllib(ct, 'thermo_set_TV', tp.tp_id, [t, v]); calllib(ct, 'thermo_set_TV', tp.tpID, [t, v]);
end end
function set.TVX(tp, input) function set.TVX(tp, input)
@ -1320,13 +1261,12 @@ classdef ThermoPhase < handle
end end
function set.UP(tp, input) function set.UP(tp, input)
checklib;
u = input{1}; u = input{1};
p = input{2}; p = input{2};
if p <= 0 if p <= 0
error('The pressure must be positive'); error('The pressure must be positive');
end end
calllib(ct, 'thermo_set_UP', tp.tp_id, [u, p]); calllib(ct, 'thermo_set_UP', tp.tpID, [u, p]);
end end
function set.UPX(tp, input) function set.UPX(tp, input)
@ -1340,13 +1280,12 @@ classdef ThermoPhase < handle
end end
function set.UV(tp, input) function set.UV(tp, input)
checklib;
u = input{1}; u = input{1};
v = input{2}; v = input{2};
if v <= 0 if v <= 0
error('The specific volume must be positive'); error('The specific volume must be positive');
end end
calllib(ct, 'thermo_set_UV', tp.tp_id, [u, v]); calllib(ct, 'thermo_set_UV', tp.tpID, [u, v]);
end end
function set.UVX(tp, input) function set.UVX(tp, input)
@ -1360,13 +1299,12 @@ classdef ThermoPhase < handle
end end
function set.VH(tp, input) function set.VH(tp, input)
checklib;
v = input{1}; v = input{1};
h = input{2}; h = input{2};
if v <= 0 if v <= 0
error('The specific volume must be positive'); error('The specific volume must be positive');
end end
calllib(ct, 'thermo_set_VH', tp.tp_id, [v, h]); calllib(ct, 'thermo_set_VH', tp.tpID, [v, h]);
end end
function set.VHX(tp, input) function set.VHX(tp, input)
@ -1381,7 +1319,6 @@ classdef ThermoPhase < handle
function tp = equilibrate(tp, xy, solver, rtol, maxsteps, ... function tp = equilibrate(tp, xy, solver, rtol, maxsteps, ...
maxiter, loglevel) maxiter, loglevel)
checklib;
% use the ChemEquil solver by default % use the ChemEquil solver by default
if nargin < 3 if nargin < 3
solver = -1; solver = -1;
@ -1398,7 +1335,7 @@ classdef ThermoPhase < handle
if nargin < 7 if nargin < 7
loglevel = 0; loglevel = 0;
end end
calllib(ct, 'thermo_equilibrate', tp.tp_id, xy, solver, ... calllib(ct, 'thermo_equilibrate', tp.tpID, xy, solver, ...
rtol, maxsteps, maxiter, loglevel); rtol, maxsteps, maxiter, loglevel);
end end

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

@ -2,7 +2,7 @@ classdef Transport < handle
properties properties
th th
tr_id trID
end end
methods methods
@ -10,7 +10,7 @@ classdef Transport < handle
function tr = Transport(tp, model, loglevel) function tr = Transport(tp, model, loglevel)
checklib; checklib;
tr.tr_id = 0; tr.trID = 0;
if nargin == 2 if nargin == 2
model = 'default' model = 'default'
end end
@ -25,23 +25,22 @@ classdef Transport < handle
else else
tr.th = tp; tr.th = tp;
if strcmp(model, 'default') if strcmp(model, 'default')
tr.tr_id = calllib(ct, 'trans_newDefault', ... tr.trID = calllib(ct, 'trans_newDefault', ...
tp.tp_id, loglevel); tp.tpID, loglevel);
else else
tr.tr_id = calllib(ct, 'trans_new', model, ... tr.trID = calllib(ct, 'trans_new', model, ...
tp.tp_id, loglevel); tp.tpID, loglevel);
end end
end end
tr.tp_id = tp.tp_id; tr.tpID = tp.tpID;
end end
%% Utility methods %% Utility methods
function tr_clear(tr) function trClear(tr)
% Delete the kernel object. % Delete the kernel object.
checklib; calllib(ct, 'trans_del', tr.trID);
calllib(ct, 'trans_del', tr.tr_id);
end end
%% Transport Methods %% Transport Methods
@ -49,11 +48,10 @@ 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; v = calllib(ct, 'trans_viscosity', tr.trID);
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
@ -64,11 +62,10 @@ classdef Transport < handle
function v = thermalConductivity(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; v = calllib(ct, 'trans_thermalConductivity', tr.trID);
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
@ -79,11 +76,10 @@ 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; v = calllib(ct, 'trans_electricalConductivity', tr.trID);
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
@ -94,81 +90,75 @@ 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.
checklib;
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(ct, 'trans_getMixDiffCoeffs', tr.tr_id, nsp, pt); calllib(ct, 'trans_getMixDiffCoeffs', tr.trID, 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.
checklib;
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(ct, 'trans_getThermalDiffCoeffs', tr.tr_id, nsp, pt); calllib(ct, 'trans_getThermalDiffCoeffs', tr.trID, 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.
checklib;
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(ct, 'trans_getBinDiffCoeffs', tr.tr_id, nsp, pt); calllib(ct, 'trans_getBinDiffCoeffs', tr.trID, 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.
checklib;
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(ct, 'trans_getMultiDiffCoeffs', tr.tr_id, nsp, pt); calllib(ct, 'trans_getMultiDiffCoeffs', tr.trID, 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.
% %
% parameter type: % :parameter type:
% parameter k: % :parameter k:
% parameter p: % :parameter p:
checklib; calllib(ct, 'trans_setParameters', tr.trID, 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.
% %
% parameter lam: % :parameter lam:
% Thermal conductivity in W/(m-K). % Thermal conductivity in W/(m-K).
checklib;
tr.setParameters(1, 0, lam); tr.setParameters(1, 0, lam);
end end

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

@ -24,7 +24,7 @@ classdef Func < handle
% object 'f' is a functor that evaluates the polynomial % object 'f' is a functor that evaluates the polynomial
% :math:'2x^2 - 3x + 1'. Then writing 'f(2)' would cause the % :math:'2x^2 - 3x + 1'. Then writing 'f(2)' would cause the
% method that evaluates the function to be invoked, and would % method that evaluates the function to be invoked, and would
% pass it the argument '2'. The return value would be 3. % pass it the argument '2'. The :return value would be 3.
% %
% The types of functors you can create in Cantera are these: % The types of functors you can create in Cantera are these:
% 1. A polynomial; % 1. A polynomial;
@ -40,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.
% %
% parameter 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'
@ -53,11 +53,11 @@ classdef Func < handle
% * 'composite' % * 'composite'
% * 'periodic' % * 'periodic'
% %
% parameter n: % :parameter n:
% Number of parameters required for the functor % Number of :parameters required for the functor
% parameter 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;
@ -70,8 +70,8 @@ classdef Func < handle
x.f2 = 0; x.f2 = 0;
x.coeffs = 0; x.coeffs = 0;
function nn = new_func(itype, n, p) function nn = newFunc(itype, n, p)
% helper function to pass the correct parameters to the C % helper function to pass the correct :parameters to the C
% library % library
if itype < 20 if itype < 20
ptr = libpointer('doublePtr', p); ptr = libpointer('doublePtr', p);
@ -89,12 +89,12 @@ classdef Func < handle
if itype > 0 if itype > 0
x.coeffs = p; x.coeffs = p;
x.id = new_func(itype, n, p); x.id = newFunc(itype, n, p);
elseif strcmp(typ, 'periodic') elseif strcmp(typ, 'periodic')
itype = 50; itype = 50;
x.f1 = n; x.f1 = n;
x.coeffs = p; x.coeffs = p;
x.id = new_func(itype, n.id, p); x.id = newFunc(itype, n.id, p);
else else
if strcmp(typ, 'sum') if strcmp(typ, 'sum')
itype = 20; itype = 20;
@ -109,7 +109,7 @@ classdef Func < handle
end end
x.f1 = n; x.f1 = n;
x.f2 = p; x.f2 = p;
x.id = new_func(itype, n.id, p.id); x.id = newFunc(itype, n.id, p.id);
end end
x.typ = typ; x.typ = typ;
@ -119,14 +119,14 @@ classdef Func < handle
function clear(f) function clear(f)
% Clear the functor from memory. % Clear the functor from memory.
checklib;
calllib(ct, '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.
% %
% parameter a: % :parameter a:
% Instance of class 'Func' % Instance of class 'Func'
disp(' '); disp(' ');
@ -158,14 +158,13 @@ classdef Func < handle
function b = subsref(a, s) function b = subsref(a, s)
% Redefine subscripted references for functors. % Redefine subscripted references for functors.
% %
% parameter a: % :parameter a:
% Instance of class 'Func' % Instance of class 'Func'
% parameter 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;
if strcmp(s.type, '()') if strcmp(s.type, '()')
ind= s.subs{:}; ind= s.subs{:};
b = zeros(1, length(ind)); b = zeros(1, length(ind));

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

@ -13,22 +13,14 @@ classdef FlowDevice < handle
function x = FlowDevice(typ) function x = FlowDevice(typ)
% Flow Device class constructor. % Flow Device class constructor.
% %
% parameter 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'.
checklib; checklib;
if nargin == 0 if nargin == 0
typ = 'MassFlowController'; error('please specify the type of flow device to be created');
end
if isa(typ, 'double')
warning(['Definition via integer type to be deprecated', ...
' after Cantera 2.5.']);
device_types = {'MassFlowController', 'PressureController', ...
'Valve'};
typ = device_types(typ);
end end
x.type = typ; x.type = typ;
@ -44,7 +36,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;
calllib(ct, 'flowdev_del', f.id); calllib(ct, 'flowdev_del', f.id);
end end
@ -53,12 +45,11 @@ 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.
% %
% parameter upstream: % :parameter upstream:
% Upsteram 'Reactor' or 'Reservoir'. % Upsteram 'Reactor' or 'Reservoir'.
% parameter downstream: % :parameter downstream:
% Downstream 'Reactor' or 'Reservoir'. % Downstream 'Reactor' or 'Reservoir'.
checklib;
if nargin == 3 if nargin == 3
if ~isa(upstream, 'Reactor') || ~isa(downstream, 'Reactor') if ~isa(upstream, 'Reactor') || ~isa(downstream, 'Reactor')
error(['Flow devices can only be installed between',... error(['Flow devices can only be installed between',...
@ -74,32 +65,21 @@ classdef FlowDevice < handle
end end
end end
function mdot = massFlowRate(f, time) function mdot = massFlowRate(f)
% Get the mass flow rate at a given time. % Get the mass flow rate.
% %
% parameter time: % :return:
% Time at which the mass flow rate is desired. % The mass flow rate through the flow device
% return:
% The mass flow rate through the flow device at the given
% time.
checklib; mdot = calllib(ct, 'flowdev_massFlowRate2', f.id);
if nargin == 1
mdot = calllib(ct, 'flowdev_massFlowRate2', f.id);
else
warning(['"Time" argument to massFlowRate is deprecated', ...
'and will be removed after Cantera 2.5.']);
mdot = calllib(ct, 'flowdev_massFlowRate', f.id, time);
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'.
% %
% parameter mf: % :parameter mf:
% Instance of class 'func'. % Instance of class 'func'.
checklib;
if strcmp(f.type, 'MassFlowController') if strcmp(f.type, 'MassFlowController')
k = calllib(ct, 'flowdev_setTimeFunction', f.id, ... k = calllib(ct, 'flowdev_setTimeFunction', f.id, ...
mf.id); mf.id);
@ -107,17 +87,16 @@ classdef FlowDevice < handle
% error(geterr); % error(geterr);
% end % end
else else
error('Mass flow rate can only be set for mass flow controllers.'); error('Time function can only be set for mass flow controllers.');
end end
end end
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.
% %
% parameter mdot: % :parameter mdot:
% Mass flow rate % Mass flow rate
checklib;
if strcmp(f.type, 'MassFlowController') if strcmp(f.type, 'MassFlowController')
k = calllib(ct, 'flowdev_setMassFlowCoeff', f.id, mdot); k = calllib(ct, 'flowdev_setMassFlowCoeff', f.id, mdot);
% if k < 0 % if k < 0
@ -134,12 +113,11 @@ classdef FlowDevice < handle
% The mass flow rate [kg/s] is computed from the expression % The mass flow rate [kg/s] is computed from the expression
% mdot = K(P_upstream - P_downstream) % mdot = K(P_upstream - P_downstream)
% 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.
% %
% parameter k: % :parameter k:
% Value fo the valve coefficient. Unit: kg/Pa-s. % Value fo the valve coefficient. Unit: kg/Pa-s.
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

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

@ -26,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.
% %
% parameter contents: % :parameter contents:
% Instance of class 'Solution' representing the contents of % Instance of class 'Solution' representing the contents of
% the reactor. % the reactor.
% parameter typ: % :parameter typ:
% Character array of reactor type. Options are: % Character array of reactor type. Options are:
% 'Reservoir' % 'Reservoir'
% 'Reactor' % 'Reactor'
@ -37,7 +37,7 @@ classdef Reactor < handle
% 'ConstPressureReactor' % 'ConstPressureReactor'
% 'IdealGasReactor' % 'IdealGasReactor'
% 'IdealGasConstPressureReactor' % 'IdealGasConstPressureReactor'
% return: % :return:
% Instance of class 'Reactor'. % Instance of class 'Reactor'.
checklib; checklib;
@ -70,16 +70,16 @@ classdef Reactor < handle
function clear(r) function clear(r)
% Clear the reactor from memory. % Clear the reactor from memory.
checklib;
calllib(ct, '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.
% %
% parameter r: % :parameter r:
% Instance of class 'Reactor'. % Instance of class 'Reactor'.
% parameter gas: % :parameter gas:
% Instance of class 'Solution'. % Instance of class 'Solution'.
r.contents = gas; r.contents = gas;
@ -95,18 +95,17 @@ 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.
% %
% parameter r: % :parameter r:
% Instance of class 'Reactor'. % Instance of class 'Reactor'.
% parameter 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;
if ~isa(t, 'ThermoPhase') if ~isa(t, 'ThermoPhase')
error('Wrong object type'); error('Wrong object type');
end end
calllib(ct, 'reactor_setThermoMgr', r.id, t.tp_id); calllib(ct, 'reactor_setThermoMgr', r.id, t.tpID);
end end
function setKineticsMgr(r, k) function setKineticsMgr(r, k)
@ -115,18 +114,17 @@ 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.
% %
% parameter r: % :parameter r:
% Instance of class 'Reactor'. % Instance of class 'Reactor'.
% parameter 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;
if ~isa(k, 'Kinetics') if ~isa(k, 'Kinetics')
error('Wrong object type'); error('Wrong object type');
end end
calllib(ct, 'reactor_setKineticsMgr', r.id, k.kin_id); calllib(ct, 'reactor_setKineticsMgr', r.id, k.kinID);
end end
%% Reactor get methods %% Reactor get methods
@ -134,91 +132,82 @@ 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;
temperature = calllib(ct, '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;
pressure = calllib(ct, '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;
rho = calllib(ct, '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;
mass = calllib(ct, '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;
volume = calllib(ct, '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;
enthalpy_mass = calllib(ct, '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;
intEnergy_mass = calllib(ct, '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.
% %
% parameter species: % :parameter species:
% String or one-based integer id of the species. % String or one-based integer id of the species.
checklib;
if ischar(species) if ischar(species)
k = r.contents.speciesIndex(species) - 1; k = r.contents.speciesIndex(species) - 1;
else k = species - 1; else k = species - 1;
@ -230,16 +219,14 @@ classdef Reactor < handle
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'.
checklib;
nsp = r.contents.nSpecies; nsp = r.contents.nSpecies;
y = zeros(1, nsp); massFractions = zeros(1, nsp);
for i = 1:nsp for i = 1:nsp
y(i) = r.massFraction(i); massFractions(i) = r.massFraction(i);
end end
end end
@ -248,20 +235,18 @@ classdef Reactor < handle
function setInitialVolume(r, v0) function setInitialVolume(r, v0)
% Set the initial reactor volume. % Set the initial reactor volume.
% %
% parameter v0: % :parameter v0:
% Initial volume in m^3. % Initial volume in m^3.
checklib;
calllib(ct, '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.
% %
% parameter MFR: % :parameter MFR:
% Mass flow rate in kg/s. % Mass flow rate in kg/s.
checklib;
calllib(ct, 'reactor_setMassFlowRate', r.id, MFR); calllib(ct, 'reactor_setMassFlowRate', r.id, MFR);
r.Mdot = MFR; r.Mdot = MFR;
end end
@ -270,21 +255,19 @@ classdef Reactor < handle
% Enable or disable changing reactor composition by reactions. % Enable or disable changing reactor composition by reactions.
% %
% If the chemistry is disabled, then the reactor composition is % If the chemistry is disabled, then the reactor composition is
% constant. The parameter should be the string "on" to enable % constant. The :parameter should be the string "on" to enable
% the species equaionts, or "off" to disable it. % the species equaionts, or "off" to disable it.
% %
% By default, Reactor objects are created with the species % By default, Reactor objects are created with the species
% 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.
% %
% parameter r: % :parameter r:
% Instance of class 'Reactor'. % Instance of class 'Reactor'.
% parameter 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.
checklib;
if strcmp(flag, 'on') if strcmp(flag, 'on')
iflag = true; iflag = true;
elseif strcmp(flag, 'off') elseif strcmp(flag, 'off')
@ -300,20 +283,18 @@ classdef Reactor < handle
% Enable or disable solving the energy equation. % Enable or disable solving the energy equation.
% %
% If the energy equation is disabled, then the reactor % If the energy equation is disabled, then the reactor
% temperature is constant. The parameter should be the string % temperature is constant. The :parameter should be the string
% "on" to enable the energy equation, or "off" to disable it. % "on" to enable the energy equation, or "off" to disable it.
% %
% 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
% %
% parameter r: % :parameter r:
% Instance of class 'Reactor'. % Instance of class 'Reactor'.
% parameter 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.
checklib;
iflag = -1; iflag = -1;
if strcmp(flag, 'on') if strcmp(flag, 'on')
iflag = 1; iflag = 1;

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

@ -17,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.
% %
% parameter 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;
@ -51,19 +51,18 @@ classdef ReactorNet < handle
function clear(r) function clear(r)
% Clear the ReactorNet object from the memory. % Clear the ReactorNet object from the memory.
checklib;
calllib(ct, '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.
% %
% parameter net: % :parameter net:
% Instance of class 'ReactorNet'. % Instance of class 'ReactorNet'.
% parameter reactor: % :parameter reactor:
% Instance of class 'Solution'. % Instance of class 'Solution'.
checklib;
calllib(ct, 'reactornet_addreactor', net.id, reactor.id); calllib(ct, 'reactornet_addreactor', net.id, reactor.id);
end end
@ -77,10 +76,9 @@ 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.
% %
% parameter tout: % :parameter tout:
% End time of the integration. Unit: s. % End time of the integration. Unit: s.
checklib;
calllib(ct, 'reactornet_advance', r.id, tout); calllib(ct, 'reactornet_advance', r.id, tout);
end end
@ -89,11 +87,10 @@ 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.
% %
% parameter 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;
calllib(ct, 'reactornet_setInitialTime', r.id, t); calllib(ct, 'reactornet_setInitialTime', r.id, t);
end end
@ -108,19 +105,17 @@ classdef ReactorNet < handle
% leads to numerical problems later. Use thismethod to set an % leads to numerical problems later. Use thismethod to set an
% upper bound on the timestep. % upper bound on the timestep.
checklib;
calllib(ct, '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.
% %
% parameter rtol: % :parameter rtol:
% Scalar relative error tolerance. % Scalar relative error tolerance.
% parameter atol: % :parameter atol:
% Scalar absolute error tolerance. % Scalar absolute error tolerance.
checklib;
calllib(ct, 'reactornet_setTolerances', r.id, rerr, aerr); calllib(ct, 'reactornet_setTolerances', r.id, rerr, aerr);
end end
@ -134,25 +129,25 @@ classdef ReactorNet < handle
% tolerance is maintained. At times when the solution is % tolerance is maintained. At times when the solution is
% rapidly changing, the time step becomes smaller to resolve % rapidly changing, the time step becomes smaller to resolve
% the solution. % the solution.
checklib;
t = calllib(ct, '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;
t = calllib(ct, '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;
t = calllib(ct, '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;
t = calllib(ct, 'reactornet_atol', r.id); t = calllib(ct, 'reactornet_atol', r.id);
end end

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

@ -1,7 +1,7 @@
classdef ReactorSurface < handle classdef ReactorSurface < handle
properties properties
surf_id surfID
area area
reactor reactor
end end
@ -21,22 +21,22 @@ 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.
% %
% parameter 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'.
% parameter 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.
% parameter 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.surf_id = calllib(ct, 'reactorsurface_new', 0); s.surfID = calllib(ct, 'reactorsurface_new', 0);
s.reactor = -1; s.reactor = -1;
% if r.id < 0 % if r.id < 0
% error(geterr); % error(geterr);
@ -68,48 +68,46 @@ classdef ReactorSurface < handle
function clear(s) function clear(s)
% Clear the ReactorSurface object from the memory. % Clear the ReactorSurface object from the memory.
checklib;
calllib(ct, 'reactorsurface_del', s.surf_id); calllib(ct, 'reactorsurface_del', s.surfID);
end end
function install(s, r) function install(s, r)
% Install a ReactorSurface in a Reactor. % Install a ReactorSurface in a Reactor.
checklib;
s.reactor = r; s.reactor = r;
calllib(ct, 'reactorsurface_install', s.surf_id, r.id); calllib(ct, 'reactorsurface_install', s.surfID, r.id);
end end
%% ReactorSurface get methods %% ReactorSurface get methods
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;
a = calllib(ct, 'reactorsurface_area', s.surf_id); a = calllib(ct, 'reactorsurface_area', s.surfID);
end end
%% ReactorSurface set methods %% ReactorSurface set methods
function set.area(s, a) function set.area(s, a)
% Set the area of a reactor surface % Set the area of a reactor surface
checklib;
calllib(ct, 'reactorsurface_setArea', s.surf_id, a); calllib(ct, 'reactorsurface_setArea', s.surfID, 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.
% parameter 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'.
checklib;
ikin = 0; ikin = 0;
if isa(kin, 'Kinetics') if isa(kin, 'Kinetics')
ikin = kin.kin_id; ikin = kin.kinID;
end end
calllib(ct, 'reactorsurface_setkinetics', s.surf_id, ikin); calllib(ct, 'reactorsurface_setkinetics', s.surfID, ikin);
end end
end end

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

@ -46,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.
% %
% parameter 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.
% parameter 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.
% parameter 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.
% parameter 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.
% parameter 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.
% parameter 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'.
% parameter 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'.
@ -137,14 +137,13 @@ classdef Wall < handle
function clear(w) function clear(w)
% Clear the Wall object from the memory. % Clear the Wall object from the memory.
checklib;
calllib(ct, 'wall_del', w.id); calllib(ct, 'wall_del', w.id);
end end
function install(w, l, r) function install(w, l, r)
% Install a wall between two reactors. % Install a wall between two reactors.
checklib;
w.left = l; w.left = l;
w.right = r; w.right = r;
calllib(ct, 'wall_install', w.id, l.id, r.id); calllib(ct, 'wall_install', w.id, l.id, r.id);
@ -152,7 +151,7 @@ classdef Wall < handle
function ok = ready(w) function ok = ready(w)
% Check whether a wall is ready. % Check whether a wall is ready.
checklib;
ok = calllib(ct, 'wall_ready', w.id); ok = calllib(ct, 'wall_ready', w.id);
end end
@ -160,37 +159,34 @@ classdef Wall < handle
function set.area(w, a) function set.area(w, a)
% Set the area of a wall. % Set the area of a wall.
checklib;
calllib(ct, 'wall_setArea', 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.
% %
% parameter r: % :parameter r:
% Thermal resistance. Unit: K*m^2/W. % Thermal resistance. Unit: K*m^2/W.
checklib;
calllib(ct, '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.
% %
% parameter u: % :parameter u:
% Heat transfer coefficient. Unit: W/(m^2-K). % Heat transfer coefficient. Unit: W/(m^2-K).
checklib;
calllib(ct, '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.
% %
% parameter k: % :parameter k:
% Expanstion rate coefficient. Unit: m/(s-Pa). % Expanstion rate coefficient. Unit: m/(s-Pa).
checklib;
calllib(ct, 'wall_setExpansionRateCoeff', w.id, k); calllib(ct, 'wall_setExpansionRateCoeff', w.id, k);
end end
@ -202,10 +198,9 @@ 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.
% %
% parameter f: % :parameter f:
% Instance of class 'Func'. Unit: W/m^2. % Instance of class 'Func'. Unit: W/m^2.
checklib;
calllib(ct, 'wall_setHeatFlux', w.id, f.id); calllib(ct, 'wall_setHeatFlux', w.id, f.id);
end end
@ -217,10 +212,9 @@ 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.
% %
% parameter f: % :parameter f:
% Instance of class 'Func'. Unit: m/s. % Instance of class 'Func'. Unit: m/s.
checklib;
calllib(ct, 'wall_setVelocity', w.id, f.id); calllib(ct, 'wall_setVelocity', w.id, f.id);
end end
@ -228,19 +222,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;
a = calllib(ct, '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;
q = calllib(ct, '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;
v = calllib(ct, 'wall_vdot', w.id, t); v = calllib(ct, 'wall_vdot', w.id, t);
end end

2
samples/matlab_experimental/ignite.m
View File

@ -11,7 +11,7 @@ function plotdata = ignite(g)
if nargin == 1 if nargin == 1
gas = g; gas = g;
else else
gas = Solution('gri30.yaml'); gas = Solution('gri30.yaml', 'gri30');
end end
% set the initial conditions % set the initial conditions

2
samples/matlab_experimental/ignite_hp.m
View File

@ -5,7 +5,7 @@ function ignite_hp(gas)
help ignite_hp help ignite_hp
if nargin == 0 if nargin == 0
gas = Solution('gri30.yaml'); gas = Solution('gri30.yaml', 'gri30');
end end
mw = gas.MolecularWeights; mw = gas.MolecularWeights;

2
samples/matlab_experimental/ignite_uv.m
View File

@ -5,7 +5,7 @@ function ignite_uv(gas)
help ignite_uv help ignite_uv
if nargin == 0 if nargin == 0
gas = Solution('gri30.yaml'); gas = Solution('gri30.yaml', 'gri30');
end end
mw = gas.MolecularWeights; mw = gas.MolecularWeights;