[1D] Avoid Jacobian update after successful solve

Calling solve() again after a successful solution no longer automatically
triggers a Jacobian update. This enables relatively efficient sensitivity
analysis by sequentially perturbing the reaction multipliers and re-solving
the system. Since the perturbed system is close to the orignal, the solution
can be found after only a few steps, even when using the original Jacobian.
This commit is contained in:
Ray Speth 2014-03-24 04:03:44 +00:00
parent 17500064c2
commit aa4dcf393c
3 changed files with 23 additions and 5 deletions

View file

@ -183,31 +183,49 @@ public:
}
void solveEnergyEqn(size_t j=npos) {
bool changed = false;
if (j == npos)
for (size_t i = 0; i < m_points; i++) {
if (!m_do_energy[i]) {
changed = true;
}
m_do_energy[i] = true;
}
else {
if (!m_do_energy[j]) {
changed = true;
}
m_do_energy[j] = true;
}
m_refiner->setActive(0, true);
m_refiner->setActive(1, true);
m_refiner->setActive(2, true);
needJacUpdate();
if (changed) {
needJacUpdate();
}
}
void fixTemperature(size_t j=npos) {
bool changed = false;
if (j == npos)
for (size_t i = 0; i < m_points; i++) {
if (m_do_energy[i]) {
changed = true;
}
m_do_energy[i] = false;
}
else {
if (m_do_energy[j]) {
changed = true;
}
m_do_energy[j] = false;
}
m_refiner->setActive(0, false);
m_refiner->setActive(1, false);
m_refiner->setActive(2, false);
needJacUpdate();
if (changed) {
needJacUpdate();
}
}
bool doSpecies(size_t k) {

View file

@ -382,7 +382,6 @@ int MultiNewton::solve(doublereal* x0, doublereal* x1,
int nJacReeval = 0;
while (1 > 0) {
// Check whether the Jacobian should be re-evaluated.
if (jac.age() > m_maxAge) {
writelog("\nMaximum Jacobian age reached ("+int2str(m_maxAge)+")\n", loglevel);
@ -427,6 +426,7 @@ int MultiNewton::solve(doublereal* x0, doublereal* x1,
// convergence
else if (m == 1) {
jac.setAge(0); // for efficient sensitivity analysis
break;
}

View file

@ -209,11 +209,11 @@ void StFlow::_finalize(const doublereal* x)
bool e = m_do_energy[0];
for (j = 0; j < m_points; j++) {
if (e || nz == 0) {
setTemperature(j, T(x, j));
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);
setTemperature(j, tt);
m_fixedtemp[j] = tt;
}
for (k = 0; k < m_nsp; k++) {
setMassFraction(j, k, Y(x, k, j));