From cd958a343ba57188bc92593ef2bc23f98a53ab3d Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Tue, 28 Jun 2016 17:35:40 -0400 Subject: [PATCH] [1D] Make better use of local variables --- src/oneD/Domain1D.cpp | 33 ++++++------ src/oneD/MultiJac.cpp | 22 ++++---- src/oneD/MultiNewton.cpp | 43 +++++++--------- src/oneD/Sim1D.cpp | 42 +++++++--------- src/oneD/StFlow.cpp | 103 ++++++++++++++++---------------------- src/oneD/boundaries1D.cpp | 95 +++++++++++++++-------------------- 6 files changed, 145 insertions(+), 193 deletions(-) diff --git a/src/oneD/Domain1D.cpp b/src/oneD/Domain1D.cpp index e1d620474..f41a0ec68 100644 --- a/src/oneD/Domain1D.cpp +++ b/src/oneD/Domain1D.cpp @@ -62,8 +62,7 @@ std::string Domain1D::componentName(size_t n) const size_t Domain1D::componentIndex(const std::string& name) const { - size_t nc = nComponents(); - for (size_t n = 0; n < nc; n++) { + for (size_t n = 0; n < nComponents(); n++) { if (name == componentName(n)) { return n; } @@ -127,8 +126,8 @@ void Domain1D::eval(size_t jg, doublereal* xg, doublereal* rg, doublereal* rsd = rg + loc(); integer* diag = mask + loc(); - size_t jmin, jmax, jpt, j, i; - jpt = jg - firstPoint(); + size_t jmin, jmax; + size_t jpt = jg - firstPoint(); if (jg == npos) { // evaluate all points jmin = 0; @@ -138,14 +137,14 @@ void Domain1D::eval(size_t jg, doublereal* xg, doublereal* rg, jmax = std::min(jpt+1,m_points-1); } - for (j = jmin; j <= jmax; j++) { + for (size_t j = jmin; j <= jmax; j++) { if (j == 0 || j == m_points - 1) { - for (i = 0; i < m_nv; i++) { + for (size_t i = 0; i < m_nv; i++) { rsd[index(i,j)] = residual(x,i,j); diag[index(i,j)] = 0; } } else { - for (i = 0; i < m_nv; i++) { + for (size_t i = 0; i < m_nv; i++) { rsd[index(i,j)] = residual(x,i,j) - timeDerivativeFlag(i)*rdt*(value(x,i,j) - prevSoln(i,j)); diag[index(i,j)] = timeDerivativeFlag(i); @@ -241,19 +240,17 @@ void Domain1D::setupGrid(size_t n, const doublereal* z) void Domain1D::showSolution(const doublereal* x) { size_t nn = m_nv/5; - size_t i, j, n; - doublereal v; - for (i = 0; i < nn; i++) { + for (size_t i = 0; i < nn; i++) { writeline('-', 79, false, true); writelog("\n z "); - for (n = 0; n < 5; n++) { + for (size_t n = 0; n < 5; n++) { writelog(" {:>10s} ", componentName(i*5 + n)); } writeline('-', 79, false, true); - for (j = 0; j < m_points; j++) { + for (size_t j = 0; j < m_points; j++) { writelog("\n {:10.4g} ", m_z[j]); - for (n = 0; n < 5; n++) { - v = value(x, i*5+n, j); + for (size_t n = 0; n < 5; n++) { + double v = value(x, i*5+n, j); writelog(" {:10.4g} ", v); } } @@ -262,14 +259,14 @@ void Domain1D::showSolution(const doublereal* x) size_t nrem = m_nv - 5*nn; writeline('-', 79, false, true); writelog("\n z "); - for (n = 0; n < nrem; n++) { + for (size_t n = 0; n < nrem; n++) { writelog(" {:>10s} ", componentName(nn*5 + n)); } writeline('-', 79, false, true); - for (j = 0; j < m_points; j++) { + for (size_t j = 0; j < m_points; j++) { writelog("\n {:10.4g} ", m_z[j]); - for (n = 0; n < nrem; n++) { - v = value(x, nn*5+n, j); + for (size_t n = 0; n < nrem; n++) { + double v = value(x, nn*5+n, j); writelog(" {:10.4g} ", v); } } diff --git a/src/oneD/MultiJac.cpp b/src/oneD/MultiJac.cpp index 9775bdd50..5faa2bd57 100644 --- a/src/oneD/MultiJac.cpp +++ b/src/oneD/MultiJac.cpp @@ -46,14 +46,14 @@ void MultiJac::eval(doublereal* x0, doublereal* resid0, doublereal rdt) m_nevals++; clock_t t0 = clock(); bfill(0.0); - size_t n, m, ipt=0, j, nv, mv, iloc; - doublereal rdx, dx, xsave; + size_t ipt=0; - for (j = 0; j < m_points; j++) { - nv = m_resid->nVars(j); - for (n = 0; n < nv; n++) { + for (size_t j = 0; j < m_points; j++) { + size_t nv = m_resid->nVars(j); + for (size_t n = 0; n < nv; n++) { // perturb x(n); preserve sign(x(n)) - xsave = x0[ipt]; + double xsave = x0[ipt]; + double dx; if (xsave >= 0) { dx = xsave*m_rtol + m_atol; } else { @@ -61,7 +61,7 @@ void MultiJac::eval(doublereal* x0, doublereal* resid0, doublereal rdt) } x0[ipt] = xsave + dx; dx = x0[ipt] - xsave; - rdx = 1.0/dx; + double rdx = 1.0/dx; // calculate perturbed residual m_resid->eval(j, x0, m_r1.data(), rdt, 0); @@ -69,9 +69,9 @@ void MultiJac::eval(doublereal* x0, doublereal* resid0, doublereal rdt) // compute nth column of Jacobian for (size_t i = j - 1; i != j+2; i++) { if (i != npos && i < m_points) { - mv = m_resid->nVars(i); - iloc = m_resid->loc(i); - for (m = 0; m < mv; m++) { + size_t mv = m_resid->nVars(i); + size_t iloc = m_resid->loc(i); + for (size_t m = 0; m < mv; m++) { value(m+iloc,ipt) = (m_r1[m+iloc] - resid0[m+iloc])*rdx; } } @@ -81,7 +81,7 @@ void MultiJac::eval(doublereal* x0, doublereal* resid0, doublereal rdt) } } - for (n = 0; n < m_size; n++) { + for (size_t n = 0; n < m_size; n++) { m_ssdiag[n] = value(n,n); } diff --git a/src/oneD/MultiNewton.cpp b/src/oneD/MultiNewton.cpp index 00efb9491..3f3e75a49 100644 --- a/src/oneD/MultiNewton.cpp +++ b/src/oneD/MultiNewton.cpp @@ -39,23 +39,21 @@ doublereal bound_step(const doublereal* x, const doublereal* step, size_t np = r.nPoints(); size_t nv = r.nComponents(); Indx index(nv, np); - doublereal above, below, val, newval; - size_t m, j; doublereal fbound = 1.0; bool wroteTitle = false; - for (m = 0; m < nv; m++) { - above = r.upperBound(m); - below = r.lowerBound(m); + for (size_t m = 0; m < nv; m++) { + double above = r.upperBound(m); + double below = r.lowerBound(m); - for (j = 0; j < np; j++) { - val = x[index(m,j)]; + for (size_t j = 0; j < np; j++) { + double val = x[index(m,j)]; if (loglevel > 0 && (val > above + 1.0e-12 || val < below - 1.0e-12)) { writelog("\nERROR: solution out of bounds.\n"); writelog("domain {:d}: {:>20s}({:d}) = {:10.3e} ({:10.3e}, {:10.3e})\n", r.domainIndex(), r.componentName(m), j, val, below, above); } - newval = val + step[index(m,j)]; + double newval = val + step[index(m,j)]; if (newval > above) { fbound = std::max(0.0, std::min(fbound, @@ -105,20 +103,19 @@ doublereal bound_step(const doublereal* x, const doublereal* step, doublereal norm_square(const doublereal* x, const doublereal* step, Domain1D& r) { - doublereal f, ewt, esum, sum = 0.0; - size_t n, j; + double sum = 0.0; doublereal f2max = 0.0; size_t nv = r.nComponents(); size_t np = r.nPoints(); - for (n = 0; n < nv; n++) { - esum = 0.0; - for (j = 0; j < np; j++) { + for (size_t n = 0; n < nv; n++) { + double esum = 0.0; + for (size_t j = 0; j < np; j++) { esum += fabs(x[nv*j + n]); } - ewt = r.rtol(n)*esum/np + r.atol(n); - for (j = 0; j < np; j++) { - f = step[nv*j + n]/ewt; + double ewt = r.rtol(n)*esum/np + r.atol(n); + for (size_t j = 0; j < np; j++) { + double f = step[nv*j + n]/ewt; sum += f*f; f2max = std::max(f*f, f2max); } @@ -153,11 +150,10 @@ void MultiNewton::resize(size_t sz) doublereal MultiNewton::norm2(const doublereal* x, const doublereal* step, OneDim& r) const { - doublereal f, sum = 0.0; + double sum = 0.0; size_t nd = r.nDomains(); for (size_t n = 0; n < nd; n++) { - f = norm_square(x + r.start(n), step + r.start(n), - r.domain(n)); + double f = norm_square(x + r.start(n), step + r.start(n), r.domain(n)); sum += f; } sum /= r.size(); @@ -167,14 +163,12 @@ doublereal MultiNewton::norm2(const doublereal* x, void MultiNewton::step(doublereal* x, doublereal* step, OneDim& r, MultiJac& jac, int loglevel) { - size_t iok; - size_t sz = r.size(); r.eval(npos, x, step); - for (size_t n = 0; n < sz; n++) { + for (size_t n = 0; n < r.size(); n++) { step[n] = -step[n]; } - iok = jac.solve(step, step); + size_t iok = jac.solve(step, step); // if iok is non-zero, then solve failed if (iok != 0) { iok--; @@ -243,10 +237,9 @@ int MultiNewton::dampStep(const doublereal* x0, const doublereal* step0, // damping coefficient starts at 1.0 doublereal damp = 1.0; - doublereal ff; size_t m; for (m = 0; m < NDAMP; m++) { - ff = fbound*damp; + double ff = fbound*damp; // step the solution by the damped step size for (size_t j = 0; j < m_n; j++) { diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp index 8b0748a8b..34b0dc99a 100644 --- a/src/oneD/Sim1D.cpp +++ b/src/oneD/Sim1D.cpp @@ -79,11 +79,10 @@ void Sim1D::setProfile(size_t dom, size_t comp, Domain1D& d = domain(dom); doublereal z0 = d.zmin(); doublereal z1 = d.zmax(); - doublereal zpt, frac, v; for (size_t n = 0; n < d.nPoints(); n++) { - zpt = d.z(n); - frac = (zpt - z0)/(z1 - z0); - v = linearInterp(frac, pos, values); + double zpt = d.z(n); + double frac = (zpt - z0)/(z1 - z0); + double v = linearInterp(frac, pos, values); setValue(dom, comp, n, v); } } @@ -138,8 +137,7 @@ void Sim1D::restore(const std::string& fname, const std::string& id, void Sim1D::setFlatProfile(size_t dom, size_t comp, doublereal v) { size_t np = domain(dom).nPoints(); - size_t n; - for (n = 0; n < np; n++) { + for (size_t n = 0; n < np; n++) { setValue(dom, comp, n, v); } } @@ -226,7 +224,6 @@ int Sim1D::newtonSolve(int loglevel) void Sim1D::solve(int loglevel, bool refine_grid) { int new_points = 1; - int nsteps; doublereal dt = m_tstep; m_nsteps = 0; int soln_number = -1; @@ -234,7 +231,7 @@ void Sim1D::solve(int loglevel, bool refine_grid) while (new_points > 0) { size_t istep = 0; - nsteps = m_steps[istep]; + size_t nsteps = m_steps[istep]; bool ok = false; if (loglevel > 0) { @@ -346,7 +343,6 @@ int Sim1D::refine(int loglevel) { int ianalyze, np = 0; vector_fp znew, xnew; - doublereal xmid, zmid; std::vector dsize; m_xlast_ss = m_x; @@ -390,14 +386,14 @@ int Sim1D::refine(int loglevel) // for this new point if (r.newPointNeeded(m) && m + 1 < npnow) { // add new point at midpoint - zmid = 0.5*(d.grid(m) + d.grid(m+1)); + double zmid = 0.5*(d.grid(m) + d.grid(m+1)); znew.push_back(zmid); np++; // for each component, linearly interpolate // the solution to this point for (size_t i = 0; i < comp; i++) { - xmid = 0.5*(value(n, i, m) + value(n, i, m+1)); + double xmid = 0.5*(value(n, i, m) + value(n, i, m+1)); xnew.push_back(xmid); } } @@ -433,14 +429,12 @@ int Sim1D::setFixedTemperature(doublereal t) { int np = 0; vector_fp znew, xnew; - doublereal xmid; - doublereal zfixed,interp_factor; + doublereal zfixed; doublereal z1 = 0.0, z2 = 0.0, t1,t2; - size_t n, m, i; size_t m1 = 0; std::vector dsize; - for (n = 0; n < nDomains(); n++) { + for (size_t n = 0; n < nDomains(); n++) { bool addnewpt=false; Domain1D& d = domain(n); size_t comp = d.nComponents(); @@ -451,7 +445,7 @@ int Sim1D::setFixedTemperature(doublereal t) size_t npnow = d.nPoints(); size_t nstart = znew.size(); if (d_free) { - for (m = 0; m < npnow-1; m++) { + for (size_t m = 0; m < npnow-1; m++) { if (value(n,2,m) == t) { zfixed = d.grid(m); d_free->m_zfixed = zfixed; @@ -475,23 +469,23 @@ int Sim1D::setFixedTemperature(doublereal t) } } - for (m = 0; m < npnow; m++) { + for (size_t m = 0; m < npnow; m++) { // add the current grid point to the new grid znew.push_back(d.grid(m)); // do the same for the solution at this point - for (i = 0; i < comp; i++) { + for (size_t i = 0; i < comp; i++) { xnew.push_back(value(n, i, m)); } if (m==m1 && addnewpt) { //add new point at zfixed znew.push_back(zfixed); np++; - interp_factor = (zfixed-z2) / (z1-z2); + double interp_factor = (zfixed-z2) / (z1-z2); // for each component, linearly interpolate // the solution to this point - for (i = 0; i < comp; i++) { - xmid = interp_factor*(value(n, i, m) - value(n, i, m+1)) + value(n,i,m+1); + for (size_t i = 0; i < comp; i++) { + double xmid = interp_factor*(value(n, i, m) - value(n, i, m+1)) + value(n,i,m+1); xnew.push_back(xmid); } } @@ -502,10 +496,10 @@ int Sim1D::setFixedTemperature(doublereal t) // At this point, the new grid znew and the new solution vector xnew have // been constructed, but the domains themselves have not yet been modified. // Now update each domain with the new grid. - size_t gridstart = 0, gridsize; - for (n = 0; n < nDomains(); n++) { + size_t gridstart = 0; + for (size_t n = 0; n < nDomains(); n++) { Domain1D& d = domain(n); - gridsize = dsize[n]; + size_t gridsize = dsize[n]; d.setupGrid(gridsize, &znew[gridstart]); gridstart += gridsize; } diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp index 366f26f1a..7ca6bec46 100644 --- a/src/oneD/StFlow.cpp +++ b/src/oneD/StFlow.cpp @@ -123,10 +123,9 @@ void StFlow::resize(size_t ncomponents, size_t points) void StFlow::setupGrid(size_t n, const doublereal* z) { resize(m_nv, n); - size_t j; m_z[0] = z[0]; - for (j = 1; j < m_points; j++) { + for (size_t j = 1; j < m_points; j++) { if (z[j] <= z[j-1]) { throw CanteraError("StFlow::setupGrid", "grid points must be monotonically increasing"); @@ -211,16 +210,14 @@ void StFlow::setGasAtMidpoint(const doublereal* x, size_t j) void StFlow::_finalize(const doublereal* x) { - size_t j; - doublereal zz, tt; size_t nz = m_zfix.size(); bool e = m_do_energy[0]; - for (j = 0; j < m_points; j++) { + for (size_t j = 0; j < m_points; j++) { if (e || nz == 0) { m_fixedtemp[j] = T(x, j); } else { - zz = (z(j) - z(0))/(z(m_points - 1) - z(0)); - tt = linearInterp(zz, m_zfix, m_tfix); + double zz = (z(j) - z(0))/(z(m_points - 1) - z(0)); + double tt = linearInterp(zz, m_zfix, m_tfix); m_fixedtemp[j] = tt; } } @@ -262,8 +259,6 @@ void StFlow::eval(size_t jg, doublereal* xg, size_t j0 = std::max(jmin, 1) - 1; size_t j1 = std::min(jmax+1,m_points-1); - size_t j, k; - // ------------ update properties ------------ updateThermo(x, j0, j1); @@ -285,7 +280,6 @@ void StFlow::eval(size_t jg, doublereal* xg, // evaluate the residual equations at all required // grid points //---------------------------------------------------- - doublereal sum, sum2, dtdzj; // calculation of qdotRadiation @@ -352,7 +346,7 @@ void StFlow::eval(size_t jg, doublereal* xg, } } - for (j = jmin; j <= jmax; j++) { + for (size_t j = jmin; j <= jmax; j++) { //---------------------------------------------- // left boundary //---------------------------------------------- @@ -377,8 +371,8 @@ void StFlow::eval(size_t jg, doublereal* xg, // The default boundary condition for species is zero flux. However, // the boundary object may modify this. - sum = 0.0; - for (k = 0; k < m_nsp; k++) { + double sum = 0.0; + for (size_t k = 0; k < m_nsp; k++) { sum += Y(x,k,0); rsd[index(c_offset_Y + k, 0)] = -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0)); @@ -408,11 +402,10 @@ void StFlow::eval(size_t jg, doublereal* xg, // = M_k\omega_k //------------------------------------------------- getWdot(x,j); - doublereal convec, diffus; - for (k = 0; k < m_nsp; k++) { - convec = rho_u(x,j)*dYdz(x,k,j); - diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1)) - /(z(j+1) - z(j-1)); + for (size_t k = 0; k < m_nsp; k++) { + double convec = rho_u(x,j)*dYdz(x,k,j); + double diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1)) + / (z(j+1) - z(j-1)); rsd[index(c_offset_Y + k, j)] = (m_wt[k]*(wdot(k,j)) - convec - diffus)/m_rho[j] @@ -434,16 +427,15 @@ void StFlow::eval(size_t jg, doublereal* xg, // heat release term const vector_fp& h_RT = m_thermo->enthalpy_RT_ref(); const vector_fp& cp_R = m_thermo->cp_R_ref(); - sum = 0.0; - sum2 = 0.0; - doublereal flxk; - for (k = 0; k < m_nsp; k++) { - flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j)); + double sum = 0.0; + double sum2 = 0.0; + for (size_t k = 0; k < m_nsp; k++) { + double flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j)); sum += wdot(k,j)*h_RT[k]; sum2 += flxk*cp_R[k]/m_wt[k]; } sum *= GasConstant * T(x,j); - dtdzj = dTdz(x,j); + double dtdzj = dTdz(x,j); sum2 *= GasConstant * dtdzj; rsd[index(c_offset_T, j)] = - m_cp[j]*rho_u(x,j)*dtdzj @@ -513,32 +505,29 @@ void StFlow::showSolution(const doublereal* x) void StFlow::updateDiffFluxes(const doublereal* x, size_t j0, size_t j1) { - size_t j, k, m; - doublereal sum, wtm, rho, dz, gradlogT; - switch (m_transport_option) { case c_Mixav_Transport: - for (j = j0; j < j1; j++) { - sum = 0.0; - wtm = m_wtm[j]; - rho = density(j); - dz = z(j+1) - z(j); - for (k = 0; k < m_nsp; k++) { + for (size_t j = j0; j < j1; j++) { + double sum = 0.0; + double wtm = m_wtm[j]; + double rho = density(j); + double dz = z(j+1) - z(j); + for (double k = 0; k < m_nsp; k++) { m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm); m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz; sum -= m_flux(k,j); } // correction flux to insure that \sum_k Y_k V_k = 0. - for (k = 0; k < m_nsp; k++) { + for (double k = 0; k < m_nsp; k++) { m_flux(k,j) += sum*Y(x,k,j); } } break; case c_Multi_Transport: - for (j = j0; j < j1; j++) { - dz = z(j+1) - z(j); - for (k = 0; k < m_nsp; k++) { + for (size_t j = j0; j < j1; j++) { + double dz = z(j+1) - z(j); + for (double k = 0; k < m_nsp; k++) { doublereal sum = 0.0; for (size_t m = 0; m < m_nsp; m++) { sum += m_wt[m] * m_multidiff[mindex(k,m,j)] * (X(x,m,j+1)-X(x,m,j)); @@ -553,10 +542,10 @@ void StFlow::updateDiffFluxes(const doublereal* x, size_t j0, size_t j1) } if (m_do_soret) { - for (m = j0; m < j1; m++) { - gradlogT = 2.0 * (T(x,m+1) - T(x,m)) / - ((T(x,m+1) + T(x,m)) * (z(m+1) - z(m))); - for (k = 0; k < m_nsp; k++) { + for (size_t m = j0; m < j1; m++) { + double gradlogT = 2.0 * (T(x,m+1) - T(x,m)) / + ((T(x,m+1) + T(x,m)) * (z(m+1) - z(m))); + for (size_t k = 0; k < m_nsp; k++) { m_flux(k,m) -= m_dthermal(k,m)*gradlogT; } } @@ -616,18 +605,15 @@ void StFlow::restore(const XML_Node& dom, doublereal* soln, int loglevel) writelog(nd["title"]+": "+nd.value()+"\n"); } - double pp = -1.0; - pp = getFloat(dom, "pressure", "pressure"); + double pp = getFloat(dom, "pressure", "pressure"); setPressure(pp); vector d = dom.child("grid_data").getChildren("floatArray"); - size_t nd = d.size(); vector_fp x; - size_t n, np = 0, j, ks, k; - string nm; + size_t np = 0; bool readgrid = false, wrote_header = false; - for (n = 0; n < nd; n++) { + for (size_t n = 0; n < d.size(); n++) { const XML_Node& fa = *d[n]; - nm = fa["title"]; + string nm = fa["title"]; if (nm == "z") { getFloatArray(fa,x,false); np = x.size(); @@ -644,9 +630,9 @@ void StFlow::restore(const XML_Node& dom, doublereal* soln, int loglevel) } debuglog("Importing datasets:\n", loglevel >= 2); - for (n = 0; n < nd; n++) { + for (size_t n = 0; n < d.size(); n++) { const XML_Node& fa = *d[n]; - nm = fa["title"]; + string nm = fa["title"]; getFloatArray(fa,x,false); if (nm == "u") { debuglog("axial velocity ", loglevel >= 2); @@ -654,7 +640,7 @@ void StFlow::restore(const XML_Node& dom, doublereal* soln, int loglevel) throw CanteraError("StFlow::restore", "axial velocity array size error"); } - for (j = 0; j < np; j++) { + for (size_t j = 0; j < np; j++) { soln[index(0,j)] = x[j]; } } else if (nm == "z") { @@ -665,7 +651,7 @@ void StFlow::restore(const XML_Node& dom, doublereal* soln, int loglevel) throw CanteraError("StFlow::restore", "radial velocity array size error"); } - for (j = 0; j < np; j++) { + for (size_t j = 0; j < np; j++) { soln[index(1,j)] = x[j]; } } else if (nm == "T") { @@ -674,7 +660,7 @@ void StFlow::restore(const XML_Node& dom, doublereal* soln, int loglevel) throw CanteraError("StFlow::restore", "temperature array size error"); } - for (j = 0; j < np; j++) { + for (size_t j = 0; j < np; j++) { soln[index(2,j)] = x[j]; } @@ -692,15 +678,15 @@ void StFlow::restore(const XML_Node& dom, doublereal* soln, int loglevel) throw CanteraError("StFlow::restore", "lambda arary size error"); } - for (j = 0; j < np; j++) { + for (size_t j = 0; j < np; j++) { soln[index(3,j)] = x[j]; } } else if (m_thermo->speciesIndex(nm) != npos) { debuglog(nm+" ", loglevel >= 2); if (x.size() == np) { - k = m_thermo->speciesIndex(nm); + size_t k = m_thermo->speciesIndex(nm); did_species[k] = 1; - for (j = 0; j < np; j++) { + for (size_t j = 0; j < np; j++) { soln[index(k+4,j)] = x[j]; } } @@ -719,7 +705,7 @@ void StFlow::restore(const XML_Node& dom, doublereal* soln, int loglevel) } if (loglevel >= 1) { - for (ks = 0; ks < nsp; ks++) { + for (size_t ks = 0; ks < nsp; ks++) { if (did_species[ks] == 0) { if (!wrote_header) { writelog("Missing data for species:\n"); @@ -770,7 +756,6 @@ void StFlow::restore(const XML_Node& dom, doublereal* soln, int loglevel) XML_Node& StFlow::save(XML_Node& o, const doublereal* const sol) { - size_t k; Array2D soln(m_nv, m_points, sol + loc()); XML_Node& flow = Domain1D::save(o, sol); flow.addAttribute("type",flowType()); @@ -797,7 +782,7 @@ XML_Node& StFlow::save(XML_Node& o, const doublereal* const sol) soln.getRow(3, x.data()); addFloatArray(gv,"L",x.size(),x.data(),"N/m^4"); - for (k = 0; k < m_nsp; k++) { + for (size_t k = 0; k < m_nsp; k++) { soln.getRow(4+k, x.data()); addFloatArray(gv,m_thermo->speciesName(k), x.size(),x.data(),"","massFraction"); diff --git a/src/oneD/boundaries1D.cpp b/src/oneD/boundaries1D.cpp index 630164be5..d62d8d5a4 100644 --- a/src/oneD/boundaries1D.cpp +++ b/src/oneD/boundaries1D.cpp @@ -177,7 +177,6 @@ void Inlet1D::eval(size_t jg, doublereal* xg, doublereal* rg, doublereal* x = xg + loc(); doublereal* r = rg + loc(); integer* diag = diagg + loc(); - doublereal* xb, *rb; // residual equations for the two local variables r[0] = m_mdot - x[0]; @@ -192,8 +191,8 @@ void Inlet1D::eval(size_t jg, doublereal* xg, doublereal* rg, // 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) { - xb = x + 2; - rb = r + 2; + double* xb = x + 2; + double* rb = r + 2; // 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 @@ -228,8 +227,7 @@ void Inlet1D::eval(size_t jg, doublereal* xg, doublereal* rg, } else { // right inlet. size_t boffset = m_flow->nComponents(); - xb = x - boffset; - rb = r - boffset; + double* rb = r - boffset; rb[1] -= m_V0; rb[2] -= x[1]; // T rb[0] += x[0]; // u @@ -353,18 +351,15 @@ void Symm1D::eval(size_t jg, doublereal* xg, doublereal* rg, integer* diagg, doublereal* x = xg + loc(); doublereal* r = rg + loc(); integer* diag = diagg + loc(); - doublereal* xb, *rb; - integer* db; r[0] = x[0]; diag[0] = 0; - size_t nc; if (m_flow_right) { - nc = m_flow_right->nComponents(); - xb = x + 1; - rb = r + 1; - db = diag + 1; + size_t nc = m_flow_right->nComponents(); + double* xb = x + 1; + double* rb = r + 1; + int* db = diag + 1; db[1] = 0; db[2] = 0; rb[1] = xb[1] - xb[1 + nc]; // zero dV/dz @@ -372,10 +367,10 @@ void Symm1D::eval(size_t jg, doublereal* xg, doublereal* rg, integer* diagg, } if (m_flow_left) { - nc = m_flow_left->nComponents(); - xb = x - nc; - rb = r - nc; - db = diag - nc; + size_t nc = m_flow_left->nComponents(); + double* xb = x - nc; + double* rb = r - nc; + int* db = diag - nc; db[1] = 0; db[2] = 0; rb[1] = xb[1] - xb[1 - nc]; // zero dV/dz @@ -441,30 +436,26 @@ void Outlet1D::eval(size_t jg, doublereal* xg, doublereal* rg, integer* diagg, doublereal* x = xg + loc(); doublereal* r = rg + loc(); integer* diag = diagg + loc(); - doublereal* xb, *rb; - integer* db; r[0] = x[0]; diag[0] = 0; - size_t nc, k; if (m_flow_right) { - nc = m_flow_right->nComponents(); - xb = x + 1; - rb = r + 1; - db = diag + 1; + size_t nc = m_flow_right->nComponents(); + double* xb = x + 1; + double* rb = r + 1; rb[0] = xb[3]; rb[2] = xb[2] - xb[2 + nc]; - for (k = 4; k < nc; k++) { + for (size_t k = 4; k < nc; k++) { rb[k] = xb[k] - xb[k + nc]; } } if (m_flow_left) { - nc = m_flow_left->nComponents(); - xb = x - nc; - rb = r - nc; - db = diag - nc; + size_t nc = m_flow_left->nComponents(); + double* xb = x - nc; + double* rb = r - nc; + int* db = diag - nc; // zero Lambda if (m_flow_left->fixed_mdot()) { @@ -473,7 +464,7 @@ void Outlet1D::eval(size_t jg, doublereal* xg, doublereal* rg, integer* diagg, rb[2] = xb[2] - xb[2 - nc]; // zero T gradient size_t kSkip = 4 + m_flow_left->rightExcessSpecies(); - for (k = 4; k < nc; k++) { + for (size_t k = 4; k < nc; k++) { if (k != kSkip) { rb[k] = xb[k] - xb[k - nc]; // zero mass fraction gradient db[k] = 0; @@ -561,19 +552,15 @@ void OutletRes1D::eval(size_t jg, doublereal* xg, doublereal* rg, doublereal* x = xg + loc(); doublereal* r = rg + loc(); integer* diag = diagg + loc(); - doublereal* xb, *rb; - integer* db; // drive dummy component to zero r[0] = x[0]; diag[0] = 0; - size_t nc, k; if (m_flow_right) { - nc = m_flow_right->nComponents(); - xb = x + 1; - rb = r + 1; - db = diag + 1; + size_t nc = m_flow_right->nComponents(); + double* xb = x + 1; + double* rb = r + 1; // this seems wrong... // zero Lambda @@ -583,16 +570,16 @@ void OutletRes1D::eval(size_t jg, doublereal* xg, doublereal* rg, rb[2] = xb[2] - xb[2 + nc]; // specified mass fractions - for (k = 4; k < nc; k++) { + for (size_t k = 4; k < nc; k++) { rb[k] = xb[k] - m_yres[k-4]; } } if (m_flow_left) { - nc = m_flow_left->nComponents(); - xb = x - nc; - rb = r - nc; - db = diag - nc; + size_t nc = m_flow_left->nComponents(); + double* xb = x - nc; + double* rb = r - nc; + int* db = diag - nc; if (!m_flow_left->fixed_mdot()) { ; @@ -601,7 +588,7 @@ void OutletRes1D::eval(size_t jg, doublereal* xg, doublereal* rg, } rb[2] = xb[2] - m_temp; // zero dT/dz size_t kSkip = m_flow_left->rightExcessSpecies(); - for (k = 4; k < nc; k++) { + for (size_t k = 4; k < nc; k++) { if (k != kSkip) { rb[k] = xb[k] - m_yres[k-4]; // fixed Y db[k] = 0; @@ -672,22 +659,20 @@ void Surf1D::eval(size_t jg, doublereal* xg, doublereal* rg, doublereal* x = xg + loc(); doublereal* r = rg + loc(); integer* diag = diagg + loc(); - doublereal* xb, *rb; r[0] = x[0] - m_temp; diag[0] = 0; - size_t nc; if (m_flow_right) { - rb = r + 1; - xb = x + 1; + double* rb = r + 1; + double* xb = x + 1; rb[2] = xb[2] - x[0]; // specified T } if (m_flow_left) { - nc = m_flow_left->nComponents(); - rb = r - nc; - xb = x - nc; + size_t nc = m_flow_left->nComponents(); + double* rb = r - nc; + double* xb = x - nc; rb[2] = xb[2] - x[0]; // specified T } } @@ -777,7 +762,6 @@ void ReactingSurf1D::eval(size_t jg, doublereal* xg, doublereal* rg, doublereal* x = xg + loc(); doublereal* r = rg + loc(); integer* diag = diagg + loc(); - doublereal* xb, *rb; // specified surface temp r[0] = x[0] - m_temp; @@ -829,16 +813,15 @@ void ReactingSurf1D::eval(size_t jg, doublereal* xg, doublereal* rg, } if (m_flow_right) { - rb = r + 1; - xb = x + 1; + double* rb = r + 1; + double* xb = x + 1; rb[2] = xb[2] - x[0]; // specified T } - size_t nc; if (m_flow_left) { - nc = m_flow_left->nComponents(); + size_t nc = m_flow_left->nComponents(); const vector_fp& mwleft = m_phase_left->molecularWeights(); - rb =r - nc; - xb = x - nc; + double* rb = r - nc; + double* xb = x - nc; rb[2] = xb[2] - x[0]; // specified T size_t nSkip = m_flow_left->rightExcessSpecies(); for (size_t nl = 0; nl < m_left_nsp; nl++) {