From 52727a955026fa432f79fe6f7d8f96fbda21cdb0 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sun, 30 Oct 2016 19:35:05 -0400 Subject: [PATCH] [1D] Eliminate unnecessary state variables from Inlet1D We do not need to solve the trivial equations 'T = Tin' and 'mdot = mdot_in'. --- include/cantera/oneD/Inlet1D.h | 7 -- interfaces/cython/cantera/test/test_onedim.py | 20 ---- src/oneD/OneDim.cpp | 2 +- src/oneD/boundaries1D.cpp | 94 ++++++------------- 4 files changed, 28 insertions(+), 95 deletions(-) diff --git a/include/cantera/oneD/Inlet1D.h b/include/cantera/oneD/Inlet1D.h index f5aa6fa88..0a631665f 100644 --- a/include/cantera/oneD/Inlet1D.h +++ b/include/cantera/oneD/Inlet1D.h @@ -80,10 +80,6 @@ public: return m_mdot; } - virtual void _getInitialSoln(doublereal* x) { - writelog("Bdry1D::_getInitialSoln called!\n"); - } - virtual void setupGrid(size_t n, const doublereal* z) {} protected: @@ -123,8 +119,6 @@ public: virtual void showSolution(const double* x); - virtual void _getInitialSoln(double* x); - virtual size_t nSpecies() { return m_nsp; } @@ -134,7 +128,6 @@ public: virtual doublereal massFraction(size_t k) { return m_yin[k]; } - virtual std::string componentName(size_t n) const; virtual void init(); virtual void eval(size_t jg, doublereal* xg, doublereal* rg, integer* diagg, doublereal rdt); diff --git a/interfaces/cython/cantera/test/test_onedim.py b/interfaces/cython/cantera/test/test_onedim.py index f960fb6bb..117bc0f98 100644 --- a/interfaces/cython/cantera/test/test_onedim.py +++ b/interfaces/cython/cantera/test/test_onedim.py @@ -115,26 +115,6 @@ class TestOnedim(utilities.CanteraTest): self.assertEqual(rtol_ss, set((5e-3, 3e-4, 7e-7))) self.assertEqual(rtol_ts, set((6e-3, 4e-4, 2e-7))) - with self.assertRaises(Exception): - left.set_steady_tolerances(default=(5e-3, 5e-5), - Y=(7e-7, 7e-9)) - - # Boundary domain - left.set_steady_tolerances(default=(5e-3, 5e-5), - temperature=(7e-7, 7e-9)) - left.set_transient_tolerances(default=(6e-3, 6e-5), - temperature=(2e-7, 2e-9)) - - atol_ss = set(left.steady_abstol()) - atol_ts = set(left.transient_abstol()) - rtol_ss = set(left.steady_reltol()) - rtol_ts = set(left.transient_reltol()) - - self.assertEqual(atol_ss, set((5e-5, 7e-9))) - self.assertEqual(atol_ts, set((6e-5, 2e-9))) - self.assertEqual(rtol_ss, set((5e-3, 7e-7))) - self.assertEqual(rtol_ts, set((6e-3, 2e-7))) - class TestFreeFlame(utilities.CanteraTest): tol_ss = [1.0e-5, 1.0e-14] # [rtol atol] for steady-state problem diff --git a/src/oneD/OneDim.cpp b/src/oneD/OneDim.cpp index cb0ae8f99..ad0a7ded4 100644 --- a/src/oneD/OneDim.cpp +++ b/src/oneD/OneDim.cpp @@ -192,7 +192,7 @@ void OneDim::resize() // bandwidth of the local block size_t bw1 = d->bandwidth(); if (bw1 == npos) { - bw1 = 2*d->nComponents() - 1; + bw1 = std::max(2*d->nComponents(), 1) - 1; } m_bw = std::max(m_bw, bw1); diff --git a/src/oneD/boundaries1D.cpp b/src/oneD/boundaries1D.cpp index 93efe40a3..01db99c81 100644 --- a/src/oneD/boundaries1D.cpp +++ b/src/oneD/boundaries1D.cpp @@ -100,12 +100,6 @@ void Inlet1D::showSolution(const double* x) writelog("\n"); } -void Inlet1D::_getInitialSoln(double* x) -{ - x[0] = m_mdot; - x[1] = m_temp; -} - void Inlet1D::setMoleFractions(const std::string& xin) { m_xstr = xin; @@ -125,25 +119,9 @@ void Inlet1D::setMoleFractions(const doublereal* xin) } } -string Inlet1D::componentName(size_t n) const -{ - switch (n) { - case 0: - return "mdot"; - case 1: - return "temperature"; - default: - break; - } - return "unknown"; -} - void Inlet1D::init() { - _init(2); - - setBounds(0, -1e5, 1e5); // mdot - setBounds(1, 200.0, 1e5); // T + _init(0); // if a flow domain is present on the left, then this must be a right inlet. // Note that an inlet object can only be a terminal object - it cannot have @@ -175,26 +153,10 @@ void Inlet1D::eval(size_t jg, doublereal* xg, doublereal* rg, return; } - // start of local part of global arrays - doublereal* x = xg + loc(); - doublereal* r = rg + loc(); - integer* diag = diagg + loc(); - - // residual equations for the two local variables - r[0] = m_mdot - x[0]; - - // Temperature - r[1] = m_temp - x[1]; - - // both are algebraic constraints - diag[0] = 0; - diag[1] = 0; - - // if it is a left inlet, then the flow solution vector - // starts 2 to the right in the global solution vector if (m_ilr == LeftInlet) { - double* xb = x + 2; - double* rb = r + 2; + // Array elements corresponding to the first point of the flow domain + double* xb = xg + m_flow->loc(); + double* rb = rg + m_flow->loc(); // The first flow residual is for u. This, however, is not modified by // the inlet, since this is set within the flow domain from the @@ -206,36 +168,36 @@ void Inlet1D::eval(size_t jg, doublereal* xg, doublereal* rg, // The third flow residual is for T, where it is set to T(0). Subtract // the local temperature to hold the flow T to the inlet T. - rb[2] -= x[1]; + rb[2] -= m_temp; - // The flow domain sets this to -rho*u. Add mdot to specify the mass - // flow rate. - rb[3] += x[0]; + if (m_flow->fixed_mdot()) { + // The flow domain sets this to -rho*u. Add mdot to specify the mass + // flow rate. + rb[3] += m_mdot; + } else { + // if the flow is a freely-propagating flame, mdot is not specified. + // Set mdot equal to rho*u, and also set lambda to zero. + m_mdot = m_flow->density(0)*xb[0]; + rb[3] = xb[3]; + } // add the convective term to the species residual equations for (size_t k = 0; k < m_nsp; k++) { if (k != m_flow_right->leftExcessSpecies()) { - rb[c_offset_Y+k] += x[0]*m_yin[k]; + rb[c_offset_Y+k] += m_mdot*m_yin[k]; } } - // if the flow is a freely-propagating flame, mdot is not specified. - // Set mdot equal to rho*u, and also set lambda to zero. - if (!m_flow->fixed_mdot()) { - m_mdot = m_flow->density(0)*xb[0]; - r[0] = m_mdot - x[0]; - rb[3] = xb[3]; - } } else { - // right inlet. - size_t boffset = m_flow->nComponents(); - double* rb = r - boffset; + // right inlet + // Array elements corresponding to the flast point in the flow domain + double* rb = rg + loc() - m_flow->nComponents(); rb[1] -= m_V0; - rb[2] -= x[1]; // T - rb[0] += x[0]; // u + rb[2] -= m_temp; // T + rb[0] += m_mdot; // u for (size_t k = 0; k < m_nsp; k++) { if (k != m_flow_left->rightExcessSpecies()) { - rb[c_offset_Y+k] += x[0]*m_yin[k]; + rb[c_offset_Y+k] += m_mdot * m_yin[k]; } } } @@ -243,12 +205,10 @@ void Inlet1D::eval(size_t jg, doublereal* xg, doublereal* rg, XML_Node& Inlet1D::save(XML_Node& o, const doublereal* const soln) { - const doublereal* s = soln + loc(); XML_Node& inlt = Domain1D::save(o, soln); inlt.addAttribute("type","inlet"); - for (size_t k = 0; k < nComponents(); k++) { - addFloat(inlt, componentName(k), s[k]); - } + addFloat(inlt, "temperature", m_temp); + addFloat(inlt, "mdot", m_mdot); for (size_t k=0; k < m_nsp; k++) { addFloat(inlt, "massFraction", m_yin[k], "", m_flow->phase().speciesName(k)); @@ -259,8 +219,8 @@ XML_Node& Inlet1D::save(XML_Node& o, const doublereal* const soln) void Inlet1D::restore(const XML_Node& dom, doublereal* soln, int loglevel) { Domain1D::restore(dom, soln, loglevel); - soln[0] = m_mdot = getFloat(dom, "mdot", "massflowrate"); - soln[1] = m_temp = getFloat(dom, "temperature", "temperature"); + m_mdot = getFloat(dom, "mdot"); + m_temp = getFloat(dom, "temperature"); m_yin.assign(m_nsp, 0.0); @@ -273,7 +233,7 @@ void Inlet1D::restore(const XML_Node& dom, doublereal* soln, int loglevel) } } } - resize(2,1); + resize(0, 1); } // ------------- Empty1D -------------