[1D] Make better use of local variables
This commit is contained in:
Ray Speth
parent
ce9a17fb76
commit
cd958a343b
6 changed files with 145 additions and 193 deletions
|
|
@ -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);
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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++) {
|
||||
|
|
|
|||
|
|
@ -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<size_t> 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<size_t> 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;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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<size_t>(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<XML_Node*> 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");
|
||||
|
|
|
|||
|
|
@ -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++) {
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue