research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Rename m_y_n_1 and m_ydot_n_1 to _trial to clarify their intent.

Browse source 
Add setPreviousTimeStep function to set m_y_nm1 and m_ydot_nm1
so that time derivative calculations can be run for time-dependent
nonlinear solves.

Change m_ydot_nm1 to a std::vector.
This commit is contained in:
Victor Brunini 2013-01-15 20:17:00 +00:00
parent a709d5f816
commit 64aef5fac7
2 changed files with 52 additions and 24 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

23
include/cantera/numerics/NonlinearSolver.h
View file

@ -544,6 +544,17 @@ namespace Cantera {
doublereal time_curr, GeneralMatrix & jac, int &num_newt_its, doublereal time_curr, GeneralMatrix & jac, int &num_newt_its,
int &num_linear_solves, int &num_backtracks, int loglevelInput); int &num_linear_solves, int &num_backtracks, int loglevelInput);
//! Set the values for the previous time step
/*!
* We set the values for the previous time step here. These are used in the nonlinear
* solve because they affect the calculation of ydot.
*
* @param y_nm1 Value of the solution vector at the previous time step
* @param ydot_nm1 Value of the solution vector derivative at the previous time step
*/
virtual void
setPreviousTimeStep(const std::vector<doublereal>& y_nm1, const std::vector<doublereal>& ydot_nm1);
private: private:
//! Set the column scales //! Set the column scales
void calcColumnScales(); void calcColumnScales();
@ -954,11 +965,14 @@ namespace Cantera {
//! Vector containing the solution at the previous time step //! Vector containing the solution at the previous time step
std::vector<doublereal> m_y_nm1; std::vector<doublereal> m_y_nm1;
//! Vector containing the solution at the previous time step //! Vector containing the solution derivative at the previous time step
std::vector<doublereal> m_y_n_1; std::vector<doublereal> m_ydot_nm1;
//! Vector containing the solution at the new point that is to be considered
std::vector<doublereal> m_y_n_trial;
//! Value of the solution time derivative at the new point that is to be considered //! Value of the solution time derivative at the new point that is to be considered
std::vector<doublereal> m_ydot_n_1; std::vector<doublereal> m_ydot_trial;
//! Value of the step to be taken in the solution //! Value of the step to be taken in the solution
std::vector<doublereal> m_step_1; std::vector<doublereal> m_step_1;
@ -1109,9 +1123,6 @@ namespace Cantera {
//! Base value of the absolute tolerance //! Base value of the absolute tolerance
doublereal atolBase_; doublereal atolBase_;
//! Pointer containing the solution derivative at the previous time step
doublereal *m_ydot_nm1;
//! absolute tolerance in the solution unknown //! absolute tolerance in the solution unknown
/*! /*!
* This is used to evaluating the weighting factor * This is used to evaluating the weighting factor

53
src/numerics/NonlinearSolver.cpp
View file

@ -109,8 +109,8 @@ namespace Cantera {
m_y_n_curr(0), m_y_n_curr(0),
m_ydot_n_curr(0), m_ydot_n_curr(0),
m_y_nm1(0), m_y_nm1(0),
m_y_n_1(0), m_y_n_trial(0),
m_ydot_n_1(0), m_ydot_trial(0),
m_colScales(0), m_colScales(0),
m_rowScales(0), m_rowScales(0),
m_rowWtScales(0), m_rowWtScales(0),
@ -192,8 +192,9 @@ namespace Cantera {
m_y_n_curr.resize(neq_, 0.0); m_y_n_curr.resize(neq_, 0.0);
m_ydot_n_curr.resize(neq_, 0.0); m_ydot_n_curr.resize(neq_, 0.0);
m_y_nm1.resize(neq_, 0.0); m_y_nm1.resize(neq_, 0.0);
m_y_n_1.resize(neq_, 0.0); m_ydot_nm1.resize(neq_, 0.0);
m_ydot_n_1.resize(neq_, 0.0); m_y_n_trial.resize(neq_, 0.0);
m_ydot_trial.resize(neq_, 0.0);
m_colScales.resize(neq_, 1.0); m_colScales.resize(neq_, 1.0);
m_rowScales.resize(neq_, 1.0); m_rowScales.resize(neq_, 1.0);
m_rowWtScales.resize(neq_, 1.0); m_rowWtScales.resize(neq_, 1.0);
@ -204,7 +205,7 @@ namespace Cantera {
atolk_.resize(neq_, atolBase_); atolk_.resize(neq_, atolBase_);
deltaX_Newton_.resize(neq_, 0.0); deltaX_Newton_.resize(neq_, 0.0);
m_step_1.resize(neq_, 0.0); m_step_1.resize(neq_, 0.0);
m_y_n_1.resize(neq_, 0.0); m_y_n_trial.resize(neq_, 0.0);
doublereal hb = std::numeric_limits<double>::max(); doublereal hb = std::numeric_limits<double>::max();
m_y_high_bounds.resize(neq_, hb); m_y_high_bounds.resize(neq_, hb);
m_y_low_bounds.resize(neq_, -hb); m_y_low_bounds.resize(neq_, -hb);
@ -232,8 +233,8 @@ namespace Cantera {
m_y_n_curr(0), m_y_n_curr(0),
m_ydot_n_curr(0), m_ydot_n_curr(0),
m_y_nm1(0), m_y_nm1(0),
m_y_n_1(0), m_y_n_trial(0),
m_ydot_n_1(0), m_ydot_trial(0),
m_step_1(0), m_step_1(0),
m_colScales(0), m_colScales(0),
m_rowScales(0), m_rowScales(0),
@ -337,8 +338,9 @@ namespace Cantera {
m_y_n_curr = right.m_y_n_curr; m_y_n_curr = right.m_y_n_curr;
m_ydot_n_curr = right.m_ydot_n_curr; m_ydot_n_curr = right.m_ydot_n_curr;
m_y_nm1 = right.m_y_nm1; m_y_nm1 = right.m_y_nm1;
m_y_n_1 = right.m_y_n_1; m_ydot_nm1 = right.m_ydot_nm1;
m_ydot_n_1 = right.m_ydot_n_1; m_y_n_trial = right.m_y_n_trial;
m_ydot_trial = right.m_ydot_trial;
m_step_1 = right.m_step_1; m_step_1 = right.m_step_1;
m_colScales = right.m_colScales; m_colScales = right.m_colScales;
m_rowScales = right.m_rowScales; m_rowScales = right.m_rowScales;
@ -3079,7 +3081,7 @@ namespace Cantera {
if (SolnType != NSOLN_TYPE_STEADY_STATE || ydot_comm) { if (SolnType != NSOLN_TYPE_STEADY_STATE || ydot_comm) {
mdp::mdp_copy_dbl_1(DATA_PTR(m_ydot_n_curr), ydot_comm, neq_); mdp::mdp_copy_dbl_1(DATA_PTR(m_ydot_n_curr), ydot_comm, neq_);
mdp::mdp_copy_dbl_1(DATA_PTR(m_ydot_n_1), ydot_comm, neq_); mdp::mdp_copy_dbl_1(DATA_PTR(m_ydot_trial), ydot_comm, neq_);
} }
// Redo the solution weights every time we enter the function // Redo the solution weights every time we enter the function
createSolnWeights(DATA_PTR(m_y_n_curr)); createSolnWeights(DATA_PTR(m_y_n_curr));
@ -3326,10 +3328,10 @@ namespace Cantera {
if (s_print_DogLeg && m_print_flag >= 4) { if (s_print_DogLeg && m_print_flag >= 4) {
printf("\t solve_nonlinear_problem(): Compare descent rates for Cauchy and Newton directions\n"); printf("\t solve_nonlinear_problem(): Compare descent rates for Cauchy and Newton directions\n");
descentComparison(time_curr, DATA_PTR(m_ydot_n_curr), DATA_PTR(m_ydot_n_1), i_numTrials); descentComparison(time_curr, DATA_PTR(m_ydot_n_curr), DATA_PTR(m_ydot_trial), i_numTrials);
} else { } else {
if (doDogLeg_) { if (doDogLeg_) {
descentComparison(time_curr, DATA_PTR(m_ydot_n_curr), DATA_PTR(m_ydot_n_1), i_numTrials); descentComparison(time_curr, DATA_PTR(m_ydot_n_curr), DATA_PTR(m_ydot_trial), i_numTrials);
} }
} }
@ -3347,7 +3349,7 @@ namespace Cantera {
printf("\t solve_nonlinear_problem(): Calculate damping along dog-leg path to ensure residual decrease\n"); printf("\t solve_nonlinear_problem(): Calculate damping along dog-leg path to ensure residual decrease\n");
} }
retnDamp = dampDogLeg(time_curr, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr), retnDamp = dampDogLeg(time_curr, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr),
m_step_1, DATA_PTR(m_y_n_1), DATA_PTR(m_ydot_n_1), stepNorm_1, stepNorm_2, jac, i_numTrials); m_step_1, DATA_PTR(m_y_n_trial), DATA_PTR(m_ydot_trial), stepNorm_1, stepNorm_2, jac, i_numTrials);
} }
#ifdef DEBUG_MODE #ifdef DEBUG_MODE
else { else {
@ -3370,7 +3372,7 @@ namespace Cantera {
*/ */
if (!doDogLeg_) { if (!doDogLeg_) {
retnDamp = dampStep(time_curr, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr), retnDamp = dampStep(time_curr, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr),
DATA_PTR(m_step_1), DATA_PTR(m_y_n_1), DATA_PTR(m_ydot_n_1), DATA_PTR(m_step_1), DATA_PTR(m_y_n_trial), DATA_PTR(m_ydot_trial),
DATA_PTR(m_wksp_2), stepNorm_2, jac, frst, i_numTrials); DATA_PTR(m_wksp_2), stepNorm_2, jac, frst, i_numTrials);
frst = false; frst = false;
num_backtracks += i_numTrials; num_backtracks += i_numTrials;
@ -3406,7 +3408,7 @@ namespace Cantera {
* Do a full residual calculation with the unlagged solution components. * Do a full residual calculation with the unlagged solution components.
* Then get the norm of the residual * Then get the norm of the residual
*/ */
info = doResidualCalc(time_curr, NSOLN_TYPE_STEADY_STATE, DATA_PTR(m_y_n_1), DATA_PTR(m_ydot_n_1)); info = doResidualCalc(time_curr, NSOLN_TYPE_STEADY_STATE, DATA_PTR(m_y_n_trial), DATA_PTR(m_ydot_trial));
if (info != 1) { if (info != 1) {
if (m_print_flag > 0) { if (m_print_flag > 0) {
printf("\t solve_nonlinear_problem(): current trial step and damping led to Residual Calc " printf("\t solve_nonlinear_problem(): current trial step and damping led to Residual Calc "
@ -3416,7 +3418,7 @@ namespace Cantera {
goto done; goto done;
} }
if (m_print_flag >= 4) { if (m_print_flag >= 4) {
m_normResid_full = residErrorNorm(DATA_PTR(m_resid), " Resulting full residual norm", 10, DATA_PTR(m_y_n_1)); m_normResid_full = residErrorNorm(DATA_PTR(m_resid), " Resulting full residual norm", 10, DATA_PTR(m_y_n_trial));
if (fabs(m_normResid_full - m_normResid_1) > 1.0E-3 * ( m_normResid_1 + m_normResid_full + 1.0E-4)) { if (fabs(m_normResid_full - m_normResid_1) > 1.0E-3 * ( m_normResid_1 + m_normResid_full + 1.0E-4)) {
if (m_print_flag >= 4) { if (m_print_flag >= 4) {
printf("\t solve_nonlinear_problem(): Full residual norm changed from %g to %g due to " printf("\t solve_nonlinear_problem(): Full residual norm changed from %g to %g due to "
@ -3441,13 +3443,13 @@ namespace Cantera {
bool m_filterIntermediate = false; bool m_filterIntermediate = false;
if (m_filterIntermediate) { if (m_filterIntermediate) {
if (retnDamp == NSOLN_RETN_CONTINUE) { if (retnDamp == NSOLN_RETN_CONTINUE) {
(void) filterNewSolution(time_n, DATA_PTR(m_y_n_1), DATA_PTR(m_ydot_n_1)); (void) filterNewSolution(time_n, DATA_PTR(m_y_n_trial), DATA_PTR(m_ydot_trial));
} }
} }
// Exchange new for curr solutions // Exchange new for curr solutions
if (retnDamp >= NSOLN_RETN_CONTINUE) { if (retnDamp >= NSOLN_RETN_CONTINUE) {
mdp::mdp_copy_dbl_1(DATA_PTR(m_y_n_curr), CONSTD_DATA_PTR(m_y_n_1), neq_); mdp::mdp_copy_dbl_1(DATA_PTR(m_y_n_curr), CONSTD_DATA_PTR(m_y_n_trial), neq_);
if (solnType_ != NSOLN_TYPE_STEADY_STATE) { if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
calc_ydot(m_order, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr)); calc_ydot(m_order, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr));
@ -3609,6 +3611,21 @@ namespace Cantera {
return retnCode; return retnCode;
} }
//==================================================================================================================== //====================================================================================================================
//! Set the values for the previous time step
/*!
* We set the values for the previous time step here. These are used in the nonlinear
* solve because they affect the calculation of ydot.
*
* @param y_nm1 Value of the solution vector at the previous time step
* @param ydot_nm1 Value of the solution vector derivative at the previous time step
*/
void NonlinearSolver::
setPreviousTimeStep(const std::vector<doublereal>& y_nm1, const std::vector<doublereal>& ydot_nm1)
{
m_y_nm1 = y_nm1;
m_ydot_nm1 = ydot_nm1;
}
//====================================================================================================================
// Print solution norm contribution // Print solution norm contribution
/* /*
* Prints out the most important entries to the update to the solution vector for the current step * Prints out the most important entries to the update to the solution vector for the current step