GreenPolyProducts/cantera
5
0
Fork 0
mirror of https://github.com/Cantera/cantera.git synced 2025-02-25 18:55:29 -06:00
Code Issues Packages Projects Releases Wiki Activity

Redo 1D class to use DomainFactory See PR #1445

Browse Source
This commit is contained in:
ssun30 2023-04-20 23:24:09 -04:00 committed by Ray Speth
parent 4f78c0f186
commit 7e6f9c999e
23 changed files with 773 additions and 905 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -1,696 +0,0 @@
classdef Domain1D < handle
% Domain1D Class ::
%
% >> d = Domain1D(a, b)
%
% :param a:
% String type of domain. Possible values are:
% - `StagnationFlow`
% - `AxisymmetricFlow`
% - `Inlet1D`
% - `Surf1D`
% - `Symm1D`
% - `Outlet1D`
% - `ReactingSurface`
% - `OutletRes`
% :param b:
% Instance of :mat:class:`Solution` (for ``a == 'StagnationFlow`` or
% ``a = 'AxisymmetricFlow'``) or :mat:class:`Interface`
% (for ``a == 'ReactingSurface'``).
% Not used for all other valid values of ``a``.
properties (SetAccess = immutable)
domainID % ID of the domain
type % Type of the domain
end
properties (SetAccess = public)
T % Boundary Temperature. Units: K.
P % Boundary Pressure. Units: Pa.
end
properties (SetAccess = protected)
% Domain index. ::
%
% >> i = d.domainIndex
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :return:
% Integer flag denoting the location of the domain,
% beginning with 1 at the left.
domainIndex
% Type of the domain. ::
%
% >> i = d.domainType
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :return:
% Integer flag denoting the domain type.
domainType
% Determines whether a domain is a flow. ::
%
% >> a = d.isFlow
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :return:
% True if the domain is a flow domain, and false otherwise.
isFlow
% Determines whether a domain is an inlet. ::
%
% >> a = d.isInlet
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :return:
% True if the domain is an inlet, and false otherwise.
isInlet
% Determines if a domain is a surface. ::
%
% >> a = d.isSurface
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :return:
% True if the domain is a surface, and false otherwise.
isSurface
% Mass flux. ::
%
% >> mdot = d.massFlux
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :return:
% The mass flux in the domain.
massFlux
% Number of components. ::
%
% >> n = d.nComponents
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :return:
% Number of variables at each grid point.
nComponents
% Get the number of grid points. ::
%
% >> n = d.nPoints
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :return:
% Integer number of grid points.
nPoints
end
methods
%% Domain1D Class Constructor.
function d = Domain1D(varargin)
% Create a :mat:class:`Domain1D` object.
ctIsLoaded;
if nargin == 1
a = varargin{1};
if strcmp(a, 'Inlet1D')
d.domainID = ctFunc('inlet_new');
elseif strcmp(a, 'Surf1D')
d.domainID = ctFunc('surf_new');
elseif strcmp(a, 'Symm1D')
d.domainID = ctFunc('symm_new');
elseif strcmp(a, 'Outlet1D')
d.domainID = ctFunc('outlet_new');
elseif strcmp(a, 'OutletRes')
d.domainID = ctFunc('outletres_new');
else
error('Not enough arguments for that job number');
end
elseif nargin == 2
a = varargin{1};
b = varargin{2};
if strcmp(a, 'StagnationFlow')
if isa(b, 'Solution')
d.domainID = ctFunc('stflow_new', b.tpID, b.kinID, b.trID, 1);
else
error('Wrong argument type. Expecting instance of class Solution.');
end
elseif strcmp(a, 'AxisymmetricFlow')
if isa(b, 'Solution')
d.domainID = ctFunc('stflow_new', b.tpID, b.kinID, b.trID, 2);
else
error('Wrong argument type. Expecting instance of class Solution.');
end
elseif strcmp(a, 'ReactingSurface')
if isa(b, 'Interface')
d.domainID = ctFunc('reactingsurf_new');
ctFunc('reactingsurf_setkineticsmgr', d.domainID, b.kinID);
else
error('Wrong argument type. Expecting instance of class Interface.');
end
else
error('Wrong object type.');
end
end
d.type = a;
end
%% Domain1D Class Destructor
function delete(d)
% Delete the :mat:class:`Domain1D` object.
ctFunc('domain_del', d.domainID);
end
%% Domain1D Utility Methods
function d = disableEnergy(d)
% Disable the energy equation. ::
%
% >> d.disableEnergy
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
disp(' ');
disp('Disabling the energy equation...');
ctFunc('stflow_solveEnergyEqn', d.domainID, 0);
end
function d = enableEnergy(d)
% Enable the energy equation. ::
%
% >> d.enableEnergy
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
disp(' ');
disp('Enabling the energy equation...');
ctFunc('stflow_solveEnergyEqn', d.domainID, 1);
end
function d = disableSoret(d)
% Disable the diffusive mass fluxes due to the Soret effect. ::
%
% >> d.disableSoret
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
disp(' ');
disp('Disabling the Soret effect...');
ctFunc('stflow_enableSoret', d.domainID, 0);
end
function d = enableSoret(d)
% Enable the diffusive mass fluxes due to the Soret effect. ::
%
% >> d.enableSoret
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
disp(' ');
disp('Disabling the Soret effect...');
ctFunc('stflow_enableSoret', d.domainID, 1);
end
%% Domain Get Methods
function b = bounds(d, component)
% Get the (lower, upper) bounds for a solution component. ::
%
% >> b = d.bounds(componoent)
%
% :param component:
% String name of the component for which the bounds are returned.
% :return:
% :math:`1\times 2` vector of the lower and upper bounds.
n = d.componentIndex(component);
lower = ctFunc('domain_lowerBound', d.domainID, n);
upper = ctFunc('domain_upperBound', d.domainID, n);
b = [lower, upper];
end
function n = componentIndex(d, name)
% Index of a component given its name. ::
%
% >>n = d.componentIndex(name)
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :param name:
% String name of the component to look up. If a numeric value
% is passed, it will be returned.
% :return:
% Index of the component, or input numeric value.
if isa(name, 'double')
n = name;
else
n = ctFunc('domain_componentIndex', d.domainID, name);
if n >= 0
n = n + 1;
end
end
if n <= 0
error('Component not found');
end
end
function s = componentName(d, index)
% Name of a component given its index. ::
%
% >> n = d.componentName(index)
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :param index:
% Integer or vector of integers of component names to get.
% :return:
% Cell array of component names.
n = length(index);
s = cell(1, n);
for i = 1:n
id = index(i) - 1;
output = ctString('domain_componentName', d.domainID, id);
s{i} = output;
end
end
function i = get.domainIndex(d)
i = ctFunc('domain_index', d.domainID);
if i >= 0
i = i + 1;
end
if i <= 0
error('Domain not found');
end
end
function i = get.domainType(d)
i = ctFunc('domain_type', d.domainID);
end
function zz = gridPoints(d, n)
% Grid points from a domain. ::
%
% >> zz = d.gridPoints(n)
%
% :param d:
% Instance of class 'Domain1D'.
% :param n:
% Optional, vector of grid points to be retrieved.
% :return:
% Vector of grid points.
if nargin == 1
np = d.nPoints;
zz = zeros(1, np);
for i = 1:np
zz(i) = ctFunc('domain_grid', d.domainID, i - 1);
end
else
m = length(n);
zz = zeros(1, m);
for i = 1:m
zz(i) = ctFunc('domain_grid', d.domainID, n(i) - 1);
end
end
end
function a = get.isFlow(d)
t = d.domainType;
% See Domain1D.h for definitions of constants.
a = t < 100;
end
function a = get.isInlet(d)
t = d.domainType;
% See Domain1D.h for definitions of constants.
a = t == 104;
end
function a = get.isSurface(d)
t = d.domainType;
% See Domain1D.h for definitions of constants.
a = t == 102;
end
function mdot = get.massFlux(d)
mdot = ctFunc('bdry_mdot', d.domainID);
end
function y = massFraction(d, k)
% Get the mass fraction of a species given its integer index. ::
%
% >> y = d.massFraction(k)
%
% 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.
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :param k:
% Integer species index.
% :return:
% Mass fraction of species.
if d.domainIndex == 0
error('No flow domain attached!')
end
if d.isInlet
y = ctFunc('bdry_massFraction', d.domainID, k - 1);
else error('Input domain must be an inlet');
end
end
function n = get.nComponents(d)
n = ctFunc('domain_nComponents', d.domainID);
end
function n = get.nPoints(d)
n = ctFunc('domain_nPoints', d.domainID);
end
function tol = tolerances(d, component)
% Return the (relative, absolute) error tolerances for a
% solution component. ::
%
% >> tol = d.tolerances(component)
%
% :param component:
% String name of the component for which the bounds are returned.
% :return:
% :math:`1\times 2` vector of the relative and absolute error tolerances.
n = d.componentIndex(component);
rerr = ctFunc('domain_rtol', d.domainID, n);
aerr = ctFunc('domain_atol', d.domainID, n);
tol = [rerr, aerr];
end
function temperature = get.T(d)
temperature = ctFunc('bdry_temperature', d.domainID);
end
function pressure = get.P(d)
pressure = calllibt(ct, 'stflow_pressure', d.domainID);
end
%% Domain Set Methods
function setBounds(d, component, lower, upper)
% Set bounds on the solution components. ::
%
% >> d.setBounds(component, lower, upper)
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :param component:
% String, component to set the bounds on.
% :param lower:
% Lower bound.
% :param upper:
n = d.componentIndex(component);
ctFunc('domain_setBounds', d.domainID, n - 1, lower, upper);
end
function setCoverageEqs(d, onoff)
% Set bounds on the solution components. ::
%
% >> d.setCoverageEqn(onoff)
%
% :param d:
% Instance of class :mat:class:`Domain1D`
% :param onoff:
% String, one of ``'on'`` or ``'yes'`` to turn solving
% the coverage equations on. One of ``'off'`` or ``'no'``
% to turn off the coverage equations.
if ~strcmp(d.type, 'ReactingSurface')
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
ctFunc('reactingsurf_enableCoverageEqs', d.domainID, ion);
end
function setFixedTempProfile(d, profile)
% Set a fixed temperature profile. ::
%
% >> d.setFixedTempProfile(profile)
%
% Set the 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.
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :param profile:
% :math:`n\times 2` or :math:`2\times n` array of ``n`` points
% at which the temperature is specified.
sz = size(profile);
if sz(1) == 2
l = length(profile(1, :));
ctFunc('stflow_setFixedTempProfile', d.domainID, ...
l, profile(1, :), l, profile(2, :));
elseif sz(2) == 2
l = length(profile(:, 1));
ctFunc('stflow_setFixedTempProfile', d.domainID, ...
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. ::
%
% >> d.setID(id)
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :param id:
% String ID to assign.
ctFunc('domain_setID', d.domainID, id);
end
function setMassFlowRate(d, mdot)
% Set the mass flow rate. ::
%
% >> d.setMassFlowRate(mdot)
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :param mdot:
% Mass flow rate.
ctFunc('bdry_setMdot', d.domainID, mdot);
end
function setMoleFractions(d, x)
% Set the mole fractions. ::
%
% >> d.setMoleFractions(x)
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :param x:
% String specifying the species and mole fractions in
% the format ``'SPEC:X,SPEC2:X2'``.
ctFunc('bdry_setMoleFractions', d.domainID, x);
end
function setSteadyTolerances(d, component, rtol, atol)
% Set the steady-state tolerances. ::
%
% >>d.setSteadyTolerances(component, rtol, atol)
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :param 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.
% :param rtol:
% Relative tolerance.
% :param atol:
% Absolute tolerance.
if strcmp(component, 'default')
nc = d.nComponents;
for ii = 1:nc
ctFunc('domain_setSteadyTolerances', ...
d.domainID, ii - 1, rtol, atol);
end
elseif iscell(component)
nc = length(component);
for ii = 1:nc
n = d.componentIndex(component{ii});
ctFunc('domain_setSteadyTolerances', d.domainID, n, rtol, atol);
end
else
n = d.componentIndex(component);
ctFunc('domain_setSteadyTolerances', d.domainID, n, rtol, atol);
end
end
function setTransientTolerances(d, component, rtol, atol)
% Set the transient tolerances. ::
%
% >> d.setTransientTolerances(component, rtol, atol)
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :param 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.
% :param rtol:
% Relative tolerance.
% :param atol:
% Absolute tolerance.
if strcmp(component, 'default')
nc = d.nComponents;
for ii = 1:nc
ctFunc('domain_setTransientTolerances', ...
d.domainID, ii - 1, rtol, atol);
end
elseif iscell(component)
nc = length(component);
for ii = 1:nc
n = d.componentIndex(component{ii});
ctFunc('domain_setTransientTolerances', ...
d.domainID, n, rtol, atol);
end
else
n = d.componentIndex(component);
ctFunc('domain_setTransientTolerances', ...
d.domainID, n, rtol, atol);
end
end
function setTransport(d, itr)
% Set the solution object used for calculating transport properties. ::
%
% >> d.setTransport(itr)
%
% :param itr:
% ID of the solution object for which transport properties
% are calculated.
ctFunc('stflow_setTransport', d.domainID, itr);
end
function setupGrid(d, grid)
% Set up the solution grid. ::
%
% >> d.setupGrid(grid)
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :param grid:
% The grid for this domain.
ctFunc('domain_setupGrid', d.domainID, numel(grid), grid);
end
function set.T(d, t)
if t <= 0
error('The temperature must be positive');
end
ctFunc('bdry_setTemperature', d.domainID, t);
end
function set.P(d, p)
if p <= 0
error('The pressure must be positive');
end
ctFunc('stflow_setPressure', d.domainID, p);
end
end
end

31
interfaces/matlab_experimental/1D/FreeFlame.m
View File

@ -1,31 +0,0 @@
classdef FreeFlame < Domain1D
% Create a freely-propagating flat flame. ::
%
% >> m = FreeFlame(gas, id)
%
% :param gas:
% Instance of class :mat:class:`Solution`.
% :param id:
% String, ID of the flow.
% :return:
% Instance of class :mat:class:`FreeFlame` representing
% a freely propagating, adiabatic flame.
methods
function m = FreeFlame(gas, id)
% Constructor
m@Domain1D('StagnationFlow', gas, 2);
if nargin == 1
m.setID('flame');
else
m.setID(id);
end
end
end
end

37
interfaces/matlab_experimental/1D/Surface.m
View File

@ -1,37 +0,0 @@
classdef Surface < Domain1D
% Create a surface domain. ::
%
% >> m = Surface(id, surface_mech)
%
% :param id:
% String ID of surface
% :param surface_mech:
% Instance of class :mat:class:`Interface` defining
% the surface reaction mechanism to be used. Optional.
% :return:
% Instance of class :mat:class:`Surface` representing a
% non-reacting or reacting surface.
methods
function m = Surface(id, surface_mech)
% Constructor
if nargin < 2
param = {'Surf1D'};
else
param = {'ReactingSurface', surface_mech};
end
m@Domain1D(param{:});
if nargin == 0
m.setID('surface');
elseif nargin == 1
m.setID(id);
end
end
end
end

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

@ -47,8 +47,6 @@ classdef Interface < handle & ThermoPhase & Kinetics
args = varargin(3:end);
s@Kinetics(t, src, id, args{:});
s.tpClear;
s.tpID = t.tpID;
end
@ -56,8 +54,6 @@ classdef Interface < handle & ThermoPhase & Kinetics
function delete(s)
% Delete :mat:class:`Interface` object.
s.tpClear;
disp('Interface class object has been deleted');
end

18
interfaces/matlab_experimental/1D/AxisymmetricFlow.m → interfaces/matlab_experimental/OneDim/AxisymmetricFlow.m
View File

@ -1,9 +1,9 @@
classdef AxisymmetricFlow < Domain1D
classdef AxisymmetricFlow < Flow1D
% Create an axisymmetric flow domain. ::
%
% >> m = AxisymmetricFlow(gas, id)
% >> m = AxisymmetricFlow(phase, id)
%
% :param gas:
% :param phase:
% Instance of class :mat:class:`Solution`.
% :param id:
% String, ID of the flow.
@ -12,17 +12,15 @@ classdef AxisymmetricFlow < Domain1D
methods
function m = AxisymmetricFlow(gas, id)
function m = AxisymmetricFlow(phase, id)
% Constructor
m@Domain1D('StagnationFlow', gas);
if nargin == 1
m.setID('flow');
else
m.setID(id);
if nargin < 2
id = 'axisymmetric-flow';
end
m@Flow1D('axisymmetric-flow', phase, id);
end
end

101
interfaces/matlab_experimental/OneDim/Boundary1D.m Normal file
View File

@ -0,0 +1,101 @@
classdef Boundary1D < Domain1D
% Create a Boundary domain. ::
%
% >> m = Boundary(phase, id)
%
% :param phase:
% Instance of class :mat:class:`Solution`.
% :param id:
% String, ID of the flow.
% :return:
% Instance of class :mat:class:`Boundary`.
properties
% Mass flux. ::
%
% >> mdot = d.massFlux
%
% :param d:
% Instance of class :mat:class:`Boundary`.
% :return:
% The mass flux in the domain. Unit: [kg/s/m^2].
massFlux
T % Boundary Temperature. Units: K.
end
methods
%% Boundary Class Constructor
function b = Boundary1D(type, phase, id)
b@Domain1D(type, phase, id);
end
%% Boundary Methods
function temperature = get.T(d)
temperature = ctFunc('bdry_temperature', d.domainID);
end
function set.T(d, t)
ctFunc('bdry_setTemperature', d.domainID, t);
end
function mdot = get.massFlux(d)
mdot = ctFunc('bdry_mdot', d.domainID);
end
function set.massFlux(d, mdot)
ctFunc('bdry_setMdot', d.domainID, mdot);
end
function y = massFraction(d, k)
% Get the mass fraction of a species given its integer index. ::
%
% >> y = d.massFraction(k)
%
% 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.
%
% :param d:
% Instance of class :mat:class:`Boundary`.
% :param k:
% Integer species index.
% :return:
% Mass fraction of species.
if d.domainIndex == 0
error('No flow domain attached!')
end
if strcmp(d.domainType, 'inlet')
y = ctFunc('bdry_massFraction', d.domainID, k - 1);
else
error('Input domain must be an inlet');
end
end
function setMoleFractions(d, x)
% Set the mole fractions. ::
%
% >> d.setMoleFractions(x)
%
% :param d:
% Instance of class :mat:class:`Boundary`.
% :param x:
% String specifying the species and mole fractions in
% the format ``'SPEC:X,SPEC2:X2'``.
ctFunc('bdry_setMoleFractions', d.domainID, x);
end
end
end

6
interfaces/matlab_experimental/1D/CounterFlowDiffusionFlame.m → interfaces/matlab_experimental/OneDim/CounterFlowDiffusionFlame.m
View File

@ -45,15 +45,15 @@ classdef CounterFlowDiffusionFlame < Sim1D
error('Oxidizer gas object must represent an ideal gas mixture.');
end
if ~left.isInlet
if ~isa(left, 'Inlet')
error('Left inlet object of wrong type.');
end
if ~flow.isFlow
if ~isa(flow, 'Flow1D')
error('Flow object of wrong type.');
end
if ~right.isInlet
if ~isa(right, 'Inlet')
error('Right inlet object of wrong type.');
end

378
interfaces/matlab_experimental/OneDim/Domain1D.m Normal file
View File

@ -0,0 +1,378 @@
classdef Domain1D < handle
% Domain1D Class ::
%
% >> d = Domain1D(a, b)
%
% :param a:
% String type of domain. Possible values are:
% - `StagnationFlow`
% - `AxisymmetricFlow`
% - `Inlet1D`
% - `Surf1D`
% - `Symm1D`
% - `Outlet1D`
% - `ReactingSurface`
% - `OutletRes`
% :param b:
% Instance of :mat:class:`Solution` (for ``a == 'StagnationFlow`` or
% ``a = 'AxisymmetricFlow'``) or :mat:class:`Interface`
% (for ``a == 'ReactingSurface'``).
% Not used for all other valid values of ``a``.
properties (SetAccess = immutable)
domainID % ID of the domain
end
properties (SetAccess = public)
% Boolean flag indicating whether the energy equation is enabled.
energyEnabled
% Boolean flag indicating whether the diffusive mass fluxes due to the Soret
% effect is enabled.
soretEnabled
% ID of the solution object used for calculating transport properties.
transport
end
properties (SetAccess = protected)
% Domain index. ::
%
% >> i = d.domainIndex
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :return:
% Integer flag denoting the location of the domain,
% beginning with 1 at the left.
domainIndex
% Type of the domain. ::
%
% >> i = d.domainType
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :return:
% Integer flag denoting the domain type.
domainType
% Number of components. ::
%
% >> n = d.nComponents
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :return:
% Number of variables at each grid point.
nComponents
% Get the number of grid points. ::
%
% >> n = d.nPoints
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :return:
% Integer number of grid points.
nPoints
end
methods
%% Domain1D Class Constructor.
function d = Domain1D(type, phase, id)
% Create a :mat:class:`Domain1D` object.
ctIsLoaded;
d.domainID = ctFunc('domain_new', type, phase.phaseID, id);
end
%% Domain1D Class Destructor
function delete(d)
% Delete the :mat:class:`Domain1D` object.
ctFunc('domain_del', d.domainID);
end
%% Domain1D Utility Methods
function set.energyEnabled(d, flag)
d.energyEnabled = flag;
ctFunc('stflow_solveEnergyEqn', d.domainID, int8(flag));
end
function set.soretEnabled(d, flag)
d.soretEnabled = flag;
ctFunc('stflow_enableSoret', d.domainID, int8(flag));
end
%% Domain Get Methods
function b = bounds(d, component)
% Get the (lower, upper) bounds for a solution component. ::
%
% >> b = d.bounds(componoent)
%
% :param component:
% String name of the component for which the bounds are returned.
% :return:
% :math:`1\times 2` vector of the lower and upper bounds.
n = d.componentIndex(component);
lower = ctFunc('domain_lowerBound', d.domainID, n);
upper = ctFunc('domain_upperBound', d.domainID, n);
b = [lower, upper];
end
function n = componentIndex(d, name)
% Index of a component given its name. ::
%
% >>n = d.componentIndex(name)
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :param name:
% String name of the component to look up. If a numeric value
% is passed, it will be returned.
% :return:
% Index of the component, or input numeric value.
if isa(name, 'double')
n = name;
else
n = ctFunc('domain_componentIndex', d.domainID, name);
if n >= 0
n = n + 1;
end
end
if n <= 0
error('Component not found');
end
end
function s = componentName(d, index)
% Name of a component given its index. ::
%
% >> n = d.componentName(index)
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :param index:
% Integer or vector of integers of component names to get.
% :return:
% Cell array of component names.
n = length(index);
s = cell(1, n);
for i = 1:n
id = index(i) - 1;
output = ctString('domain_componentName', d.domainID, id);
s{i} = output;
end
end
function i = get.domainIndex(d)
i = ctFunc('domain_index', d.domainID) + 1;
if i <= 0
error('Domain not found');
end
end
function str = get.domainType(d)
str = ctString('domain_type3', d.domainID);
end
function zz = gridPoints(d, n)
% Grid points from a domain. ::
%
% >> zz = d.gridPoints(n)
%
% :param d:
% Instance of class 'Domain1D'.
% :param n:
% Optional, vector of grid points to be retrieved.
% :return:
% Vector of grid points.
if nargin == 1
np = d.nPoints;
zz = zeros(1, np);
for i = 1:np
zz(i) = ctFunc('domain_grid', d.domainID, i - 1);
end
else
m = length(n);
zz = zeros(1, m);
for i = 1:m
zz(i) = ctFunc('domain_grid', d.domainID, n(i) - 1);
end
end
end
function n = get.nComponents(d)
n = ctFunc('domain_nComponents', d.domainID);
end
function n = get.nPoints(d)
n = ctFunc('domain_nPoints', d.domainID);
end
function tol = tolerances(d, component)
% Return the (relative, absolute) error tolerances for a
% solution component. ::
%
% >> tol = d.tolerances(component)
%
% :param component:
% String name of the component for which the bounds are returned.
% :return:
% :math:`1\times 2` vector of the relative and absolute error tolerances.
n = d.componentIndex(component);
rerr = ctFunc('domain_rtol', d.domainID, n);
aerr = ctFunc('domain_atol', d.domainID, n);
tol = [rerr, aerr];
end
%% Domain Set Methods
function setBounds(d, component, lower, upper)
% Set bounds on the solution components. ::
%
% >> d.setBounds(component, lower, upper)
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :param component:
% String, component to set the bounds on.
% :param lower:
% Lower bound.
% :param upper:
n = d.componentIndex(component);
ctFunc('domain_setBounds', d.domainID, n - 1, lower, upper);
end
function setSteadyTolerances(d, component, rtol, atol)
% Set the steady-state tolerances. ::
%
% >>d.setSteadyTolerances(component, rtol, atol)
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :param 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.
% :param rtol:
% Relative tolerance.
% :param atol:
% Absolute tolerance.
if strcmp(component, 'default')
nc = d.nComponents;
for ii = 1:nc
ctFunc('domain_setSteadyTolerances', ...
d.domainID, ii - 1, rtol, atol);
end
elseif iscell(component)
nc = length(component);
for ii = 1:nc
n = d.componentIndex(component{ii});
ctFunc('domain_setSteadyTolerances', d.domainID, n, rtol, atol);
end
else
n = d.componentIndex(component);
ctFunc('domain_setSteadyTolerances', d.domainID, n, rtol, atol);
end
end
function setTransientTolerances(d, component, rtol, atol)
% Set the transient tolerances. ::
%
% >> d.setTransientTolerances(component, rtol, atol)
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :param 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.
% :param rtol:
% Relative tolerance.
% :param atol:
% Absolute tolerance.
if strcmp(component, 'default')
nc = d.nComponents;
for ii = 1:nc
ctFunc('domain_setTransientTolerances', ...
d.domainID, ii - 1, rtol, atol);
end
elseif iscell(component)
nc = length(component);
for ii = 1:nc
n = d.componentIndex(component{ii});
ctFunc('domain_setTransientTolerances', ...
d.domainID, n, rtol, atol);
end
else
n = d.componentIndex(component);
ctFunc('domain_setTransientTolerances', ...
d.domainID, n, rtol, atol);
end
end
function set.transport(d, itr)
ctFunc('stflow_setTransport', d.domainID, itr);
end
function setupGrid(d, grid)
% Set up the solution grid. ::
%
% >> d.setupGrid(grid)
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :param grid:
% The grid for this domain.
ctFunc('domain_setupGrid', d.domainID, numel(grid), grid);
end
end
end

70
interfaces/matlab_experimental/OneDim/Flow1D.m Normal file
View File

@ -0,0 +1,70 @@
classdef Flow1D < Domain1D
% Create a Flow domain. ::
%
% >> m = Flow1D(phase, id)
%
% :param phase:
% Instance of class :mat:class:`Solution`.
% :param id:
% String, ID of the flow.
% :return:
% Instance of class :mat:class:`Flow1D`.
properties
P % Flow Pressure. Units: Pa.
end
methods
%% Flow1D Class Constructor
function f = Flow1D(type, phase, id)
f@Domain1D(type, phase, id);
end
%% Flow1D Class Methods
function pressure = get.P(d)
pressure = ctFunc('stflow_pressure', d.domainID);
end
function set.P(d, p)
ctFunc('stflow_setPressure', d.domainID, p);
end
function setFixedTempProfile(d, profile)
% Set a fixed temperature profile. ::
%
% >> d.setFixedTempProfile(profile)
%
% Set the 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.
%
% :param d:
% Instance of class :mat:class:`Domain1D`.
% :param profile:
% :math:`n\times 2` or :math:`2\times n` array of ``n`` points
% at which the temperature is specified.
sz = size(profile);
if sz(1) == 2
l = length(profile(1, :));
ctFunc('stflow_setFixedTempProfile', d.domainID, ...
l, profile(1, :), l, profile(2, :));
elseif sz(2) == 2
l = length(profile(:, 1));
ctFunc('stflow_setFixedTempProfile', d.domainID, ...
l, profile(:, 1), l, profile(:, 2));
else
error('Wrong temperature profile array shape.');
end
end
end
end

28
interfaces/matlab_experimental/OneDim/FreeFlow.m Normal file
View File

@ -0,0 +1,28 @@
classdef FreeFlow < Flow1D
% Create an free flow domain. ::
%
% >> m = FreeFlow(phase, id)
%
% :param phase:
% Instance of class :mat:class:`Solution`.
% :param id:
% String, ID of the flow.
% :return:
% Instance of class :mat:class:`FreeFlow`.
methods
function m = FreeFlow(phase, id)
% Constructor
if nargin < 2
id = 'free-flow';
end
m@Flow1D('free-flow', phase, id);
end
end
end

17
interfaces/matlab_experimental/1D/Inlet.m → interfaces/matlab_experimental/OneDim/Inlet.m
View File

@ -1,4 +1,4 @@
classdef Inlet < Domain1D
classdef Inlet < Boundary1D
% Create an inlet domain. ::
%
% >> m = Inlet(id)
@ -6,6 +6,8 @@ classdef Inlet < Domain1D
% Note that an inlet can only be a terminal domain - it must be
% either the leftmost or rightmost domain in a stack.
%
% :param phase:
% Instance of class :mat:class:`Solution`.
% :param id:
% String name of the inlet.
% :return:
@ -13,17 +15,14 @@ classdef Inlet < Domain1D
methods
function m = Inlet(id)
function m = Inlet(phase, id)
% Constructor
m@Domain1D('Inlet1D');
if nargin == 0
m.setID('inlet');
else
m.setID(id);
if nargin < 2
id = 'inlet';
end
m@Boundary1D('inlet', phase, id);
end
end

15
interfaces/matlab_experimental/1D/Outlet.m → interfaces/matlab_experimental/OneDim/Outlet.m
View File

@ -1,8 +1,10 @@
classdef Outlet < Domain1D
classdef Outlet < Boundary1D
% Create an outlet domain. ::
%
% >> m = Outlet(id)
%
% :param phase:
% Instance of class :mat:class:`Solution`.
% :param id:
% String ID of the outlet.
% :return:
@ -10,17 +12,14 @@ classdef Outlet < Domain1D
methods
function m = Outlet(id)
function m = Outlet(phase, id)
% Constructor
m@Domain1D('Outlet1D');
if nargin == 0
m.setID('outlet');
else
m.setID(id);
if nargin < 2
id = 'outlet';
end
m@Boundary1D('outlet', phase, id);
end
end

16
interfaces/matlab_experimental/1D/OutletRes.m → interfaces/matlab_experimental/OneDim/OutletRes.m
View File

@ -1,8 +1,10 @@
classdef OutletRes < Domain1D
classdef OutletRes < Boundary1D
% Create an outlet reservoir domain. ::
%
% >> m = OutletRes(id)
%
% :param phase:
% Instance of class :mat:class:`Solution`.
% :param id:
% String ID of the outlet reservoir.
% :return:
@ -10,17 +12,15 @@ classdef OutletRes < Domain1D
methods
function m = OutletRes(id)
function m = OutletRes(phase, id)
% Constructor
m@Domain1D('OutletRes');
if nargin == 0
m.setID('outletres');
else
m.setID(id);
if nargin < 2
id = 'outlet-reservoir';
end
m@Boundary1D('outlet-reservoir', phase, id);
end
end

55
interfaces/matlab_experimental/OneDim/ReactingSurface.m Normal file
View File

@ -0,0 +1,55 @@
classdef ReactingSurface < Boundary1D
% Create a reacting surface domain. ::
%
% >> m = Surface(surface_mech, id)
%
% :param surface_mech:
% Instance of class :mat:class:`Interface` defining
% the surface reaction mechanism to be used.
% :param id:
% String ID of the reacting surface.
% :return:
% Instance of class :mat:class:`ReactingSurface`.
properties
% Set bounds on the solution components. ::
%
% >> d.coverageEnabled = flag
%
% :param d:
% Instance of class :mat:class:`Surface`
% :param flag:
% Boolean flag indicating whether coverage equations are enabled.
coverageEnabled
end
methods
%% ReactingSurface Class Constructor
function s = ReactingSurface(surface_mech, id)
if nargin < 2
id = 'reacting-surface';
end
if ~isa(surface_mech, 'Interface')
error('Wrong argument type. Expecting an instance of Interface class.');
end
s@Boundary1D('reacting-surface', surface_mech, id);
ctFunc('reactingsurf_setkineticsmgr', s.domainID, surface_mech.kinID);
s.coverageEnabled = false;
end
%% ReactingSurface Class Methods
function set.coverageEnabled(d, flag)
d.coverageEnabled = flag;
ctFunc('reactingsurf_enableCoverageEqs', d.domainID, int8(flag));
end
end
end

52
interfaces/matlab_experimental/1D/Sim1D.m → interfaces/matlab_experimental/OneDim/Sim1D.m
View File

@ -34,11 +34,6 @@ classdef Sim1D < handle
s.stID = -1;
s.domains = domains;
if nargin ~= 1
help(Sim1D);
error('Wrong number of parameters.');
end
nd = length(domains);
ids = zeros(1, nd);
@ -70,28 +65,6 @@ classdef Sim1D < handle
ctFunc('sim1D_showSolution', s.stID, fname);
end
function plotSolution(s, domain, component)
% Plot a specified solution component. ::
%
% >> s.plotSolution(domain, component)
%
% :param s:
% Instance of class :mat:class:`Sim1D`.
% :param domain:
% Name of domain from which the component should be
% retrieved.
% :param 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 restore(s, fname, id)
% Restore a previously-saved solution. ::
%
@ -111,10 +84,10 @@ classdef Sim1D < handle
ctFunc('sim1D_restore', s.stID, fname, id)
end
function saveSoln(s, fname, id, desc)
function save(s, fname, id, desc)
% Save a solution to a file. ::
%
% >> s.saveSoln(fname, id, desc)
% >> s.save(fname, id, desc)
%
% The output file is in a format that can be used by :mat:class:`restore`
%
@ -127,24 +100,19 @@ classdef Sim1D < handle
% :param desc:
% Description to be written to the output file.
if nargin == 1
fname = 'soln.yaml';
id = 'solution';
desc = '--';
elseif nargin == 2
id = 'solution';
desc = '--';
if nargin < 3
error('Not enough input arguments');
elseif nargin == 3
desc = '--';
desc = '';
end
ctFunc('sim1D_save', s.stID, fname, id, desc);
end
function x = solution(s, domain, component)
function x = getSolution(s, domain, component)
% Get a solution component in one domain. ::
%
% >> s.solution(domain, component)
% >> s.getSolution(domain, component)
%
% :param s:
% Instance of class :mat:class:`Sim1D`.
@ -269,10 +237,10 @@ classdef Sim1D < handle
z = d.gridPoints;
end
function r = resid(s, domain, rdt, count)
function r = residual(s, domain, rdt, count)
% Evaluate the multi-domain residual function. ::
%
% >> r = s.resid(domain, rdt, count)
% >> r = s.residual(domain, rdt, count)
%
% :param s:
% Instance of class :mat:class:`Sim1D`.
@ -500,7 +468,7 @@ classdef Sim1D < handle
% taken first time the steady-state solution attempted.
% If this failed, two time steps would be taken.
ctFunc('sim1D_TimeStep', s.stID, stepsize, length(steps), steps);
ctFunc('sim1D_setTimeStep', s.stID, stepsize, length(steps), steps);
end
function setValue(s, n, comp, localPoints, v)

30
interfaces/matlab_experimental/OneDim/Surface.m Normal file
View File

@ -0,0 +1,30 @@
classdef Surface < Bondary1D
% Create a surface domain. ::
%
% >> m = Surface(id)
%
% :param phase:
% Instance of class :mat:class:`Solution`.
% :param id:
% String ID of surface
% :return:
% Instance of class :mat:class:`Surface` representing a
% non-reacting surface.
methods
%% Surface Class Constructor
function m = Surface(phase, id)
if nargin < 2
id= 'surface';
end
m@Boundary1D('surface', phase, id);
end
end
end

16
interfaces/matlab_experimental/1D/SymmPlane.m → interfaces/matlab_experimental/OneDim/SymmPlane.m
View File

@ -1,8 +1,10 @@
classdef SymmPlane < Domain1D
classdef SymmPlane < Boundary1D
% Create a symmetry plane domain. ::
%
% >> m = SymmPlane(id)
%
% :param phase:
% Instance of class :mat:class:`Solution`.
% :param id:
% String ID of the symmetry plane.
% :return:
@ -10,17 +12,15 @@ classdef SymmPlane < Domain1D
methods
function m = SymmPlane(id)
function m = SymmPlane(phase, id)
% Constructor
m@Domain1D('Symm1D');
if nargin == 0
m.setID('symmetry_plane');
else
m.setID(id);
if nargin < 2
id = 'symmetry-plane';
end
m@Boundary1D('symmetry-plane', phase, id);
end
end

42
samples/matlab_experimental/catcomb.m
View File

@ -106,11 +106,11 @@ flow.setTransientTolerances('default', tol_ts{:});
% specified. This object provides the inlet boundary conditions for
% the flow equations.
inlt = Inlet('inlet');
inlt = Inlet(gas, 'inlet');
% set the inlet parameters. Start with comp1 (hydrogen/air)
inlt.T = tinlet;
inlt.setMassFlowRate(mdot);
inlt.massFlux = mdot;
inlt.setMoleFractions(comp1);
%% create the surface
@ -121,7 +121,7 @@ inlt.setMoleFractions(comp1);
% equation set, and used to compute the surface production rates of
% the gas-phase species.
surf = Surface('surface', surf_phase);
surf = ReactingSurface(surf_phase, 'surface');
surf.T = tsurf;
%% create the stack
@ -144,20 +144,18 @@ for k = 1:gas.nSpecies
stack.setProfile(2, names{k}, [0, 1; y, y]);
end
stack
%setTimeStep(fl, 1.0e-5, [1, 3, 6, 12]);
%setMaxJacAge(fl, 4, 5);
stack.setTimeStep(1.0e-5, [1, 3, 6, 12]);
stack.setMaxJacAge(4, 5);
%% Solution
% start with the energy equation on
flow.enableEnergy;
flow.energyEnabled = true;
% disable the surface coverage equations, and turn off all gas and
% surface chemistry
surf.setCoverageEqs('off');
surf.coverageEnabled = false;
surf_phase.setMultiplier(0.0);
gas.setMultiplier(0.0);
@ -167,7 +165,7 @@ stack.solve(1, refine_grid);
% now turn on the surface coverage equations, and turn the
% chemistry on slowly
surf.setCoverageEqs('on');
surf.coverageEnabled = true;
for iter = 1:6
mult = 10.0^(iter - 6);
@ -186,12 +184,6 @@ stack.setRefineCriteria(2, 100.0, 0.15, 0.2);
% solve the problem for the final time
stack.solve(loglevel, refine_grid);
% show the solution
stack
% save the solution
stack.saveSoln('catcomb.xml', 'energy', ['solution with energy equation']);
%% Show statistics
stack.writeStats;
@ -204,37 +196,37 @@ disp(e);
clf;
subplot(3, 3, 1);
stack.plotSolution('flow', 'T');
plotSolution(stack, 'flow', 'T');
title('Temperature [K]');
subplot(3, 3, 2);
stack.plotSolution('flow', 'velocity');
plotSolution(stack, 'flow', 'velocity');
title('Axial Velocity [m/s]');
subplot(3, 3, 3);
stack.plotSolution('flow', 'spread_rate');
plotSolution(stack, 'flow', 'spread_rate');
title('Radial Velocity / Radius [1/s]');
subplot(3, 3, 4);
stack.plotSolution('flow', 'CH4');
plotSolution(stack, 'flow', 'CH4');
title('CH4 Mass Fraction');
subplot(3, 3, 5);
stack.plotSolution('flow', 'O2');
plotSolution(stack, 'flow', 'O2');
title('O2 Mass Fraction');
subplot(3, 3, 6);
stack.plotSolution('flow', 'CO');
plotSolution(stack, 'flow', 'CO');
title('CO Mass Fraction');
subplot(3, 3, 7);
stack.plotSolution('flow', 'CO2');
plotSolution(stack, 'flow', 'CO2');
title('CO2 Mass Fraction');
subplot(3, 3, 8);
stack.plotSolution('flow', 'H2O');
plotSolution(stack, 'flow', 'H2O');
title('H2O Mass Fraction');
subplot(3, 3, 9);
stack.plotSolution('flow', 'H2');
plotSolution(stack, 'flow', 'H2');
title('H2 Mass Fraction');

17
samples/matlab_experimental/diff_flame.m
View File

@ -22,7 +22,7 @@ p = OneAtm; % Pressure
tin = 300.0; % Inlet temperature
mdot_o = 0.72; % Air mass flux, kg/m^2/s
mdot_f = 0.24; % Fuel mass flux, kg/m^2/s
transport = 'Mix'; % Transport model
transport = 'mixture-averaged'; % Transport model
% NOTE: Transport model needed if mechanism file does not have transport
% properties.
@ -65,15 +65,15 @@ f.setTransientTolerances('default', tol_ts{:});
% Specify the temperature, mass flux, and composition correspondingly.
% Set the oxidizer inlet.
inlet_o = Inlet('air_inlet');
inlet_o = Inlet(ox, 'air_inlet');
inlet_o.T = tin;
inlet_o.setMassFlowRate(mdot_o);
inlet_o.massFlux = mdot_o;
inlet_o.setMoleFractions(oxcomp);
% Set the fuel inlet.
inlet_f = Inlet('fuel_inlet');
inlet_f = Inlet(fuel, 'fuel_inlet');
inlet_f.T = tin;
inlet_f.setMassFlowRate(mdot_f);
inlet_f.massFlux = mdot_f;
inlet_f.setMoleFractions(fuelcomp);
%% Create the flame object.
@ -102,10 +102,9 @@ fl.solve(loglevel, 0);
% on cantera.org in the Matlab User's Guide and can be accessed by
% help setRefineCriteria
f.enableEnergy;
f.energyEnabled = true;
fl.setRefineCriteria(2, 200.0, 0.1, 0.2);
fl.solve(loglevel, refine_grid);
fl.saveSoln('c2h6.xml', 'energy', ['solution with energy equation']);
%% Show statistics of solution and elapsed time.
@ -120,11 +119,11 @@ disp(e);
z = fl.grid('flow'); % Get grid points of flow
spec = fuel.speciesNames; % Get species names in gas
T = fl.solution('flow', 'T'); % Get temperature solution
T = fl.getSolution('flow', 'T'); % Get temperature solution
for i = 1:length(spec)
% Get mass fraction of all species from solution
y(i, :) = fl.solution('flow', spec{i});
y(i, :) = fl.getSolution('flow', spec{i});
end
j = fuel.speciesIndex('O2'); % Get index of O2 in gas object

8
samples/matlab_experimental/flame.m
View File

@ -10,19 +10,19 @@ function f = flame(gas, left, flow, right)
error('gas object must represent an ideal gas mixture.');
end
if ~left.isInlet
if ~isa(left, 'Inlet')
error('burner object of wrong type.');
end
if ~flow.isFlow
if ~isa(flow, 'Flow1D')
error('flow object of wrong type.');
end
flametype = 0;
if right.isSurface
if isa(right, 'ReactingSurface')
flametype = 1;
elseif right.isInlet
elseif isa(right, 'Inlet')
flametype = 3;
end

9
samples/matlab_experimental/flame1.m
View File

@ -59,9 +59,9 @@ f.setTransientTolerances('default', tol_ts{:});
%
% The burner is an Inlet object. The temperature, mass flux,
% and composition (relative molar) may be specified.
burner = Inlet('burner');
burner = Inlet(gas, 'burner');
burner.T = tburner;
burner.setMassFlowRate(mdot);
burner.massFlux = mdot;
burner.setMoleFractions(comp);
%% Create the outlet
@ -71,7 +71,7 @@ burner.setMoleFractions(comp);
% conditions for the temperature and mass fractions, and zero
% radial velocity and radial pressure gradient.
s = Outlet('out');
s = Outlet(gas, 'out');
%% Create the flame object
%
@ -93,10 +93,9 @@ fl.solve(loglevel, refine_grid);
% temperature profile. We also tighten the grid refinement
% criteria to get an accurate final solution.
f.enableEnergy;
f.energyEnabled = true;
fl.setRefineCriteria(2, 200.0, 0.05, 0.1);
fl.solve(1, 1);
fl.saveSoln('h2fl.xml', 'energy', ['solution with energy equation']);
%% Show statistics

11
samples/matlab_experimental/flame2.m
View File

@ -58,16 +58,16 @@ f.setTransientTolerances('default', tol_ts{:});
% The temperature, mass flux, and composition (relative molar) may be
% specified.
inlet_o = Inlet('air_inlet');
inlet_o = Inlet(gas, 'air_inlet');
inlet_o.T = tin;
inlet_o.setMassFlowRate(mdot_o);
inlet_o.massFlux = mdot_o;
inlet_o.setMoleFractions(comp1);
%% Create the fuel inlet
%
inlet_f = Inlet('fuel_inlet');
inlet_f = Inlet(gas, 'fuel_inlet');
inlet_f.T = tin;
inlet_f.setMassFlowRate(mdot_f);
inlet_f.massFlux = mdot_f;
inlet_f.setMoleFractions(comp2);
%% Create the flame object
@ -90,10 +90,9 @@ fl.solve(loglevel, refine_grid);
% temperature profile. We also tighten the grid refinement
% criteria to get an accurate final solution.
f.enableEnergy;
f.energyEnabled = true;
fl.setRefineCriteria(2, 200.0, 0.1, 0.1);
fl.solve(loglevel, refine_grid);
fl.saveSoln('c2h6.xml', 'energy', ['solution with energy equation']);
%% Show statistics

21
samples/matlab_experimental/plotSolution.m Normal file
View File

@ -0,0 +1,21 @@
function plotSolution(s, domain, component)
% Plot a specified solution component. ::
%
% >> plotSolution(s, domain, component)
%
% :param s:
% Instance of class :mat:class:`Sim1D`.
% :param domain:
% Name of domain from which the component should be
% retrieved.
% :param component:
% Name of the component to be plotted.
n = s.stackIndex(domain);
d = s.domains{n};
z = d.gridPoints;
x = s.getSolution(domain, component);
plot(z, x);
xlabel('z (m)');
ylabel(component);
end