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

[Matlab] Initial version of 'experimental' Matlab toolbox

Browse Source 
New toolbox replacing old-style classes with 'classdef' classes,
and using 'calllib' to interact directly with the Cantera C interface
rather than the 'ctmethods' MEX file.
This commit is contained in:
Su Sun
2021-11-01 15:20:01 -04:00
committed by Ray Speth
parent 63313b08d5
commit 03ca9302b0
23 changed files with 2660 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
interfaces/matlab_experimental/Base/Air.m Normal file
View File

@@ -0,0 +1,4 @@
function gas = Air()
% Create an object representing air.
gas = Solution('air.yaml', 'air');
end

4
interfaces/matlab_experimental/Base/CarbonDioxide.m Normal file
View File

@@ -0,0 +1,4 @@
function c = CarbonDioxide()
% Return an object representing carbon dioxide.
h = Solution('liquidvapor.yaml', 'carbon-dioxide');
end

4
interfaces/matlab_experimental/Base/HFC134a.m Normal file
View File

@@ -0,0 +1,4 @@
function h = HFC134a()
% Return an object representing refrigerant HFC134a.
h = Solution('liquidvapor.yaml', 'HFC-134a');
end

4
interfaces/matlab_experimental/Base/Heptane.m Normal file
View File

@@ -0,0 +1,4 @@
function h = Heptane
% Return an object representing n-heptane.
h = Solution('liquidvapor.yaml', 'heptane');
end

4
interfaces/matlab_experimental/Base/Hydrogen.m Normal file
View File

@@ -0,0 +1,4 @@
function h = Hydrogen()
% Return an object representing hydrogen.
h = Solution('liquidvapor.yaml', 'hydrogen');
end

130
interfaces/matlab_experimental/Base/Interface.m Normal file
View File

@@ -0,0 +1,130 @@
classdef Interface < handle
properties
th
kin
end
properties(Constant = true)
lib = 'cantera_shared'
end
methods
%% Interface class constructor
function s = Interface(src, id, p1, p2, p3, p4)
checklib;
t = ThermoPhase(src, id);
if nargin == 2
k = Kinetics(t, src, id);
elseif nargin == 3
k = Kinetics(t, src, id, p1);
elseif nargin == 4
k = Kinetics(t, src, id, p1, p2);
elseif nargin == 5
k = Kinetics(t, src, id, p1, p2, p3);
elseif nargin == 6
k = Kinetics(t, src, id, p1, p2, p3, p4);
end
s.kin = k;
s.th = t;
end
%% Interface methods
function c = coverages(s)
% Get the surface coverages of the species on an interface.
%
% :return:
% If no output value is assigned, a bar graph will be
% plotted. Otherwise, a vector of length "n_surf_species"
% will be returned.
checklib;
surf_id = s.th.tr_id;
nsp = s.th.nSpecies;
xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx);
calllib(s.lib, 'surf_getCoverages', surf_id, xx);
c = pt.Value;
if nargout == 0
figure
set(gcf, 'Name', 'Coverages')
bar(c);
colormap(summer);
nm = speciesNames(s);
set(gca, 'XTickLabel', nm);
xlabel('Species Name');
ylabel('Coverage');
title('Surface Species Coverages');
end
end
function c = concentrations(s)
% Get the concentrations of the species on an interface.
%
% :return:
% If no output value is assigned, a bar graph will be
% plotted. Otherwise, a vector of length "n_surf_species"
% will be returned.
checklib;
surf_id = s.th.tr_id;
nsp = s.th.nSpecies;
xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx);
calllib(s.lib, 'surf_getConcentrations', surf_id, xx);
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
function setCoverages(s, cov, norm)
% Set surface coverages of the species on an interface.
%
% :param cov:
% Coverage of the species. "Cov" can be either a vector of
% length "n_surf_species", or a string in the format
% "Species:Coverage, Species:Coverage".
checklib;
if nargin == 3 && strcmp(norm, 'nonorm')
norm_flag = 0;
else
norm_flag = 1;
end
surf_id = s.th.tr_id;
nsp = s.th.nSpecies;
[m, n] = size(cov);
if isa(cov, 'double')
sz = length(cov);
if sz == nsp
if ((m == nsp && n == 1) || (m == 1 & n == nsp))
calllib(s.lib, 'surf_setCoverages', surf_id, cov, norm_flag);
else error('wrong size for coverage array');
end
else
error('wrong size for coverage array');
end
elseif isa(cov, 'char')
calllib(s.lib, 'surf_setCoveragesByName', surf_id, cov);
end
end
end
end

577
interfaces/matlab_experimental/Base/Kinetics.m Normal file
View File

@@ -0,0 +1,577 @@
classdef Kinetics < handle
properties
kin_owner
kin_id
Kc % equilibrium constant
Kf % forward reaction rate
Kr % reverse reaction rate
dH % enthalpy of reaction
dS % entropy of reaction
dG % gibbs free energy of reaction
end
properties(Constant = true)
lib = 'cantera_shared'
end
methods
%% Kinetics class constructor
function kin = Kinetics(ph, src, id, n1, n2, n3, n4)
checklib;
% indices for bulk phases in a heterogeneous mechanism
inb1 = -1;
inb2 = -1;
inb3 = -1;
inb4 = -1;
if nargin == 2
id = '-';
end
kin.kin_owner = 1;
% get the integer indices used to find the stored objects
% representing the phases participating in the mechanism
iph = ph.tp_id;
if nargin > 3
inb1 = n1.tp_id;
if nargin > 4
inb2 = n2.tp_id;
if nargin > 5
inb3 = n3.tp_id;
if nargin > 6
inb4 = n4.tp_id;
end
end
end
end
kin.kin_id = calllib(kin.lib, 'kin_newFromFile', src, id, ...
iph, inb1, inb2, inb3, inb4);
end
%% Utility methods
function kin_clear(kin)
% Delete the kernel object
checklib;
calllib(kin.lib, 'kin_del', kin.kin_id);
end
%% Get scalar attributes
function n = nReactions(kin)
% Get the number of reactions.
%
% :return:
% Integer number of reactions
checklib;
n = calllib(kin.lib, 'kin_nReactions', kin.id);
end
function n = multiplier(kin, irxn)
% Get the multiplier for reaction rate of progress.
%
% :param irxn:
% Integer reaction number for which the multiplier is
% desired.
% :return:
% Multiplier of the rate of progress of reaction irxn.
checklib;
n = calllib(kin.lib, 'kin_multiplier', kin.id, irxn-1);
end
function n = nSpecies(kin)
% Get the total number of species.
%
% :return:
% Integer total number of species.
checklib;
n = calllib(kin.lib, 'kin_nSpecies', kin.id);
end
function n = isReversible(kin, i)
% Get an array of flags indicating reversibility of a reaction.
%
% :param i:
% Integer reaction number.
% :return:
% 1 if reaction number i is reversible. 0 if irreversible.
checklib;
n = calllib(kin.lib, 'kin_isReversible', kin.id, i);
end
function n = stoich_r(kin, species, rxns)
% Get the reactant stoichiometric coefficients.
%
% :param species:
% Species indices for which reactant stoichiometric
% coefficients should be retrieved. Optional argument; if
% specified, "rxns" must be specified as well.
% :param rxns:
% Reaction indicies for which reactant stoichiometric
% coefficients should be retrieved. Optional argument; if
% specified, "species" must be specified as well.
% :return:
% Returns a sparse matrix of all reactant stoichiometric
% coefficients. The matrix elements "nu(k, i)" is the
% stoichiometric coefficient of species k as a reactant in
% reaction i. If "species" and "rxns" are specified, the
% matrix will contain only entries for the specified species
% and reactions. For example, "stoich_p(a, 3, [1, 3, 5,
% 7])" returns a sparse matrix containing only the
% coefficients for specis 3 in reactions 1, 3, 5, and 7.
checklib;
nsp = kin.nSpecies;
nr = kin.nReactions;
temp = sparse(nsp, nr);
if nargin == 1
krange = 1:nsp;
irange = 1:nr;
elseif nargin == 3
krange = species;
irange = rxns;
else error('stoich_r requires 1 or 3 arguments.')
end
for k = range
for i = irange
t = calllib(kin.lib, 'kin_reactantStoichCoeff', ...
kin.id, k-1, i-1);
if t ~= 0.0
temp(k, i) = t;
end
end
end
n = temp;
end
function n = stoich_p(kin, species, rxns)
% Get the product stoichiometric coefficients.
%
% :param species:
% Species indices for which product stoichiometric
% coefficients should be retrieved. Optional argument; if
% specified, "rxns" must be specified as well.
% :param rxns:
% Reaction indicies for which product stoichiometric
% coefficients should be retrieved. Optional argument; if
% specified, "species" must be specified as well.
% :return:
% Returns a sparse matrix of all product stoichiometric
% coefficients. The matrix elements "nu(k, i)" is the
% stoichiometric coefficient of species k as a product in
% reaction i. If "species" and "rxns" are specified, the
% matrix will contain only entries for the specified species
% and reactions. For example, "stoich_p(a, 3, [1, 3, 5,
% 7])" returns a sparse matrix containing only the
% coefficients for specis 3 in reactions 1, 3, 5, and 7.
checklib;
nsp = kin.nSpecies;
nr = kin.nReactions;
temp = sparse(nsp, nr);
if nargin == 1
krange = 1:nsp;
irange = 1:nr;
elseif nargin == 3
krange = species;
irange = rxns;
else error('stoich_p requires 1 or 3 arguments.')
end
for k = krange
for i = irange
t = calllib(kin.lib, 'kin_productStoichCoeff', ...
kin.id, k-1, i-1);
if t ~= 0.0
temp(k, i) = t;
end
end
end
n = temp;
end
function n = stoich_net(kin, species, rxns)
% Get the net stoichiometric coefficients.
%
% :param species:
% Species indices for which net stoichiometric coefficients
% should be retrieved. Optional argument; if specified,
% "rxns" must be specified as well.
% :param rxns:
% Reaction indicies for which net stoichiometric
% coefficients should be retrieved. Optional argument; if
% specified, "species" must be specified as well.
% :return:
% Returns a sparse matrix of all net stoichiometric
% coefficients. The matrix elements "nu(k, i)" is the
% stoichiometric coefficient of species k as a net in
% reaction i. If "species" and "rxns" are specified, the
% matrix will contain only entries for the specified species
% and reactions. For example, "stoich_net(a, 3, [1, 3, 5,
% 7])" returns a sparse matrix containing only the
% coefficients for specis 3 in reactions 1, 3, 5, and 7.
checklib;
if nargin == 1
nu = stoich_p(kin)-stoich_r(kin);
elseif nargin == 3
nu = stoich_p(a, species, rxns) - stoich_r (a, species, rxns);
else error('stoich_net requires 1 or 3 arguments.');
end
end
%% Get reaction array attributes
function q = rop_f(kin)
% Get the forward rates of progress for all reactions.
%
% :return:
% Returns a column vector of the forward rates of progress
% for all reactions. If this function is called without
% argument, a bar graph is produced.
checklib;
nr = kin.nReactions;
xx = zeros(1, nr);
pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getFwdRateOfProgress', kin.id, nr, pt);
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
function q = rop_r(kin)
% Get the reverse rates of progress for all reactions.
%
% :return:
% Returns a column vector of the reverse rates of progress
% for all reactions. If this function is called without
% argument, a bar graph is produced.
checklib;
nr = kin.nReactions;
xx = zeros(1, nr);
pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getRevRateOfProgress', kin.id, nr, pt);
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
function q = rop(kin)
% Get the forward and reverse rates of progress.
%
% :return:
% Returns an I x 2 array of reaction rates of progress,
% where I is the number of reactions. The first column
% contains the forward rates of progress, and the second
% column the reverse rates. If this function is called
% without arguments, a bar graph is produced.
checklib;
f = rop_f(kin);
r = rop_r(kin);
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
function q = rop_net(kin)
% Get the net rates of progress for all reactions.
%
% :return:
% Returns a column vector of the net rates of progress
% for all reactions. If this function is called without
% argument, a bar graph is produced.
checklib;
nr = kin.nReactions;
xx = zeros(1, nr);
pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getNetRateOfProgress', kin.id, nr, pt);
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
function k = get.Kc(kin)
% Get the equilibrium constants for all reactions.
%
% :return:
% Returns a column vector of the equilibrium constants for
% all reactions. The vector has an entry for every reaction,
% whether reversible or not, but non-zero values occur only
% for the reversible reactions. If the output is not
% assigned to a variable, a bar graph is produced.
checklib;
nr = kin.nReactions;
xx = zeros(1, nr);
pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getEquilibriumConstants', kin.id, nr, pt);
k = pt.Value;
if nargout == 0
figure
set(gcf, 'Name', 'Equilibrium Constants')
bar(q)
xlabel('Reaction Number')
ylabel('log_{10} Kc [kmol,m, s]')
title('Equilibrium Constants')
end
end
function k = get.Kf(kin)
% Get the forward reaction rate constants for all reactions.
%
% :return:
% Returns a column vector of the forward rates constants of
% all of the reactions.
checklib;
nr = kin.nReactions;
xx = zeros(1, nr);
pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getFwdRateConstants', kin.id, nr, pt);
k = pt.Value;
end
function k = get.Kr(kin)
% Get the reverse reaction rate constants for all reactions.
%
% :return:
% Returns a column vector of the reverse rates constants of
% all of the reactions.
checklib;
nr = kin.nReactions;
xx = zeros(1, nr);
pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getRevRateConstants', kin.id, 1, nr, pt);
k = pt.Value;
end
function enthalpy = get.dH(kin)
% Get the enthalpy of reaction for each reaction.
%
% :return:
% Returns a vector of the enthalpy of reaction for each
% reaction. Unit: J/kmol.
checklib;
nr = kin.nReactions;
xx = zeros(1, nr);
pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getDelta', kin.id, 0, nr, pt);
enthalpy = pt.Value;
end
function entropy = get.dS(kin)
% Get the entropy of reaction for each reaction.
%
% :return:
% Returns a vector of the entropy of reaction for each
% reaction. Unit: J/kmol-K.
checklib;
nr = kin.nReactions;
xx = zeros(1, nr);
pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getDelta', kin.id, 2, nr, pt);
entropy = pt.Value;
end
function gibbs = get.dG(kin)
% Get the Gibbs free energy of reaction for each reaction.
%
% :return:
% Returns a vector of the Gibbs free energy of reaction for
% each reaction. Unit: J/kmol.
checklib;
nr = kin.nReactions;
xx = zeros(1, nr);
pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getDelta', kin.id, 1, nr, pt);
gibbs = pt.Value;
end
function cdot = creationRates(kin)
% Get the chemical reaction rates.
%
% :return:
% Returns a vector of the creation rates of all species. If
% the output is not assigned to a variable, a bar graph is
% produced. Unit: kmol/m^3-s.
checklib;
nsp = kin.nSpecies;
xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getCreationRates', kin.id, nsp, pt);
cdot = pt.Value;
if nargout == 0
figure
set(gcf, 'Name', 'Creation Rates')
bar(q)
xlabel('Species Number')
ylabel('Creation Rate [kmol/m^3-s]')
title('Species Chemical Reaction Rates')
end
end
function ddot = destructionRates(kin)
% Get the chemical destruction rates.
%
% :return:
% Returns a vector of the destruction rates of all species.
% If the output is not assigned to a variable, a bar graph
% is produced. Unit: kmol/m^3-s.
checklib;
nsp = kin.nSpecies;
xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getDestructionRates', kin.id, nsp, pt);
ddot = pt.Value;
if nargout == 0
figure
set(gcf, 'Name', 'Destruction Rates')
bar(q)
xlabel('Species Number')
ylabel('Destruction Rate [kmol/m^3-s]')
title('Species Chemical Reaction Rates')
end
end
function wdot = netProdRates(kin)
% Get the net chemical production rates for all species.
%
% :return:
% Returns a vector of the net production (creation-destruction)
% rates of all species. If the output is not assigned to a
% variable, a bar graph is produced. Unit: kmol/m^3-s.
checklib;
nsp = kin.nSpecies;
xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getNetProductionRates', kin.id, nsp, pt);
wdot = pt.Value;
if nargout == 0
figure
set(gcf, 'Name', 'Production Rates')
bar(q)
xlabel('Species Number')
ylabel('Net Production Rate [kmol/m^3-s]')
title('Species Net Chemical Reaction Rates')
end
end
function ydot = massProdRates(kin)
% Get the mass production rates of the species.
%
% :return:
% Returns a vector of the mass production rates.
checklib;
nsp = kin.nSpecies;
xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx);
calllib(kin.lib, 'kin_getSourceTerms', kin.id, nsp, pt);
ydot = pt.Value;
end
% function e = reactionEqn(kin, irxn)
% % Get the reaction equation of a reaction
% %
% % :param irxn:
% % Optional. Integer or vector of reaction numbers.
% % :return:
% % String or cell arrray of strings of the reaction
% % equations.
% Need to resolve string printing issues first!
%
% if nargin == 1
% nr = kin.nReactions;
% irxn = (1:m)';
% elseif nargin == 2
% if isa(irxn, 'double')
% [m, n] = size(irxn);
% else
% error('reaction numbers must be numeric');
% end
% end
% end
%% Set attributes
function n = setMultiplier(kin, irxn, v)
% Set the multiplier for the reaction rate of progress.
%
% :param irxn:
% Integer of vector reaction numbers for which the
% multiplier should be set. Optional.
% :param v:
% Value by which the reaction rate of progress should be
% multiplied.
if nargin == 2
v = irxn;
nr = kin.nReactions;
irxn = (1:nr)';
n = 1;
elseif nargin == 3
[nr, n] = size(irxn);
else
error('setMultiplier requires 2 or 3 arguments.')
end
for i = 1:nr
for j = 1:n
calllib(kin.lib, 'kin_setMultiplier', kin.id, ...
irxn(i, j)-1, v);
end
end
end
function advanceCoeverages(kin, dt)
% Advance the surface coveages forward in time
%
% :param dt:
% Time interval by which the coverages should be advanced.
calllib(kin.lib, 'kin_advanceCoverages', kin.id, dt);
end
end
end

4
interfaces/matlab_experimental/Base/Methane.m Normal file
View File

@@ -0,0 +1,4 @@
function m = Methane()
% Return an object representing methane.
h = Solution('liquidvapor.yaml', 'methane');
end

331
interfaces/matlab_experimental/Base/Mixture.m Normal file
View File

@@ -0,0 +1,331 @@
classdef Mixture < handle
properties
mixindex
phases
T
P
end
properties(Constant = true)
lib = 'cantera_shared'
end
methods
%% Mixture class constructor
function m = Mixture(phases)
% To construct a mixture, supply a cell array of phases and mole
% numbers:
% >> gas = Solution('gas.cti');
% >> graphite = Solution('graphite.cti');
% >> mix = Mixture({gas, 1.0; graphite, 0.1});
%
% Phases may also be added later using the addPhase method:
% >> water = Solution('water.cti');
% >> addPhase(mix, water, 3.0);
%
% Note that the objects representing each phase compute only
% the intensive state of the phase - they do not store any
% information on the amount of this phase. Mixture objects, on
% the other hand, represent full extensive state.
%
% Mixture objects are 'lightweight'in the sense that they do
% not store parameters needed to compute thermodynamic or
% kinetic properties of the phases. These are contained in the
% ('heavyweight') phase objects. Multiple mixture objects are
% constructed using the same set of phase objects. Each one
% store its own state information locally, and syncrhonizes the
% phase objects whenever itrequires phase properties.
%
% :param phases:
% Cell array of phases and mole numbers.
% :return:
% Instance of class 'Mixture'.
checklib;
if nargin > 1
error('Mixture: wrong number of arguments');
end
% Create an empty mixture.
m.mixindex = calllib(m.lib, 'mix_new');
m.phases = phases;
% If phases are supplied, add them
if nargin == 1
if ~isa(phases, 'cell')
error('Enter phasesas a cell array.');
end
% First column contains the phase objects, and the second
% column contains the mole numbers of each phase.
[np, nc] = size(phases);
if nc ~= 2
error('Cell array of phases should have each phase on a new row');
end
for n = 1:np
addPhase(m, phases{n, 1}, phases{n, 2});
end
m.T = phases{n, 1}.T;
m.P = phases{n, 1}.P;
end
end
%% Utility methods
function display(m)
% Display the state of the mixture on the terminal.
calllib(m.lib, 'mix_updatePhases', m.mixindex);
[np, nc] = size(m.phases);
for n = 1:np
s = [sprintf('\n******************* Phase %d', n) ...
sprintf(' ******************************\n\n Moles: %12.6g', ...
phaseMoles(m, n))];
disp(s);
display(m.phases{n, 1});
end
end
function mix_clear(m)
% Delete the MultiPhase object.
checklib;
calllib(m.lib, 'mix_del', m.mixindex);
end
%% Mixture Get methods
function addPhase(m, phase, moles)
% Add a phase to the mixture
%
% :param m:
% Instance of class 'Mixture' to which phases is added.
% :param phase:
% Instance of class 'ThermoPhase' which should be added.
% :param moles:
% Number of moles of the phase to be added. Unit: kmol.
checklib;
if ~isa(phase, 'ThermoPhase')
error('Phase object of wrong type.');
end
if ~isa(moles, 'numeric')
error('Number of moles must be numeric');
end
if moles < 0.0
error('Negative moles');
end
iok = calllib(m.lib, 'mix_addPhase', m.mixindex, phase.tp_id, ...
moles);
if iok < 0
error('Error adding phase');
end
end
function temperature = get.T(m)
% Get the temperature of the mixture.
%
% :return:
% Temperature in K.
checklib;
temperature = calllib(m.lib, 'mix_temperature', m.mixindex);
end
function pressure = get.P(m)
% Get the pressure of themixture.
%
% :return:
% Pressure in Pa.
checklib;
pressure = calllib(m.lib, 'mix_pressure', m.mixindex);
end
function n = nPhases(m)
% Get the number of phases in the mixture.
checklib;
n = calllib(m.lib, 'mix_nPhases', m.mixindex);
end
function n = nElements(m)
% Get the number of elements in the mixture.
checklib;
n = calllib(m.lib, 'mix_nElements', m.mixindex);
end
function n = nSpecies(m)
% Get the number of species in the mixture.
checklib;
n = calllib(m.lib, 'mix_nSpecies', m.mixindex);
end
function n = elementIndex(m, name)
% Get the index of element 'name'.
%
% Note: In keeping with the conventions used by Matlab, the
% indices start from 1 instead of 0 as in Cantera C++ and
% Python interfaces.
checklib;
n = calllib(m.lib, 'mix_elementIndex', m.mixindex, name) + 1;
end
function n = speciesIndex(m, k, p)
% Get the index of element 'name'.
%
% Note: In keeping with the conventions used by Matlab, the
% indices start from 1 instead of 0 as in Cantera C++ and
% Python interfaces.
checklib;
n = calllib(m.lib, 'mix_speciesIndex', m.mixindex, k-1, p-1) + 1;
% check back on this one!
end
function moles = phaseMoles(m, n)
% Get the number of moles of a phase in a mixture.
%
% :param n:
% Integer phase number in the input.
% :return:
% Moles of phase number 'n'. Unit: kmol.
checklib;
if nargin == 2
moles = calllib(m.lib, 'mix_phaseMoles', m.mixindex, n);
elseif nargin == 1
np = m.nPhases;
moles = zeros(1, np);
for i = 1:np
moles(i) = calllib(m.lib, 'mix_phaseMoles', ...
m.mixindex, i);
end
else error('wrong number of arguments');
end
end
function mu = chemPotentials(m)
% Get the chemical potentials of species in the mixture.
%
% :return:
% Vector of chemical potentials. Unit: J/kmol.
checklib;
nsp = m.nSpecies;
xx = zeros(1, nsp);
ptr = libpointer('doublePtr', xx);
calllib(m.lib, 'mix_getChemPotential', m.mixindex, nsp, ptr);
mu = ptr.Value;
end
%% Mixture Set methods
function m = set.T(m, temp)
% Set the mixture temperature.
%
% :param temp:
% Temperature to set. Unit: K.
checklib;
calllib(m.lib, 'mix_setTemperature', m.mixindex, temp);
end
function m = set.P(m, pressure)
% Set the mixture pressure.
%
% :param pressure:
% Pressure to set. Unit: Pa.
checklib;
calllib(m.lib, 'mix_setPressure', m.mixindex, pressure);
end
function setPhaseMoles(m, n, moles)
% Set the number of moles of phase n in the mixture.
%
% :param n:
% Phase number.
% :param moles:
% Number of moles to set. Unit: kmol.
checklib;
calllib(m.lib, 'mix_setPhaseMoles', m.mixindex, n-1, moles);
end
function setSpeciesMoles(m, moles)
% Set the moles of the species. The moles may be specified as a
% string or as a vector. If a vector is used, it must be
% dimensioned at least as large as the total number of species
% in the mixture. Note that the species may belong to any
% phase, and unspecified species are set to zero.
%
% :param moles:
% Vector or string specifying the moles of species.
checklib;
calllib(m.lib, 'mix_setMolesByName', m.mixindex, moles);
% check back on this one!
end
function r = equilibrate(m, XY, err, maxsteps, maxiter, loglevel)
% Set the mixture to a state of chemical equilibrium.
%
% This method uses a version of the VCS algorithm to find the
% composition that minimizes the total Gibbs free energy of the
% mixture, subject to element conservation constraints. For a
% description of the theory, see Smith and Missen, "Chemical
% Reaction Equilibrium". The VCS algorithm is implemented in
% Cantera kernel class MultiPhaseEquil.
%
% The VCS algorithm solves for the equilibrium composition for
% specified temperature and pressure. If any other property
% pair other than 'TP' is specified, then an outer iteration
% loop is used to adjust T and/or P so that the specified
% property values are obtained:
% >> equilibrate(mix, 'TP);
% >> equilibrate('TP', 1.0e-6, 500);
%
% :param XY:
% Two-letter string specifying the two properties to hold
% fixed. Currently 'TP', 'HP', 'TV', and 'SP' have been
% implemented. Default: 'TP'.
% :param err:
% Error tolerance. Iteration will continue until delta_Mu/RT
% is less than this value for each reaction. Default:
% 1.0e-9.
% :param maxsteps:
% Maximum number of steps to take while solving the
% equilibrium problem for specified T and P. Default: 1000.
% :param maxiter:
% Maximum number of temperature and/or pressure iterations.
% This is only relevant if a property pair other than (T,
% P)is specified. Default: 200.
% :param loglevel:
% Set to a value > 0 to write diagnostic output. Larger
% values generate more detailed information.
% :return:
% The error in the solution.
checklib;
if nargin < 6
loglevel = 0;
end
if nargin < 5
maxiter = 200;
end
if nargin < 4
maxsteps = 1000;
end
if nargin < 3
err = 1.0e-9;
end
if nargin < 2
XY = 'TP'
end
r = calllib(m.lib, 'mix_equilibrate', m.mixindex, XY, err, ...
maxsteps, maxiter, loglevel);
end
end
end

4
interfaces/matlab_experimental/Base/Nitrogen.m Normal file
View File

@@ -0,0 +1,4 @@
function n = Nitrogen()
% Return an object representing nitrogen.
h = Solution('liquidvapor.yaml', 'nitrogen');
end

4
interfaces/matlab_experimental/Base/Oxygen.m Normal file
View File

@@ -0,0 +1,4 @@
function o = Oxygen()
% Return an object representing oxygen.
h = Solution('liquidvapor.yaml', 'oxygen');
end

37
interfaces/matlab_experimental/Base/Solution.m Normal file
View File

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

893
interfaces/matlab_experimental/Base/ThermoPhase.m Normal file
View File

@@ -0,0 +1,893 @@
classdef ThermoPhase < handle
properties
tp_owner
tp_id
T % temperature
P % pressure
D % density
X % mole fractions
Y % mass fractions
H % enthalpy
S % entropy
U % internal energy
% V % specific volume not added yet
basis
end
properties (Constant = true)
lib = 'cantera_shared'
end
properties (Dependent)
DP
DPX
DPY
HP
HPX
HPY
SP
SPX
SPY
% SVX
% SVY
TD
TDX
TDY
TP
TPX
TPY
% UV
% UVX
% UVY
end
methods
%% ThermoPhase class constructor
function tp = ThermoPhase(src, id)
checklib;
if nargin > 2
error('ThermoPhase expects 1 or 2 input arguments.');
end
if nargin == 1
id = '-';
end
tp.tp_owner = 1;
tp.tp_id = calllib(tp.lib, 'thermo_newFromFile', src, id);
tp.basis = 'molar';
end
%% Utility methods
function display(tp, threshold)
% Display thermo properties
if nargin < 2 || ~isnumeric(threshold)
threshold = 1e-14;
end
calllib(tp.lib, 'thermo_print', tp.tp_id, 1, threshold);
end
function tp_clear(tp)
% Delete the kernel object.
checklib;
calllib(tp.lib, 'thermo_del', tp.tp_id);
end
function tp = set.basis(tp, b)
% Set basis of thermodynamic properties
% Default is molar
if strcmp(b, 'mole') || strcmp(b, 'molar') ...
|| strcmp(b, 'Mole') || strcmp(b, 'Molar')
tp.basis = 'molar'
elseif strcmp(b, 'mass') || strcmp(b, 'Mass')
tp.basis = 'mass'
else error("Basis must be mass or molar")
end
end
%% PhaseGet single methods
function temperature = get.T(tp)
% Get the temperature.
%
% :return:
% Double temperature. Unit: K
checklib;
temperature = calllib(tp.lib, 'thermo_temperature', tp.tp_id);
end
function pressure = get.P(tp)
% Get the pressure.
%
% :return:
% Double pressure. Unit: Pa
checklib;
pressure = calllib(tp.lib, 'thermo_pressure', tp.tp_id);
end
function density = get.D(tp)
% Get the density depending on the basis.
%
% :return:
% Density depending on the basis. Units:
% kmol/m^3 (molar) kg/m^3 (mass).
checklib;
if strcmp(tp.basis, 'molar')
density = calllib(tp.lib, 'thermo_density', tp.tp_id);
elseif strcmp(tp.basis, 'mass')
density = calllib(tp.lib, 'thermo_molarDensity', tp.tp_id);
else error("basis not specified");
end
end
function mmw = meanMolecularWeight(tp)
% Get the mean molecular weight.
%
% :return:
% Double mean molecular weight. Unit: kg/kmol
checklib;
mmw = calllib(tp.lib, 'thermo_meanMolecularWeight', tp.tp_id);
end
function nel = nElements(tp)
% Get the number of elements.
%
% :return:
% Integer number of elements in the phase.
checklib;
nel = calllib(tp.lib, 'thermo_nElements', tp.tp_id);
end
function nsp = nSpecies(tp)
% Get the number of species.
%
% :return:
% Integer number of species in the phase.
checklib;
nsp = calllib(tp.lib, 'thermo_nSpecies', tp.tp_id);
end
function k = speciesIndex(tp, name)
% Get the index of a species given the name.
% The index is an integer assigned to each species in sequence
% as it is read in from the input file.
%
% Note: In keeping with the conventions used by Matlab, the
% indices start from 1 instead of 0 as in Cantera C++ and
% Python interfaces.
%
% :param name:
% String or cell array of species whose index is requested.
% :return:
% Integer number of species in the phase.
checklib;
if iscell(name)
[m, n] = size(name);
k = zeros(m, n);
for i = 1:m
for j = 1:n
k(i, j) = calllib(tp.lib, 'thermo_speciesIndex', ...
tp.tp_id, name{i, j}) + 1;
if k(i, j) > 1e3
warning(['Species ', name{i, j}, ...
' does not exist in the phase']);
k(i, j) = -1;
end
end
end
elseif ischar(name)
k = calllib(tp.lib, 'thermo_speciesIndex', ...
tp.tp_id, name) + 1;
if k > 1e3
warning(['Species ', name, ...
' does not exist in the phase.']);
k = -1;
end
else
error('Name must be either a string or cell array of strings.')
end
end
function k = elementIndex(tp, name)
% Get the index of an element given the name.
% The index is an integer assigned to each element in sequence
% as it is read in from the input file.
%
% Note: In keeping with the conventions used by Matlab, the
% indices start from 1 instead of 0 as in Cantera C++ and
% Python interfaces.
%
% :param name:
% String or cell array of elements whose index is requested
% :return:
% Integer number of elements in the phase.
checklib;
if iscell(name)
[m, n] = size(name);
k = zeros(m, n);
for i = 1:m
for j = 1:n
k(i, j) = calllib(tp.lib, 'thermo_elementIndex', ...
tp.tp_id, name{i, j}) + 1;
if k(i, j) > 1e3
warning(['Element ', name{i, j}, ...
' does not exist in the phase']);
k(i, j) = -1;
end
end
end
elseif ischar(name)
k = calllib(tp.lib, 'thermo_elementIndex', ...
tp.tp_id, name) + 1;
if k > 1e3
warning(['Element ', name, ...
' does not exist in the phase']);
k = -1;
end
else
error('name must be either a string or cell array of strings')
end
end
function n = nAtoms(tp, species, element)
% Get the number of atoms of an element in a species.
%
% :param k:
% String species name or integer species number.
% :param m:
% String element name or integer element number.
% :return:
% Integer number of atoms of the element in the species.
if nargin == 3
if ischar(species)
k = tp.speciesIndex(species);
else k = species;
end
if k < 0
n = -1;
return
end
if ischar(element)
m = tp.elementIndex(element);
else m = element;
end
if m < 0
n = -1;
return
end
n = calllib(tp.lib, 'thermo_nAtoms', tp.tp_id, k-1, m-1);
else
error('Two input arguments required.')
end
end
function moleFractions = get.X(tp)
% Get the mole fractions of all species.
%
% :return:
% Vector of species mole fractions.
checklib;
nsp = tp.nSpecies;
xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx);
calllib(tp.lib, 'thermo_getMoleFractions', ...
tp.tp_id, nsp, pt);
moleFractions = pt.Value;
% if no output argument is specified, a bar plot is produced.
if nargout == 0
figure
set(gcf, 'Name', 'Mole Fractions')
bar(moleFractions)
xlabel('Species Number')
ylabel('Mole Fraction')
title('Species Mole Fractions')
end
end
function x = moleFraction(tp, species)
% Get the mole fraction of one or a list of species.
%
% :param species:
% String or cell array of species whose mole fraction is
% requested.
% :return:
% Scalar or vector of species mole fractions.
xarray = tp.X;
if isa(species, 'char')
k = tp.speciesIndex(tp, species);
if k > 0
x = xarray(k);
else error("species not found.");
end
elseif isa(species, 'cell')
n = length(species);
x = zeros(1, n);
for j = 1:n
k = tp.speciesIndex(tp, species{j});
if k > 0
x(j) = xarray(k);
else error("species not found.");
end
end
end
end
function massFractions = get.Y(tp)
% Get the mass fractions of all species.
%
% :return:
% Vector of species mass fractions.
checklib;
nsp = tp.nSpecies;
yy = zeros(1, nsp);
pt = libpointer('doublePtr', yy);
calllib(tp.lib, 'thermo_getMassFractions', ...
tp.tp_id, nsp, pt);
massFractions = pt.Value;
% If no output argument is specified, a bar plot is produced.
if nargout == 0
figure
set(gcf, 'Name', 'Mole Fractions')
bar(massFractions)
xlabel('Species Number')
ylabel('Mole Fraction')
title('Species Mole Fractions')
end
end
function y = massFraction(tp, species)
% Get the mass fraction of one or a list of species.
%
% :param species:
% String or cell array of species whose mass fraction is
% requested.
% :return:
% Scalar or vector of species mass fractions.
yarray = tp.Y;
if isa(species, 'char')
k = tp.speciesIndex(tp, species);
if k > 0
y = yarray(k);
else error("Error: species not found.")
end
elseif isa(species, 'cell')
n = length(species);
y = zeros(1, n);
for j = 1:n
k = tp.speciesIndex(tp, species{j});
if k > 0
y(j) = yarray(k);
else error("Error: species not found.")
end
end
end
end
function mw = MolecularWeights(tp)
% Get the array of molecular weights of all species.
%
% :return:
% Vector of species molecular weights. Unit: kg/kmol
checklib;
nsp = tp.nSpecies;
yy = zeros(1, nsp);
pt = libpointer('doublePtr', yy);
calllib(tp.lib, 'thermo_getMolecularWeights', ...
tp.tp_id, nsp, pt);
mw = pt.Value;
end
function e = charges(tp)
% Get the array of species charges.
%
% :return:
% Vector of species charges. Unit: elem. charge
checklib;
nsp = tp.nSpecies;
yy = zeros(1, nsp);
pt = libpointer('doublePtr', yy);
calllib(tp.lib, 'thermo_getCharges', ...
tp.tp_id, nsp, pt);
e = pt.Value;
end
function amu = atomicMasses(tp)
% Get the atomic masses of the elements.
%
% :return:
% Vector of element atomic masses. Unit: kg/kmol
checklib;
nel = tp.nElements;
aa = zeros(1, nel);
pt = libpointer('doublePtr', aa);
calllib(tp.lib, 'thermo_getAtomicWeights', ...
tp.tp_id, nel, pt);
amu = pt.Value;
end
function nm = speciesName(tp, k) %need to solve string ptr issue
% Get the name of a species given the index.
%
% :param k:
% Scalar or array of integer species index.
% :return:
% Cell array of strings species name.
[m, n] = size(k);
nm = cell(m, n);
for i = 1:m
for j = 1:n
ksp = k(i, j) - 1;
buflen = calllib(tp.lib, 'thermo_getSpeciesName', ...
tp.tp_id, ksp, 0, '');
if buflen > 0
aa = char(zeros(1, buflen));
pt = libpointer('voidPtr', int8(aa));
pt2 = libpointer('cstring', aa);
out_buf = calllib(tp.lib, 'thermo_getSpeciesName', ...
tp.tp_id, ksp, buflen, pt);
out2_buf = calllib(tp.lib, 'thermo_getSpeciesName', ...
tp.tp_id, ksp, buflen, pt2);
out3_buf = calllib(tp.lib, 'thermo_getSpeciesName', ...
tp.tp_id, ksp, buflen, aa);
nm = pt2.Value;
end
end
end
end
%% ThermoGet single methods
function c = cv(tp)
% Get the specific heat at constant volume.
%
% :return:
% Specific heat of the mixture at constant volume depending
% on the basis. Units: J/kmol-K (molar) J/kg-K (mass).
checklib;
if strcmp(tp.basis, 'molar')
c = calllib(tp.lib, 'thermo_cv_mole', tp.tp_id);
elseif strcmp(tp.basis, 'mass')
c = calllib(tp.lib, 'thermo_cv_mass', tp.tp_id);
else error("basis not specified");
end
end
function c = cp(tp)
% Get the specific heat at constant pressure.
%
% :return:
% Specific heat of the mixture at constant pressure depending
% on the basis. Units: J/kmol-K (molar) J/kg-K (mass).
checklib;
if strcmp(tp.basis, 'molar')
c = calllib(tp.lib, 'thermo_cp_mole', tp.tp_id);
elseif strcmp(tp.basis, 'mass')
c = calllib(tp.lib, 'thermo_cp_mass', tp.tp_id);
else error("basis not specified");
end
end
function enthalpy = get.H(tp)
% Get the enthalpy.
%
% :return:
% Enthalpy of the mixture depending on the basis.
% Units: J/kmol (molar) J/kg (mass).
checklib;
if strcmp(tp.basis, 'molar')
enthalpy = calllib(tp.lib, 'thermo_enthalpy_mole', tp.tp_id);
elseif strcmp(tp.basis, 'mass')
enthalpy = calllib(tp.lib, 'thermo_enthalpy_mass', tp.tp_id);
else error("basis not specified");
end
end
function entropy = get.S(tp)
% Get the entropy.
%
% :return:
% Entropy of the mixture depending on the basis.
% Units: J/kmol-K (molar) J/kg-K (mass).
checklib;
if strcmp(tp.basis, 'molar')
entropy = calllib(tp.lib, 'thermo_entropy_mole', tp.tp_id);
elseif strcmp(tp.basis, 'mass')
entropy = calllib(tp.lib, 'thermo_entropy_mass', tp.tp_id);
else error("basis not specified");
end
end
function intEnergy = get.U(tp)
% Get the internal energy.
%
% :return:
% Internal energy of the mixture depending on the basis.
% Units: J/kmol (molar) J/kg (mass).
checklib;
if strcmp(tp.basis, 'molar')
intEnergy = calllib(tp.lib, 'thermo_intEnergy_mole', tp.tp_id);
elseif strcmp(tp.basis, 'mass')
intEnergy = calllib(tp.lib, 'thermo_intEnergy_mass', tp.tp_id);
else error("basis not specified");
end
end
function gibbs = G(tp)
% Get the Gibss free energy.
%
% :return:
% Gibbs free energy of the mixture depending on the basis.
% Units: J/kmol (molar) J/kg (mass).
checklib;
if strcmp(tp.basis, 'molar')
gibbs = calllib(tp.lib, 'thermo_gibbs_mole', tp.tp_id);
elseif strcmp(tp.basis, 'mass')
gibbs = calllib(tp.lib, 'thermo_gibbs_mass', tp.tp_id);
else error("basis not specified");
end
end
function t = minTemp(tp)
% Get the minimum temperature of the parameter fits.
%
% :return:
% Vector of minimum temperatures of all species.
checklib;
t = calllib(tp.lib, 'thermo_minTemp', tp.tp_id, -1);
end
function t = maxTemp(tp)
% Get the maximum temperature of the parameter fits.
%
% :return:
% Vector of maximum temperatures of all species.
checklib;
t = calllib(tp.lib, 'thermo_maxTemp', tp.tp_id, -1);
end
function p = refPressure(tp)
% Get the reference pressure.
%
% :return:
% Reference pressure. Unit: Pa.
checklib;
p = calllib(tp.lib, 'thermo_refPressure', tp.tp_id, -1);
end
function t = critTemperature(tp)
% Get the critical temperature.
%
% :return:
% Critical temperature. Unit: K.
checklib;
t = calllib(tp.lib, 'thermo_critTemperature', tp.tp_id);
end
function p = critPressure(tp)
% Get the critical pressure.
%
% :return:
% Critical temperature. Unit: Pa.
checklib;
p = calllib(tp.lib, 'thermo_critPressure', tp.tp_id);
end
function d = critDensity(tp)
% Get the critical density.
%
% :return:
% Critical density. Unit: K.
checklib;
d = calllib(tp.lib, 'thermo_critDensity', tp.tp_id);
end
function v = vaporFraction(tp)
% Get the vapor fractions.
%
% :return:
% Vapor fraction.
checklib;
v = calllib(tp.lib, 'thermo_vaporFraction', tp.tp_id);
end
function t = T_sat(tp, p)
% Get the saturation temperature for a given pressure.
%
% :param p:
% Pressure. Unit: Pa.
% :return:
% Saturation temperature for pressure p. Unit: K.
checklib;
t = calllib(tp.lib, 'thermo_satTemperature', tp.tp_id, p);
end
function p = P_sat(tp, t)
% Get the saturation pressure for a given temperature.
%
% :param t:
% Temperature. Unit: K.
% :return:
% Saturation pressure for temperature t. Unit: Pa.
checklib;
p = calllib(tp.lib, 'thermo_satPressure', tp.tp_id, t);
end
function v = electricPotential(tp)
% Get the electric potential
%
% :return:
% Electric potential of the phase. Unit: V.
checklib;
v = calllib(tp.lib, 'thermo_electricPotential', tp.tp_id);
end
function b = isothermalCompressibility(tp)
% Get the isothermal compressibility
%
% :return:
% Isothermal compressibility. Unit: 1/Pa.
checklib;
b = calllib(tp.lib, 'thermo_isothermalCompressibility', tp.tp_id);
end
function a = thermalExpansionCoeff(tp)
% Get the thermal expansion coefficient.
%
% :return:
% Thermal expansion coefficient. Unit: 1/K.
checklib;
a = calllib(tp.lib, 'thermo_thermalExpansionCoeff', tp.tp_id);
end
function mu = chemical_potentials(tp)
% Get the chemical potentials of the species.
%
% :return:
% Vector of species chemical potentials. Unit: J/kmol.
checklib;
nsp = tp.nSpecies;
xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx);
calllib(tp.lib, 'thermo_chemPotentials', ...
tp.tp_id, nsp, pt);
moleFractions = pt.Value;
end
%% ThermoGet multi methods
function output = get.DP(tp)
% Get density and pressure depending on the basis.
% :return:
% Density. Unit: kmol/m^3 (molar) kg/m^3 (mass)
% Pressure. Unit: Pa
output = [tp.D, tp.P];
end
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];
end
function output = get.DPY(tp)
% Get density, pressure, and mass fractions depending on the basis.
% :return:
% Density. Unit: kmol/m^3 (molar) kg/m^3 (mass).
% Pressure. Unit: Pa.
% Mass fractions of all species.
output = [tp.D, tp.P, tp.Y];
end
function output = get.HP(tp)
% Get enthalpy and pressure depending on the basis.
% :return:
% Enthalpy. Unit: J/kmol (molar) J/kg (mass).
% Pressure. Unit: Pa.
output = [tp.H, tp.P];
end
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];
end
function output = get.HPY(tp)
% Get enthalpy, pressure, and mole fractions depending on the basis.
% :return:
% Enthalpy. Unit: J/kmol (molar) J/kg (mass).
% Pressure. Unit: Pa.
% Mass fractions of all species.
output = [tp.H, tp.P, tp.Y];
end
function output = get.SP(tp)
% Get entropy, pressure, and mole fractions depending on the basis.
% :return:
% Entropy. Unit: J/kmol-K (molar) J/kg-K (mass).
% Pressure. Unit: Pa.
output = [tp.S, tp.P];
end
function output = get.SPX(tp)
% Get entropy, pressure, and mole fractions depending on the basis.
% :return:
% Entropy. Unit: J/kmol-K (molar) J/kg-K (mass).
% Pressure. Unit: Pa.
% Mole fractions of all species.
output = [tp.S, tp.P, tp.X];
end
function output = get.SPY(tp)
% Get entropy, pressure, and mass fractions depending on the basis.
% :return:
% Entropy. Unit: J/kmol-K (molar) J/kg-K (mass).
% Pressure. Unit: Pa.
% Mass fractions of all species.
output = [tp.H, tp.P, tp.Y];
end
function output = get.TD(tp)
% Get temperature and density depending on the basis.
% :return:
% Temperature. Unit: K.
% Density. Unit: kmol/m^3 (molar) kg/m^3 (mass).
output = [tp.T, tp.D];
end
function output = get.TDX(tp)
% Get entropy, pressure, and mole fractions depending on the basis.
% :return:
% Temperature. Unit: K.
% Density. Unit: kmol/m^3 (molar) kg/m^3 (mass).
% Mole fractions of all species.
output = [tp.T, tp.D, tp.X];
end
function output = get.TDY(tp)
% Get entropy, pressure, and mass fractions depending on the basis.
% :return:
% Temperature. Unit: K.
% Density. Unit: kmol/m^3 (molar) kg/m^3 (mass).
% Mass fractions of all species.
output = [tp.T, tp.D, tp.Y];
end
function output = get.TP(tp)
% Get temperature and density depending on the basis.
% :return:
% Temperature. Unit: K.
% Pressure. Unit: Pa.
output = [tp.T, tp.P];
end
function output = get.TPX(tp)
% Get entropy, pressure, and mole fractions depending on the basis.
% :return:
% Temperature. Unit: K.
% Pressure. Unit: Pa.
% Mole fractions of all species.
output = [tp.T, tp.P, tp.X];
end
function output = get.TPY(tp)
% Get entropy, pressure, and mass fractions depending on the basis.
% :return:
% Temperature. Unit: K.
% Pressure. Unit: Pa.
% Mass fractions of all species.
output = [tp.T, tp.P, tp.Y];
end
%% PhaseSet Single Methods
function tp = set.T(tp, temperature)
checklib;
if temperature <= 0
error('The temperature must be positive');
end
calllib(tp.lib, 'thermo_setTemperature', tp.tp_id, temperature);
end
function tp = set.P(tp, pressure)
checklib;
if pressure <= 0
error('The pressure must be positive');
end
calllib(tp.lib, 'thermo_setPressure', tp.tp_id, pressure);
end
function tp = set.D(tp, density)
checklib;
if density <= 0
error('The density must be positive');
end
calllib(tp.lib, 'thermo_setDensity', tp.tp_id, density);
end
% function tp = set.X(tp,
function tp = equilibrate(tp, xy, solver, rtol, maxsteps, ...
maxiter, loglevel)
checklib;
% use the ChemEquil solver by default
if nargin < 3
solver = -1;
end
if nargin < 4
rtol = 1.0e-9;
end
if nargin < 5
maxsteps = 1000;
end
if nargin < 6
maxiter = 100;
end
if nargin < 7
loglevel = 0;
end
calllib(tp.lib, 'thermo_equilibrate', tp.tp_id, xy, solver, ...
rtol, maxsteps, maxiter, loglevel);
end
end
end

178
interfaces/matlab_experimental/Base/Transport.m Normal file
View File

@@ -0,0 +1,178 @@
classdef Transport < handle
properties
th
tr_id
end
properties(Constant = true)
lib = 'cantera_shared'
end
methods
%% Transport class constructor
function tr = Transport(tp, model, loglevel)
checklib;
tr.tr_id = 0;
if nargin == 2
model = 'default'
end
if nargin < 3
loglevel = 4;
end
if ~isa(tp, 'ThermoPhase')
error(['The first argument must be an ', ...
'instance of class ThermoPhase']);
else
tr.th = tp;
if strcmp(model, 'default')
tr.tr_id = calllib(tr.lib, 'trans_newDefault', ...
tp.tp_id, loglevel);
else
tr.tr_id = calllib(tr.lib, 'trans_new', model, ...
tp.tp_id, loglevel);
end
end
end
%% Utility methods
function tr_clear(tr)
% Delete the kernel object.
checklib;
calllib(tr.lib, 'trans_del', tr.tr_id);
end
%% Transport Methods
function v = viscosity(tr)
% Get the dynamic viscosity.
%
% :return:
% Double dynamic viscosity. Unit: Pa*s.
checklib;
v = calllib(tr.lib, 'trans_viscosity', tr.tr_id);
if v == -1.0
error(geterr);
elseif v < 0.0
error('exception raised');
end
end
function v = thermaoConductivity(tr)
% Get the thermal conductivity.
%
% :return:
% Double thermal conductivity. Unit: W/m-K.
checklib;
v = calllib(tr.lib, 'trans_thermalConductivity', tr.tr_id);
if v == -1.0
error(geterr);
elseif v < 0.0
error('exception raised');
end
end
function v = electricalConductivity(tr)
% Get the electrical conductivity.
%
% :return:
% Double electrical conductivity. Unit: S/m.
checklib;
v = calllib(tr.lib, 'trans_electricalConductivity', tr.tr_id);
if v == -1.0
error(geterr);
elseif v < 0.0
error('exception raised');
end
end
function v = mixDiffCoeffs(tr)
% Get the mixture-averaged diffusion coefficients.
%
% :return:
% Vector of length nSpecies with the mixture-averaged
% diffusion coefficients. Unit: m^2/s.
checklib;
nsp = tr.th.nSpecies;
xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx);
calllib(tr.lib, 'trans_getMixDiffCoeffs', tr.tr_id, nsp, pt);
v = pt.Value;
end
function v = thermalDiffCoeffs(tr)
% Get the thermal diffusion coefficients.
%
% :return:
% Vector of length nSpecies with the thermal diffusion
% coefficients.
checklib;
nsp = tr.th.nSpecies;
xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx);
calllib(tr.lib, 'trans_getThermalDiffCoeffs', tr.tr_id, nsp, pt);
v = pt.Value;
end
function v = binDiffCoeffs(tr)
% Get the binary diffusion coefficients.
%
% :return:
% A matrix of binary diffusion coefficients. The matrix is
% symmetric: d(i, j) = d(j, i). Unit: m^2/s.
checklib;
nsp = tr.th.nSpecies;
xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx);
calllib(tr.lib, 'trans_getBinDiffCoeffs', tr.tr_id, nsp, pt);
v = pt.Value;
end
function v = multiDiffCoeffs(tr)
% Get the multicomponent diffusion coefficients.
%
% :return:
% Vector of length nSpecies with the multicomponent
% diffusion coefficients. Unit: m^2/s.
checklib;
nsp = tr.th.nSpecies;
xx = zeros(1, nsp);
pt = libpointer('doublePtr', xx);
calllib(tr.lib, 'trans_getMultiDiffCoeffs', tr.tr_id, nsp, pt);
v = pt.Value;
end
function setParameters(tr, type, k, p)
% Set the parameters.
%
% :param type:
% :param k:
% :param p:
checklib;
calllib(tr.lib, 'trans_setParameters', tr.tr_id, type, k, p);
end
function setThermalConductivity(tr, lam)
% Set the thermal conductivity.
%
% :param lam:
% Thermal conductivity in W/(m-K).
checklib;
setParameters(tr, 1, 0, lam);
end
end
end

4
interfaces/matlab_experimental/Base/Water.m Normal file
View File

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

275
interfaces/matlab_experimental/Func/Func.m Normal file
View File

@@ -0,0 +1,275 @@
classdef Func < handle
properties
f1
f2
coeffs
index
typ
end
properties(Constant = true)
lib = 'cantera_shared'
end
methods
%% Functor class constructor
function x = Func(typ, n, p)
% Functor class constructor.
%
% A functor is an object that behaves like a function. Cantera
% defines a set of functors to use to create arbitrary
% functions to specify things like heat fluxes, piston speeds,
% etc., in reactor network simulations.
%
% The main feature of a functor class is that it overloads the
% '()' operator to evaluate the function. For example, suppose
% object 'f' is a functor that evaluates the polynomial
% :math:'2x^2 - 3x + 1'. Then writing 'f(2)' would cause the
% method that evaluates the function to be invoked, and would
% pass it the argument '2'. The return value would be 3.
%
% The types of functors you can create in Cantera are these:
% 1. A polynomial;
% 2. A Fourier series;
% 3. A sum of Arrhenius terms;
% 4. A Gaussian.
%
% You can also create composite functors by adding,
% multiplying, or dividing these basic functors or other
% composite functors.
%
% Note: this MATLAB class shadows the underlying C++ Cantera
% class 'Func1'. See the Cantera C++ documentation for more
% details.
%
% :param typ:
% String indicating type of functor to create. Possible
% values are:
% * 'polynomial'
% * 'fourier'
% * 'gaussian'
% * 'arrhenius'
% * 'sum'
% * 'diff'
% * 'ratio'
% * 'composite'
% * 'periodic'
%
% :param n:
% Number of parameters required for the functor
% :param p:
% Vector of parameters
% :return:
% Instance of class :mat:func:`Func`
checklib;
if ~isa(typ, 'char')
error('Function type must be a string');
end
x.f1 = 0;
x.f2 = 0;
x.coeffs = 0;
function nn = new_func(itype, n, p)
% helper function to pass the correct parameters to the C
% library
if itype < 20
ptr = libpointer('doublePtr', p);
[m, n] = size(p);
lenp = m * n;
nn = calllib(x.lib, 'func_new', type, n, lenp, ptr);
elseif itype < 45
m = n;
nn = calllib(x.lib, 'func_new', type, n, m, 0);
else
ptr = libpointer('doublePtr', p);
nn = calllib(x.lib, 'func_new', type, n, 0, ptr);
end
end
if itype > 0
x.coeffs = p;
x.index = new_func(itype, n, p);
elseif strcmp(typ, 'periodic')
itype = 50;
x.f1 = n;
x.coeffs = p;
x.index = new_func(itype, n.index, p);
else
if strcmp(typ, 'sum')
itype = 20;
elseif strcmp(typ, 'diff')
itype = 25;
elseif strcmp(typ, 'prod')
itype = 30;
elseif strcmp(typ, 'ratio')
itype = 40;
elseif strcmp(typ, 'composite')
itype = 60;
end
x.f1 = n;
x.f2 = p;
x.index = new_func(itype, n.index, p.index);
end
x.typ = typ;
end
%% Utility methods
function func_clear(f)
% Clear the specified flow device from memory.
checklib;
calllib(f.lib, 'func_del', f.index);
end
function display(a)
% Display the equation of the input function on the terminal.
%
% :param a:
% Instance of class 'Func'
disp(' ');
disp([inputname(1),' = '])
disp(' ');
disp([' ' char(a)])
disp(' ');
end
%% Functor methods
function r = plus(a, b)
% Get a functor representing the sum two input functors 'a' and
% 'b'.
r = Func('sum', a, b);
end
function r = rdivide(a, b)
% Get a functor that is the ratio of the input functors 'a' and
% 'b'.
r = Func('ratio', a, b);
end
function r = times(a, b)
% Get a functor that multiplies two functors 'a' and 'b'
r = Func('prod', a, b);
end
function b = subsref(a, s)
% Redefine subscripted references for functors.
%
% :param a:
% Instance of class 'Func'
% :param s:
% Value at which the function should be evaluated.
% :return:
% The value of the function evaluated at 's'.
checklib;
if strcmp(s.type, '()')
ind= s.subs{:};
b = zeros(1, length(ind));
for k = 1:length(ind)
b(k) = calllib(a.lib, 'func_value', a.index, ind(k));
end
else error('Specify value for x as p(x)');
end
end
function s = char(p)
% Get the formatted string to display the function.
if strcmp(p.typ,'sum')
s = ['(' char(p.f1) ') + (' char(p.f2) ')'];
elseif strcmp(p.typ,'diff')
s = ['(' char(p.f1) ') - (' char(p.f2) ')'];
elseif strcmp(p.typ,'prod')
s = ['(' char(p.f1) ') * (' char(p.f2) ')'];
elseif strcmp(p.typ,'ratio')
s = ['(' char(p.f1) ') / (' char(p.f2) ')'];
elseif all(p.coeffs == 0)
s = '0';
else
if strcmp(p.typ,'polynomial')
d = length(p.coeffs) - 1;
s = [];
nn = 0;
for b = p.coeffs;
cc(d+1-nn) = b;
nn = nn + 1;
end
for a = cc;
if a ~= 0;
if ~isempty(s)
if a > 0
s = [s ' + '];
else
s = [s ' - '];
a = -a;
end
end
if a ~= 1 || d == 0
s = [s num2str(a)];
if d > 0
s = [s '*'];
end
end
if d >= 2
s = [s 'x^' int2str(d)];
elseif d == 1
s = [s 'x'];
end
end
d = d - 1;
end
elseif strcmp(p.typ, 'gaussian')
s = ['Gaussian(' num2str(p.coeffs(1)) ',' ...
num2str(p.coeffs(2)) ',' ...
num2str(p.coeffs(3)) ')'];
elseif strcmp(p.typ, 'fourier')
c = reshape(p.coeffs, [], 2);
Ao = c(1, 1);
w = c(1, 2);
A = c(2:end, 1);
B = c(2:end, 2);
N = size(c, 1) - 1;
if Ao ~= 0
s = num2str(Ao/2);
else
s = '';
end
for n=1:N
if A(n) ~= 0
if A(n) < 0
prefix = ' - ';
elseif s
prefix = ' + ';
else
prefix = '';
end
s = [s prefix num2str(abs(A(n))), ...
'*cos(' num2str(n*w) '*x)'];
end
if B(n) ~= 0
if B(n) < 0
prefix = ' - ';
elseif s
prefix = ' + ';
else
prefix = '';
end
s = [s prefix num2str(abs(B(n))),...
'*sin(' num2str(n*w) '*x)'];
end
end
else
s = ['*** char not yet implemented for' p.typ ' ***'];
end
end
end
end
end

15
interfaces/matlab_experimental/LoadCantera.m Normal file
View File

@@ -0,0 +1,15 @@
function Load_Cantera
addpath('Class', 'Utility', 'PresetObjects');
if ~libisloaded('cantera_shared')
load('Utility/cantera_root.mat');
[~,warnings] = loadlibrary([cantera_root '/lib/cantera_shared.dll'], ...
[cantera_root '/include/cantera/clib/ctmatlab.h'], ...
'includepath', [cantera_root '/include'], ...
'addheader','ct','addheader','ctfunc', ...
'addheader','ctmultiphase','addheader', ...
'ctonedim','addheader','ctreactor', ...
'addheader','ctrpath','addheader','ctsurf', ...
'addheader','ctxml');
end
disp('Cantera is ready for use');
end

157
interfaces/matlab_experimental/Reactor/FlowDevice.m Normal file
View File

@@ -0,0 +1,157 @@
classdef FlowDevice < handle
properties
type
index
upstream
downstream
end
properties(Constant = true)
lib = 'cantera_shared'
end
methods
%% FlowDevice class constructor
function x = FlowDevice(typ)
% Flow Device class constructor.
%
% :param typ:
% Type of flow device to be created. Type =
% 'MassFlowController', 'PressureController' or 'Valve'.
checklib;
if nargin == 0
typ = 'MassFlowController';
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
x.type = typ;
x.index = calllib(x.lib, 'flowdev_new2', typ);
% if x.index < 0
% error(geterr);
% end
x.upstream = -1;
x.downstream = -1;
end
%% Utility methods
function flow_clear(f)
% Clear the specified flow device from memory.
checklib;
calllib(f.lib, 'flowdev_del', f.index);
end
%% FlowDevice methods
function install(f, upstream, downstream)
% Install a flow device between reactors or reservoirs.
%
% :param upstream:
% Upsteram 'Reactor' or 'Reservoir'.
% :param downstream:
% Downstream 'Reactor' or 'Reservoir'.
checklib;
if nargin == 3
if ~isa(upstream, 'Reactor') || ~isa(downstream, 'Reactor')
error(['Flow devices can only be installed between',...
'reactors or reservoirs']);
end
i = upstream.id;
j = downstream.id;
ok = calllib(f.lib, 'flowdev_install', f.index, i, j);
% if ok < 0
% error(geterr)
% end
else error('install requires 3 arguments');
end
end
function mdot = massFlowRate(f, time)
% Get the mass flow rate at a given time.
%
% :param time:
% Time at which the mass flow rate is desired.
% :return:
% The mass flow rate through the flow device at the given
% time.
checklib;
if nargin == 1
mdot = calllib(f.lib, 'flowdev_massFlowRate2', f.id);
else
warning(['"Time" argument to massFlowRate is deprecated', ...
'and will be removed after Cantera 2.5.']);
mdot = calllib(f.lib, 'flowdev_massFlowRate', f.id, time);
end
end
function setFunction(f, mf)
% Set the time function with class 'func'.
%
% :param mf:
% Instance of class 'func'.
checklib;
if strcmp(f.type, 'MassFlowController')
k = calllib(f.lib, 'flowdev_setTimeFunction', f.index, ...
mf.id);
% if k < 0
% error(geterr);
% end
else
error('Mass flow rate can only be set for mass flow controllers.');
end
end
function setMassFlowRate(f, mdot)
% Set the mass flow rate to a constant value.
%
% :param mdot:
% Mass flow rate
checklib;
if strcmp(f.type, 'MassFlowController')
k = calllib(f.lib, 'flowdev_setMassFlowCoeff', f.index, mdot);
% if k < 0
% error(geterr);
% end
else
error('Mass flow rate can only be set for mass flow controllers.');
end
end
function setValveCoeff(f, k)
% Set the valve coefficient 'K'.
%
% The mass flow rate [kg/s] is computed from the expression
% mdot = K(P_upstream - P_downstream)
% as long as this produces a positive value. If thsi expression
% is negative, zero is returned.
%
% :param k:
% Value fo the valve coefficient. Unit: kg/Pa-s.
checklib;
if ~strcmp(f.type, 'Valve')
error('Valve coefficient can only be set for valves.');
end
ok = calllib(f.lib, 'flowdev_setValveCoeff', f.index, k);
% if k < 0
% error(geterr);
% end
end
end
end

3
interfaces/matlab_experimental/SetCanteraPath.m Normal file
View File

@@ -0,0 +1,3 @@
function SetCanteraPath(cantera_root)
save('Utility/cantera_root.mat','cantera_root');
end

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

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

3
interfaces/matlab_experimental/Utility/checklib.m Normal file
View File

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

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

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

12
interfaces/matlab_experimental/cleanup.m Normal file
View File

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