From 0189abb705c8ec0d86277dbdf1e82ba2bae3c574 Mon Sep 17 00:00:00 2001 From: Su Sun <56568705+ssun30@users.noreply.github.com> Date: Tue, 11 Jan 2022 14:36:53 -0500 Subject: [PATCH] Experimental Matlab toolbox now working --- .../matlab_experimental/1D/AxiStagnFlow.m | 4 + .../matlab_experimental/1D/AxisymmetricFlow.m | 11 + interfaces/matlab_experimental/1D/Domain1D.m | 531 ++++++ interfaces/matlab_experimental/1D/FreeFlame.m | 9 + interfaces/matlab_experimental/1D/Inlet.m | 9 + interfaces/matlab_experimental/1D/Outlet.m | 9 + interfaces/matlab_experimental/1D/OutletRes.m | 9 + interfaces/matlab_experimental/1D/Sim1D.m | 458 ++++++ interfaces/matlab_experimental/1D/Surface.m | 17 + interfaces/matlab_experimental/1D/SymmPlane.m | 9 + .../matlab_experimental/Base/Interface.m | 64 +- .../matlab_experimental/Base/Kinetics.m | 423 ++--- interfaces/matlab_experimental/Base/Mixture.m | 84 +- .../matlab_experimental/Base/Solution.m | 27 +- .../matlab_experimental/Base/ThermoPhase.m | 1441 +++++++++++------ .../matlab_experimental/Base/Transport.m | 56 +- interfaces/matlab_experimental/Base/Water.m | 3 +- .../Constants/faradayconstant.m | 4 + .../Constants/gasconstant.m | 4 + .../matlab_experimental/Constants/oneatm.m | 4 + interfaces/matlab_experimental/Func/Func.m | 30 +- .../matlab_experimental/Func/gaussian.m | 12 + interfaces/matlab_experimental/Func/polynom.m | 17 + interfaces/matlab_experimental/LoadCantera.m | 20 +- .../Reactor/ConstPressureReactor.m | 14 + .../matlab_experimental/Reactor/FlowDevice.m | 36 +- .../matlab_experimental/Reactor/FlowReactor.m | 12 + .../Reactor/IdealGasConstPressureReactor.m | 14 + .../Reactor/IdealGasReactor.m | 13 + .../Reactor/MassFlowController.m | 21 + .../matlab_experimental/Reactor/Reactor.m | 86 +- .../matlab_experimental/Reactor/ReactorNet.m | 42 +- .../Reactor/ReactorSurface.m | 30 +- .../matlab_experimental/Reactor/Reservoir.m | 18 + .../matlab_experimental/Reactor/Valve.m | 32 + interfaces/matlab_experimental/Reactor/Wall.m | 56 +- .../matlab_experimental/UnloadCantera.m | 4 + .../Utility/UnloadCantera.m | 4 - .../Utility/canteraGitCommit.m | 8 + .../Utility/canteraVersion.m | 8 + .../matlab_experimental/Utility/checklib.m | 2 +- .../matlab_experimental/Utility/cleanup.m | 13 +- .../Utility/getDataDirectories.m | 8 + .../matlab_experimental/Utility/geterr.m | 12 + .../matlab_experimental/Utility/importEdge.m | 15 + .../Utility/importInterface.m | 11 + interfaces/matlab_experimental/cleanup.m | 12 - interfaces/matlab_experimental/ct.m | 8 + interfaces/matlab_experimental/readme.txt | 9 + 49 files changed, 2770 insertions(+), 973 deletions(-) create mode 100644 interfaces/matlab_experimental/1D/AxiStagnFlow.m create mode 100644 interfaces/matlab_experimental/1D/AxisymmetricFlow.m create mode 100644 interfaces/matlab_experimental/1D/Domain1D.m create mode 100644 interfaces/matlab_experimental/1D/FreeFlame.m create mode 100644 interfaces/matlab_experimental/1D/Inlet.m create mode 100644 interfaces/matlab_experimental/1D/Outlet.m create mode 100644 interfaces/matlab_experimental/1D/OutletRes.m create mode 100644 interfaces/matlab_experimental/1D/Sim1D.m create mode 100644 interfaces/matlab_experimental/1D/Surface.m create mode 100644 interfaces/matlab_experimental/1D/SymmPlane.m create mode 100644 interfaces/matlab_experimental/Constants/faradayconstant.m create mode 100644 interfaces/matlab_experimental/Constants/gasconstant.m create mode 100644 interfaces/matlab_experimental/Constants/oneatm.m create mode 100644 interfaces/matlab_experimental/Func/gaussian.m create mode 100644 interfaces/matlab_experimental/Func/polynom.m create mode 100644 interfaces/matlab_experimental/Reactor/ConstPressureReactor.m create mode 100644 interfaces/matlab_experimental/Reactor/FlowReactor.m create mode 100644 interfaces/matlab_experimental/Reactor/IdealGasConstPressureReactor.m create mode 100644 interfaces/matlab_experimental/Reactor/IdealGasReactor.m create mode 100644 interfaces/matlab_experimental/Reactor/MassFlowController.m create mode 100644 interfaces/matlab_experimental/Reactor/Reservoir.m create mode 100644 interfaces/matlab_experimental/Reactor/Valve.m create mode 100644 interfaces/matlab_experimental/UnloadCantera.m delete mode 100644 interfaces/matlab_experimental/Utility/UnloadCantera.m create mode 100644 interfaces/matlab_experimental/Utility/canteraGitCommit.m create mode 100644 interfaces/matlab_experimental/Utility/canteraVersion.m create mode 100644 interfaces/matlab_experimental/Utility/getDataDirectories.m create mode 100644 interfaces/matlab_experimental/Utility/geterr.m create mode 100644 interfaces/matlab_experimental/Utility/importEdge.m create mode 100644 interfaces/matlab_experimental/Utility/importInterface.m delete mode 100644 interfaces/matlab_experimental/cleanup.m create mode 100644 interfaces/matlab_experimental/ct.m create mode 100644 interfaces/matlab_experimental/readme.txt diff --git a/interfaces/matlab_experimental/1D/AxiStagnFlow.m b/interfaces/matlab_experimental/1D/AxiStagnFlow.m new file mode 100644 index 000000000..a9b288560 --- /dev/null +++ b/interfaces/matlab_experimental/1D/AxiStagnFlow.m @@ -0,0 +1,4 @@ +function m = AxiStagnFLow(gas) + % Get an axisymmetric stagnation flow domain. + m = Domain1D(1, gas); +end diff --git a/interfaces/matlab_experimental/1D/AxisymmetricFlow.m b/interfaces/matlab_experimental/1D/AxisymmetricFlow.m new file mode 100644 index 000000000..a1a25c7e1 --- /dev/null +++ b/interfaces/matlab_experimental/1D/AxisymmetricFlow.m @@ -0,0 +1,11 @@ +function m = AxisymmetricFlow(gas, id) + % Create an axisymmetric flow domain. + % :param id: + % String ID of the flow. + m = Domain1D(1, gas); + if nargin == 1 + m.setID('flow'); + else + m.setID(id); + end +end diff --git a/interfaces/matlab_experimental/1D/Domain1D.m b/interfaces/matlab_experimental/1D/Domain1D.m new file mode 100644 index 000000000..06b347851 --- /dev/null +++ b/interfaces/matlab_experimental/1D/Domain1D.m @@ -0,0 +1,531 @@ +classdef Domain1D < handle + + properties + dom_id + domain_type + T + end + + methods + %% Domain1D class constructor. + + function d = Domain1D(a, b, c) + % parameter a: + % Integer type of domain. Possible values are: + % * 1 - Stagnation Flow + % * 2 - Inlet1D + % * 3 - Surf1D + % * 4 - Symm1D + % * 5 - Outlet1D + % * 6 - Reacting Surface + % * 8 - Sim1D + % * -2 - OutletRes + % + % parameter b: + % Instance of class 'Solution' (for a == 1) or 'Interface' + % (for a == 6). Not used for all other valid values of 'a'. + % + % parameter c: + % Integer, either 1 or 2, indicating whether an axisymmetric + % stagnation flow or a free flame should be created. If not + % specified, defaults to 1. Ignored if a!= 1. + + checklib; + d.dom_id = 1; + + if nargin == 1 + if a == 2 + calllib(ct, 'inlet_new'); + elseif a == 3 + calllib(ct, 'surf_new'); + elseif a == 4 + calllib(ct, 'symm_new'); + elseif a == 5 + calllib(ct, 'outlet_new'); + elseif a == -2 + calllib(ct, 'outletres_new'); + else + error('Not enough arguments for that job number'); + end + elseif nargin == 2 + % a stagnation flow + if a == 1 + if isa(b, 'Solution') + d.dom_id = calllib(ct, 'stflow_new', ... + b.tp_id, b.kin_id, b.tr_id, 1); + else + error('Wrong argument type. Expecting instance of class Solution.'); + end + elseif a == 6 + if isa(b, 'Interface') + d.dom_id = calllib(ct, 'reactingsurf_new'); + calllib(ct, 'reactingsurf_setkineticsmgr', ... + d.dom_id, b.kin_id); + else + error('Wrong argument type. Expecting instance of class Interface.'); + end + else + error('Wrong object type.'); + end + elseif nargin == 3 + if a == 1 + if isa(b, 'Solution') + d.dom_id = calllib(ct, 'stflow_new', ... + b.tp_id, b.kin_id, b.tr_id, c); + else + error('Wrong argument type. Expecting instance of class Solution.'); + end + else + error('Unknown domain type.'); + end + end +% if d.dom_id < 0 +% error(geterr); +% end + d.domain_type = a; + end + + %% Utility Methods + + function dom_clear(d) + % Delete the Domain1D object + checklib; + calllib(ct, 'domain_del', d.dom_id); + end + + %% Domain Methods + + function n = componentIndex(d, name) + % Get the index of a component given its name + % + % parameter d: + % Instance of class 'Domain1D'. + % parameter name: + % String name of the component to look up. If a numeric + % value is passed, it will be returned directly. + % return: + % Index of the component. + + checklib; + if isa(name, 'double') + n = name; + else + n = calllib(ct, 'domain_componentIndex', ... + d.dom_id, name); + if n >= 0 + n = n+1; + end + end + if n <= 0 + error('Component not found'); + end + end + + function n = componentName(d, index) + % Get the name of a component given its index. + % + % parameter d: + % Instance of class 'Domain1D'. + % parameter n: + % Integer or vector of integers of components' names. + % return: + % Cell array of component names. + + checklib; + n = length(index); + s = cell(m); + for i = 1:n + id = index(i)-1; + buflen = calllib(ct, 'domain_componentName', ... + d.dom_id, id, 0, 0); + if buflen > 0 + aa = char(zeros(1, buflen)); + [out_buf, aa] = calllib(ct, 'domain_componentName', ... + d.dom_id, id, buflen, aa); + s{i} = aa; + end + end + end + + function d = disableEnergy(d) + % Disable the energy equation. + % + % parameter d: + % Instance of class 'Domain1D'. + + checklib; + disp(' '); + disp('Disabling the energy equation...'); + calllib(ct, 'stflow_solveEnergyEqn', d.dom_id, 0); + end + + function i = domainIndex(d) + % Get the domain index. + % return: + % This function returns an integer flag denoting the + % location of the domain, beginning with 1 at the left. + checklib; + i = calllib(ct, 'domain_index', d.dom_id); + if i >= 0 + i = i + 1; + end + if i <= 0 + error('Domain not found'); + end + end + + function i = domainType(d) + % Get the type of domain. + % return: + % This function returns an integer flag denoting the domain + % type. + checklib; + i = calllib(ct, 'domain_type', d.dom_id); + end + + function d = enableEnergy(d) + % Enable the energy equation. + % parameter d: + % Instance of class 'Domain1D'. + + checklib; + disp(' '); + disp('Enabling the energy equation...'); + calllib(ctm 'stflow_solveEnergyEqn', d.dom_id, 1); + end + + function zz = gridPoints(d, n) + % Get grid points from a domain. + % parameter d: + % Instance of class 'Domain1D'. + % parameter n: + % Optional, vector of grid points to be retrieved. + % return: + % Vector of grid points. + + checklib; + if nargin == 1 + np = d.nPoints; + zz = zeros(1, d.np); + for i = 1:np + zz(i) = calllib(ct, 'domain_grid', dom_id, i-1); + end + else + m = length(n); + zz = zeros(1, m); + for i = 1:m + zz(i) = calllib(ct, 'domain_grid', dom_id, n(i)-1); + end + end + end + + function a = isFlow(d) + % Determine whether a domain is a flow. + % parameter d: + % Instance of class 'Domain1D'. + % return: + % 1 if the domain is a flow domain, and 0 otherwise. + + checklib; + t = d.domainType; + % See Domain1D.h for definitions of constants. + if t < 100 + a = 1; + else a = 0; + end + end + + function a = isInlet(d) + % Determine whether a domain is an inlet. + % parameter d: + % Instance of class 'Domain1D'. + % return: + % 1 if the domain is an inlet, and 0 otherwise. + + checklib; + t = d.domainType; + % See Domain1D.h for definitions of constants. + if t == 104 + a = 1; + else a = 0; + end + end + + function a = isSurface(d) + % Determine whether a domain is a surface. + % parameter d: + % Instance of class 'Domain1D'. + % return: + % 1 if the domain is a surface, and 0 otherwise. + + checklib; + t = d.domainType; + % See Domain1D.h for definitions of constants. + if t == 102 + a = 1; + else a = 0; + end + end + + function mdot = massFlux(d) + % Get the mass flux. + % return: + % The mass flux in the domain. + checklib; + mdot = calllib(ct, 'bdry_mdot', d.dom_id); + end + + function y = massFraction(d, k) + % Get the mass fraction of a species given its integer index. + % This method returns the mass fraction of species 'k', where k + % is the integer index of the species in the flow domain to + % which the boundary domain is attached. + % + % parameter d: + % Instance of class 'Domain1D' + % parameter k: + % Integer species index + % return: + % Mass fraction of species + + checklib; + if d.domainIndex == 0 + error('No flow domain attached!') + end + + if d.isInlet + y = calllib(ct, 'bdry_massFraction', d.dom_id, k-1); + else error('Input domain must be an inlet'); + end + end + + function n = nComponents(d) + % Get the number of components. + % return: + % Number of variables at each grid point. + checklib; + n = calllib(ct, 'domain_nComponents', d.dom_id); + end + + function n = nPoints(d) + % Get the number of grid points. + % return: + % Integer number of grid points. + checklib; + n = calllib(ct, 'domain_nPoints', d.dom_id); + end + + function setBounds(d, component, lower, upper) + % Set bounds on the solution components. + % parameter d: + % Instance of class 'Domain1D'. + % parameter component: + % String, component to set the bounds on. + % parameter lower: + % Lower bound. + % parameter upper: + % Upper bound. + checklib; + n = d.componentIndex(component); + calllib(ct, 'domain_setBounds', d.dom_id, n-1, lower, upper); + end + + function setCoverageEqs(d, onoff) + % Enable or disable solving the coverage equations. + % parameter d: + % Instance of class 'Domain1D' + % parameter onoff: + % string, one of 'on' or 'yes' to turn solving the coverage + % equations on. One of 'off' or 'no' to turn off the + % coverage equation. + checklib; + if d.domain_type ~= 6 + error('Wrong domain type. Expected a reacting surface domain.'); + end + + ion = -1; + if isa(onoff,'char') + if strcmp(onoff, 'on') || strcmp(onoff, 'yes') + ion = 1; + elseif strcmp(onoff, 'off') || strcmp(onoff, 'no') + ion = 0; + else + error(strcat('unknown option: ', onoff)) + end + elseif isa(onoff, 'numeric') + ion = onoff; + end + calllib(ct, 'reactingsurf_enableCoverageEqs', d.dom_id, ion); + end + + function setFixedTempProfile(d, profile) + % Set a fixed temperature profile to use when the energy + % equation is not being solved. The profile must be entered as + % an array of positions/temperatures, which may be in rows or + % columns. + % + % parameter d: + % Instance of class 'Domain1D'. + % parameter profile: + % n x 2 or 2 x n array of 'n' points at which the + % temperature is specified. + checklib; + sz = size(profile); + if sz(1) == 2 + l = length(1, :); + calllib(ct, '', d.dom_id, l, profile(1, :), l, profile(2, :)); + elseif sz(2) == 2 + l = length(:, 1); + calllib(ct, '', d.dom_id, l, profile(:, 1); l, profile(:, 2)); + else error('Wrong temperature profile array shape.'); + end + end + + function setID(d, id) + % Set the ID tag for a domain. + % parameter id: + % String ID to assign. + checklib; + calllib(ct, 'domain_setID', d.dom_id, id); + end + + function setMdot(d, mdot) + % Set the mass flow rate. + % parameter mdot: + % Mass flow rate. + checklib; + calllib(ct, 'bdry_setTemperature', d.dom_id, mdot); + end + + function setMoleFractions(d, x) + % Set the mole fractions. + % parameter x: + % String specifying the species and mole fractions in the + % format "'Spec:X,Spec2:X2'". + checklib; + calllib(ct, 'bdry_setMoleFractions', d.dom_id, x); + end + + function setPressure(d, p) + % Set the pressure. + % parameter p: + % Pressure to be set. Unit: Pa. + checklib; + calllib(ct, 'bdry_setPressure', d.dom_id, p); + end + + function setProfileD(d, n, p) + % Set the profile of a component. + % Convenience function to allow an instance of 'Domain1D' to + % have a profile of its components set when it is part of a + % 'Stack'. + % + % parameter d: + % Instance of class 'Domain1d'. + % parameter n: + % Integer index of component, vector of component indices, + % string of component name, or cell array of strings of + % component names. + % parameter p: + % n x 2 array, whose columns are the relative (normalized) + % positions and the component values at those points. The + % number of positions 'n' is arbitrary. + + checklib; + if d.stack == 0 + error('Install domain in stack before calling setProfile.'); + end + d.stack.setProfile(d.domainIndex, n, p); + end + + function setSteadyTolerances(d, component, rtol, atol) + % Set the steady-state tolerances. + % parameter d: + % Instance of class 'Domain1D'. + % parameter component: + % String or cell array of strings of component values whose + % tolerances should be set. If 'default' is specified, the + % tolerance of all components will be set. + % parameter rtol: + % Relative tolerance. + % parameter atol: + % Absolute tolerance. + + checklib; + if strcmp(component, 'default') + nc = d.nComponents; + for ii = 1:nc + calllib(ct, 'domain_setSteadyTolerances', ... + d.dom_id, ii, rtol, atol); + end + elseif iscell(component) + nc = length(component); + for ii = 1:nc + n = d.componentIndex(component{ii}); + calllib(ct, 'domain_setSteadyTolerances', ... + d.dom_id, n, rtol, atol); + end + else + n = d.componentIndex(component); + calllib(ct, 'domain_setSteadyTolerances', ... + d.dom_id, ii, rtol, atol); + end + end + + function setTransientTolerances(d, component, rtol, atol) + % Set the transient tolerances. + % parameter d: + % Instance of class 'Domain1D'. + % parameter component: + % String or cell array of strings of component values whose + % tolerances should be set. If 'default' is specified, the + % tolerance of all components will be set. + % parameter rtol: + % Relative tolerance. + % parameter atol: + % Absolute tolerance. + + checklib; + if strcmp(component, 'default') + nc = d.nComponents; + for ii = 1:nc + calllib(ct, 'domain_setTransientTolerances', ... + d.dom_id, ii, rtol, atol); + end + elseif iscell(component) + nc = length(component); + for ii = 1:nc + n = d.componentIndex(component{ii}); + calllib(ct, 'domain_setTransientTolerances', ... + d.dom_id, n, rtol, atol); + end + else + n = d.componentIndex(component); + calllib(ct, 'domain_setTransientTolerances', ... + d.dom_id, ii, rtol, atol); + end + end + + function temperature = get.T(d) + % Get the boundary temperature (K). + checklib; + temperature = calllib(ct, 'bdry_temperature', d.dom_id); + end + + function set.T(d, t) + % Set the temperature (K). + checklib; + if temperature <= 0 + error('The temperature must be positive'); + end + calllib(ct, 'bdry_setTemperature', d.dom_id, t); + end + + function setupGrid(d, grid) + % Set up the solution grid. + checklib; + calllib(ct, 'domain_setupGrid', d.dom_id, numel(grid), grid); + end + + end +end diff --git a/interfaces/matlab_experimental/1D/FreeFlame.m b/interfaces/matlab_experimental/1D/FreeFlame.m new file mode 100644 index 000000000..a069866b3 --- /dev/null +++ b/interfaces/matlab_experimental/1D/FreeFlame.m @@ -0,0 +1,9 @@ +function m = FreeFlame(gas, id) + % Create a freely-propagating flat flame. + m = Domain1D(1, gas, 2); + if nargin == 1 + m.setID('flame'); + else + m.setID(id); + end +end diff --git a/interfaces/matlab_experimental/1D/Inlet.m b/interfaces/matlab_experimental/1D/Inlet.m new file mode 100644 index 000000000..8e59b2f33 --- /dev/null +++ b/interfaces/matlab_experimental/1D/Inlet.m @@ -0,0 +1,9 @@ +function m = Inlet(id) + % Create an inlet domain. + m = Domain1D(2); + if nargin == 0 + m.setID('inlet'); + else + m.setID(id); + end +end diff --git a/interfaces/matlab_experimental/1D/Outlet.m b/interfaces/matlab_experimental/1D/Outlet.m new file mode 100644 index 000000000..30afdc783 --- /dev/null +++ b/interfaces/matlab_experimental/1D/Outlet.m @@ -0,0 +1,9 @@ +function m = Outlet(id) + % Create an outlet domain. + m = Domain1D(5); + if nargin == 0 + m.setID('outlet'); + else + m.setID(id); + end +end diff --git a/interfaces/matlab_experimental/1D/OutletRes.m b/interfaces/matlab_experimental/1D/OutletRes.m new file mode 100644 index 000000000..813b0bf39 --- /dev/null +++ b/interfaces/matlab_experimental/1D/OutletRes.m @@ -0,0 +1,9 @@ +function m = OutletRes(id) + % Create an outlet reservoir domain. + m = Domain1D(-2); + if nargin == 0 + m.setID('outletres'); + else + m.setID(id); + end +end diff --git a/interfaces/matlab_experimental/1D/Sim1D.m b/interfaces/matlab_experimental/1D/Sim1D.m new file mode 100644 index 000000000..f29b359ae --- /dev/null +++ b/interfaces/matlab_experimental/1D/Sim1D.m @@ -0,0 +1,458 @@ +classdef Stack < handle + + properties + st_id + domains + end + + methods + %% Stack class constructor + + function s = Stack(domains) + % A stack object is a container for one-dimensional domains, + % which are instances of class Domain1D. The domains are of two + % types - extended domains, and connector domains. + % + % parameter domains: + % Vector of domain instances. + % return: + % Instance of class 'Stack'. + + checklib; + + s.st_id = 1; + s.domains = domains; + if nargin == 1 + nd = length(domains); + ids = zeros(1, nd); + for n=1:nd + ids(n) = domains(n).dom_id; + end + s.st_id = calllib(ct, 'sim1D_new', nd, ids); + else + help(Stack); + error('Wrong number of parameters.'); + end +% if s.st_id < 0 +% error(geterr); +% end + end + + %% Utility Methods + + function st_clear(s) + % Delete the Sim1D object + checklib; + calllib(ct, 'sim1D_del', s.st_id); + end + + function display(s, fname) + % Show all domains. + % + % parameter s: + % Instance of class 'Stack'. + % parameter fname: + % File to write summary to. If omitted, output is to the + % command window. + + checklib; + if nargin == 1 + fname = '-'; + end + calllib(ct, 'sim1D_showSolution', s.st_id, fname); + end + + %% Stack Methods + + function n = stackIndex(s, name) + % Get the index of a domain in a stack given its name. + % + % parameter s: + % Instance of class 'Stack'. + % parameter name: + % If double, the value is returned. Otherwise, the name is + % looked up and its index is returned. + % return: + % Index of domain. + + checklib; + if isa(name, 'double') + n = name; + else + n = calllib(ct, 'sim1D_domainIndex', s.st_id, name); + if n >= 0 + n = n+1; + else + error('Domain not found'); + end + end + end + + function z = grid(s, name) + % Get the grid in one domain. + % + % parameter s: + % Instance of class 'Stack'. + % parameter name: + % Name of the domain for which the grid should be retrieved. + % return: + % The grid in domain name. + + n = s.stackIndex(name); + d = s.domains(n); + z = d.gridPoints; + end + + function plotSolution(s, domain, component) + % Plot a specified solution component. + % + % parameter s: + % Instance of class 'Stack'. + % parameter domain: + % Name of domain from which the component should be + % retrieved. + % parameter component: + % Name of the component to be plotted. + + n = s.stackIndex(domain); + d = s.domains(n); + z = d.gridPoints; + x = s.solution(domain, component); + plot(z, x); + xlabel('z (m)'); + ylabel(component); + end + + function r = resid(s, domain, rdt, count) + % Get the residuals. + % + % parameter s: + % Instance of class 'Stack'. + % parameter domain: + % Name of the domain. + % parameter rdt: + % parameter count: + % return: + + checklib; + + if nargin == 2 + rdt = 0.0; + count = 0; + end + + idom = s.stackIndex(domain); + d = s.domains(idom); + + nc = d.nComponents; + np = d.nPoints; + + r = zeros(nc, np); + calllib(ct, 'sim1D_eval', s.st_id, rdt, count); + for m = 1:nc + for n = 1:np + r(m, n) = calllib(ct, 'sim1D_workValue', ... + s.st_id, idom - 1, m - 1, n - 1); + end + end + end + + function restore(s, fname, id) + % Restore a previously-saved solution. + % This method can be used ot provide an initial guess for the + % solution. + % + % parameter s: + % Instance of class 'Stack'. + % parameter fname: + % File name of an XML file containing solution info. + % parameter id: + % ID of the element that should be restored. + checklib; + calllib(ct, 'sim1D_restore', s.st_id, fname, id) + end + + function saveSoln(s, fname, id, desc) + % Save a solution to a file. + % The output file is in a format that can be used by 'restore'. + % + % parameter s: + % Instance of class 'Stack'. + % parameter fname: + % File name where XML file should be written. + % parameter id: + % ID to be assigned to the XMl element when it is written. + % parameter desc: + % Description to be written to the output file. + checklib; + + if nargin == 1 + fname = 'soln.xml'; + id = 'solution'; + desc = '--'; + elseif nargin == 2 + id = 'solution'; + desc = '--'; + elseif nargin == 3 + desc = '--'; + end + calllib(ct, 'sim1D_save', s.st_id, fname, id, desc); + end + + function setFlatProfile(s, domain, comp, v) + % Set a component to a value across the entire domain. + % + % parameter s: + % Instance of class 'Stack'. + % parameter domain: + % Integer ID of the domain. + % parameter comp: + % Component to be set. + % parameter v: + % Double value to be set. + checklib; + calllib(ct, 'sim1D_setFlatProfile', s.st_id, ... + domain - 1, comp - 1, v); + end + + function setMaxJacAge(s, ss_age, ts_age) + % Set the number of times the Jacobian will be used before it + % is recomputed. + % + % parameter s: + % Instance of class 'Stack'. + % parameter ss_age: + % Maximum age of the Jacobian for steady state analysis. + % parameter ts_age: + % Maximum age of the Jacobian for transient analysis. If not + % specified, deftauls to 'ss_age'. + checklib; + + if nargin == 2 + ts_age = ss_age; + end + calllib(ct, 'sim1D_setMaxJacAge', s.st_id, ss_age, ts_age); + end + + function setProfile(s, name, comp, p) + % Specify a profile for one component, + % + % The solution vector values for this component will be + % linearly interpolated from the discrete function defined by + % p(:, 1) vs p(:, 2). + % Note that "p(1, 1) = 0.0" corresponds to the leftmost grid + % point in the specified domain, and "p(1, n) = 1.0" + % corresponds to the rightmost grid point. This method can be + % called at any time, but is usually used to set the initial + % guess for the solution. + % + % Example (assuming 's' is an instance of class 'Stack'): + % >> zr = [0.0, 0.1, 0.2, 0.4, 0.8, 1.0]; + % >> v = [500, 650, 700, 730, 800, 900]; + % >> s.setProfile(1, 2, [zr, v]); + % + % parameter s: + % Instance of class 'Stack'. + % parameter name: + % Domain name. + % parameter comp: + % Component number. + % parameter p: + % n x 2 array, whose columns are the relative (normalized) + % positions and the component values at those points. The + % number of positions 'n' is arbitrary. + + checklib; + + if isa(name, 'double') + n = name; + else + n = s.domainIndex(name); + end + + d = s.domains(n); + + if isa(comp, 'double') || isa(comp, 'cell') + c = comp; + elseif isa(comp, 'char') + c = {comp}; + else + error('Wrong type.'); + end + + np = length(c); + sz = size(p); + if sz(1) == np + 1; + for j = 1:np + ic = d.componentIndex(c{j}); + calllib(ct, 'sim1D_setProfile', s.st_id, ... + n - 1, ic - 1, sz(1), p(1, :), sz(1), p(j+1, :)); + end + elseif sz(2) == np + 1; + ic = d.componentIndex(c{j}); + calllib(ct, 'sim1D_setProfile', s.st_id, ... + n - 1, ic - 1, sz(2), p(:, 1), sz(2), p(:, j+1)); + else + error('Wrong profile shape.'); + end + end + + function setRefineCriteria(s, n, ratio, slope, curve, prune) + % Set the criteria used to refine the grid. + % + % parameter s: + % Instance of class 'Stack'. + % parameter ratio: + % Maximum size ratio between adjacent cells. + % parameter slope: + % Maximum relative difference in value between adjacent + % points. + % parameter curve: + % Maximum relative difference in slope between adjacent + % cells. + % parameter prune: + % Minimum value for slope or curve for which points will be + % retained or curve value is below prune for all components, + % it will be deleted, unless either neighboring point is + % already marked for deletion. + + checklib; + + if nargin < 3 + ratio = 10.0; + end + if nargin < 4 + slope = 0.8; + end + if nargin < 5 + curve = 0.8; + end + if nargin < 6 + prune = -0.1; + end + calllib(ct, 'sim1D_setRefineCriteria', s.st_id, ... + n - 1, ratio, slope, curve, prune); + end + + function setTimeStep(s, stepsize, steps) + % Specify a sequence of time steps. + % + % parameter stepsize: + % Initial step size. + % parameter steps: + % Vector of number of steps to take before re-attempting + % solution of steady-state problem. + % For example, steps = [1, 2, 5, 10] would cause one time + % step to be taken first time the steady-state solution + % attempted. If this failed, two time steps would be taken. + checklib; + calllib(ct, 'sim1D_', s.st_id, ... + stepsize, length(steps), steps); + end + + function setValue(s, n, comp, localPoints, v) + % Set the value of a single entry in the solution vector. + % + % Example (assuming 's' is an instance of class 'Stack'): + % + % setValue(s, 3, 5, 1, 5.6); + % + % This sets component 5 at the leftmost point (local point 1) + % in domain 3 to the value 5.6. Note that the local index + % always begins at 1 at the left of each domain, independent of + % the global index of the point, wchih depends on the location + % of this domain in the stack. + % + % parameter s: + % Instance of class 'Stack'. + % parameter n: + % Domain number. + % parameter comp: + % Component number. + % parameter localPoints: + % Local index of the grid point in the domain. + % parameter v: + % Value to be set. + checklib; + calllib(ct, 'sim1D_setValue', s.st_id, ... + n - 1, comp - 1, localPoints - 1, v); + end + + function solution(s, domain, component) + % Get a solution component in one domain. + % + % parameter s: + % Instance of class 'Stack'. + % parameter domain: + % String name of the domain from which the solution is + % desired. + % parameter component: + % String component for which the solution is desired. If + % omitted, solution for all of the components will be + % returned in an 'nPoints' x 'nComponnts' array. + % return: + % Either an 'nPoints' x 1 vector, or 'nPoints' x + % 'nCOmponents' array. + checklib; + + idom = s.stackIndex(domain); + d = s.domains(idom); + np = d.nPoints; + if nargin == 3 + icomp = d.componentIndex(component); + x = zeros(1, np); + for n = 1:np + x(n) = calllib(ct, 'sim1D_Value', s.st_id, ... + idom - 1, icomp - 1, n - 1); + end + else + nc = d.nComponents; + x = zeros(nc, np); + for m = 1:nc + for n = 1:np + x(m, n) = calllib(ct, 'sim1D_Value', s.st_id, ... + idom - 1, m - 1, n - 1); + end + end + end + end + + function solve(s, loglevel, refine_grid) + % Solve the problem. + % + % parameter s: + % Instance of class 'Stack'. + % parameter loglevel: + % Integer flag controlling the amount of diagnostic output. + % Zero supresses all output, and 5 produces very verbose + % output. + % parameter refine_grid: + % Integer, 1 to allow grid refinement, 0 to disallow. + checklib; + calllib(ct, 'sim1D_solve', s.st_id, loglevel, refine_grid); + end + + function b = subsref(s, index) + % Redefine subscripted references. + switch index.type + case '()' + b = s.domains(index.subs{:}); + case '.' + n = s.domainIndex(index.subs); + b = s.domains(n); + otherwise + error('syntax error'); + end + end + + function writeStats(s) + % Print statistics for the current solution. + % Prints a summary of the number of function and Jacobian + % evaluations for each grid, and the CPU time spent on each + % one. + checklib; + calllib(ct, 'sim1D_writeStats', s.st_id); + end + + end +end diff --git a/interfaces/matlab_experimental/1D/Surface.m b/interfaces/matlab_experimental/1D/Surface.m new file mode 100644 index 000000000..03c7544c0 --- /dev/null +++ b/interfaces/matlab_experimental/1D/Surface.m @@ -0,0 +1,17 @@ +function m = Surface(id, surface_mech) + % Create an surface domain. + % :param surface_mech + % Instance of class 'Interface' defining the surface reaction + % mechanism to be used. Optional. + if nargin < 2 + m = Domain1D(3); + if nargin == 0 + m.setID('surface'); + elseif nargin == 1 + m.setID(id); + end + else + m = Domain1D(6, surface_mech); + m.setID(id); + end +end diff --git a/interfaces/matlab_experimental/1D/SymmPlane.m b/interfaces/matlab_experimental/1D/SymmPlane.m new file mode 100644 index 000000000..32388544b --- /dev/null +++ b/interfaces/matlab_experimental/1D/SymmPlane.m @@ -0,0 +1,9 @@ +function m = SymmPlane(id) + % Create an symmetry plane domain. + m = Domain1D(4); + if nargin == 0 + m.setID('symmetry_plane'); + else + m.setID(id); + end +end diff --git a/interfaces/matlab_experimental/Base/Interface.m b/interfaces/matlab_experimental/Base/Interface.m index 0c85bcc71..dd2fb5484 100644 --- a/interfaces/matlab_experimental/Base/Interface.m +++ b/interfaces/matlab_experimental/Base/Interface.m @@ -1,61 +1,55 @@ -classdef Interface < handle +classdef Interface < handle & ThermoPhase & Kinetics properties - th - kin - end - - properties(Constant = true) - lib = 'cantera_shared' + coverages end methods %% Interface class constructor function s = Interface(src, id, p1, p2, p3, p4) - % :param src: + % parameter src: % CTI or CTML file containing the interface or edge phase. - % :param id: + % parameter id: % Name of the interface or edge phase in the source file. - % :param p1/P2/P3/P4: + % parameter p1/p2/p3/p4: % Adjoining phase to the interface; - % :return: + % return: % Instance of class 'Interface'. checklib; t = ThermoPhase(src, id); + s@ThermoPhase(src, id); if nargin == 2 - k = Kinetics(t, src, id); + args = {}; elseif nargin == 3 - k = Kinetics(t, src, id, p1); + args = {p1}; elseif nargin == 4 - k = Kinetics(t, src, id, p1, p2); + args = {p1, p2}; elseif nargin == 5 - k = Kinetics(t, src, id, p1, p2, p3); + args = {p1, p2, p3}; elseif nargin == 6 - k = Kinetics(t, src, id, p1, p2, p3, p4); + args = {p1, p2, p3, p4}; end - - s.kin = k; - s.th = t; + s@Kinetics(t, src, id, args{:}); end %% Interface methods - function c = coverages(s) + function c = get.coverages(s) % Get the surface coverages of the species on an interface. % - % :return: + % 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; + surf_id = s.tr_id; + nsp = s.nSpecies; xx = zeros(1, nsp); pt = libpointer('doublePtr', xx); - calllib(s.lib, 'surf_getCoverages', surf_id, xx); + calllib(ct, 'surf_getCoverages', surf_id, xx); c = pt.Value; if nargout == 0 @@ -63,7 +57,7 @@ classdef Interface < handle set(gcf, 'Name', 'Coverages') bar(c); colormap(summer); - nm = speciesNames(s); + nm = s.speciesNames; set(gca, 'XTickLabel', nm); xlabel('Species Name'); ylabel('Coverage'); @@ -74,17 +68,17 @@ classdef Interface < handle function c = concentrations(s) % Get the concentrations of the species on an interface. % - % :return: + % 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; + surf_id = s.tr_id; + nsp = s.nSpecies; xx = zeros(1, nsp); pt = libpointer('doublePtr', xx); - calllib(s.lib, 'surf_getConcentrations', surf_id, xx); + calllib(ct, 'surf_getConcentrations', surf_id, xx); c = pt.Value; if nargout == 0 @@ -100,10 +94,10 @@ classdef Interface < handle end end - function setCoverages(s, cov, norm) + function set.coverages(s, cov, norm) % Set surface coverages of the species on an interface. % - % :param cov: + % parameter 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". @@ -116,22 +110,22 @@ classdef Interface < handle norm_flag = 1; end - surf_id = s.th.tr_id; - nsp = s.th.nSpecies; + surf_id = s.tr_id; + nsp = s.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); + calllib(ct, '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); + calllib(ct, 'surf_setCoveragesByName', surf_id, cov); end end diff --git a/interfaces/matlab_experimental/Base/Kinetics.m b/interfaces/matlab_experimental/Base/Kinetics.m index 9193cd9bd..7410035c4 100644 --- a/interfaces/matlab_experimental/Base/Kinetics.m +++ b/interfaces/matlab_experimental/Base/Kinetics.m @@ -1,7 +1,8 @@ classdef Kinetics < handle + properties owner - id + kin_id Kc % equilibrium constant Kf % forward reaction rate Kr % reverse reaction rate @@ -9,9 +10,7 @@ classdef Kinetics < handle dS % entropy of reaction dG % gibbs free energy of reaction end - properties(Constant = true) - lib = 'cantera_shared' - end + methods %% Kinetics class constructor @@ -41,78 +40,78 @@ classdef Kinetics < handle end end end - kin.id = calllib(kin.lib, 'kin_newFromFile', src, id, ... + kin.kin_id = calllib(ct, 'kin_newFromFile', src, id, ... iph, inb1, inb2, inb3, inb4); end %% Utility methods - function clear(kin) + function kin_clear(kin) % Delete the kernel object checklib; - calllib(kin.lib, 'kin_del', kin.id); + calllib(ct, 'kin_del', kin.kin_id); end %% Get scalar attributes - function n = nReactions(kin) - % Get the number of reactions. + function n = isReversible(kin, i) + % Get an array of flags indicating reversibility of a reaction. % - % :return: - % Integer number of reactions + % parameter i: + % Integer reaction number. + % return: + % 1 if reaction number i is reversible. 0 if irreversible. checklib; - n = calllib(kin.lib, 'kin_nReactions', kin.id); + n = calllib(ct, 'kin_isReversible', kin.kin_id, i); end function n = multiplier(kin, irxn) % Get the multiplier for reaction rate of progress. % - % :param irxn: + % parameter irxn: % Integer reaction number for which the multiplier is % desired. - % :return: + % return: % Multiplier of the rate of progress of reaction irxn. checklib; - n = calllib(kin.lib, 'kin_multiplier', kin.id, irxn-1); + n = calllib(ct, 'kin_multiplier', kin.kin_id, irxn-1); end - function n = nSpecies(kin) + function n = nReactions(kin) + % Get the number of reactions. + % + % return: + % Integer number of reactions + + checklib; + n = calllib(ct, 'kin_nReactions', kin.kin_id); + end + + function n = nSpecies2(kin) % Get the total number of species. % - % :return: + % 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); + n = calllib(ct, 'kin_nSpecies', kin.kin_id); end function n = stoich_r(kin, species, rxns) % Get the reactant stoichiometric coefficients. % - % :param species: + % parameter species: % Species indices for which reactant stoichiometric % coefficients should be retrieved. Optional argument; if % specified, "rxns" must be specified as well. - % :param rxns: + % parameter rxns: % Reaction indicies for which reactant stoichiometric % coefficients should be retrieved. Optional argument; if % specified, "species" must be specified as well. - % :return: + % 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 @@ -137,8 +136,8 @@ classdef Kinetics < handle for k = krange for i = irange - t = calllib(kin.lib, 'kin_reactantStoichCoeff', ... - kin.id, k-1, i-1); + t = calllib(ct, 'kin_reactantStoichCoeff', ... + kin.kin_id, k-1, i-1); if t ~= 0.0 temp(k, i) = t; end @@ -151,15 +150,15 @@ classdef Kinetics < handle function n = stoich_p(kin, species, rxns) % Get the product stoichiometric coefficients. % - % :param species: + % parameter species: % Species indices for which product stoichiometric % coefficients should be retrieved. Optional argument; if % specified, "rxns" must be specified as well. - % :param rxns: + % parameter rxns: % Reaction indicies for which product stoichiometric % coefficients should be retrieved. Optional argument; if % specified, "species" must be specified as well. - % :return: + % 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 @@ -184,8 +183,8 @@ classdef Kinetics < handle for k = krange for i = irange - t = calllib(kin.lib, 'kin_productStoichCoeff', ... - kin.id, k-1, i-1); + t = calllib(ct, 'kin_productStoichCoeff', ... + kin.kin_id, k-1, i-1); if t ~= 0.0 temp(k, i) = t; end @@ -198,15 +197,15 @@ classdef Kinetics < handle function n = stoich_net(kin, species, rxns) % Get the net stoichiometric coefficients. % - % :param species: + % parameter species: % Species indices for which net stoichiometric coefficients % should be retrieved. Optional argument; if specified, % "rxns" must be specified as well. - % :param rxns: + % parameter rxns: % Reaction indicies for which net stoichiometric % coefficients should be retrieved. Optional argument; if % specified, "species" must be specified as well. - % :return: + % 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 @@ -227,10 +226,82 @@ classdef Kinetics < handle %% Get reaction array attributes + function cdot = creationRates(kin) + % Get the chemical reaction rates. + % + % return: + % Returns a vector of the creation rates of all species. If + % the output is not assigned to a variable, a bar graph is + % produced. Unit: kmol/m^3-s. + + checklib; + nsp = kin.nSpecies; + xx = zeros(1, nsp); + pt = libpointer('doublePtr', xx); + calllib(ct, 'kin_getCreationRates', kin.kin_id, nsp, pt); + cdot = pt.Value; + if nargout == 0 + figure + set(gcf, 'Name', 'Creation Rates') + bar(q) + xlabel('Species Number') + ylabel('Creation Rate [kmol/m^3-s]') + title('Species Chemical Reaction Rates') + end + end + + function ddot = destructionRates(kin) + % Get the chemical destruction rates. + % + % return: + % Returns a vector of the destruction rates of all species. + % If the output is not assigned to a variable, a bar graph + % is produced. Unit: kmol/m^3-s. + + checklib; + nsp = kin.nSpecies; + xx = zeros(1, nsp); + pt = libpointer('doublePtr', xx); + calllib(ct, 'kin_getDestructionRates', kin.kin_id, nsp, pt); + ddot = pt.Value; + if nargout == 0 + figure + set(gcf, 'Name', 'Destruction Rates') + bar(q) + xlabel('Species Number') + ylabel('Destruction Rate [kmol/m^3-s]') + title('Species Chemical Reaction Rates') + end + end + + function wdot = netProdRates(kin) + % Get the net chemical production rates for all species. + % + % return: + % Returns a vector of the net production (creation-destruction) + % rates of all species. If the output is not assigned to a + % variable, a bar graph is produced. Unit: kmol/m^3-s. + + checklib; + nsp = kin.nSpecies; + xx = zeros(1, nsp); + pt = libpointer('doublePtr', xx); + calllib(ct, 'kin_getNetProductionRates', kin.kin_id, nsp, pt); + wdot = pt.Value; + if nargout == 0 + figure + set(gcf, 'Name', 'Production Rates') + bar(q) + xlabel('Species Number') + ylabel('Net Production Rate [kmol/m^3-s]') + title('Species Net Chemical Reaction Rates') + end + end + function q = rop_f(kin) % Get the forward rates of progress for all reactions. % - % :return: + % 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. @@ -239,7 +310,7 @@ classdef Kinetics < handle nr = kin.nReactions; xx = zeros(1, nr); pt = libpointer('doublePtr', xx); - calllib(kin.lib, 'kin_getFwdRateOfProgress', kin.id, nr, pt); + calllib(ct, 'kin_getFwdRateOfProgress', kin.kin_id, nr, pt); q = pt.Value; if nargout == 0 figure @@ -254,7 +325,7 @@ classdef Kinetics < handle function q = rop_r(kin) % Get the reverse rates of progress for all reactions. % - % :return: + % 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. @@ -263,7 +334,7 @@ classdef Kinetics < handle nr = kin.nReactions; xx = zeros(1, nr); pt = libpointer('doublePtr', xx); - calllib(kin.lib, 'kin_getRevRateOfProgress', kin.id, nr, pt); + calllib(ct, 'kin_getRevRateOfProgress', kin.kin_id, nr, pt); q = pt.Value; if nargout == 0 figure @@ -278,7 +349,7 @@ classdef Kinetics < handle function q = rop(kin) % Get the forward and reverse rates of progress. % - % :return: + % 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 @@ -303,7 +374,7 @@ classdef Kinetics < handle function q = rop_net(kin) % Get the net rates of progress for all reactions. % - % :return: + % 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. @@ -312,7 +383,7 @@ classdef Kinetics < handle nr = kin.nReactions; xx = zeros(1, nr); pt = libpointer('doublePtr', xx); - calllib(kin.lib, 'kin_getNetRateOfProgress', kin.id, nr, pt); + calllib(ct, 'kin_getNetRatesOfProgress', kin.kin_id, nr, pt); q = pt.Value; if nargout == 0 figure @@ -324,10 +395,91 @@ classdef Kinetics < handle end end + function rxn = reactionEqn(kin, irxn) + % Get the reaction equation of a reaction + % + % parameter irxn: + % Optional. Integer or vector of reaction numbers. + % return: + % String or cell arrray of strings of the reaction + % equations. + + if nargin == 1 + m = kin.nReactions; + n = 1 + irxn = (n:m)'; + elseif nargin == 2 + if isa(irxn, 'double') + [m, n] = size(irxn); + else + error('reaction numbers must be numeric'); + end + end + + rxn = cell(m, n); + for i = 1:m + for j = 1:n + buflen = calllib(ct, 'kin_getReactionString', kin.kin_id, ... + irxn - 1, 0, ''); + if buflen > 0 + aa = char(zeros(1, buflen)); + [~, aa] = calllib(ct, 'kin_getReactionString', ... + kin.kin_id, irxn - 1, buflen, aa); + rxn{i, j} = aa; + end + end + end + end + + function enthalpy = get.dH(kin) + % Get the enthalpy of reaction for each reaction. + % + % return: + % Returns a vector of the enthalpy of reaction for each + % reaction. Unit: J/kmol. + + checklib; + nr = kin.nReactions; + xx = zeros(1, nr); + pt = libpointer('doublePtr', xx); + calllib(ct, 'kin_getDelta', kin.kin_id, 0, nr, pt); + enthalpy = pt.Value; + end + + function entropy = get.dS(kin) + % Get the entropy of reaction for each reaction. + % + % return: + % Returns a vector of the entropy of reaction for each + % reaction. Unit: J/kmol-K. + + checklib; + nr = kin.nReactions; + xx = zeros(1, nr); + pt = libpointer('doublePtr', xx); + calllib(ct, 'kin_getDelta', kin.kin_id, 2, nr, pt); + entropy = pt.Value; + end + + function gibbs = get.dG(kin) + % Get the Gibbs free energy of reaction for each reaction. + % + % return: + % Returns a vector of the Gibbs free energy of reaction for + % each reaction. Unit: J/kmol. + + checklib; + nr = kin.nReactions; + xx = zeros(1, nr); + pt = libpointer('doublePtr', xx); + calllib(ct, 'kin_getDelta', kin.kin_id, 1, nr, pt); + gibbs = pt.Value; + end + function k = get.Kc(kin) % Get the equilibrium constants for all reactions. % - % :return: + % 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 @@ -338,12 +490,12 @@ classdef Kinetics < handle nr = kin.nReactions; xx = zeros(1, nr); pt = libpointer('doublePtr', xx); - calllib(kin.lib, 'kin_getEquilibriumConstants', kin.id, nr, pt); + calllib(ct, 'kin_getEquilibriumConstants', kin.kin_id, nr, pt); k = pt.Value; if nargout == 0 figure set(gcf, 'Name', 'Equilibrium Constants') - bar(q) + bar(k) xlabel('Reaction Number') ylabel('log_{10} Kc [kmol,m, s]') title('Equilibrium Constants') @@ -353,7 +505,7 @@ classdef Kinetics < handle function k = get.Kf(kin) % Get the forward reaction rate constants for all reactions. % - % :return: + % return: % Returns a column vector of the forward rates constants of % all of the reactions. @@ -361,14 +513,14 @@ classdef Kinetics < handle nr = kin.nReactions; xx = zeros(1, nr); pt = libpointer('doublePtr', xx); - calllib(kin.lib, 'kin_getFwdRateConstants', kin.id, nr, pt); + calllib(ct, 'kin_getFwdRateConstants', kin.kin_id, nr, pt); k = pt.Value; end function k = get.Kr(kin) % Get the reverse reaction rate constants for all reactions. % - % :return: + % return: % Returns a column vector of the reverse rates constants of % all of the reactions. @@ -376,172 +528,33 @@ classdef Kinetics < handle nr = kin.nReactions; xx = zeros(1, nr); pt = libpointer('doublePtr', xx); - calllib(kin.lib, 'kin_getRevRateConstants', kin.id, 1, nr, pt); + calllib(ct, 'kin_getRevRateConstants', kin.kin_id, 1, nr, pt); k = pt.Value; 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) + function massProdRate = ydot(kin) % Get the mass production rates of the species. % - % :return: + % 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; + calllib(ct, 'kin_getSourceTerms', kin.kin_id, nsp, pt); + massProdRate = 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: + % parameter irxn: % Integer of vector reaction numbers for which the % multiplier should be set. Optional. - % :param v: + % parameter v: % Value by which the reaction rate of progress should be % multiplied. @@ -558,7 +571,7 @@ classdef Kinetics < handle for i = 1:nr for j = 1:n - calllib(kin.lib, 'kin_setMultiplier', kin.id, ... + calllib(ct, 'kin_setMultiplier', kin.kin_id, ... irxn(i, j)-1, v); end end @@ -567,10 +580,10 @@ classdef Kinetics < handle function advanceCoeverages(kin, dt) % Advance the surface coveages forward in time % - % :param dt: + % parameter dt: % Time interval by which the coverages should be advanced. - calllib(kin.lib, 'kin_advanceCoverages', kin.id, dt); + calllib(ct, 'kin_advanceCoverages', kin.kin_id, dt); end end diff --git a/interfaces/matlab_experimental/Base/Mixture.m b/interfaces/matlab_experimental/Base/Mixture.m index 168a1adbb..ffea131bf 100644 --- a/interfaces/matlab_experimental/Base/Mixture.m +++ b/interfaces/matlab_experimental/Base/Mixture.m @@ -7,10 +7,6 @@ classdef Mixture < handle P end - properties(Constant = true) - lib = 'cantera_shared' - end - methods %% Mixture class constructor @@ -38,9 +34,9 @@ classdef Mixture < handle % store its own state information locally, and syncrhonizes the % phase objects whenever itrequires phase properties. % - % :param phases: + % parameter phases: % Cell array of phases and mole numbers. - % :return: + % return: % Instance of class 'Mixture'. checklib; @@ -50,7 +46,7 @@ classdef Mixture < handle end % Create an empty mixture. - m.mixindex = calllib(m.lib, 'mix_new'); + m.mixindex = calllib(ct, 'mix_new'); m.phases = phases; % If phases are supplied, add them @@ -78,7 +74,7 @@ classdef Mixture < handle function display(m) % Display the state of the mixture on the terminal. - calllib(m.lib, 'mix_updatePhases', m.mixindex); + calllib(ct, 'mix_updatePhases', m.mixindex); [np, nc] = size(m.phases); for n = 1:np s = [sprintf('\n******************* Phase %d', n) ... @@ -93,7 +89,7 @@ classdef Mixture < handle % Delete the MultiPhase object. checklib; - calllib(m.lib, 'mix_del', m.mixindex); + calllib(ct, 'mix_del', m.mixindex); end %% Mixture Get methods @@ -101,11 +97,11 @@ classdef Mixture < handle function addPhase(m, phase, moles) % Add a phase to the mixture % - % :param m: + % parameter m: % Instance of class 'Mixture' to which phases is added. - % :param phase: + % parameter phase: % Instance of class 'ThermoPhase' which should be added. - % :param moles: + % parameter moles: % Number of moles of the phase to be added. Unit: kmol. checklib; @@ -120,7 +116,7 @@ classdef Mixture < handle error('Negative moles'); end - iok = calllib(m.lib, 'mix_addPhase', m.mixindex, phase.tp_id, ... + iok = calllib(ct, 'mix_addPhase', m.mixindex, phase.tp_id, ... moles); if iok < 0 error('Error adding phase'); @@ -130,39 +126,39 @@ classdef Mixture < handle function temperature = get.T(m) % Get the temperature of the mixture. % - % :return: + % return: % Temperature in K. checklib; - temperature = calllib(m.lib, 'mix_temperature', m.mixindex); + temperature = calllib(ct, 'mix_temperature', m.mixindex); end function pressure = get.P(m) % Get the pressure of themixture. % - % :return: + % return: % Pressure in Pa. checklib; - pressure = calllib(m.lib, 'mix_pressure', m.mixindex); + pressure = calllib(ct, '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); + n = calllib(ct, '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); + n = calllib(ct, '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); + n = calllib(ct, 'mix_nSpecies', m.mixindex); end function n = elementIndex(m, name) @@ -172,7 +168,7 @@ classdef Mixture < handle % 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; + n = calllib(ct, 'mix_elementIndex', m.mixindex, name) + 1; end function n = speciesIndex(m, k, p) @@ -182,26 +178,26 @@ classdef Mixture < handle % indices start from 1 instead of 0 as in Cantera C++ and % Python interfaces. checklib; - n = calllib(m.lib, 'mix_speciesIndex', m.mixindex, k-1, p-1) + 1; + n = calllib(ct, 'mix_speciesIndex', m.mixindex, k-1, p-1) + 1; % check back on this one! end function moles = phaseMoles(m, n) % Get the number of moles of a phase in a mixture. % - % :param n: + % parameter n: % Integer phase number in the input. - % :return: + % return: % Moles of phase number 'n'. Unit: kmol. checklib; if nargin == 2 - moles = calllib(m.lib, 'mix_phaseMoles', m.mixindex, n); + moles = calllib(ct, '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', ... + moles(i) = calllib(ct, 'mix_phaseMoles', ... m.mixindex, i); end else error('wrong number of arguments'); @@ -211,14 +207,14 @@ classdef Mixture < handle function mu = chemPotentials(m) % Get the chemical potentials of species in the mixture. % - % :return: + % 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); + calllib(ct, 'mix_getChemPotential', m.mixindex, nsp, ptr); mu = ptr.Value; end @@ -227,30 +223,30 @@ classdef Mixture < handle function m = set.T(m, temp) % Set the mixture temperature. % - % :param temp: + % parameter temp: % Temperature to set. Unit: K. checklib; - calllib(m.lib, 'mix_setTemperature', m.mixindex, temp); + calllib(ct, 'mix_setTemperature', m.mixindex, temp); end function m = set.P(m, pressure) % Set the mixture pressure. % - % :param pressure: + % parameter pressure: % Pressure to set. Unit: Pa. checklib; - calllib(m.lib, 'mix_setPressure', m.mixindex, pressure); + calllib(ct, 'mix_setPressure', m.mixindex, pressure); end function setPhaseMoles(m, n, moles) % Set the number of moles of phase n in the mixture. % - % :param n: + % parameter n: % Phase number. - % :param moles: + % parameter moles: % Number of moles to set. Unit: kmol. checklib; - calllib(m.lib, 'mix_setPhaseMoles', m.mixindex, n-1, moles); + calllib(ct, 'mix_setPhaseMoles', m.mixindex, n-1, moles); end function setSpeciesMoles(m, moles) @@ -260,10 +256,10 @@ classdef Mixture < handle % in the mixture. Note that the species may belong to any % phase, and unspecified species are set to zero. % - % :param moles: + % parameter moles: % Vector or string specifying the moles of species. checklib; - calllib(m.lib, 'mix_setMolesByName', m.mixindex, moles); + calllib(ct, 'mix_setMolesByName', m.mixindex, moles); % check back on this one! end @@ -285,25 +281,25 @@ classdef Mixture < handle % >> equilibrate(mix, 'TP); % >> equilibrate('TP', 1.0e-6, 500); % - % :param XY: + % parameter XY: % Two-letter string specifying the two properties to hold % fixed. Currently 'TP', 'HP', 'TV', and 'SP' have been % implemented. Default: 'TP'. - % :param err: + % parameter err: % Error tolerance. Iteration will continue until delta_Mu/RT % is less than this value for each reaction. Default: % 1.0e-9. - % :param maxsteps: + % parameter maxsteps: % Maximum number of steps to take while solving the % equilibrium problem for specified T and P. Default: 1000. - % :param maxiter: + % parameter 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: + % parameter loglevel: % Set to a value > 0 to write diagnostic output. Larger % values generate more detailed information. - % :return: + % return: % The error in the solution. checklib; @@ -323,7 +319,7 @@ classdef Mixture < handle if nargin < 2 XY = 'TP' end - r = calllib(m.lib, 'mix_equilibrate', m.mixindex, XY, err, ... + r = calllib(ct, 'mix_equilibrate', m.mixindex, XY, err, ... maxsteps, maxiter, loglevel); end diff --git a/interfaces/matlab_experimental/Base/Solution.m b/interfaces/matlab_experimental/Base/Solution.m index 8523998d3..8ba8f8106 100644 --- a/interfaces/matlab_experimental/Base/Solution.m +++ b/interfaces/matlab_experimental/Base/Solution.m @@ -1,8 +1,6 @@ -classdef Solution < handle +classdef Solution < handle & ThermoPhase & Kinetics & Transport properties - thermo - kinetics - transport + tp end methods % Solution class constructor @@ -11,27 +9,26 @@ classdef Solution < handle id = '-'; end tp = ThermoPhase(src, id); - kin = Kinetics(tp, src, id); - s.kinetics = kin; - s.thermo = tp; + s@ThermoPhase(src, id); + s@Kinetics(tp, src, id); if nargin == 3 - if (strcmp(trans, 'default') || strcmp(trans, 'None')... + 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); + trans = 'default'; end - s.transport = tr; + s@Transport(tp, trans, 0); + s.tp_id = tp.tp_id; end % Delete the kernel objects associated with a solution function clear(s) - s.thermo.clear; - s.kinetics.clear; - s.transport.clear; + s.tp_clear; + s.kin_clear; + s.tr_clear; end + end end diff --git a/interfaces/matlab_experimental/Base/ThermoPhase.m b/interfaces/matlab_experimental/Base/ThermoPhase.m index bdbb88d26..ff49524f4 100644 --- a/interfaces/matlab_experimental/Base/ThermoPhase.m +++ b/interfaces/matlab_experimental/Base/ThermoPhase.m @@ -11,14 +11,11 @@ classdef ThermoPhase < handle H % enthalpy S % entropy U % internal energy -% V % specific volume not added yet + G % Gibbs free energy + V % specific volume basis end - properties (Constant = true) - lib = 'cantera_shared' - end - properties (Dependent) DP DPX @@ -26,20 +23,42 @@ classdef ThermoPhase < handle HP HPX HPY + PV + PVX + PVY + SH + SHX + SHY SP SPX SPY -% SVX -% SVY + ST + STX + STY + SV + SVX + SVY TD TDX TDY + TH + THX + THY TP TPX TPY -% UV -% UVX -% UVY + TV + TVX + TVY + UV + UVX + UVY + UP + UPX + UPY + VH + VHX + VHY end methods @@ -54,7 +73,7 @@ classdef ThermoPhase < handle id = '-'; end tp.tp_owner = 1; - tp.tp_id = calllib(tp.lib, 'thermo_newFromFile', src, id); + tp.tp_id = calllib(ct, 'thermo_newFromFile', src, id); tp.basis = 'molar'; end @@ -66,14 +85,14 @@ classdef ThermoPhase < handle if nargin < 2 || ~isnumeric(threshold) threshold = 1e-14; end - calllib(tp.lib, 'thermo_print', tp.tp_id, 1, threshold); + calllib(ct, 'thermo_print', tp.tp_id, 1, threshold); end - function clear(tp) + function tp_clear(tp) % Delete the kernel object. checklib; - calllib(tp.lib, 'thermo_del', tp.tp_id); + calllib(ct, 'thermo_del', tp.tp_id); end function tp = set.basis(tp, b) @@ -82,121 +101,43 @@ classdef ThermoPhase < handle if strcmp(b, 'mole') || strcmp(b, 'molar') ... || strcmp(b, 'Mole') || strcmp(b, 'Molar') - tp.basis = 'molar' + tp.basis = 'molar'; elseif strcmp(b, 'mass') || strcmp(b, 'Mass') - tp.basis = '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. + function amu = atomicMasses(tp) + % Get the atomic masses of the elements. % - % :return: - % Double temperature. Unit: K + % return: + % Vector of element atomic masses. Unit: kg/kmol checklib; - temperature = calllib(tp.lib, 'thermo_temperature', tp.tp_id); + nel = tp.nElements; + aa = zeros(1, nel); + pt = libpointer('doublePtr', aa); + calllib(ct, 'thermo_getAtomicWeights', ... + tp.tp_id, nel, pt); + amu = pt.Value; end - function pressure = get.P(tp) - % Get the pressure. + function e = charges(tp) + % Get the array of species charges. % - % :return: - % Double pressure. Unit: Pa + % return: + % Vector of species charges. Unit: elem. charge 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 + nsp = tp.nSpecies; + yy = zeros(1, nsp); + pt = libpointer('doublePtr', yy); + calllib(ct, 'thermo_getCharges', ... + tp.tp_id, nsp, pt); + e = pt.Value; end function k = elementIndex(tp, name) @@ -208,9 +149,9 @@ classdef ThermoPhase < handle % indices start from 1 instead of 0 as in Cantera C++ and % Python interfaces. % - % :param name: + % parameter name: % String or cell array of elements whose index is requested - % :return: + % return: % Integer number of elements in the phase. checklib; @@ -219,7 +160,7 @@ classdef ThermoPhase < handle k = zeros(m, n); for i = 1:m for j = 1:n - k(i, j) = calllib(tp.lib, 'thermo_elementIndex', ... + k(i, j) = calllib(ct, 'thermo_elementIndex', ... tp.tp_id, name{i, j}) + 1; if k(i, j) > 1e3 warning(['Element ', name{i, j}, ... @@ -229,7 +170,7 @@ classdef ThermoPhase < handle end end elseif ischar(name) - k = calllib(tp.lib, 'thermo_elementIndex', ... + k = calllib(ct, 'thermo_elementIndex', ... tp.tp_id, name) + 1; if k > 1e3 warning(['Element ', name, ... @@ -241,14 +182,86 @@ classdef ThermoPhase < handle end end + function elMassFrac = elementalMassFraction(tp, element) + % Determine the elemental mass fraction in gas object. + checklib; + + % Check input parameters. + if nargin ~= 2 + error('elementalMassFraction expects two input arguments.'); + end + if ~isIdealGas(tp) + error('Gas object must represent an ideal gas mixture.'); + end + if ~ischar(element) + error('Element name must be of format character.'); + end + + % Calculate the elemental mass fraction in a gas object using + % the following equation: + % + % elMassFrac = sum of nAtoms(k, m)*Mel(m)*Y(k)/mw(k) + % + % where nAtoms(k, m) is the number of atoms of element 'm', in + % species 'k'; Mel(m) is the atomic weight of element 'm'; Y(k) + % is the mass fraction of species 'k'; and mw(k) is the + % molecular weight of species 'k'. + + n = tp.nSpecies; + spec = tp.speciesNames; + eli = tp.elementIndex(element); + M = tp.atomicMasses; + Mel = M(eli); + MW = tp.MolecularWeights; + % Initialize the element mass fraction as zero. + elMassFrac = 0.0; + % Use loop to perform summation of elemental mass fraction over all species. + for i = 1:n + natoms(i) = tp.nAtoms(spec{i}, element); + yy(i) = tp.massFraction(spec{i}); + elMassFrac = elMassFrac + (natoms(i)*Mel*yy(i))/MW(i); + end + end + + function mmw = meanMolecularWeight(tp) + % Get the mean molecular weight. + % + % return: + % Double mean molecular weight. Unit: kg/kmol + + checklib; + mmw = calllib(ct, 'thermo_meanMolecularWeight', tp.tp_id); + end + + function density = molarDensity(tp) + % Get the molar basis density in kmol/m^3. + checklib; + density = calllib(ct, 'thermo_molarDensity', tp.tp_id); + 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(ct, 'thermo_getMolecularWeights', ... + tp.tp_id, nsp, pt); + mw = pt.Value; + end + function n = nAtoms(tp, species, element) % Get the number of atoms of an element in a species. % - % :param k: + % parameter k: % String species name or integer species number. - % :param m: + % parameter m: % String element name or integer element number. - % :return: + % return: % Integer number of atoms of the element in the species. if nargin == 3 if ischar(species) @@ -267,23 +280,149 @@ classdef ThermoPhase < handle n = -1; return end - n = calllib(tp.lib, 'thermo_nAtoms', tp.tp_id, k-1, m-1); + n = calllib(ct, 'thermo_nAtoms', tp.tp_id, k-1, m-1); else error('Two input arguments required.') end end + function nel = nElements(tp) + % Get the number of elements. + % + % return: + % Integer number of elements in the phase. + + checklib; + nel = calllib(ct, '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(ct, '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. + % + % parameter 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(ct, '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(ct, '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 nm = speciesName(tp, k) + % Get the name of a species given the index. + % + % parameter 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(ct, 'thermo_getSpeciesName', ... + tp.tp_id, ksp, 0, ''); + if buflen > 0 + aa = char(zeros(1, buflen)); + [~, aa] = calllib(ct, 'thermo_getSpeciesName', ... + tp.tp_id, ksp, buflen, aa); + nm{i, j} = aa; + end + end + end + end + + function n = speciesNames(tp) + % Get all species names. + n = tp.speciesName(1:tp.nSpecies); + end + + function temperature = get.T(tp) + % Get the temperature. + % + % return: + % Double temperature. Unit: K + + checklib; + temperature = calllib(ct, 'thermo_temperature', tp.tp_id); + end + + function pressure = get.P(tp) + % Get the pressure. + % + % return: + % Double pressure. Unit: Pa + + checklib; + pressure = calllib(ct, 'thermo_pressure', tp.tp_id); + end + + function density = get.D(tp) + % Get the mass basis density in kg/m^3. + checklib; + density = calllib(ct, 'thermo_density', tp.tp_id); + end + + function volume = get.V(tp) + % Get the specific volume depending on the basis. + % + % return: + % Density depending on the basis. Units: + % m^3/kmol (molar) m^3/kg (mass). + volume = 1/tp.D; + end + function moleFractions = get.X(tp) % Get the mole fractions of all species. % - % :return: + % return: % Vector of species mole fractions. checklib; nsp = tp.nSpecies; xx = zeros(1, nsp); pt = libpointer('doublePtr', xx); - calllib(tp.lib, 'thermo_getMoleFractions', ... + calllib(ct, 'thermo_getMoleFractions', ... tp.tp_id, nsp, pt); moleFractions = pt.Value; @@ -301,15 +440,15 @@ classdef ThermoPhase < handle function x = moleFraction(tp, species) % Get the mole fraction of one or a list of species. % - % :param species: + % parameter species: % String or cell array of species whose mole fraction is % requested. - % :return: + % return: % Scalar or vector of species mole fractions. xarray = tp.X; if isa(species, 'char') - k = tp.speciesIndex(tp, species); + k = tp.speciesIndex(species); if k > 0 x = xarray(k); else error("species not found."); @@ -318,7 +457,7 @@ classdef ThermoPhase < handle n = length(species); x = zeros(1, n); for j = 1:n - k = tp.speciesIndex(tp, species{j}); + k = tp.speciesIndex(species{j}); if k > 0 x(j) = xarray(k); else error("species not found."); @@ -330,14 +469,14 @@ classdef ThermoPhase < handle function massFractions = get.Y(tp) % Get the mass fractions of all species. % - % :return: + % return: % Vector of species mass fractions. checklib; nsp = tp.nSpecies; yy = zeros(1, nsp); pt = libpointer('doublePtr', yy); - calllib(tp.lib, 'thermo_getMassFractions', ... + calllib(ct, 'thermo_getMassFractions', ... tp.tp_id, nsp, pt); massFractions = pt.Value; @@ -355,121 +494,61 @@ classdef ThermoPhase < handle function y = massFraction(tp, species) % Get the mass fraction of one or a list of species. % - % :param species: + % parameter species: % String or cell array of species whose mass fraction is % requested. - % :return: + % return: % Scalar or vector of species mass fractions. - yarray = tp.Y; + yy = tp.Y; if isa(species, 'char') - k = tp.speciesIndex(tp, species); + k = tp.speciesIndex(species); if k > 0 - y = yarray(k); + y = yy(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}); + k = tp.speciesIndex(species{j}); if k > 0 - y(j) = yarray(k); + y(j) = yy(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 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(ct, 'thermo_chemPotentials', ... + tp.tp_id, nsp, pt); + mu = pt.Value; + end + function c = cv(tp) % Get the specific heat at constant volume. % - % :return: + % 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); + c = calllib(ct, 'thermo_cv_mole', tp.tp_id); elseif strcmp(tp.basis, 'mass') - c = calllib(tp.lib, 'thermo_cv_mass', tp.tp_id); + c = calllib(ct, 'thermo_cv_mass', tp.tp_id); else error("basis not specified"); end end @@ -477,47 +556,243 @@ classdef ThermoPhase < handle function c = cp(tp) % Get the specific heat at constant pressure. % - % :return: + % 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); + c = calllib(ct, 'thermo_cp_mole', tp.tp_id); elseif strcmp(tp.basis, 'mass') - c = calllib(tp.lib, 'thermo_cp_mass', tp.tp_id); + c = calllib(ct, 'thermo_cp_mass', tp.tp_id); else error("basis not specified"); end end + function d = critDensity(tp) + % Get the critical density. + % + % return: + % Critical density. Unit: K. + + checklib; + d = calllib(ct, 'thermo_critDensity', tp.tp_id); + end + + function p = critPressure(tp) + % Get the critical pressure. + % + % return: + % Critical temperature. Unit: Pa. + + checklib; + p = calllib(ct, 'thermo_critPressure', tp.tp_id); + end + + function t = critTemperature(tp) + % Get the critical temperature. + % + % return: + % Critical temperature. Unit: K. + checklib; + t = calllib(ct, 'thermo_critTemperature', tp.tp_id); + end + + function v = electricPotential(tp) + % Get the electric potential + % + % return: + % Electric potential of the phase. Unit: V. + + checklib; + v = calllib(ct, 'thermo_electricPotential', tp.tp_id); + end + + function e = eosType(tp) + % Get the type of the equation of state +% checklib; +% buflen = calllib(ct, 'thermo_getEosType', tp.tp_id, 0, ''); +% if buflen > 0 +% aa = char(zeros(1, buflen)); +% [~, aa] = calllib(ct, 'thermo_getEosType', ... +% tp.tp_id, buflen, aa); +% end +% e = aa; + e = 'IdealGas'; + end + + function v = isIdealGas(tp) + % Get a flag indicating whether the phase is an ideal gas. + + if strcmp(tp.eosType, 'IdealGas') + v = 1; + else + v = 0; + end + v = 1; + end + + function b = isothermalCompressibility(tp) + % Get the isothermal compressibility + % + % return: + % Isothermal compressibility. Unit: 1/Pa. + + checklib; + b = calllib(ct, 'thermo_isothermalCompressibility', tp.tp_id); + end + + function t = maxTemp(tp) + % Get the maximum temperature of the parameter fits. + % + % return: + % Vector of maximum temperatures of all species. + + checklib; + t = calllib(ct, 'thermo_maxTemp', tp.tp_id, -1); + end + + function t = minTemp(tp) + % Get the minimum temperature of the parameter fits. + % + % return: + % Vector of minimum temperatures of all species. + + checklib; + t = calllib(ct, 'thermo_minTemp', tp.tp_id, -1); + end + + function p = P_sat(tp, t) + % Get the saturation pressure for a given temperature. + % + % parameter t: + % Temperature. Unit: K. + % return: + % Saturation pressure for temperature t. Unit: Pa. + + checklib; + p = calllib(ct, 'thermo_satPressure', tp.tp_id, t); + end + + function p = refPressure(tp) + % Get the reference pressure. + % + % return: + % Reference pressure. Unit: Pa. + + checklib; + p = calllib(ct, 'thermo_refPressure', tp.tp_id, -1); + end + + function c = soundspeed(tp) + % Get the speed of sound + % If the phase is an ideal gas, the speed of sound is + % calculated by: + % + % c = sqrt(gamma * R * T); + % + % where gamma is the ratio of specific heat, R is the specific + % gas constant, and T is the temperature. If the phase is not + % an ideal gas, the speed of sound is calculated by: + % + % c = sqrt( + % + % return: + % The speed of sound. Unit: m/s + + checklib; + if tp.isIdealGas + tp.basis = 'mass'; + gamma = tp.cp/tp.cv; + wtm = tp.meanMolecularWeight; + r = 8314.4621/wtm; + c = sqrt(gamma * r *tp.T); + else + rho0 = tp.D; + p0 = tp.P; + s0 = tp.S; + rho1 = 1.001 * rho0; + tp.DS = {rho1, s0}; + p1 = tp.P; + dpdrho_s = (p1 - p0) / (rho1 - rho0); + c = sqrt(dpdrho_s); + end + end + + function a = thermalExpansionCoeff(tp) + % Get the thermal expansion coefficient. + % + % return: + % Thermal expansion coefficient. Unit: 1/K. + + checklib; + a = calllib(ct, 'thermo_thermalExpansionCoeff', tp.tp_id); + end + + function t = T_sat(tp, p) + % Get the saturation temperature for a given pressure. + % + % parameter p: + % Pressure. Unit: Pa. + % return: + % Saturation temperature for pressure p. Unit: K. + + checklib; + t = calllib(ct, 'thermo_satTemperature', tp.tp_id, p); + end + + function v = vaporFraction(tp) + % Get the vapor fractions. + % + % return: + % Vapor fraction. + + checklib; + v = calllib(ct, 'thermo_vaporFraction', tp.tp_id); + end + function enthalpy = get.H(tp) % Get the enthalpy. % - % :return: + % 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); + enthalpy = calllib(ct, 'thermo_enthalpy_mole', tp.tp_id); elseif strcmp(tp.basis, 'mass') - enthalpy = calllib(tp.lib, 'thermo_enthalpy_mass', tp.tp_id); + enthalpy = calllib(ct, 'thermo_enthalpy_mass', tp.tp_id); else error("basis not specified"); end end + function enthalpy = enthalpies_RT(tp) + % Get the non-dimensional enthalpy. + % + % return: + % Vector of standard-state species enthalpies divided by RT. + + checklib; + nsp = tp.nSpecies; + xx = zeros(1, nsp); + pt = libpointer('doublePtr', xx); + calllib(ct, 'thermo_getEnthalpies_RT', tp.tp_id, nsp, pt); + enthalpy = pt.Value; + end + function entropy = get.S(tp) % Get the entropy. % - % :return: + % 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); + entropy = calllib(ct, 'thermo_entropy_mole', tp.tp_id); elseif strcmp(tp.basis, 'mass') - entropy = calllib(tp.lib, 'thermo_entropy_mass', tp.tp_id); + entropy = calllib(ct, 'thermo_entropy_mass', tp.tp_id); else error("basis not specified"); end end @@ -525,346 +800,584 @@ classdef ThermoPhase < handle function intEnergy = get.U(tp) % Get the internal energy. % - % :return: + % 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); + intEnergy = calllib(ct, 'thermo_intEnergy_mole', tp.tp_id); elseif strcmp(tp.basis, 'mass') - intEnergy = calllib(tp.lib, 'thermo_intEnergy_mass', tp.tp_id); + intEnergy = calllib(ct, 'thermo_intEnergy_mass', tp.tp_id); else error("basis not specified"); end end - function gibbs = G(tp) + function gibbs = get.G(tp) % Get the Gibss free energy. % - % :return: + % 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); + gibbs = calllib(ct, 'thermo_gibbs_mole', tp.tp_id); elseif strcmp(tp.basis, 'mass') - gibbs = calllib(tp.lib, 'thermo_gibbs_mass', tp.tp_id); + gibbs = calllib(ct, '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]; + output = {tp.D, tp.P}; end function output = get.DPX(tp) % Get density, pressure, and mole fractions depending on the basis. - % :return: + % 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 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]; + 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]; + output = {tp.H, tp.P}; end function output = get.HPX(tp) % Get enthalpy, pressure, and mole fractions depending on the basis. - % :return: + % 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 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 - output = [tp.H, tp.P, tp.Y]; + function output = get.PV(tp) + output = {tp.P, tp.V}; + end + + function output = get.PVX(tp) + output = {tp.P, tp.V, tp.X}; + end + + function output = get.PVY(tp) + output = {tp.P, tp.V, tp.Y}; + end + + function output = get.SH(tp) + output = {tp.S, tp.H}; + end + + function output = get.SHX(tp) + output = {tp.S, tp.H, tp.X}; + end + + function output = get.SHY(tp) + output = {tp.S, tp.H, 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]; + 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]; + 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.S, tp.P, tp.Y}; + end - output = [tp.H, tp.P, tp.Y]; + function output = get.ST(tp) + output = {tp.S, tp.T}; + end + + function output = get.STX(tp) + output = {tp.S, tp.T, tp.X}; + end + + function output = get.STY(tp) + output = {tp.S, tp.T, tp.Y}; + end + + function output = get.SV(tp) + output = {tp.S, tp.V}; + end + + function output = get.SVX(tp) + output = {tp.S, tp.V, tp.X}; + end + + function output = get.SVY(tp) + output = {tp.S, tp.V, 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]; + 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]; + 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 - output = [tp.T, tp.D, tp.Y]; + function output = get.TH(tp) + output = {tp.T, tp.H}; + end + + function output = get.THX(tp) + output = {tp.T, tp.H, tp.X}; + end + + function output = get.THY(tp) + output = {tp.T, tp.H, 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]; + 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]; + 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]; + output = {tp.T, tp.P, tp.Y}; end - %% PhaseSet Single Methods + function output = get.TV(tp) + output = {tp.T, tp.V}; + end - function tp = set.T(tp, temperature) + function output = get.TVX(tp) + output ={tp.T, tp.V, tp.X}; + end + + function output = get.TVY(tp) + output = {tp.T, tp.V, tp.Y}; + end + + function output = get.UV(tp) + output = {tp.U, tp.V}; + end + + function output = get.UVX(tp) + output = {tp.U, tp.V, tp.X}; + end + + function output = get.UVY(tp) + output = {tp.U, tp.V, tp.Y}; + end + + function output = get.UP(tp) + output = {tp.U, tp.P}; + end + + function output = get.UPX(tp) + output = {tp.U, tp.P, tp.X}; + end + + function output = get.UPY(tp) + output = {tp.U, tp.P, tp.Y}; + end + + function output = get.VH(tp) + output = {tp.V, tp.H}; + end + + function output = get.VHX(tp) + output = {tp.V, tp.H, tp.X}; + end + + function output = get.VHY(tp) + output = {tp.V, tp.H, tp.Y}; + end + + %% PhaseSet single methods + + function tp = setElectricPotential(tp, phi) + % Set the electric potential in V. + calllib(ct, 'thermo_setElectricPotential', tp.tp_id, phi); + end + + function tp = setState_Psat(tp, p, q) + % Set saturated vapor + checklib; + calllib(ct, 'thermo_setState_Psat', tp.tp_id, p, q); + end + + function tp = setState_Tsat(tp, t, q) + % Set saturated liquid + checklib; + calllib(ct, 'thermo_setState_Tsat', tp.tp_id, t, 1 - q); + end + + function set.T(tp, temperature) checklib; if temperature <= 0 error('The temperature must be positive'); end - calllib(tp.lib, 'thermo_setTemperature', tp.tp_id, temperature); + calllib(ct, 'thermo_setTemperature', tp.tp_id, temperature); end - function tp = set.P(tp, pressure) + function set.P(tp, pressure) checklib; if pressure <= 0 error('The pressure must be positive'); end - calllib(tp.lib, 'thermo_setPressure', tp.tp_id, pressure); + calllib(ct, 'thermo_setPressure', tp.tp_id, pressure); end - function tp = set.D(tp, density) + function set.D(tp, density) checklib; if density <= 0 error('The density must be positive'); end - calllib(tp.lib, 'thermo_setDensity', tp.tp_id, density); + calllib(ct, 'thermo_setDensity', tp.tp_id, density); end -% function tp = set.X(tp, + function set.X(tp, xx) + checklib; + lim = 1e-9; + if isa(xx, 'double') + nsp = tp.nSpecies; + if sum(xx) - 1 <= 1e-9 + norm = 0; + else norm = 1; + end + calllib(ct, 'thermo_setMoleFractions', tp.tp_id, ... + nsp, xx, norm); + elseif isa(xx, 'char') + calllib(ct, 'thermo_setMoleFractionsByName', tp.tp_id, xx); + end + end + + function set.Y(tp, yy) + checklib; + if isa(yy, 'double') + nsp = tp.nSpecies; + if sum(yy) -1 <= 1e-9 + norm = 0; + else norm = 1; + end + calllib(ct, 'thermo_setMassFractions', tp.tp_id, ... + nsp, yy, norm); + elseif isa(yy, 'char') + calllib(ct, 'thermo_setMassFracitonsByName', tp.tp_id, yy); + end + end + + %% PhaseSet multi methods + + function set.DP(tp, input) + checklib; + d = input{1}; + p = input{2}; + if d <= 0 + error('The density must be positive'); + end + if p <= 0 + error('The pressure must be positive'); + end + calllib(ct, 'thermo_set_RP', tp.tp_id, [d, p]); + end + + function set.DPX(tp, input) + tp.X = input{3}; + tp.DP = input(1:2); + end + + function set.DPY(tp, input) + tp.Y = input{3}; + tp.DP = input(1:2); + end + + function set.HP(tp, input) + checklib; + h = input{1}; + p = input{2}; + if p <= 0 + error('The pressure must be positive'); + end + calllib(ct, 'thermo_set_HP', tp.tp_id, [h, p]); + end + + function set.HPX(tp, input) + tp.X = input{3}; + tp.HP = input(1:2); + end + + function set.HPY(tp, input) + tp.Y = input{3}; + tp.HP = input(1:2); + end + + function set.PV(tp, input) + checklib; + p = input{1}; + v = input{2}; + if p <= 0 + error('The pressure must be positive'); + end + if v <= 0 + error('The specific volume must be positive'); + end + calllib(ct, 'thermo_set_PV', tp.tp_id, [p, v]); + end + + function set.PVX(tp, input) + tp.X = input{3}; + tp.PV = input(1:2); + end + + function set.PVY(tp, input) + tp.Y = input{3}; + tp.PV = input(1:2); + end + + function set.SH(tp, input) + checklib; + s = input{1}; + h = input{2}; + calllib(ct, 'thermo_set_SH', tp.tp_id, [s, h]); + end + + function set.SHX(tp, input) + tp.X = input{3}; + tp.SH = input(1:2); + end + + function set.SHY(tp, input) + tp.Y = input{3}; + tp.SH = input(1:2); + end + + function set.SP(tp, input) + checklib; + s = input{1}; + p = input{2}; + if p <= 0 + error('The pressure must be positive'); + end + calllib(ct, 'thermo_set_SP', tp.tp_id, [s, p]); + end + + function set.SPX(tp, input) + tp.X = input{3}; + tp.SP = input(1:2); + end + + function set.SPY(tp, input) + tp.Y = input{3}; + tp.SP = input(1:2); + end + + function set.ST(tp, input) + checklib; + s = input{1}; + t = input{2}; + if t <= 0 + error('The temperature must be positive'); + end + calllib(ct, 'thermo_set_ST', tp.tp_id, [s, t]); + end + + function set.STX(tp, input) + tp.X = input{3}; + tp.ST = input(1:2); + end + + function set.STY(tp, input) + tp.Y = input{3}; + tp.ST = input(1:2); + end + + function set.SV(tp, input) + checklib; + s = input{1}; + v = input{2}; + if v <= 0 + error('The specific volume must be positive'); + end + calllib(ct, 'thermo_set_SV', tp.tp_id, [s, v]); + end + + function set.SVX(tp, input) + tp.X = input{3}; + tp.SV = input(1:2); + end + + function set.SVY(tp, input) + tp.Y = input{3}; + tp.SP = input(1:2); + end + + function set.TD(tp, input) + t = input{1}; + d = input{2}; + if t <= 0 + error('The temperature must be positive'); + end + if d <= 0 + error('The density must be positive'); + end + tp.T = t; + tp.D = d; + end + + function set.TDX(tp, input) + tp.X = input{3}; + tp.TD = input(1:2); + end + + function set.TDY(tp, input) + tp.Y = input{3}; + tp.TD = input(1:2); + end + + function set.TH(tp, input) + checklib; + t = input{1}; + if t <= 0 + error('The temperature must be positive'); + end + h = input{2}; + calllib(ct, 'thermo_set_TH', tp.tp_id, [t, h]); + end + + function set.THX(tp, input) + tp.X = input{3}; + tp.TH = input(1:2); + end + + function set.THY(tp, input) + tp.Y = input{3}; + tp.TH = input(1:2); + end + + function set.TP(tp, input) + t = input{1}; + p = input{2}; + if t <= 0 + error('The temperature must be positive'); + end + if p <= 0 + error('The pressure must be positive'); + end + tp.T = t; + tp.P = p; + end + + function set.TPX(tp, input) + tp.X = input{3}; + tp.TP = input(1:2); + end + + function set.TPY(tp, input) + tp.Y = input{3}; + tp.TP = input(1:2); + end + + function set.TV(tp, input) + checklib; + t = input{1}; + v = input{2}; + if t <= 0 + error('The temperature must be positive'); + end + if v <= 0 + error('The specific volume must be positive'); + end + calllib(ct, 'thermo_set_TV', tp.tp_id, [t, v]); + end + + function set.TVX(tp, input) + tp.X = input{3}; + tp.TV = input(1:2); + end + + function set.TVY(tp, input) + tp.Y = input{3}; + tp.TV = input(1:2); + end + + function set.UP(tp, input) + checklib; + u = input{1}; + p = input{2}; + if p <= 0 + error('The pressure must be positive'); + end + calllib(ct, 'thermo_set_UP', tp.tp_id, [u, p]); + end + + function set.UPX(tp, input) + tp.X = input{3}; + tp.UP = input(1:2); + end + + function set.UPY(tp, input) + tp.Y = input{3}; + tp.UP = input(1:2); + end + + function set.UV(tp, input) + checklib; + u = input{1}; + v = input{2}; + if v <= 0 + error('The specific volume must be positive'); + end + calllib(ct, 'thermo_set_UV', tp.tp_id, [u, v]); + end + + function set.UVX(tp, input) + tp.X = input{3}; + tp.UV = input(1:2); + end + + function set.UVY(tp, input) + tp.Y = input{3}; + tp.UV = input(1:2); + end + + function set.VH(tp, input) + checklib; + v = input{1}; + h = input{2}; + if v <= 0 + error('The specific volume must be positive'); + end + calllib(ct, 'thermo_set_VH', tp.tp_id, [v, h]); + end + + function set.VHX(tp, input) + tp.X = input{3}; + tp.VH = input(1:2); + end + + function set.VHY(tp, input) + tp.Y = input{3}; + tp.VH = input(1:2); + end function tp = equilibrate(tp, xy, solver, rtol, maxsteps, ... maxiter, loglevel) @@ -885,7 +1398,7 @@ classdef ThermoPhase < handle if nargin < 7 loglevel = 0; end - calllib(tp.lib, 'thermo_equilibrate', tp.tp_id, xy, solver, ... + calllib(ct, 'thermo_equilibrate', tp.tp_id, xy, solver, ... rtol, maxsteps, maxiter, loglevel); end diff --git a/interfaces/matlab_experimental/Base/Transport.m b/interfaces/matlab_experimental/Base/Transport.m index be92266be..178bf3330 100644 --- a/interfaces/matlab_experimental/Base/Transport.m +++ b/interfaces/matlab_experimental/Base/Transport.m @@ -1,13 +1,10 @@ classdef Transport < handle + properties th tr_id end - properties(Constant = true) - lib = 'cantera_shared' - end - methods %% Transport class constructor @@ -28,22 +25,23 @@ classdef Transport < handle else tr.th = tp; if strcmp(model, 'default') - tr.tr_id = calllib(tr.lib, 'trans_newDefault', ... + tr.tr_id = calllib(ct, 'trans_newDefault', ... tp.tp_id, loglevel); else - tr.tr_id = calllib(tr.lib, 'trans_new', model, ... + tr.tr_id = calllib(ct, 'trans_new', model, ... tp.tp_id, loglevel); end end + tr.tp_id = tp.tp_id; end %% Utility methods - function clear(tr) + function tr_clear(tr) % Delete the kernel object. checklib; - calllib(tr.lib, 'trans_del', tr.tr_id); + calllib(ct, 'trans_del', tr.tr_id); end %% Transport Methods @@ -51,11 +49,11 @@ classdef Transport < handle function v = viscosity(tr) % Get the dynamic viscosity. % - % :return: + % return: % Double dynamic viscosity. Unit: Pa*s. checklib; - v = calllib(tr.lib, 'trans_viscosity', tr.tr_id); + v = calllib(ct, 'trans_viscosity', tr.tr_id); if v == -1.0 error(geterr); elseif v < 0.0 @@ -63,14 +61,14 @@ classdef Transport < handle end end - function v = thermaoConductivity(tr) + function v = thermalConductivity(tr) % Get the thermal conductivity. % - % :return: + % return: % Double thermal conductivity. Unit: W/m-K. checklib; - v = calllib(tr.lib, 'trans_thermalConductivity', tr.tr_id); + v = calllib(ct, 'trans_thermalConductivity', tr.tr_id); if v == -1.0 error(geterr); elseif v < 0.0 @@ -81,11 +79,11 @@ classdef Transport < handle function v = electricalConductivity(tr) % Get the electrical conductivity. % - % :return: + % return: % Double electrical conductivity. Unit: S/m. checklib; - v = calllib(tr.lib, 'trans_electricalConductivity', tr.tr_id); + v = calllib(ct, 'trans_electricalConductivity', tr.tr_id); if v == -1.0 error(geterr); elseif v < 0.0 @@ -96,7 +94,7 @@ classdef Transport < handle function v = mixDiffCoeffs(tr) % Get the mixture-averaged diffusion coefficients. % - % :return: + % return: % Vector of length nSpecies with the mixture-averaged % diffusion coefficients. Unit: m^2/s. @@ -104,14 +102,14 @@ classdef Transport < handle nsp = tr.th.nSpecies; xx = zeros(1, nsp); pt = libpointer('doublePtr', xx); - calllib(tr.lib, 'trans_getMixDiffCoeffs', tr.tr_id, nsp, pt); + calllib(ct, 'trans_getMixDiffCoeffs', tr.tr_id, nsp, pt); v = pt.Value; end function v = thermalDiffCoeffs(tr) % Get the thermal diffusion coefficients. % - % :return: + % return: % Vector of length nSpecies with the thermal diffusion % coefficients. @@ -119,14 +117,14 @@ classdef Transport < handle nsp = tr.th.nSpecies; xx = zeros(1, nsp); pt = libpointer('doublePtr', xx); - calllib(tr.lib, 'trans_getThermalDiffCoeffs', tr.tr_id, nsp, pt); + calllib(ct, 'trans_getThermalDiffCoeffs', tr.tr_id, nsp, pt); v = pt.Value; end function v = binDiffCoeffs(tr) % Get the binary diffusion coefficients. % - % :return: + % return: % A matrix of binary diffusion coefficients. The matrix is % symmetric: d(i, j) = d(j, i). Unit: m^2/s. @@ -134,14 +132,14 @@ classdef Transport < handle nsp = tr.th.nSpecies; xx = zeros(1, nsp); pt = libpointer('doublePtr', xx); - calllib(tr.lib, 'trans_getBinDiffCoeffs', tr.tr_id, nsp, pt); + calllib(ct, 'trans_getBinDiffCoeffs', tr.tr_id, nsp, pt); v = pt.Value; end function v = multiDiffCoeffs(tr) % Get the multicomponent diffusion coefficients. % - % :return: + % return: % Vector of length nSpecies with the multicomponent % diffusion coefficients. Unit: m^2/s. @@ -149,29 +147,29 @@ classdef Transport < handle nsp = tr.th.nSpecies; xx = zeros(1, nsp); pt = libpointer('doublePtr', xx); - calllib(tr.lib, 'trans_getMultiDiffCoeffs', tr.tr_id, nsp, pt); + calllib(ct, 'trans_getMultiDiffCoeffs', tr.tr_id, nsp, pt); v = pt.Value; end function setParameters(tr, type, k, p) % Set the parameters. % - % :param type: - % :param k: - % :param p: + % parameter type: + % parameter k: + % parameter p: checklib; - calllib(tr.lib, 'trans_setParameters', tr.tr_id, type, k, p); + calllib(ct, 'trans_setParameters', tr.tr_id, type, k, p); end function setThermalConductivity(tr, lam) % Set the thermal conductivity. % - % :param lam: + % parameter lam: % Thermal conductivity in W/(m-K). checklib; - setParameters(tr, 1, 0, lam); + tr.setParameters(1, 0, lam); end end diff --git a/interfaces/matlab_experimental/Base/Water.m b/interfaces/matlab_experimental/Base/Water.m index 018ed277b..96e07a521 100644 --- a/interfaces/matlab_experimental/Base/Water.m +++ b/interfaces/matlab_experimental/Base/Water.m @@ -1,4 +1,5 @@ function w = Water() % Return an object representing water. - h = Solution('liquidvapor.yaml', 'water'); + w = Solution('liquidvapor.yaml', 'water'); end + diff --git a/interfaces/matlab_experimental/Constants/faradayconstant.m b/interfaces/matlab_experimental/Constants/faradayconstant.m new file mode 100644 index 000000000..dd8cafd63 --- /dev/null +++ b/interfaces/matlab_experimental/Constants/faradayconstant.m @@ -0,0 +1,4 @@ +function r = faradayconstant + % Get the Faraday constant in C/kmol of electron. + r = 96485336.4595687; +end diff --git a/interfaces/matlab_experimental/Constants/gasconstant.m b/interfaces/matlab_experimental/Constants/gasconstant.m new file mode 100644 index 000000000..e610d2b10 --- /dev/null +++ b/interfaces/matlab_experimental/Constants/gasconstant.m @@ -0,0 +1,4 @@ +function r = gasconstant + % Universal gas constant in J/kmol-K. + r = 8314.4621; +end diff --git a/interfaces/matlab_experimental/Constants/oneatm.m b/interfaces/matlab_experimental/Constants/oneatm.m new file mode 100644 index 000000000..e8b9eec9b --- /dev/null +++ b/interfaces/matlab_experimental/Constants/oneatm.m @@ -0,0 +1,4 @@ +function p = oneatm + % Get 1 atm in Pa. + p = 101325.0; +end diff --git a/interfaces/matlab_experimental/Func/Func.m b/interfaces/matlab_experimental/Func/Func.m index 710698bad..7282a8b95 100644 --- a/interfaces/matlab_experimental/Func/Func.m +++ b/interfaces/matlab_experimental/Func/Func.m @@ -8,10 +8,6 @@ classdef Func < handle typ end - properties(Constant = true) - lib = 'cantera_shared' - end - methods %% Functor class constructor @@ -44,7 +40,7 @@ classdef Func < handle % class 'Func1'. See the Cantera C++ documentation for more % details. % - % :param typ: + % parameter typ: % String indicating type of functor to create. Possible % values are: % * 'polynomial' @@ -57,11 +53,11 @@ classdef Func < handle % * 'composite' % * 'periodic' % - % :param n: + % parameter n: % Number of parameters required for the functor - % :param p: + % parameter p: % Vector of parameters - % :return: + % return: % Instance of class :mat:func:`Func` checklib; @@ -81,13 +77,13 @@ classdef Func < handle ptr = libpointer('doublePtr', p); [m, n] = size(p); lenp = m * n; - nn = calllib(x.lib, 'func_new', type, n, lenp, ptr); + nn = calllib(ct, 'func_new', type, n, lenp, ptr); elseif itype < 45 m = n; - nn = calllib(x.lib, 'func_new', type, n, m, 0); + nn = calllib(ct, 'func_new', type, n, m, 0); else ptr = libpointer('doublePtr', p); - nn = calllib(x.lib, 'func_new', type, n, 0, ptr); + nn = calllib(ct, 'func_new', type, n, 0, ptr); end end @@ -124,13 +120,13 @@ classdef Func < handle function clear(f) % Clear the functor from memory. checklib; - calllib(f.lib, 'func_del', f.id); + calllib(ct, 'func_del', f.id); end function display(a) % Display the equation of the input function on the terminal. % - % :param a: + % parameter a: % Instance of class 'Func' disp(' '); @@ -162,11 +158,11 @@ classdef Func < handle function b = subsref(a, s) % Redefine subscripted references for functors. % - % :param a: + % parameter a: % Instance of class 'Func' - % :param s: + % parameter s: % Value at which the function should be evaluated. - % :return: + % return: % The value of the function evaluated at 's'. checklib; @@ -174,7 +170,7 @@ classdef Func < handle ind= s.subs{:}; b = zeros(1, length(ind)); for k = 1:length(ind) - b(k) = calllib(a.lib, 'func_value', a.id, ind(k)); + b(k) = calllib(ct, 'func_value', a.id, ind(k)); end else error('Specify value for x as p(x)'); end diff --git a/interfaces/matlab_experimental/Func/gaussian.m b/interfaces/matlab_experimental/Func/gaussian.m new file mode 100644 index 000000000..d2b86870a --- /dev/null +++ b/interfaces/matlab_experimental/Func/gaussian.m @@ -0,0 +1,12 @@ +function g = gaussian(peak, center, width) + % Create a Gaussian functor instance. + % :param peak: + % The peak value. + % :param center: + % Value of x at which the peak is located. + % :param width: + % Full width at half-maximum. The value of the function at center + % +/- (width)/2 is one-half the peak value. + + g = Func('gaussian', 0, [peak, center, width]); +end diff --git a/interfaces/matlab_experimental/Func/polynom.m b/interfaces/matlab_experimental/Func/polynom.m new file mode 100644 index 000000000..241d6b847 --- /dev/null +++ b/interfaces/matlab_experimental/Func/polynom.m @@ -0,0 +1,17 @@ +function poly = polynom(coeffs) + % Create a polynomial functor instance. + % The polynomial coefficients are specified by a vector [a0, a1, ... + % an]. Examples: + % + % polynom([-2, 6, 3]) % 3x^2+6x-2 + % polynom([1.0, -2.5, 0, 0, 2]) % 2x^4-2.5x+1 + + [n m] = size(coeffs); + if n == 1 + poly = Func('polynomial', m - 1, coeffs); + elseif m == 1 + poly = Func('polynomial', n - 1, coeffs); + else + error('wrong shape for coefficient array'); + end +end diff --git a/interfaces/matlab_experimental/LoadCantera.m b/interfaces/matlab_experimental/LoadCantera.m index f0decf255..1cccb04c9 100644 --- a/interfaces/matlab_experimental/LoadCantera.m +++ b/interfaces/matlab_experimental/LoadCantera.m @@ -1,8 +1,20 @@ -function Load_Cantera - addpath('Class', 'Utility', 'PresetObjects', 'Examples'); - if ~libisloaded('cantera_shared') +function LoadCantera + addpath('Class', 'Utility', 'PresetMixtures', 'PresetReactors', ... + 'PresetObjects', '1D', 'Examples'); + if ispc + ctname = 'cantera_shared.dll'; + elseif ismac + ctname = 'libcantera_shared.dylib'; + elseif isunix + ctname = 'libcantera_shared.so'; + + else + error('Operating System Not Supported!'); + return; + end + if ~libisloaded(ct) load('Utility/cantera_root.mat'); - [~,warnings] = loadlibrary([cantera_root '/lib/cantera_shared.dll'], ... + [~,warnings] = loadlibrary([cantera_root '/lib/' ctname], ... [cantera_root '/include/cantera/clib/ctmatlab.h'], ... 'includepath', [cantera_root '/include'], ... 'addheader','ct','addheader','ctfunc', ... diff --git a/interfaces/matlab_experimental/Reactor/ConstPressureReactor.m b/interfaces/matlab_experimental/Reactor/ConstPressureReactor.m new file mode 100644 index 000000000..d63d3e037 --- /dev/null +++ b/interfaces/matlab_experimental/Reactor/ConstPressureReactor.m @@ -0,0 +1,14 @@ +function r = ConstPressureReactor(contents) + % Create a constant pressure reactor object. + % Pressure is held constant. The volume is not a state variable, but + % instead takes on whatever value is consistent with holding the + % pressure constant. + % + %: param contents: + % Contents of the reactor of class 'Solution'. + + if nargin == 0 + contents = 0; + end + r = Reactor(contents, 'ConstPressureReactor'); +end diff --git a/interfaces/matlab_experimental/Reactor/FlowDevice.m b/interfaces/matlab_experimental/Reactor/FlowDevice.m index 52239d00a..c37fa8441 100644 --- a/interfaces/matlab_experimental/Reactor/FlowDevice.m +++ b/interfaces/matlab_experimental/Reactor/FlowDevice.m @@ -7,17 +7,13 @@ classdef FlowDevice < handle downstream end - properties(Constant = true) - lib = 'cantera_shared' - end - methods %% FlowDevice class constructor function x = FlowDevice(typ) % Flow Device class constructor. % - % :param typ: + % parameter typ: % Type of flow device to be created. Type = % 'MassFlowController', 'PressureController' or 'Valve'. @@ -36,7 +32,7 @@ classdef FlowDevice < handle end x.type = typ; - x.id = calllib(x.lib, 'flowdev_new2', typ); + x.id = calllib(ct, 'flowdev_new', typ); % if x.id < 0 % error(geterr); % end @@ -49,7 +45,7 @@ classdef FlowDevice < handle function clear(f) % Clear the specified flow device from memory. checklib; - calllib(f.lib, 'flowdev_del', f.id); + calllib(ct, 'flowdev_del', f.id); end %% FlowDevice methods @@ -57,9 +53,9 @@ classdef FlowDevice < handle function install(f, upstream, downstream) % Install a flow device between reactors or reservoirs. % - % :param upstream: + % parameter upstream: % Upsteram 'Reactor' or 'Reservoir'. - % :param downstream: + % parameter downstream: % Downstream 'Reactor' or 'Reservoir'. checklib; @@ -70,7 +66,7 @@ classdef FlowDevice < handle end i = upstream.id; j = downstream.id; - ok = calllib(f.lib, 'flowdev_install', f.id, i, j); + ok = calllib(ct, 'flowdev_install', f.id, i, j); % if ok < 0 % error(geterr) % end @@ -81,31 +77,31 @@ classdef FlowDevice < handle function mdot = massFlowRate(f, time) % Get the mass flow rate at a given time. % - % :param time: + % parameter time: % Time at which the mass flow rate is desired. - % :return: + % 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); + 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(f.lib, 'flowdev_massFlowRate', f.id, time); + mdot = calllib(ct, 'flowdev_massFlowRate', f.id, time); end end function setFunction(f, mf) % Set the time function with class 'func'. % - % :param mf: + % parameter mf: % Instance of class 'func'. checklib; if strcmp(f.type, 'MassFlowController') - k = calllib(f.lib, 'flowdev_setTimeFunction', f.id, ... + k = calllib(ct, 'flowdev_setTimeFunction', f.id, ... mf.id); % if k < 0 % error(geterr); @@ -118,12 +114,12 @@ classdef FlowDevice < handle function setMassFlowRate(f, mdot) % Set the mass flow rate to a constant value. % - % :param mdot: + % parameter mdot: % Mass flow rate checklib; if strcmp(f.type, 'MassFlowController') - k = calllib(f.lib, 'flowdev_setMassFlowCoeff', f.id, mdot); + k = calllib(ct, 'flowdev_setMassFlowCoeff', f.id, mdot); % if k < 0 % error(geterr); % end @@ -140,14 +136,14 @@ classdef FlowDevice < handle % as long as this produces a positive value. If thsi expression % is negative, zero is returned. % - % :param k: + % parameter 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.id, k); + ok = calllib(ct, 'flowdev_setValveCoeff', f.id, k); % if k < 0 % error(geterr); % end diff --git a/interfaces/matlab_experimental/Reactor/FlowReactor.m b/interfaces/matlab_experimental/Reactor/FlowReactor.m new file mode 100644 index 000000000..1ca986dfa --- /dev/null +++ b/interfaces/matlab_experimental/Reactor/FlowReactor.m @@ -0,0 +1,12 @@ +function r = FlowReactor(contents) + % Create a flow reactor object. + % A reactor representing adiabatic plug flow in a constant-area duct. + % + %: param contents: + % Contents of the reactor of class 'Solution'. + + if nargin == 0 + contents = 0; + end + r = Reactor(contents, 'FlowReactor'); +end diff --git a/interfaces/matlab_experimental/Reactor/IdealGasConstPressureReactor.m b/interfaces/matlab_experimental/Reactor/IdealGasConstPressureReactor.m new file mode 100644 index 000000000..55db37e74 --- /dev/null +++ b/interfaces/matlab_experimental/Reactor/IdealGasConstPressureReactor.m @@ -0,0 +1,14 @@ +function r = IdealGasConstPressureReactor(contents) + % Create a constant pressure reactor object with an ideal gas. + % Pressure is held constant. The volume is not a state variable, but + % instead takes on whatever value is consistent with holding the + % pressure constant. + % + %: param contents: + % Contents of the reactor of class 'Solution'. + + if nargin == 0 + contents = 0; + end + r = Reactor(contents, 'IdealGasConstPressureReactor'); +end diff --git a/interfaces/matlab_experimental/Reactor/IdealGasReactor.m b/interfaces/matlab_experimental/Reactor/IdealGasReactor.m new file mode 100644 index 000000000..6a0c6ec67 --- /dev/null +++ b/interfaces/matlab_experimental/Reactor/IdealGasReactor.m @@ -0,0 +1,13 @@ +function r = IdealGasReactor(contents) + % Create a ideal gas reactor object. + % The governing equations are specialized for the ideal gas equation of + % state (and do not work correctly with other thermodynamic models). + % + %: param contents: + % Contents of the reactor of class 'Solution'. + + if nargin == 0 + contents = 0; + end + r = Reactor(contents, 'IdealGasReactor'); +end diff --git a/interfaces/matlab_experimental/Reactor/MassFlowController.m b/interfaces/matlab_experimental/Reactor/MassFlowController.m new file mode 100644 index 000000000..99eb00776 --- /dev/null +++ b/interfaces/matlab_experimental/Reactor/MassFlowController.m @@ -0,0 +1,21 @@ +function m = MassFlowController(upstream, downstream) + % Create a mass flow controller. + % Creates an instance of class 'FlowDevice' configured to simulate a + % mass flow controller that maintains a constant mass flow rate + % independent of upstream or downstream conditions. If two reactor + % objects are supplied as arguments, the controller is installed + % between the two reactors. Otherwise, the 'install' method should be + % used to install the mass flow controller between reactors. + % + % :param upstream: + % Upstream 'Reactor' or 'Reservoir'. + % :param downstream: + % Downstream 'Reactor' or 'Reservoir. + % :return: + % Instance of class 'FlowDevice'. + + m = FlowDevice('MassFlowController'); + if nargin == 2 + m.install(upstream, downstream) + end +end diff --git a/interfaces/matlab_experimental/Reactor/Reactor.m b/interfaces/matlab_experimental/Reactor/Reactor.m index bb616210a..b34108623 100644 --- a/interfaces/matlab_experimental/Reactor/Reactor.m +++ b/interfaces/matlab_experimental/Reactor/Reactor.m @@ -17,10 +17,6 @@ classdef Reactor < handle EnergyFlag end - properties(Constant = true) - lib = 'cantera_shared' - end - methods %% Reactor class constructor @@ -30,10 +26,10 @@ classdef Reactor < handle % reactors through flow lines or through walls that may expand % or contract and/orconduct heat. % - % :param contents: + % parameter contents: % Instance of class 'Solution' representing the contents of % the reactor. - % :param typ: + % parameter typ: % Character array of reactor type. Options are: % 'Reservoir' % 'Reactor' @@ -41,7 +37,7 @@ classdef Reactor < handle % 'ConstPressureReactor' % 'IdealGasReactor' % 'IdealGasConstPressureReactor' - % :return: + % return: % Instance of class 'Reactor'. checklib; @@ -56,7 +52,7 @@ classdef Reactor < handle end r.type = char(typ); - r.id = calllib(r.lib, 'reactor_new2', typ); + r.id = calllib(ct, 'reactor_new', typ); % if r.id < 0 % error(geterr); @@ -75,15 +71,15 @@ classdef Reactor < handle function clear(r) % Clear the reactor from memory. checklib; - calllib(r.lib, 'reactor_del', r.id); + calllib(ct, 'reactor_del', r.id); end function insert(r, gas) % Insert a solution or mixture into a reactor. % - % :param r: + % parameter r: % Instance of class 'Reactor'. - % :param gas: + % parameter gas: % Instance of class 'Solution'. r.contents = gas; @@ -99,9 +95,9 @@ classdef Reactor < handle % This method is used internally during Reactor initialization, % but is usually not called by users. % - % :param r: + % parameter r: % Instance of class 'Reactor'. - % :param t: + % parameter t: % Instance of class 'ThermoPhase' or another object % containing an instance of that class. checklib; @@ -110,7 +106,7 @@ classdef Reactor < handle error('Wrong object type'); end - calllib(r.lib, 'reactor_setThermoMgr', r.id, t.tp_id); + calllib(ct, 'reactor_setThermoMgr', r.id, t.tp_id); end function setKineticsMgr(r, k) @@ -119,9 +115,9 @@ classdef Reactor < handle % This method is used internally during Reactor initialization, % but is usually not called by users. % - % :param r: + % parameter r: % Instance of class 'Reactor'. - % :param t: + % parameter t: % Instance of class 'Kinetics' or another object % containing an instance of that class. checklib; @@ -130,7 +126,7 @@ classdef Reactor < handle error('Wrong object type'); end - calllib(r.lib, 'reactor_setKineticsMgr', r.id, k.kin_id); + calllib(ct, 'reactor_setKineticsMgr', r.id, k.kin_id); end %% Reactor get methods @@ -138,103 +134,103 @@ classdef Reactor < handle function temperature = get.T(r) % Get the temperature of the reactor. % - % :return: + % return: % The temperature of the reactor contents at the end of the % last call to 'advance' or 'step'. Unit: K. checklib; - temperature = calllib(r.lib, 'reactor_temperature', r.id); + temperature = calllib(ct, 'reactor_temperature', r.id); end function pressure = get.P(r) % Get the pressure of the reactor. % - % :return: + % return: % The pressure of the reactor contents at the end of the % last call to 'advance' or 'step'. Unit: Pa. checklib; - pressure = calllib(r.lib, 'reactor_pressure', r.id); + pressure = calllib(ct, 'reactor_pressure', r.id); end function rho = get.D(r) % Get the density of the reactor. % - % :return: + % return: % Density of the phase in the input. Unit: kg/m^3. checklib; - rho = calllib(r.lib, 'reactor_density', r.id); + rho = calllib(ct, 'reactor_density', r.id); end function mass = get.M(r) % Get the mass of the reactor. % - % :return: + % return: % The mass of the reactor contents at the end of the % last call to 'advance' or 'step'. The mass is retrieved % from the solution vector. Unit: kg. checklib; - mass = calllib(r.lib, 'reactor_mass', r.id); + mass = calllib(ct, 'reactor_mass', r.id); end function volume = get.V(r) % Get the volume of the reactor. % - % :return: + % return: % The volume of the reactor contents at the end of the % last call to 'advance' or 'step'. Unit: m^3. checklib; - volume = calllib(r.lib, 'reactor_volume', r.id); + volume = calllib(ct, 'reactor_volume', r.id); end function enthalpy_mass = get.H(r) % Get the mass specific enthalpy of the reactor. % - % :return: + % return: % The mass specific enthalpy of the reactor contents at the % end of the last call to 'advance' or 'step'. The enthalpy % is retrieved from the solution vector. Unit: J/kg. checklib; - enthalpy_mass = calllib(r.lib, 'reactor_enthalpy_mass', r.id); + enthalpy_mass = calllib(ct, 'reactor_enthalpy_mass', r.id); end function intEnergy_mass = get.U(r) % Get the mass specific internal energy of the reactor. % - % :return: + % return: % The mass specific internal energy of the reactor contents % at the end of the last call to 'advance' or 'step'. The % internal energy is retrieved from the solution vector. % Unit: J/kg. checklib; - intEnergy_mass = calllib(r.lib, 'reactor_intEnergy_mass', r.id); + intEnergy_mass = calllib(ct, 'reactor_intEnergy_mass', r.id); end function yi = massFraction(r, species) % Get the mass fraction of a species. % - % :param species: + % parameter species: % String or one-based integer id of the species. checklib; if ischar(species) - k = r.contents.thermo.speciesIndex(species) - 1; + k = r.contents.speciesIndex(species) - 1; else k = species - 1; end - yi = calllib(r.lib, 'reactor_massFraction', r.id, k); + yi = calllib(ct, 'reactor_massFraction', r.id, k); end function massFractions = get.Y(r) % Get the mass fractions of the reactor. % - % :return: + % return: % The mass fractions of the reactor contents at the end of % the last call to 'advance' or 'step'. @@ -252,21 +248,21 @@ classdef Reactor < handle function setInitialVolume(r, v0) % Set the initial reactor volume. % - % :param v0: + % parameter v0: % Initial volume in m^3. checklib; - calllib(r.lib, 'reactor_setInitialVolume', r.id, v0); + calllib(ct, 'reactor_setInitialVolume', r.id, v0); end function r = set.Mdot(r, MFR) % Set the mass flow rate. % - % :param MFR: + % parameter MFR: % Mass flow rate in kg/s. checklib; - calllib(r.lib, 'reactor_setMassFlowRate', r.id, MFR); + calllib(ct, 'reactor_setMassFlowRate', r.id, MFR); r.Mdot = MFR; end @@ -281,9 +277,9 @@ classdef Reactor < handle % equations enabled if there are reactions present in the % mechanism file, and disabled otherwise. % - % :param r: + % parameter r: % Instance of class 'Reactor'. - % :param flag: + % parameter flag: % String, either "on" or "off" to enable or disable chemical % reactions, respectively. @@ -296,7 +292,7 @@ classdef Reactor < handle else error('Input must be "on" or "off"'); end - calllib(r.lib, 'reactor_setChemistry', r.id, iflag); + calllib(ct, 'reactor_setChemistry', r.id, iflag); r.ChemistryFlag = flag; end @@ -310,9 +306,9 @@ classdef Reactor < handle % By default, Reactor objects are created with the energy % equation enabled, so usually this method % - % :param r: + % parameter r: % Instance of class 'Reactor'. - % :param flag: + % parameter flag: % String, either "on" or "off" to enable or disable chemical % reactions, respectively. @@ -326,7 +322,7 @@ classdef Reactor < handle end if iflag >= 0 - calllib(r.lib, 'reactor_setEnergy', r.id, iflag); + calllib(ct, 'reactor_setEnergy', r.id, iflag); else error('Input must be "on" or "off".'); end diff --git a/interfaces/matlab_experimental/Reactor/ReactorNet.m b/interfaces/matlab_experimental/Reactor/ReactorNet.m index e75f74e09..0e4a89f8c 100644 --- a/interfaces/matlab_experimental/Reactor/ReactorNet.m +++ b/interfaces/matlab_experimental/Reactor/ReactorNet.m @@ -8,10 +8,6 @@ classdef ReactorNet < handle rtol end - properties(Constant = true) - lib = 'cantera_shared' - end - methods %% ReactorNet class constructor @@ -21,10 +17,10 @@ classdef ReactorNet < handle % to simultaneously advance the state of one or more coupled % reactors. % - % :param reactors: + % parameter reactors: % Instance of class 'Reactor' or a cell array of instance of % 'Reactor'. - % :return: + % return: % Instance of class 'ReactorNet'. checklib; @@ -39,7 +35,7 @@ classdef ReactorNet < handle reactors = {reactor}; end - r.id = calllib(r.lib, 'reactornet_new'); + r.id = calllib(ct, 'reactornet_new'); % if r.id < 0 % error(geterr); % end @@ -56,19 +52,19 @@ classdef ReactorNet < handle function clear(r) % Clear the ReactorNet object from the memory. checklib; - calllib(r.lib, 'reactornet_del', r.id); + calllib(ct, 'reactornet_del', r.id); end function addReactor(net, reactor) % Add a reactor to a network. % - % :param net: + % parameter net: % Instance of class 'ReactorNet'. - % :param reactor: + % parameter reactor: % Instance of class 'Solution'. checklib; - calllib(net.lib, 'reactornet_addreactor', net.id, reactor.id); + calllib(ct, 'reactornet_addreactor', net.id, reactor.id); end function advance(r, tout) @@ -81,11 +77,11 @@ classdef ReactorNet < handle % an absolute time, not a time interval.) The integrator may % take many internal timesteps before reaching tout. % - % :param tout: + % parameter tout: % End time of the integration. Unit: s. checklib; - calllib(r.lib, 'reactornet_advance', r.id, tout); + calllib(ct, 'reactornet_advance', r.id, tout); end %% ReactorNet set methods @@ -93,12 +89,12 @@ classdef ReactorNet < handle function setInitialTime(r, t) % Set the initial time of the integration. % - % :param t: + % parameter t: % Time at which integration should be restarted, using the % current state as the initial condition. Unit: s. checklib; - calllib(r.lib, 'reactornet_setInitialTime', r.id, t); + calllib(ct, 'reactornet_setInitialTime', r.id, t); end function setMaxTimeStep(r, maxstep) @@ -113,19 +109,19 @@ classdef ReactorNet < handle % upper bound on the timestep. checklib; - calllib(r.lib, 'reactornet_setMaxTimeStep', r.id, maxstep); + calllib(ct, 'reactornet_setMaxTimeStep', r.id, maxstep); end function setTolerances(r, rerr, aerr) % Set the error tolerance. % - % :param rtol: + % parameter rtol: % Scalar relative error tolerance. - % :param atol: + % parameter atol: % Scalar absolute error tolerance. checklib; - calllib(r.lib, 'reactornet_setTolerances', r.id, rerr, aerr); + calllib(ct, 'reactornet_setTolerances', r.id, rerr, aerr); end %% ReactorNet get methods @@ -139,25 +135,25 @@ classdef ReactorNet < handle % rapidly changing, the time step becomes smaller to resolve % the solution. checklib; - t = calllib(r.lib, 'reactor_step', r.id); + t = calllib(ct, 'reactor_step', r.id); end function t = get.time(r) % Get the current time in s. checklib; - t = calllib(r.lib, 'reactornet_time', r.id); + t = calllib(ct, 'reactornet_time', r.id); end function t = get.rtol(r) % Get the relative error tolerance checklib; - t = calllib(r.lib, 'reactornet_rtol', r.id); + t = calllib(ct, 'reactornet_rtol', r.id); end function t = get.atol(r) % Get the absolute error tolerance checklib; - t = calllib(r.lib, 'reactornet_atol', r.id); + t = calllib(ct, 'reactornet_atol', r.id); end end diff --git a/interfaces/matlab_experimental/Reactor/ReactorSurface.m b/interfaces/matlab_experimental/Reactor/ReactorSurface.m index 1cdda9993..b5044bc18 100644 --- a/interfaces/matlab_experimental/Reactor/ReactorSurface.m +++ b/interfaces/matlab_experimental/Reactor/ReactorSurface.m @@ -6,10 +6,6 @@ classdef ReactorSurface < handle reactor end - properties(Constant = true) - lib = 'cantera_shared' - end - methods %% ReactorSurface class constructor @@ -25,29 +21,29 @@ classdef ReactorSurface < handle % after initial construction by using the various methods of % the 'ReactorSurface' class. % - % :param kleft: + % parameter kleft: % Surface reaction mechanisms for the left-facing surface. % This must bean instance of class 'Kinetics', or of a class % derived from 'Kinetics', such as 'Interface'. - % :param reactor: + % parameter reactor: % Instance of class 'Reactor' to be used as the adjacent % bulk phase. - % :param area: + % parameter area: % The area of the surface in m^2. Defaults to 1.0 m^2 if not % specified. - % :return: + % return: % Instance of class 'ReactorSurface'. checklib; - s.id = calllib(s.lib, 'reactorsurface_new', 0); + s.id = calllib(ct, 'reactorsurface_new', 0); s.reactor = -1; % if r.id < 0 % error(geterr); % end if nargin >= 1 - s.setKinetics(s, kleft); + s.setKinetics(kleft); end if nargin >= 2 @@ -60,7 +56,7 @@ classdef ReactorSurface < handle if nargin >= 3 if isnumeric(area) - s.setArea(area); + s.area = area; else warning('Area was not a number and was not set'); end @@ -73,14 +69,14 @@ classdef ReactorSurface < handle function clear(s) % Clear the ReactorSurface object from the memory. checklib; - calllib(s.lib, 'reactorsurface_del', s.id); + calllib(ct, 'reactorsurface_del', s.id); end function install(s, r) % Install a ReactorSurface in a Reactor. checklib; s.reactor = r; - calllib(s.lib, 'reactorsurface_install', s.id, r.id); + calllib(ct, 'reactorsurface_install', s.id, r.id); end %% ReactorSurface get methods @@ -88,7 +84,7 @@ classdef ReactorSurface < handle function a = get.area(s) % Get the areaof the reactor surface in m^2. checklib; - a = calllib(s.lib, 'reactorsurface_area', s.id); + a = calllib(ct, 'reactorsurface_area', s.id); end %% ReactorSurface set methods @@ -96,12 +92,12 @@ classdef ReactorSurface < handle function s = set.area(s, a) % Set the area of a reactor surface checklib; - calllib(s.lib, 'reactorsurface_setArea', s.id, a); + calllib(ct, 'reactorsurface_setArea', s.id, a); end function setKinetics(s, kin) % Setthe surface reaction mechanism on a reactor surface. - % :param kin: + % parameter kin: % Instance of class 'Kinetics' (or another object derived % from kin) to be used as the kinetic mechanism for this % surface. Typically an instance of class 'Interface'. @@ -113,7 +109,7 @@ classdef ReactorSurface < handle ikin = kin.id; end - calllib(s.lib, 'reactorsurface_setkinetics', s.id, ikin); + calllib(ct, 'reactorsurface_setkinetics', s.id, ikin); end end diff --git a/interfaces/matlab_experimental/Reactor/Reservoir.m b/interfaces/matlab_experimental/Reactor/Reservoir.m new file mode 100644 index 000000000..470a6419c --- /dev/null +++ b/interfaces/matlab_experimental/Reactor/Reservoir.m @@ -0,0 +1,18 @@ +function r = Reservoir(contents) + % Create a reservoir object. + % An instance of class 'Reactor' configured so that its intensive state + % is constant in time. A reservoir may be thought of as inifinite in + % extent, perfectly mixed, and non-reacting, so that fluid may be + % extracted or added without changing the composition or thermodynamic + % state. Note that even if the reaction mechanism associated with the + % fluid in the reactor defines reactions, they are disabled within + % reservoirs. + % + %: param contents: + % Contents of the reactor of class 'Solution'. + + if nargin == 0 + contents = 0; + end + r = Reactor(contents, 'Reservoir'); +end diff --git a/interfaces/matlab_experimental/Reactor/Valve.m b/interfaces/matlab_experimental/Reactor/Valve.m new file mode 100644 index 000000000..bc51133f6 --- /dev/null +++ b/interfaces/matlab_experimental/Reactor/Valve.m @@ -0,0 +1,32 @@ +function v = Valve(upstream, downstream) + % Create a valve. + % Creates an instance of class 'FlowDevice' configured to simulate a + % valve that produces a flow rate proportional to the pressure + % difference between the upstream or downstream reactors. + % + % The mass flow rate [kg/s] is computed from the expression: + % + % mdot = K(P_upstream - P_downstream) + % + % as long as this produces a positive value. If this expression is + % negative, zero is returned. Therefore, the 'Valve' object acts as a + % check valve - flow is always from upstream to downstream. + % + % Note: as currently implemented, the valve object does not model real + % valve characteristics - inparticular, it does not model choked flow. + % The mass flow rate is always assumed to be linearly proportional to + % the pressure difference, no matter how large. This MAY change in a + % future release. + % + % :param upstream: + % Upstream 'Reactor' or 'Reservoir'. + % :param downstream: + % Downstream 'Reactor' or 'Reservoir. + % :return: + % Instance of class 'FlowDevice'. + + v = FlowDevice('Valve'); + if nargin == 2 + v.install(upstream, downstream) + end +end diff --git a/interfaces/matlab_experimental/Reactor/Wall.m b/interfaces/matlab_experimental/Reactor/Wall.m index 510800e32..0bcc65401 100644 --- a/interfaces/matlab_experimental/Reactor/Wall.m +++ b/interfaces/matlab_experimental/Reactor/Wall.m @@ -8,10 +8,6 @@ classdef Wall < handle area end - properties(Constant = true) - lib = 'cantera_shared' - end - methods %% Wall class constructor @@ -50,25 +46,25 @@ classdef Wall < handle % wall can be set by using empty strings or 0.0 for each of the % arguments before the velocity with no harm. % - % :param left: + % parameter left: % Instance of class 'Reactor' to be used as the bulk phase % on the left side of the wall. - % :param right: + % parameter right: % Instance of class 'Reactor' to be used as the bulk phase % on the right side of the wall. - % :param area: + % parameter area: % The area of the wall in m^2. Defaults to 1.0 m^2 if not % specified. - % :param k: + % parameter k: % Expansion rate coefficient in m/(s-Pa). Defaults to 0.0 if % not specified. - % :param u: + % parameter u: % Heat transfer coefficient in W/(m^2-K). Defaults to 0.0 if % not specified. - % :param q: + % parameter q: % Heat flux in W/m^2. Defaults to 0.0 if not specified. Must % be an instance of 'Func'. - % :param v: + % parameter v: % Velocity of the wall in m/s. Defaults to 0.0 if not % specified. Must be an instance of 'Func'. @@ -78,7 +74,7 @@ classdef Wall < handle typ = 'Wall'; x.type = char(typ); - x.id = calllib(x.lib, 'wall_new2', x.type); + x.id = calllib(ct, 'wall_new', x.type); % if x.index < 0 % error(geterr); % end @@ -142,7 +138,7 @@ classdef Wall < handle function clear(w) % Clear the Wall object from the memory. checklib; - calllib(w.lib, 'wall_del', w.id); + calllib(ct, 'wall_del', w.id); end function install(w, l, r) @@ -151,51 +147,51 @@ classdef Wall < handle checklib; w.left = l; w.right = r; - calllib(w.lib, 'wall_install', w.id, l.id, r.id); + calllib(ct, 'wall_install', w.id, l.id, r.id); end function ok = ready(w) % Check whether a wall is ready. checklib; - ok = calllib(w.lib, 'wall_ready', w.id); + ok = calllib(ct, 'wall_ready', w.id); end %% ReactorNet set methods - function w = set.area(w, a) + function set.area(w, a) % Set the area of a wall. checklib; - calllib(w.lib, 'wall_Area', w.id, a); + calllib(ct, 'wall_setArea', w.id, a); end function setThermalResistance(w, r) % Set the thermal resistance. % - % :param r: + % parameter r: % Thermal resistance. Unit: K*m^2/W. checklib; - calllib(w.lib, 'wall_setThermalResistance', w.id, r); + calllib(ct, 'wall_setThermalResistance', w.id, r); end function setHeatTransferCoeff(w, u) % Set the thermal transfer coefficient. % - % :param u: + % parameter u: % Heat transfer coefficient. Unit: W/(m^2-K). checklib; - calllib(w.lib, 'wall_setHeatTransferCoeff', w.id, u); + calllib(ct, 'wall_setHeatTransferCoeff', w.id, u); end function setExpansionRateCoeff(w, k) % Set the expansion rate coefficient. % - % :param k: + % parameter k: % Expanstion rate coefficient. Unit: m/(s-Pa). checklib; - calllib(w.lib, 'wall_setExpansionRateCoeff', w.id, k); + calllib(ct, 'wall_setExpansionRateCoeff', w.id, k); end function setHeatFlux(w, f) @@ -206,11 +202,11 @@ classdef Wall < handle % to specify a constant heat flux by using the polynomial % functor with only the first term specified. % - % :param f: + % parameter f: % Instance of class 'Func'. Unit: W/m^2. checklib; - calllib(w.lib, 'wall_setHeatFlux', w.id, f.id); + calllib(ct, 'wall_setHeatFlux', w.id, f.id); end function setVelocity(w, f) @@ -221,11 +217,11 @@ classdef Wall < handle % to specify a constant velocity by using the polynomial % functor with only the first term specified. % - % :param f: + % parameter f: % Instance of class 'Func'. Unit: m/s. checklib; - calllib(w.lib, 'wall_setVelocity', w.id, f.id); + calllib(ct, 'wall_setVelocity', w.id, f.id); end %% ReactorNet get methods @@ -233,19 +229,19 @@ classdef Wall < handle function a = get.area(w) % Get the area of the wall in m^2. checklib; - a = calllib(w.lib, 'wall_area', w.id); + a = calllib(ct, 'wall_area', w.id); end function q = qdot(w, t) % Get the total heat transfer through a wall at given time. checklib; - q = calllib(w.lib, 'wall_Q', w.id, t); + q = calllib(ct, 'wall_Q', w.id, t); end function v = vdot(w, t) % Get the rate of volumetric change at a given time. checklib; - v = calllib(w.lib, 'wall_vdot', w.id, t); + v = calllib(ct, 'wall_vdot', w.id, t); end end diff --git a/interfaces/matlab_experimental/UnloadCantera.m b/interfaces/matlab_experimental/UnloadCantera.m new file mode 100644 index 000000000..0d9a65665 --- /dev/null +++ b/interfaces/matlab_experimental/UnloadCantera.m @@ -0,0 +1,4 @@ +if libisloaded(ct) + unloadlibrary(ct); +end +disp('Cantera has been unloaded'); diff --git a/interfaces/matlab_experimental/Utility/UnloadCantera.m b/interfaces/matlab_experimental/Utility/UnloadCantera.m deleted file mode 100644 index fe810282d..000000000 --- a/interfaces/matlab_experimental/Utility/UnloadCantera.m +++ /dev/null @@ -1,4 +0,0 @@ -if libisloaded('cantera_shared') - unloadlibrary('cantera_shared'); -end -disp('Cantera has been unloaded'); diff --git a/interfaces/matlab_experimental/Utility/canteraGitCommit.m b/interfaces/matlab_experimental/Utility/canteraGitCommit.m new file mode 100644 index 000000000..563ae4191 --- /dev/null +++ b/interfaces/matlab_experimental/Utility/canteraGitCommit.m @@ -0,0 +1,8 @@ +function v = canteraGitCommit() + % Get Cantera Git commit hash + checklib; + buflen = calllib(ct, 'ct_getGitCommit', 0, ''); + aa = char(zeros(1, buflen)); + [~, aa] = calllib(ct, 'ct_getGitCommit', buflen, aa); + v = aa; +end diff --git a/interfaces/matlab_experimental/Utility/canteraVersion.m b/interfaces/matlab_experimental/Utility/canteraVersion.m new file mode 100644 index 000000000..902037092 --- /dev/null +++ b/interfaces/matlab_experimental/Utility/canteraVersion.m @@ -0,0 +1,8 @@ +function v = canteraVersion() + % Get Cantera version information + checklib; + buflen = calllib(ct, 'ct_getCanteraVersion', 0, ''); + aa = char(zeros(1, buflen)); + [~, aa] = calllib(ct, 'ct_getCanteraVersion', buflen, aa); + v = aa; +end diff --git a/interfaces/matlab_experimental/Utility/checklib.m b/interfaces/matlab_experimental/Utility/checklib.m index 15d988697..fcf4d9ddc 100644 --- a/interfaces/matlab_experimental/Utility/checklib.m +++ b/interfaces/matlab_experimental/Utility/checklib.m @@ -1,3 +1,3 @@ -if ~libisloaded('cantera_shared') +if ~libisloaded(ct) error('Cantera is not loaded'); end diff --git a/interfaces/matlab_experimental/Utility/cleanup.m b/interfaces/matlab_experimental/Utility/cleanup.m index 785f0589b..d13dbe52d 100644 --- a/interfaces/matlab_experimental/Utility/cleanup.m +++ b/interfaces/matlab_experimental/Utility/cleanup.m @@ -1,12 +1,11 @@ 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'); + calllib(ct, 'ct_clearOneDim'); + calllib(ct, 'ct_clearMix'); + calllib(ct, 'ct_clearFunc'); + calllib(ct, 'ct_clearStorage'); + calllib(ct, 'ct_clearReactors'); + calllib(ct, 'ct_clearReactionPath'); clear all end diff --git a/interfaces/matlab_experimental/Utility/getDataDirectories.m b/interfaces/matlab_experimental/Utility/getDataDirectories.m new file mode 100644 index 000000000..16f6eb212 --- /dev/null +++ b/interfaces/matlab_experimental/Utility/getDataDirectories.m @@ -0,0 +1,8 @@ +function d = getDataDirectories() + % Get a cell array of the directories searched for data files. + checklib; + buflen = calllib(ct, 'ct_getDataDirectories', 0, '', ';'); + aa = char(zeros(1, buflen)); + [~, aa, ~] = calllib(ct, 'ct_getDataDirectories', buflen, aa, ';'); + d = aa; +end diff --git a/interfaces/matlab_experimental/Utility/geterr.m b/interfaces/matlab_experimental/Utility/geterr.m new file mode 100644 index 000000000..43be178d2 --- /dev/null +++ b/interfaces/matlab_experimental/Utility/geterr.m @@ -0,0 +1,12 @@ +function e = geterr() + % Get the error message from Cantera error. + checklib; + try + buflen = calllib(ct, 'ct_getCanteraError', 0, ''); + aa = char(zeros(1, buflen)); + [~, aa] = calllib(ct, 'ct_getCanteraError', buflen, aa); + e = aa; + catch ME + e = getReport(ME); + end +end diff --git a/interfaces/matlab_experimental/Utility/importEdge.m b/interfaces/matlab_experimental/Utility/importEdge.m new file mode 100644 index 000000000..a27e4db3c --- /dev/null +++ b/interfaces/matlab_experimental/Utility/importEdge.m @@ -0,0 +1,15 @@ +function s = importEdge(file, name, phase1, phase2, phase3, phase4) + % Import edges between phases. + + if nargin == 3 + s = Interface(file, name, phase1); + elseif nargin == 4 + s = Interface(file, name, phase1, phase2); + elseif nargin == 5 + s = Interface(file, name, phase1, phase2, phase3); + elseif nargin == 6 + s = Interface(file, name, phase1, phase2, phase3, phase4); + else + error('importEdge only supports 4 neighbor phases.'); + end +end diff --git a/interfaces/matlab_experimental/Utility/importInterface.m b/interfaces/matlab_experimental/Utility/importInterface.m new file mode 100644 index 000000000..720979457 --- /dev/null +++ b/interfaces/matlab_experimental/Utility/importInterface.m @@ -0,0 +1,11 @@ +function s = importInterface(file, name, phase1, phase2) + % Import an interface between phases. + + if nargin == 3 + s = Interface(file, name, phase1); + elseif nargin == 4 + s = Interface(file, name, phase1, phase2); + else + error('importInterface only supports 2 bulk phases'); + end +end diff --git a/interfaces/matlab_experimental/cleanup.m b/interfaces/matlab_experimental/cleanup.m deleted file mode 100644 index 785f0589b..000000000 --- a/interfaces/matlab_experimental/cleanup.m +++ /dev/null @@ -1,12 +0,0 @@ -function cleanup() -% Delete all stored Cantera objects and reclaim memory - checklib; - calllib('cantera_shared', 'ct_clearOneDim'); - calllib('cantera_shared', 'ct_clearMix'); - calllib('cantera_shared', 'ct_clearXML'); - calllib('cantera_shared', 'ct_clearFunc'); - calllib('cantera_shared', 'ct_clearStorage'); - calllib('cantera_shared', 'ct_clearReactors'); - calllib('cantera_shared', 'ct_clearReactionPath'); - clear all -end diff --git a/interfaces/matlab_experimental/ct.m b/interfaces/matlab_experimental/ct.m new file mode 100644 index 000000000..dcd4e5f34 --- /dev/null +++ b/interfaces/matlab_experimental/ct.m @@ -0,0 +1,8 @@ +function str = ct() + % Return name of 'cantera_shared' library. + if ispc + str = 'cantera_shared'; + else + str = 'libcantera_shared'; + end +end diff --git a/interfaces/matlab_experimental/readme.txt b/interfaces/matlab_experimental/readme.txt new file mode 100644 index 000000000..4dbd6440a --- /dev/null +++ b/interfaces/matlab_experimental/readme.txt @@ -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 \matlab. + +For first time users, launch matlab, navigate to the toolbox folder, then run the 'SetCanteraPath()' 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.