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

More fixes to new MatlabToolbox

Browse Source
This commit is contained in:
ssun30 2022-04-22 15:38:12 -04:00 committed by Ray Speth
parent 71fc640f56
commit 2318dd9c21
15 changed files with 382 additions and 103 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -2,8 +2,9 @@ classdef Domain1D < handle
properties
domainID
domainType
type
T
P
end
methods
@ -89,17 +90,64 @@ classdef Domain1D < handle
% if d.domainID < 0
% error(geterr);
% end
d.domainType = a;
d.type = a;
end
%% Utility Methods
function domClear(d)
function clear(d)
% Delete the Domain1D object
calllib(ct, 'domain_del', d.domainID);
end
%% Domain Methods
function d = disableEnergy(d)
% Disable the energy equation.
disp(' ');
disp('Disabling the energy equation...');
calllib(ct, 'stflow_solveEnergyEqn', d.domainID, 0);
end
function d = enableEnergy(d)
% Enable the energy equation.
disp(' ');
disp('Enabling the energy equation...');
calllib(ct, 'stflow_solveEnergyEqn', d.domainID, 1);
end
function d = disableSoret(d)
% Disable the diffusive mass fluxes due to the Soret effect.
disp(' ');
disp('Disabling the Soret effect...');
calllib(ct, 'stflow_enableSoret', d.domainID, 0);
end
function d = enableSoret(d)
% Enable the diffusive mass fluxes due to the Soret effect.
disp(' ');
disp('Disabling the Soret effect...');
calllib(ct, 'stflow_enableSoret', d.domainID, 1);
end
%% Domain Get Methods
function b = bounds(d, component)
% Return the (lower, upper) bounds for a solution component.
%
% :param component:
% String name of the component for which the bounds are
% returned.
% :return:
% 1x2 Vector of the lower and upper bounds.
n = d.componentIndex(component);
lower = calllib(ct, 'domain_lowerBound', d.domainID, n);
upper = calllib(ct, 'domain_upperBound', d.domainID, n);
b = [lower, upper];
end
function n = componentIndex(d, name)
% Get the index of a component given its name
@ -151,17 +199,6 @@ classdef Domain1D < handle
end
end
function d = disableEnergy(d)
% Disable the energy equation.
%
% :parameter d:
% Instance of class 'Domain1D'.
disp(' ');
disp('Disabling the energy equation...');
calllib(ct, 'stflow_solveEnergyEqn', d.domainID, 0);
end
function i = domainIndex(d)
% Get the domain index.
% :return:
@ -183,17 +220,7 @@ classdef Domain1D < handle
% This function :returns an integer flag denoting the domain
% type.
i = calllib(ct, 'domainType', d.domainID);
end
function d = enableEnergy(d)
% Enable the energy equation.
% :parameter d:
% Instance of class 'Domain1D'.
disp(' ');
disp('Enabling the energy equation...');
calllib(ct, 'stflow_solveEnergyEqn', d.domainID, 1);
i = calllib(ct, 'domain_type', d.domainID);
end
function zz = gridPoints(d, n)
@ -312,6 +339,54 @@ classdef Domain1D < handle
n = calllib(ct, 'domain_nPoints', d.domainID);
end
function tol = tolerances(d, component)
% Return the (relative, absolute) error tolerances for a
% solution component.
%
% :param component:
% String name of the component for which the bounds are
% returned.
% :return:
% 1x2 Vector of the relative and absolute error tolerances.
n = d.componentIndex(component);
rerr = calllib(ct, 'domain_rtol', d.domainID, n);
aerr = calllib(ct, 'domain_atol', d.domainID, n);
tol = [rerr, aerr];
end
function temperature = get.T(d)
% Get the boundary temperature (K).
temperature = calllib(ct, 'bdry_temperature', d.domainID);
end
function pressure = get.P(d)
% Get the pressure (Pa).
pressure = calllibt(ct, 'stflow_pressure', d.domainID);
end
%% Domain Set Methods
function set.T(d, t)
% Set the temperature (K).
if t <= 0
error('The temperature must be positive');
end
calllib(ct, 'bdry_setTemperature', d.domainID, t);
end
function set.P(d, p)
% Set the pressure (K).
if p <= 0
error('The pressure must be positive');
end
calllib(ct, 'stflow_setPressure', d.domainID, p);
end
function setBounds(d, component, lower, upper)
% Set bounds on the solution components.
% :parameter d:
@ -405,14 +480,6 @@ classdef Domain1D < handle
calllib(ct, 'bdry_setMoleFractions', d.domainID, x);
end
function setPressure(d, p)
% Set the pressure.
% :parameter p:
% Pressure to be set. Unit: Pa.
calllib(ct, 'stflow_setPressure', d.domainID, p);
end
function setProfileD(d, n, p)
% Set the profile of a component.
% Convenience function to allow an instance of 'Domain1D' to
@ -502,19 +569,15 @@ classdef Domain1D < handle
end
end
function temperature = get.T(d)
% Get the boundary temperature (K).
function setTransport(d, itr)
% Set the solution object used for calculating transport
% properties.
%
% :param itr:
% ID of the solution object for which transport properties
% are calculated.
temperature = calllib(ct, 'bdry_temperature', d.domainID);
end
function set.T(d, t)
% Set the temperature (K).
if t <= 0
error('The temperature must be positive');
end
calllib(ct, 'bdry_setTemperature', d.domainID, t);
calllib(ct, 'stflow_setTransport', d.domainID, itr);
end
function setupGrid(d, grid)

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

@ -40,7 +40,7 @@ classdef Stack < handle
%% Utility Methods
function st_clear(s)
function clear(s)
% Delete the Sim1D object
calllib(ct, 'sim1D_del', s.stID);
@ -86,6 +86,12 @@ classdef Stack < handle
end
end
function getInitialSoln(s)
% Get the initial solution.
calllib(ct, 'sim1D_getInitialSoln', s.stID);
end
function z = grid(s, name)
% Get the grid in one domain.
%
@ -194,6 +200,19 @@ classdef Stack < handle
calllib(ct, 'sim1D_save', s.stID, fname, id, desc);
end
function setFixedTemperature(s, T)
% Set the temperature used to fix the spatial location of a
% freely propagating flame.
%
% :parameter T:
% Double Temperature to be set. Unit: K.
if T <= 0
error('temperature must be positive');
end
calllib(ct, 'sim1D_setFixedTemperature', s.stID, T);
end
function setFlatProfile(s, domain, comp, v)
% Set a component to a value across the entire domain.
%
@ -210,6 +229,17 @@ classdef Stack < handle
domain - 1, comp - 1, v);
end
function setGridMin(s, domain, gridmin)
% Set the minimum grid spacing on domain.
%
% :parameter domain:
% Integer ID of the domain.
% :parameter gridmin:
% Double minimum grid spacing.
calllib(ct, 'sim1D_setGridMin', s.stID, domain-1, gridmin);
end
function setMaxJacAge(s, ss_age, ts_age)
% Set the number of times the Jacobian will be used before it
% is recomputed.
@ -406,7 +436,7 @@ classdef Stack < handle
end
end
function solve(s, loglevel, refine_grid)
function solve(s, loglevel, refineGrid)
% Solve the problem.
%
% :parameter s:
@ -418,7 +448,7 @@ classdef Stack < handle
% :parameter refine_grid:
% Integer, 1 to allow grid refinement, 0 to disallow.
calllib(ct, 'sim1D_solve', s.stID, loglevel, refine_grid);
calllib(ct, 'sim1D_solve', s.stID, loglevel, refineGrid);
end
function writeStats(s)

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

@ -2,6 +2,7 @@ classdef Interface < handle & ThermoPhase & Kinetics
properties
coverages
siteDensity
end
methods
@ -51,6 +52,17 @@ classdef Interface < handle & ThermoPhase & Kinetics
c = pt.Value;
end
function d = get.siteDensity(s)
% Get the site density.
%
% :return:
% Double site density. Unit: kmol/m^2 for surface phases,
% kmol/m for edge phases.
surfID = s.tpID;
d = calllibt(ct, 'surf_siteDensity', surfID);
end
function c = concentrations(s)
% Get the concentrations of the species on an interface.
%
@ -79,7 +91,7 @@ classdef Interface < handle & ThermoPhase & Kinetics
norm_flag = 1;
end
surfID = s.tr_id;
surfID = s.tpID;
nsp = s.nSpecies;
[m, n] = size(cov);
@ -98,5 +110,17 @@ classdef Interface < handle & ThermoPhase & Kinetics
end
end
function set.siteDensity(s, d)
% Set surface site densities.
%
% :parameter density:
% Double site density. Unit: kmol/m^2 for surface phases,
% kmol/m for edge phases.
surfID = s.tpID;
calllib(ct, 'surf_setSiteDensity', surfID, d);
end
end
end

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

@ -364,6 +364,20 @@ classdef Kinetics < handle
enthalpy = pt.Value;
end
function enthalpy = get.dHss(kin)
% Get the standard state enthalpy of reaction for each reaction.
%
% :return:
% A vector of the standard state enthalpy of reaction for each reaction.
% Unit: J/kmol.
nr = kin.nReactions;
xx = zeros(1, nr);
pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getDelta', kin.kinID, 3, nr, pt);
enthalpy = pt.Value;
end
function entropy = get.dS(kin)
% Get the entropy of reaction for each reaction.
%
@ -378,6 +392,20 @@ classdef Kinetics < handle
entropy = pt.Value;
end
function entropy = get.dSss(kin)
% Get the standard state entropy of reaction for each reaction.
%
% :return:
% A vector of the standard state entropy of reaction for each reaction.
% Unit: J/kmol-K.
nr = kin.nReactions;
xx = zeros(1, nr);
pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getDelta', kin.kinID, 5, nr, pt);
entropy = pt.Value;
end
function gibbs = get.dG(kin)
% Get the Gibbs free energy of reaction for each reaction.
%
@ -392,6 +420,20 @@ classdef Kinetics < handle
gibbs = pt.Value;
end
function gibbs = get.dGss(kin)
% Get the standard state Gibbs free energy of reaction for each reaction.
%
% :return:
% A vector of the standard state Gibbs free energy of reaction for each
% reaction. Unit: J/kmol.
nr = kin.nReactions;
xx = zeros(1, nr);
pt = libpointer('doublePtr', xx);
calllib(ct, 'kin_getDelta', kin.kinID, 4, nr, pt);
gibbs = pt.Value;
end
function k = get.Kc(kin)
% Get the equilibrium constants for all reactions.
%

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

@ -138,10 +138,15 @@ classdef Mixture < handle
pressure = calllib(ct, 'mix_pressure', m.mixID);
end
function n = nPhases(m)
% Get the number of phases in the mixture.
function n = nAtoms(m, k, e)
% Get the number of atoms of element e in species k.
%
% Note: In keeping with the conventions used by Matlab, the
% indices start from 1 instead of 0 as in Cantera C++ and
% Python interfaces.
n = calllib(ct, 'mix_nPhases', m.mixID);
n = calllib(ct, 'mix_nPhases', m.mixID, k-1, e-1);
% Check back on this one!
end
function n = nElements(m)
@ -150,6 +155,12 @@ classdef Mixture < handle
n = calllib(ct, 'mix_nElements', m.mixID);
end
function n = nPhases(m)
% Get the number of phases in the mixture.
n = calllib(ct, 'mix_nPhases', m.mixID);
end
function n = nSpecies(m)
% Get the number of species in the mixture.
@ -167,23 +178,47 @@ classdef Mixture < handle
end
function n = speciesIndex(m, k, p)
% Get the index of element 'name'.
% Get the index of species k in phase p.
%
% :param k:
% Double
% Note: In keeping with the conventions used by Matlab, the
% indices start from 1 instead of 0 as in Cantera C++ and
% Python interfaces.
n = calllib(ct, 'mix_speciesIndex', m.mixID, k-1, p-1) + 1;
% check back on this one!
% Check back on this one!
end
function moles = elementMoles(m, e)
% Get the number of moles of an element in a mixture.
%
% :parameter e:
% Integer element number.
% :return:
% Moles of element number 'e'. If input 'e' is empty, return
% moles of every element in the mixture. Unit: kmol.
if nargin == 2
moles = calllib(ct, 'mix_elementMoles', m.mixID, e)
elseif nargin == 1
nel = m.nElements;
moles = zeros(1, nel);
for i = 1:nel
moles(i) = calllib(ct, 'mix_elementMoles', m.mixID, i);
end
else error('wrong number of arguments');
end
end
function moles = phaseMoles(m, n)
% Get the number of moles of a phase in a mixture.
%
% :parameter n:
% Integer phase number in the input.
% Integer phase number.
% :return:
% Moles of phase number 'n'. Unit: kmol.
% Moles of phase number 'n'. If input 'n' is empty, return
% moles of phase element in the mixture. Unit: kmol.
if nargin == 2
moles = calllib(ct, 'mix_phaseMoles', m.mixID, n);
@ -198,6 +233,28 @@ classdef Mixture < handle
end
end
function moles = speciesMoles(m, k)
% Get the number of moles of a species in a mixture.
%
% :parameter k:
% Integer species number.
% :return:
% Moles of species number 'k'. If input 'k' is empty, return
% moles of every species in the mixture. Unit: kmol.
if nargin == 2
moles = calllib(ct, 'mix_speciesMoles', m.mixID, k);
elseif nargin == 1
nsp = m.nSpecies;
moles = zeros(1, nsp);
for i = 1:nsp
moles(i) = calllib(ct, 'mix_speciesMoles', ...
m.mixID, i);
end
else error('wrong number of arguments');
end
end
function mu = chemPotentials(m)
% Get the chemical potentials of species in the mixture.
%
@ -252,8 +309,14 @@ classdef Mixture < handle
% :parameter moles:
% Vector or string specifying the moles of species.
calllib(ct, 'mix_setMolesByName', m.mixID, moles);
% check back on this one!
if isa(moles, 'double')
l = length(moles);
calllib(ct, 'mix_setMoles', m.mixID, l, moles);
elseif isa(moles, 'string')
calllib(ct, 'mix_setMolesByName', m.mixID, moles);
else
error('The input must be a vector or string!');
end
end
function r = equilibrate(m, XY, err, maxsteps, maxiter, loglevel)

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

@ -60,8 +60,6 @@ classdef Func < handle
% :return:
% Instance of class :mat:func:`Func`
checklib;
if ~isa(typ, 'char')
error('Function type must be a string');
end
@ -69,24 +67,33 @@ classdef Func < handle
x.f1 = 0;
x.f2 = 0;
x.coeffs = 0;
itype = -1;
function nn = newFunc(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(ct, 'func_new', type, n, lenp, ptr);
nn = calllib(ct, 'func_new', itype, n, lenp, p);
elseif itype < 45
m = n;
nn = calllib(ct, 'func_new', type, n, m, 0);
nn = calllib(ct, 'func_new', itype, n, m, 0);
else
ptr = libpointer('doublePtr', p);
nn = calllib(ct, 'func_new', type, n, 0, ptr);
nn = calllib(ct, 'func_new', itype, n, 0, p);
end
end
if strcmp(typ, 'polynomial')
itype = 2;
elseif strcmp(typ, 'fourier')
itype = 1;
elseif strcmp(typ, 'arrhenius')
itype = 3;
elseif strcmp(typ, 'gaussian')
itype = 4;
end
if itype > 0
x.coeffs = p;
x.id = newFunc(itype, n, p);
@ -136,25 +143,6 @@ classdef Func < handle
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.
%

4
interfaces/matlab_experimental/LoadCantera.m
View File

@ -1,6 +1,6 @@
function LoadCantera
addpath('Class', 'Utility', 'PresetMixtures', 'PresetReactors', ...
'PresetObjects', '1D');
addpath('Class', 'Utility', 'PresetFlowDevices', 'PresetFunctors', ...
'PresetMixtures', 'PresetReactors', 'PresetObjects', '1D');
if ispc
ctname = 'cantera_shared.dll';
elseif ismac

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

@ -107,6 +107,17 @@ classdef FlowDevice < handle
end
end
function setMaster(f, d)
% Set the Master flow device used to compute this device's mass
% flow rate.
if strcmp(f.type, 'PressureController')
k = calllib(ct, 'flowdev_setMaster', f.id, d);
else
error('Master flow device can only be set for pressure controllers.');
end
end
function setValveCoeff(f, k)
% Set the valve coefficient 'K'.
%

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

@ -74,6 +74,17 @@ classdef Reactor < handle
calllib(ct, 'reactor_del', r.id);
end
function addSensitivityReaction(r, m)
% Specifies that the sensitivity of the state variables with
% respect to reaction m should be computed. The reactor must be
% part of a network first.
%
% :parameter m:
% Index number of reaction.
calllib(ct, 'reactor_addSensitivityReaction', r.id, m);
end
function insert(r, gas)
% Insert a solution or mixture into a reactor.
%

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

@ -4,15 +4,13 @@ classdef ReactorNet < handle
id
time
dt
atol
rtol
end
methods
%% ReactorNet class constructor
function r = ReactorNet(reactors)
% A 'ReactorNet' object is a container that holds one ormore
% A 'ReactorNet' object is a container that holds one or more
% 'Reactor' objects in a network. 'ReactorNet' objects are used
% to simultaneously advance the state of one or more coupled
% reactors.
@ -108,12 +106,23 @@ classdef ReactorNet < handle
calllib(ct, 'reactornet_setMaxTimeStep', r.id, maxstep);
end
function setSensitivityTolerances(r, rerr, aerr)
% Set the error tolerance for sensitivity analysis.
%
% :parameter rerr:
% Scalar relative error tolerance.
% :parameter aerr:
% Scalar absolute error tolerance.
calllib(ct, 'reactornet_setSensitivityTolerances', r.id, rerr, aerr);
end
function setTolerances(r, rerr, aerr)
% Set the error tolerance.
%
% :parameter rtol:
% :parameter rerr:
% Scalar relative error tolerance.
% :parameter atol:
% :parameter aerr:
% Scalar absolute error tolerance.
calllib(ct, 'reactornet_setTolerances', r.id, rerr, aerr);
@ -139,17 +148,35 @@ classdef ReactorNet < handle
t = calllib(ct, 'reactornet_time', r.id);
end
function t = get.rtol(r)
% Get the relative error tolerance
t = calllib(ct, 'reactornet_rtol', r.id);
end
function t = get.atol(r)
% Get the absolute error tolerance
function t = atol(r)
% Get the absolute error tolerance.
t = calllib(ct, 'reactornet_atol', r.id);
end
function t = rtol(r)
% Get the relative error tolerance.
t = calllib(ct, 'reactornet_rtol', r.id);
end
function s = sensitivity(r, component, p, rxtr)
% Get the sensitivity of the solution variable c in reactor
% rxtr with respect to the parameter p.
%
% :param component:
% String name of variable.
% :param p:
% Integer sensitivity parameter.
% :param rxtr:
% Instance of class 'reactor'.
if isa(component, 'string')
calllib(ct, 'reactornet_sensitivity', r.id, component, ...
p, rxtr.id);
end
% Check back on this one to add cases for component type integer.
end
end
end

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

@ -79,6 +79,17 @@ classdef ReactorSurface < handle
calllib(ct, 'reactorsurface_install', s.surfID, r.id);
end
function addSensitivityReaction(s, r)
% Specifies that the sensitivity of the state variables with
% respect to reaction m should be computed. The surface must be
% installed on a reactor and part of a network first.
%
% :parameter m:
% Index number of reaction.
calllib(ct, 'reactorsurface_addSensitivityReaction', s.surfID, r);
end
%% ReactorSurface get methods
function a = get.area(s)

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

@ -181,6 +181,15 @@ classdef Wall < handle
calllib(ct, 'wall_setHeatTransferCoeff', w.id, u);
end
function setEmissivity(w, epsilon)
% Set the emissivity.
%
% :param epsilon:
% Nondimensional emissivity.
calllib(ct, 'wall_setEmissivity', w.id, epsilon);
end
function setExpansionRateCoeff(w, k)
% Set the expansion rate coefficient.
%

2
samples/matlab_experimental/diff_flame.m
View File

@ -53,7 +53,7 @@ ox.TPX = {tin, p, oxcomp};
% prior, same for the tolerances.
f = AxisymmetricFlow(fuel, 'flow');
f.setPressure(p);
f.P = p;
f.setupGrid(initial_grid);
f.setSteadyTolerances('default', tol_ss{:});
f.setTransientTolerances('default', tol_ts{:});

2
samples/matlab_experimental/flame1.m
View File

@ -50,7 +50,7 @@ gas.TPX = {tburner, p, comp};
%% Create the flow object
f = AxisymmetricFlow(gas, 'flow');
f.setPressure(p);
f.P = p;
f.setupGrid(initial_grid);
f.setSteadyTolerances('default', tol_ss{:});
f.setTransientTolerances('default', tol_ts{:});

2
samples/matlab_experimental/flame2.m
View File

@ -49,7 +49,7 @@ gas.TPX = {tin, p, comp2};
%% Create the flow object
%
f = AxisymmetricFlow(gas,'flow');
f.setPressure(p);
f.P = p;
f.setupGrid(initial_grid);
f.setSteadyTolerances('default', tol_ss{:});
f.setTransientTolerances('default', tol_ts{:});