Fix MATLAB interface samples

interfaces/matlab_experimental/Base/Interface.m

@@ -84,7 +84,7 @@
         end
 
         function c = get.concentrations(s)
-            surfID = s.tr_id;
+            surfID = s.tpID;
             nsp = s.nSpecies;
             xx = zeros(1, nsp);
             pt = libpointer('doublePtr', xx);

interfaces/matlab_experimental/Base/ThermoPhase.m

@@ -66,6 +66,8 @@
         %     Scalar double mean molecular weight. Units: kg/kmol.
         meanMolecularWeight
 
+        massDensity % Mass basis density. Units: kg/m^3.
+
         molarDensity % Molar basis density. Units: kmol/m^3.
 
         molecularWeights % Molecular weights of the species. Units: kg/kmol.
@@ -903,6 +905,10 @@ function display(tp)
             mmw = ctFunc('thermo_meanMolecularWeight', tp.tpID);
         end
 
+        function density = get.massDensity(tp)
+            density = ctFunc('thermo_density', tp.tpID);
+        end
+
         function density = get.molarDensity(tp)
             density = ctFunc('thermo_molarDensity', tp.tpID);
         end

interfaces/matlab_experimental/OneDim/Sim1D.m

@@ -58,14 +58,9 @@ function delete(s)
 
         %% Sim1D Utility Methods
 
-        function display(s, fname)
+        function display(s)
             % Show all domains.
-
-            if nargin == 1
-                fname = '-';
-            end
-
-            ctFunc('sim1D_show', s.stID, fname);
+            ctFunc('sim1D_show', s.stID);
         end
 
         function restore(s, fname, id)

interfaces/matlab_experimental/Reactor/ReactorSurface.m

@@ -40,7 +40,7 @@
                 name = '(none)';
             end
 
-            s.surfID = ctFunc('reactor_new', 'ReactorSurface', surf.solnID, name);
+            s@Reactor(surf, 'ReactorSurface', name)
             ctFunc('reactorsurface_install', s.id, reactor.id);
         end
 

samples/matlab_experimental/conhp.m

@@ -6,14 +6,14 @@
     % It assumes that the ``gas`` object represents a reacting ideal gas mixture.
 
     % Set the state of the gas, based on the current solution vector.
+    gas.basis = 'mass';
     gas.Y = y(2:end);
     gas.TP = {y(1), gas.P};
     nsp = gas.nSpecies;
 
     % energy equation
     wdot = gas.netProdRates;
     H = gas.partialMolarEnthalpies';
-    gas.basis = 'mass';
     tdot =- 1 / (gas.D * gas.cp) .* wdot * H;
 
     % set up column vector for dydt

samples/matlab_experimental/conuv.m

@@ -6,14 +6,14 @@
     % It assumes that the ``gas`` object represents a reacting ideal gas mixture.
 
     % Set the state of the gas, based on the current solution vector.
+    gas.basis = 'mass';
     gas.Y = y(2:end);
     gas.TD = {y(1), gas.D};
     nsp = gas.nSpecies;
 
     % energy equation
     wdot = gas.netProdRates;
     U = gas.partialMolarIntEnergies';
-    gas.basis = 'mass';
     tdot =- 1 / (gas.D * gas.cv) .* wdot * U;
 
     % set up column vector for dydt

samples/matlab_experimental/diamond_cvd.m

@@ -69,7 +69,7 @@
     r = surf_phase.netProdRates;
     carbon_dot = r(iC);
     mdot = mw * carbon_dot;
-    rate = mdot / dbulk.D;
+    rate = mdot / dbulk.massDensity;
     xx = [xx; x(ih)];
     rr = [rr; rate * 1.0e6 * 3600.0];
     cov = [cov; surf_phase.coverages];

samples/matlab_experimental/diff_flame.m

@@ -104,7 +104,7 @@
 % ``help setRefineCriteria``.
 
 f.energyEnabled = true;
-fl.setRefineCriteria(2, 200.0, 0.1, 0.2);
+fl.setRefineCriteria(2, 4, 0.2, 0.3, 0.04);
 fl.solve(loglevel, refine_grid);
 
 %% Show statistics of solution and elapsed time

samples/matlab_experimental/equil.m

@@ -20,7 +20,7 @@ function equil(g)
 
     nsp = gas.nSpecies;
 
-    % find methane, nitrogen, and oxygen indices
+    %% find methane, nitrogen, and oxygen indices
     ich4 = gas.speciesIndex('CH4');
     io2 = gas.speciesIndex('O2');
     in2 = gas.speciesIndex('N2');
@@ -41,7 +41,7 @@ function equil(g)
         xeq(:, i) = gas.X;
     end
 
-    % make plots
+    %% make plots
     clf;
     subplot(1, 2, 1);
     plot(phi, tad);

samples/matlab_experimental/flame.m

@@ -32,7 +32,7 @@
     f = Sim1D({left flow right});
 
     % set default initial profiles.
-    rho0 = gas.D;
+    rho0 = gas.massDensity;
 
     % find the adiabatic flame temperature and corresponding
     % equilibrium composition

samples/matlab_experimental/ignite.m

@@ -25,7 +25,7 @@
     gas.TPX = {1001.0, OneAtm, 'H2:2,O2:1,N2:4'};
     gas.basis = 'mass';
     y0 = [gas.U
-          1.0 / gas.D
+          1.0 / gas.massDensity
           gas.Y'];
 
     time_interval = [0 0.001];

samples/matlab_experimental/isentropic.m

@@ -1,5 +1,5 @@
 function isentropic(g)
-    %% Isentropic, adiabatic flow
+    %% Adiabatic, isentropic flow
     %
     % In this example, the area ratio vs. Mach number curve is computed for a
     % hydrogen/nitrogen gas mixture.

samples/matlab_experimental/periodic_cstr.m

@@ -30,7 +30,7 @@
     tic
     help periodic_cstr
 
-    % create the gas mixture
+    %% create the gas mixture
     gas = Solution('h2o2.yaml', 'ohmech');
 
     % pressure = 60 Torr, T = 770 K
@@ -39,26 +39,27 @@
 
     gas.TPX = {300, p, 'H2:2, O2:1'};
 
-    % create an upstream reservoir that will supply the reactor.  The
-    % temperature, pressure, and composition of the upstream reservoir are
+    %%
+    % Create an upstream reservoir that will supply the reactor.  
+    % The temperature, pressure, and composition of the upstream reservoir are
     % set to those of the 'gas' object at the time the reservoir is
     % created.
     upstream = Reservoir(gas);
-
+    %%
     % Now set the gas to the initial temperature of the reactor, and create
     % the reactor object.
     gas.TP = {t, p};
     cstr = IdealGasReactor(gas);
-
+    %%
     % Set its volume to 10 cm^3. In this problem, the reactor volume is
     % fixed, so the initial volume is the volume at all later times.
     cstr.V = 10.0 * 1.0e-6;
-
+    %%
     % We need to have heat loss to see the oscillations. Create a
     % reservoir to represent the environment, and initialize its
     % temperature to the reactor temperature.
     env = Reservoir(gas);
-
+    %%
     % Create a heat-conducting wall between the reactor and the
     % environment. Set its area, and its overall heat transfer
     % coefficient. Larger U causes the reactor to be closer to isothermal.
@@ -67,27 +68,27 @@
     w = Wall(cstr, env);
     w.area = 1.0;
     w.heatTransferCoeff = 0.02;
-
+    %%
     % Connect the upstream reservoir to the reactor with a mass flow
     % controller (constant mdot). Set the mass flow rate to 1.25 sccm.
     sccm = 1.25;
     vdot = sccm * 1.0e-6/60.0 * ((OneAtm / gas.P) * (gas.T / 273.15)); % m^3/s
-    mdot = gas.D * vdot; % kg/s
+    mdot = gas.massDensity * vdot; % kg/s
     mfc = MassFlowController(upstream, cstr);
     mfc.massFlowRate = mdot;
-
+    %%
     % now create a downstream reservoir to exhaust into.
     downstream = Reservoir(gas);
-
+    %%
     % connect the reactor to the downstream reservoir with a valve, and
     % set the coefficient sufficiently large to keep the reactor pressure
     % close to the downstream pressure of 60 Torr.
     v = Valve(cstr, downstream);
     v.valveCoeff = 1.0e-9;
-
+    %%
     % create the network
     network = ReactorNet({cstr});
-
+    %%
     % now integrate in time
     tme = 0.0;
     dt = 0.1;
@@ -104,11 +105,19 @@
         y(3, n) = cstr.massFraction('H2O');
     end
 
+    for i = 2:length(tm)
+        if abs(y(3, i-1) - y(3, i)) > (y(3, i-1) * 5)
+            disp(tm(i))
+        end
+    end
+    %% plot
     clf
     figure(1)
     plot(tm, y)
     legend('H2', 'O2', 'H2O')
     title('Mass Fractions')
+    ylabel('Mass Fractions')
+    xlabel('Time (s)')
 
     toc
 end

samples/matlab_experimental/plug_flow_reactor.m

@@ -99,7 +99,7 @@
 
 T_calc(1) = gas_calc.T;
 Y_calc(1, :) = gas_calc.Y;
-rho_calc(1) = gas_calc.D;
+rho_calc(1) = gas_calc.massDensity;
 
 for i = 2:length(x_calc)
 

samples/matlab_experimental/prandtl1.m

@@ -56,7 +56,7 @@ function prandtl1(g)
 
     disp(['CPU time = ' num2str(cputime - t0)]);
 
-    % plot results
+    %% plot results
 
     clf;
     %figure(1);

samples/matlab_experimental/prandtl2.m

@@ -55,7 +55,7 @@ function prandtl2(g)
 
     disp(['CPU time = ' num2str(cputime - t0)]);
 
-    % plot results
+    %% plot results
 
     clf;
     subplot(2, 2, 1);

samples/matlab_experimental/rankine.m

@@ -10,35 +10,39 @@
 
 help rankine
 
-% Initialize parameters
+%% Initialize parameters
 eta_pump = 0.6;
 eta_turbine = 0.8;
 p_max = 8.0 * OneAtm;
 t1 = 300.0;
 
-% create an object representing water
+%% create an object representing water
 w = Water;
 
+%%
 % start with saturated liquid water at t1
 basis = 'mass';
 w.TQ = {t1, 0};
 h1 = w.H;
 p1 = w.P;
 
+%%
 % pump it to p2
 pump_work = pump(w, p_max, eta_pump);
 h2 = w.H;
 
+%%
 % heat to saturated vapor
 w.PQ = {p_max, 1.0};
 h3 = w.H;
 
 heat_added = h3 - h2;
 
+%%
 % expand adiabatically back to the initial pressure
 turbine_work = expand(w, p1, eta_turbine);
 
-% compute the efficiency
+%% compute the efficiency
 efficiency = (turbine_work - pump_work) / heat_added;
 disp(sprintf('efficiency = %.2f%%', efficiency*100));
 

samples/matlab_experimental/reactor1.m

@@ -20,7 +20,7 @@ function reactor1(g)
     end
 
     P = 101325.0;
-    % set the initial conditions
+    %% set the initial conditions
     gas.TP = {1001.0, P};
     nsp = gas.nSpecies;
     xx = zeros(nsp, 1);
@@ -29,26 +29,29 @@ function reactor1(g)
     xx(48) = 0.573;
     gas.X = xx;
 
-    % create a reactor, and insert the gas
+    %% create a reactor, and insert the gas
     r = IdealGasReactor(gas);
 
-    % create a reservoir to represent the environment
+    %% create a reservoir to represent the environment
     a = Solution('air.yaml', 'air', 'none');
     a.TP = {a.T, P};
     env = Reservoir(a);
 
+    %%
     % Define a wall between the reactor and the environment and
     % make it flexible, so that the pressure in the reactor is held
     % at the environment pressure.
     w = Wall(r, env);
 
+    %%
     % set expansion parameter. dV/dt = KA(P_1 - P_2)
     w.expansionRateCoeff = 1.0e6;
 
+    %%
     % set wall area
     w.area = 1.0;
 
-    % create a reactor network and insert the reactor:
+    %% create a reactor network and insert the reactor:
     network = ReactorNet({r});
 
     nSteps = 100;

samples/matlab_experimental/reactor2.m

@@ -19,13 +19,13 @@ function reactor2(g)
         gas = Solution('gri30.yaml', 'gri30', 'none');
     end
 
-    % set the initial conditions
+    %% set the initial conditions
     gas.TPX = {1001.0, OneAtm, 'H2:2,O2:1,N2:4'};
 
-    % create a reactor, and insert the gas
+    %% create a reactor, and insert the gas
     r = IdealGasReactor(gas);
 
-    % create a reactor network and insert the reactor
+    %% create a reactor network and insert the reactor
     network = ReactorNet({r});
 
     nSteps = 100;