Various minor fixes to the Matlab toolbox

doc/sphinx/matlab/onedim.rst

@@ -78,6 +78,10 @@ Boundary
 ^^^^^^^^
 .. autoclass:: Boundary
 
+Flow
+^^^^
+.. autoclass:: Flow
+
 Sim1D
 ^^^^^
 .. autoclass:: Sim1D(domains)

doc/sphinx/matlab/zeroD.rst

@@ -80,6 +80,10 @@ MassFlowController
 ^^^^^^^^^^^^^^^^^^
 .. autoclass:: MassFlowController(upstream, downstream[, name])
 
+PressureController
+^^^^^^^^^^^^^^^^^^
+.. autoclass:: PressureController(upstream, downstream[, name])
+
 Valve
 ^^^^^
 .. autoclass:: Valve(upstream, downstream[, name])

doc/sphinx/reference/releasenotes/v4.0.md

@@ -53,8 +53,31 @@ have been updated to return `std::span` instead. For example,
 `const std::vector<double> Phase::molecularWeights()` has been replaced with
 `span<const double> Phase::molecularWeights()`.
 
+### CLib
+
+- The generated CLib function `reactornet_sensitivity`, which wrapped
+  `ReactorNet::sensitivity(const string& component, size_t p, int reactor)`, has been
+  renamed to `reactornet_sensitivityByName`. The name `reactornet_sensitivity` now wraps
+  the index-based overload `ReactorNet::sensitivity(size_t k, size_t p)`, where `k` is
+  an index into the global state vector of the network.
+
 ## New features
 
+- Map between reactor state-vector indices and component names from the CLib and
+  MATLAB interfaces. The generated CLib gains `reactor_neq`,
+  `reactor_componentName` and `reactor_componentIndex`, wrapping
+  `ReactorBase::neq()`, `ReactorBase::componentName(size_t)` and
+  `ReactorBase::componentIndex(const string&)`, plus
+  `reactornet_globalComponentIndex`, wrapping
+  `ReactorNet::globalComponentIndex(const string&, size_t)`; these join the
+  existing `reactornet_neq` and `reactornet_componentName`. The MATLAB toolbox
+  gains the corresponding `ReactorBase.nVars`, `ReactorBase.componentIndex`,
+  `ReactorBase.componentName`, `ReactorNet.nVars`, `ReactorNet.componentName`
+  and `ReactorNet.globalComponentIndex`, all using 1-based indices. Note that
+  `ReactorNet.componentName` returns a name prefixed with that of the reactor it
+  belongs to, for example `reactor1: temperature`, while
+  `ReactorNet.globalComponentIndex` takes an unprefixed name resolved within a
+  given reactor.
 - Make an analytic Jacobian the default for 1D flame simulations. The species
   mass-fraction Jacobian columns at interior grid points are evaluated
   analytically from kinetics concentration-derivative functions instead of by

interfaces/clib/SConscript

@@ -30,7 +30,7 @@ env.Alias("doxygen_tags", tags)
 # Generated file names can be anticipated
 generated_files = []
 auto_path = Path(Dir("#interfaces/sourcegen/src/sourcegen/headers").abspath)
-yaml_files = auto_path.glob("ct*.yaml")
+yaml_files = sorted(auto_path.glob("ct*.yaml"))
 for yaml_file in yaml_files:
     base_name = yaml_file.stem
     generated_files.extend([

interfaces/julia/src/reactornet.jl

@@ -100,10 +100,16 @@ end
 
 """
     sensitivity(net, component, p, reactor) -> Float64
+    sensitivity(net, k, p) -> Float64
 
-Normalized sensitivity of `component` (e.g. `"temperature"` or a species name)
-in `reactor` with respect to sensitivity parameter `p` (1-based).  `reactor` may
-be a [`Reactor`](@ref) belonging to the network or its 1-based position.
+Normalized sensitivity with respect to sensitivity parameter `p` (1-based).
+
+In the first form, `component` is the name of a state variable (for example
+`"temperature"` or a species name) resolved within `reactor`, which may be a
+[`Reactor`](@ref) belonging to the network or its 1-based position.
+
+In the second form, `k` is the 1-based index of the component in the global
+state vector of the network, as ordered by [`component_names`](@ref).
 """
 function sensitivity(net::ReactorNet, component::AbstractString, p::Integer,
                      reactor::Reactor)
@@ -114,8 +120,14 @@ end
 
 function sensitivity(net::ReactorNet, component::AbstractString, p::Integer,
                      reactor::Integer)
-    return checkd(LibCantera.reactornet_sensitivity(net.handle, component,
-                                                    Int32(p - 1), Int32(reactor - 1)))
+    return checkd(LibCantera.reactornet_sensitivityByName(net.handle, component,
+                                                          Int32(p - 1),
+                                                          Int32(reactor - 1)))
+end
+
+function sensitivity(net::ReactorNet, k::Integer, p::Integer)
+    return checkd(LibCantera.reactornet_sensitivity(net.handle, Int32(k - 1),
+                                                    Int32(p - 1)))
 end
 
 """

interfaces/julia/test/test_connectors.jl

@@ -100,6 +100,12 @@ end
     advance!(net, 2e-4)
     s = sensitivity(net, "temperature", 1, r)
     @test isfinite(s)
+    # the same value is reachable through the global state-vector index
+    k = findfirst(endswith(": temperature"), component_names(net))
+    @test k !== nothing
+    @test sensitivity(net, k, 1) == s
+    @test sensitivity(net, "temperature", 1, 1) == s
+    @test_throws CanteraError sensitivity(net, "spam", 1, 1)
     @test rtol(net) > 0
     @test atol(net) > 0
 end

interfaces/matlab/+ct/+impl/errorCode.m

@@ -1,3 +1,10 @@
 function err = errorCode()
-    err = [-1, -2, -999.999, double(intmax('uint64'))];
+    % Values returned by CLib functions to signal that an exception was raised.
+    %
+    % `ERR` (-999) and `DERR` (-999.999) are the sentinels defined by
+    % `clib_defs.h`; -1 and -2 are returned in their place by generated functions
+    % whose normal return value is a length or a handle. `intmax('uint64')` is the
+    % C++ `npos`.
+
+    err = [-1, -2, -999, -999.999, double(intmax('uint64'))];
 end

interfaces/matlab/+ct/+oneD/ReactingSurface.m

@@ -10,12 +10,8 @@
     %     String ID of the reacting surface.
 
     properties
-        % Set bounds on the solution components. ::
-        %
-        %     >> d.coverageEnabled = flag
-        %
-        % :param flag:
-        %     Boolean flag indicating whether coverage equations are enabled.
+        % Boolean flag that controls whether or not the surface coverage equations
+        % are solved as part of the 1D problem.
         coverageEnabled
     end
 

interfaces/matlab/+ct/+zeroD/FlowDevice.m

@@ -12,12 +12,13 @@
     % pressure difference equals the difference in pressure between the
     % upstream and downstream reactors.
     %
-    % See also: :mat:class:`ct.zeroD.MassFlowController`, :mat:class:`ct.zeroD.Valve`
+    % See also: :mat:class:`ct.zeroD.MassFlowController`,
+    % :mat:class:`ct.zeroD.PressureController`, :mat:class:`ct.zeroD.Valve`
     %
     % :param typ:
     %     Type of :mat:class:`ct.zeroD.FlowDevice` to be created. ``typ='MassFlowController'``
     %     for :mat:class:`ct.zeroD.MassFlowController`,  ``typ='PressureController'`` for
-    %     :mat:class:`ct.PressureController`, and ``typ='Valve'`` for
+    %     :mat:class:`ct.zeroD.PressureController`, and ``typ='Valve'`` for
     %     :mat:class:`ct.zeroD.Valve`.
     % :param upstream:
     %     Upstream reactor or reservoir.
@@ -53,6 +54,14 @@
         % as long as this produces a positive value.  If this expression is
         % negative, zero is returned.
         valveCoeff
+
+        % The coefficient defining the behavior of this
+        % :mat:class:`ct.zeroD.FlowDevice`, with a meaning that depends on the type
+        % of the device: the mass flow rate [kg/s] for a
+        % :mat:class:`ct.zeroD.MassFlowController`, the valve coefficient [kg/s/Pa]
+        % for a :mat:class:`ct.zeroD.Valve`, and the pressure coefficient [kg/s/Pa]
+        % for a :mat:class:`ct.zeroD.PressureController`.
+        deviceCoefficient
     end
 
     methods
@@ -77,6 +86,10 @@
             mdot = ct.impl.call('mFlowdev_massFlowRate', obj.id);
         end
 
+        function c = get.deviceCoefficient(obj)
+            c = ct.impl.call('mFlowdev_deviceCoefficient', obj.id);
+        end
+
         %% FlowDevice Set Methods
 
         function set.massFlowRate(obj, mdot)
@@ -104,8 +117,12 @@ function setPrimary(obj, d)
             % :param d:
             %     Instance of class :mat:class:`ct.zeroD.FlowDevice`.
 
+            if ~isa(d, 'ct.zeroD.FlowDevice')
+                error('Primary device must be an instance of ct.zeroD.FlowDevice.');
+            end
+
             if strcmp(obj.type, 'PressureController')
-                ct.impl.call('mFlowdev_setPrimary', obj.id, d);
+                ct.impl.call('mFlowdev_setPrimary', obj.id, d.id);
             else
                 error('Primary flow device can only be set for pressure controllers.');
             end
@@ -118,7 +135,11 @@ function setPrimary(obj, d)
                 error('Valve coefficient can only be set for valves.');
             end
 
-            ct.impl.call('mFlowdev_setDeviceCoefficient', obj.id, k);
+            obj.deviceCoefficient = k;
+        end
+
+        function set.deviceCoefficient(obj, c)
+            ct.impl.call('mFlowdev_setDeviceCoefficient', obj.id, c);
         end
 
     end

interfaces/matlab/+ct/+zeroD/FlowReactor.m

@@ -24,7 +24,13 @@
 
     properties (SetAccess = public)
 
-        massFlowRate % Mass flow rate [kg/s].
+        % Mass flow rate through the reactor [kg/s].
+        %
+        % Setting the mass flow rate sets the flow speed based on the current
+        % density of the reactor contents and the reactor area. The value read
+        % back is computed from the flow speed, density and area at the end of
+        % the last call to "advance" or "step".
+        massFlowRate
 
     end
 
@@ -43,15 +49,14 @@
 
         %% FlowReactor Get Methods
 
-        function flag = get.massFlowRate(obj)
+        function rate = get.massFlowRate(obj)
             rate = ct.impl.call('mReactor_massFlowRate', obj.id);
         end
 
         %% FlowReactor Set Methods
 
         function set.massFlowRate(obj, MFR)
             ct.impl.call('mReactor_setMassFlowRate', obj.id, MFR);
-            obj.massFlowRate = MFR;
         end
 
     end

interfaces/matlab/+ct/+zeroD/PressureController.m

@@ -0,0 +1,68 @@
+classdef PressureController < ct.zeroD.FlowDevice
+    % Create a pressure controller. ::
+    %
+    %     >> p = ct.zeroD.PressureController(upstream, downstream, name)
+    %
+    % Creates an instance of class :mat:class:`ct.zeroD.FlowDevice` configured to
+    % simulate a pressure controller. A :mat:class:`ct.zeroD.PressureController` is
+    % designed to be used in conjunction with another primary flow device, typically a
+    % :mat:class:`ct.zeroD.MassFlowController`. The primary flow device is installed on
+    % the inlet of the reactor, and the corresponding
+    % :mat:class:`ct.zeroD.PressureController` is installed on the outlet of the
+    % reactor. The mass flow rate through the :mat:class:`ct.zeroD.PressureController`
+    % is equal to the mass flow rate of the primary device, plus a correction that
+    % depends on the pressure difference:
+    %
+    % .. math:: \dot{m} = \dot{m}_{primary} + K_v(P_{upstream} - P_{downstream})
+    %
+    % as long as this produces a positive value. If this expression is negative, zero
+    % is returned. The primary device is set using
+    % :mat:meth:`ct.zeroD.FlowDevice.setPrimary` and must be set before the mass flow
+    % rate can be evaluated.
+    %
+    % see also: :mat:class:`ct.zeroD.FlowDevice`,
+    % :mat:class:`ct.zeroD.MassFlowController`, :mat:class:`ct.zeroD.Valve`
+    %
+    % :param upstream:
+    %     Upstream :mat:class:`ct.zeroD.ReactorBase`.
+    % :param downstream:
+    %     Downstream :mat:class:`ct.zeroD.ReactorBase`.
+    % :param name:
+    %     Flow device name (optional; default is ``(none)``).
+
+    properties (SetAccess = public)
+
+        % Proportionality constant :math:`K_v` in kg/s/Pa between the pressure drop
+        % and the mass flow rate correction added to the primary device's mass flow
+        % rate.
+        pressureCoeff
+
+    end
+
+    methods
+
+        function obj = PressureController(upstream, downstream, name)
+            arguments
+                upstream (1,1) ct.zeroD.ReactorBase
+                downstream (1,1) ct.zeroD.ReactorBase
+                name (1,1) string = "(none)"
+            end
+
+            obj@ct.zeroD.FlowDevice('PressureController', upstream, downstream, name);
+        end
+
+        %% PressureController Get Methods
+
+        function k = get.pressureCoeff(obj)
+            k = obj.deviceCoefficient;
+        end
+
+        %% PressureController Set Methods
+
+        function set.pressureCoeff(obj, k)
+            obj.deviceCoefficient = k;
+        end
+
+    end
+
+end