diff --git a/Cantera/src/numerics/DAE_Solver.h b/Cantera/src/numerics/DAE_Solver.h index d87f3b9a8..22105df57 100644 --- a/Cantera/src/numerics/DAE_Solver.h +++ b/Cantera/src/numerics/DAE_Solver.h @@ -26,6 +26,12 @@ namespace Cantera { + /** + * @defgroup numerics Numerical Utilities within Cantera + * + * + */ + class Jacobian { public: Jacobian(){} diff --git a/Cantera/src/numerics/Makefile.in b/Cantera/src/numerics/Makefile.in index d1c5cf93e..c11673501 100644 --- a/Cantera/src/numerics/Makefile.in +++ b/Cantera/src/numerics/Makefile.in @@ -34,13 +34,13 @@ CXX_FLAGS = @CXXFLAGS@ $(LOCAL_DEFS) $(CXX_OPT) $(PIC_FLAG) $(DEBUG_FLAG) NUMERICS_OBJ = DenseMatrix.o funcs.o Func1.o \ ODE_integrators.o BandMatrix.o DAE_solvers.o \ - funcs.o sort.o SquareMatrix.o + funcs.o sort.o SquareMatrix.o ResidJacEval.o NonlinearSolve.o NUMERICS_H = ArrayViewer.h DenseMatrix.h \ funcs.h ctlapack.h Func1.h FuncEval.h \ polyfit.h\ BandMatrix.h Integrator.h DAE_Solver.h ResidEval.h sort.h \ - SquareMatrix.h + SquareMatrix.h ResidJacEval.h NonlinearSolve.h ifeq ($(use_sundials), 1) ODEPACKAGE_H = CVodesIntegrator.h diff --git a/Cantera/src/numerics/NonlinearSolver.cpp b/Cantera/src/numerics/NonlinearSolver.cpp new file mode 100644 index 000000000..0b9a9599f --- /dev/null +++ b/Cantera/src/numerics/NonlinearSolver.cpp @@ -0,0 +1,840 @@ +/** + * + * @file NonlinearSolver.cpp + * + * Damped Newton solver for 1D multi-domain problems + */ + +/* + * $Author$ + * $Date$ + * $Revision$ + */ +/* + * Copywrite 2004 Sandia Corporation. Under the terms of Contract + * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government + * retains certain rights in this software. + * See file License.txt for licensing information. + */ + + + +#include "SquareMatrix.h" +#include "NonlinearSolver.h" +//#include "md_timer.h" +#include "clockWC.h" +#include "vec_functions.h" +#include +extern double second(); +#include "mdp_allo.h" +extern void print_line(const char *, int); + +#include + +#include +#include + +#ifndef MAX +#define MAX(x,y) (( (x) > (y) ) ? (x) : (y)) +#define MIN(x,y) (( (x) < (y) ) ? (x) : (y)) +#endif + +using namespace std; + +namespace Cantera { + + + + //----------------------------------------------------------- + // Constants + //----------------------------------------------------------- + + const double DampFactor = 4; + const int NDAMP = 10; + + //----------------------------------------------------------- + // Static Functions + //----------------------------------------------------------- + + static void print_line(const char *str, int n) { + for (int i = 0; i < n; i++) { + printf("%s", str); + } + printf("\n"); + } + + // Default constructor + /* + * @param func Residual and jacobian evaluator function object + */ + NonlinearSolver::NonlinearSolver(ResidJacEval *func) : + m_func(func), + neq_(0), + delta_t_n(-1.0), + m_nfe(0), + m_colScaling(0), + m_rowScaling(0), + m_numTotalLinearSolves(0), + m_numTotalNewtIts(0), + m_min_newt_its(0), + filterNewstep(0), + time_n(0.0), + m_matrixConditioning(0), + m_order(1) + { + neq_ = m_func->nEquations(); + } + + NonlinearSolver::NonlinearSolver(const NonlinearSolver &right) { + *this =operator=(right); + } + + + NonlinearSolver::~NonlinearSolver() { + } + + NonlinearSolver& NonlinearSolver::operator=(const NonlinearSolver &right) { + if (this == &right) { + return *this; + } + // rely on the ResidJacEval duplMyselfAsresidJacEval() function to + // create a deep copy + m_func = right.m_func->duplMyselfAsResidJacEval(); + neq_ = right.neq_; + delta_t_n = right.delta_t_n; + m_nfe = right.m_nfe; + m_colScaling = right.m_colScaling; + m_rowScaling = right.m_rowScaling; + m_numTotalLinearSolves = right.m_numTotalLinearSolves; + m_numTotalNewtIts = right.m_numTotalNewtIts; + m_min_newt_its = right.m_min_newt_its; + filterNewstep = right.filterNewstep; + time_n = right.time_n; + m_matrixConditioning = right.m_matrixConditioning; + m_order = right.m_order; + + return *this; + } + + /** + * L2 Norm of a delta in the solution + * + * The second argument has a default of false. However, + * if true, then a table of the largest values is printed + * out to standard output. + */ + double NonlinearSolver::soln_error_norm(const double * const delta_y, + bool printLargest) + { + int i; + double sum_norm = 0.0, error; + for (i = 0; i < neq_; i++) { + error = delta_y[i] / m_ewt[i]; + sum_norm += (error * error); + } + sum_norm = sqrt(sum_norm / neq_); + if (printLargest) { + const int num_entries = 8; + double dmax1, normContrib; + int j; + int *imax = mdp::mdp_alloc_int_1(num_entries, -1); + printf("\t\tPrintout of Largest Contributors to norm " + "of value (%g)\n", sum_norm); + printf("\t\t I ysoln deltaY weightY " + "Error_Norm**2\n"); + printf("\t\t "); print_line("-", 80); + for (int jnum = 0; jnum < num_entries; jnum++) { + dmax1 = -1.0; + for (i = 0; i < neq_; i++) { + bool used = false; + for (j = 0; j < jnum; j++) { + if (imax[j] == i) used = true; + } + if (!used) { + error = delta_y[i] / m_ewt[i]; + normContrib = sqrt(error * error); + if (normContrib > dmax1) { + imax[jnum] = i; + dmax1 = normContrib; + } + } + } + i = imax[jnum]; + if (i >= 0) { + printf("\t\t %4d %12.4e %12.4e %12.4e %12.4e\n", + i, m_y_n[i], delta_y[i], m_ewt[i], dmax1); + } + } + printf("\t\t "); print_line("-", 80); + mdp::mdp_safe_free((void **) &imax); + } + return sum_norm; + } + + + /** + * setColumnScales(): + * + * Set the column scaling vector at the current time + */ + void NonlinearSolver::setColumnScales() { + m_func->calcSolnScales(time_n, DATA_PTR(m_y_n), DATA_PTR(m_y_nm1), + DATA_PTR(m_colScales)); + } + + + void NonlinearSolver::doResidualCalc(const double time_curr, const int typeCalc, + const double * const y_curr, + const double * const ydot_curr, double* const residual, + int loglevel) + { + + + // Calculate the current residual + // Put the current residual into the vector, delta_y[] + // We need to pull this out of this function and carry it in. + m_func->evalResidNJ(time_curr, delta_t_n, y_curr, ydot_curr, residual); + m_nfe++; + } + + + // Compute the undamped Newton step + /* + * Compute the undamped Newton step. The residual function is + * evaluated at the current time, t_n, at the current values of the + * solution vector, m_y_n, and the solution time derivative, m_ydot_n. + * The Jacobian is not recomputed. + * + * A factored jacobian is reused, if available. If a factored jacobian + * is not available, then the jacobian is factored. Before factoring, + * the jacobian is row and column-scaled. Column scaling is not + * recomputed. The row scales are recomputed here, after column + * scaling has been implemented. + */ + void NonlinearSolver::doNewtonSolve(const double time_curr, const double * const y_curr, + const double * const ydot_curr, double* const delta_y, + SquareMatrix& jac, int loglevel) + { + int irow, jcol; + + + //! multiply the residual by -1 + for (int n = 0; n < neq_; n++) { + delta_y[n] = -delta_y[n]; + } + + + /* + * Column scaling -> We scale the columns of the Jacobian + * by the nominal important change in the solution vector + */ + if (m_colScaling) { + if (!jac.m_factored) { + /* + * Go get new scales -> Took this out of this inner loop. + * Needs to be done at a larger scale. + */ + // setColumnScales(); + + /* + * Scale the new Jacobian + */ + double *jptr = &(*(jac.begin())); + for (jcol = 0; jcol < neq_; jcol++) { + for (irow = 0; irow < neq_; irow++) { + *jptr *= m_colScales[jcol]; + jptr++; + } + } + } + } + + // if (m_matrixConditioning) { + // if (jac.m_factored) { + // m_func->matrixConditioning(0, neq_, delta_y); + // } else { + //double *jptr = &(*(jac.begin())); + // m_func->matrixConditioning(jptr, neq_, delta_y); + // } + //} + + /* + * row sum scaling -> Note, this is an unequivical success + * at keeping the small numbers well balanced and + * nonnegative. + */ + if (m_rowScaling) { + if (! jac.m_factored) { + /* + * Ok, this is ugly. jac.begin() returns an vector iterator + * to the first data location. + * Then &(*()) reverts it to a double *. + */ + double *jptr = &(*(jac.begin())); + for (irow = 0; irow < neq_; irow++) { + m_rowScales[irow] = 0.0; + } + for (jcol = 0; jcol < neq_; jcol++) { + for (irow = 0; irow < neq_; irow++) { + m_rowScales[irow] += fabs(*jptr); + jptr++; + } + } + + jptr = &(*(jac.begin())); + for (jcol = 0; jcol < neq_; jcol++) { + for (irow = 0; irow < neq_; irow++) { + *jptr /= m_rowScales[irow]; + jptr++; + } + } + } + for (irow = 0; irow < neq_; irow++) { + delta_y[irow] /= m_rowScales[irow]; + } + } + + + /* + * Solve the system -> This also involves inverting the + * matrix + */ + (void) jac.solve(delta_y); + + + /* + * reverse the column scaling if there was any. + */ + if (m_colScaling) { + for (irow = 0; irow < neq_; irow++) { + delta_y[irow] *= m_colScales[irow]; + } + } + +#ifdef DEBUG_JAC + if (printJacContributions) { + for (int iNum = 0; iNum < numRows; iNum++) { + if (iNum > 0) focusRow++; + double dsum = 0.0; + vector_fp& Jdata = jacBack.data(); + double dRow = Jdata[neq_ * focusRow + focusRow]; + printf("\n Details on delta_Y for row %d \n", focusRow); + printf(" Value before = %15.5e, delta = %15.5e," + "value after = %15.5e\n", y_curr[focusRow], + delta_y[focusRow], + y_curr[focusRow] + delta_y[focusRow]); + if (!freshJac) { + printf(" Old Jacobian\n"); + } + printf(" col delta_y aij " + "contrib \n"); + printf("--------------------------------------------------" + "---------------------------------------------\n"); + printf(" Res(%d) %15.5e %15.5e %15.5e (Res = %g)\n", + focusRow, delta_y[focusRow], + dRow, RRow[iNum] / dRow, RRow[iNum]); + dsum += RRow[iNum] / dRow; + for (int ii = 0; ii < neq_; ii++) { + if (ii != focusRow) { + double aij = Jdata[neq_ * ii + focusRow]; + double contrib = aij * delta_y[ii] * (-1.0) / dRow; + dsum += contrib; + if (fabs(contrib) > Pcutoff) { + printf("%6d %15.5e %15.5e %15.5e\n", ii, + delta_y[ii] , aij, contrib); + } + } + } + printf("--------------------------------------------------" + "---------------------------------------------\n"); + printf(" %15.5e %15.5e\n", + delta_y[focusRow], dsum); + } + } + +#endif + + m_numTotalLinearSolves++; + } + + /************************************************************************** + * + * boundStep(): + * + * Return the factor by which the undamped Newton step 'step0' + * must be multiplied in order to keep all solution components in + * all domains between their specified lower and upper bounds. + * Other bounds may be applied here as well. + * + * Currently the bounds are hard coded into this routine: + * + * Minimum value for all variables: - 0.01 * m_ewt[i] + * Maximum value = none. + * + * Thus, this means that all solution components are expected + * to be numerical greater than zero in the limit of time step + * truncation errors going to zero. + * + * Delta bounds: The idea behind these is that the Jacobian + * couldn't possibly be representative if the + * variable is changed by a lot. (true for + * nonlinear systems, false for linear systems) + * Maximum increase in variable in any one newton iteration: + * factor of 2 + * Maximum decrease in variable in any one newton iteration: + * factor of 5 + */ + double NonlinearSolver::boundStep(const double* y, + const double* step0, int loglevel) { + int i, i_lower = -1, i_fbounds, ifbd = 0, i_fbd = 0; + double fbound = 1.0, f_lowbounds = 1.0, f_delta_bounds = 1.0; + double ff, y_new, ff_alt; + for (i = 0; i < neq_; i++) { + y_new = y[i] + step0[i]; + if ((y_new < (-0.01 * m_ewt[i])) && y[i] >= 0.0) { + ff = 0.9 * (y[i] / (y[i] - y_new)); + if (ff < f_lowbounds) { + f_lowbounds = ff; + i_lower = i; + } + } + /** + * Now do a delta bounds + * Increase variables by a factor of 2 only + * decrease variables by a factor of 5 only + */ + ff = 1.0; + if ((fabs(y_new) > 2.0 * fabs(y[i])) && + (fabs(y_new-y[i]) > m_ewt[i])) { + ff = fabs(y[i]/(y_new - y[i])); + ff_alt = fabs(m_ewt[i] / (y_new - y[i])); + ff = MAX(ff, ff_alt); + ifbd = 1; + } + if ((fabs(5.0 * y_new) < fabs(y[i])) && + (fabs(y_new - y[i]) > m_ewt[i])) { + ff = y[i]/(y_new-y[i]) * (1.0 - 5.0)/5.0; + ff_alt = fabs(m_ewt[i] / (y_new - y[i])); + ff = MAX(ff, ff_alt); + ifbd = 0; + } + if (ff < f_delta_bounds) { + f_delta_bounds = ff; + i_fbounds = i; + i_fbd = ifbd; + } + f_delta_bounds = MIN(f_delta_bounds, ff); + } + fbound = MIN(f_lowbounds, f_delta_bounds); + /* + * Report on any corrections + */ + if (loglevel > 1) { + if (fbound != 1.0) { + if (f_lowbounds < f_delta_bounds) { + printf("\t\tboundStep: Variable %d causing lower bounds " + "damping of %g\n", + i_lower, f_lowbounds); + } else { + if (ifbd) { + printf("\t\tboundStep: Decrease of Variable %d causing " + "delta damping of %g\n", + i_fbd, f_delta_bounds); + } else { + printf("\t\tboundStep: Increase of variable %d causing" + "delta damping of %g\n", + i_fbd, f_delta_bounds); + } + } + } + } + //return fbound; + return 1.0; + } + + /************************************************************************** + * + * dampStep(): + * + * On entry, step0 must contain an undamped Newton step for the + * solution x0. This method attempts to find a damping coefficient + * such that the next undamped step would have a norm smaller than + * that of step0. If successful, the new solution after taking the + * damped step is returned in y1, and the undamped step at y1 is + * returned in step1. + */ + int NonlinearSolver::dampStep(double time_curr, const double* y0, + const double *ydot0, const double* step0, + double* y1, double* ydot1, double* step1, + double& s1, SquareMatrix& jac, + int& loglevel, bool writetitle, + int& num_backtracks) { + + + // Compute the weighted norm of the undamped step size step0 + double s0 = soln_error_norm(step0); + + // Compute the multiplier to keep all components in bounds + // A value of one indicates that there is no limitation + // on the current step size in the nonlinear method due to + // bounds constraints (either negative values of delta + // bounds constraints. + double fbound = boundStep(y0, step0, loglevel); + + // if fbound is very small, then y0 is already close to the + // boundary and step0 points out of the allowed domain. In + // this case, the Newton algorithm fails, so return an error + // condition. + if (fbound < 1.e-10) { + if (loglevel > 1) printf("\t\t\tdampStep: At limits.\n"); + return -3; + } + + //-------------------------------------------- + // Attempt damped step + //-------------------------------------------- + + // damping coefficient starts at 1.0 + double damp = 1.0; + int j, m; + double ff; + num_backtracks = 0; + for (m = 0; m < NDAMP; m++) { + + ff = fbound*damp; + + // step the solution by the damped step size + /* + * Whenever we update the solution, we must also always + * update the time derivative. + */ + for (j = 0; j < neq_; j++) { + y1[j] = y0[j] + ff*step0[j]; + // HKM setting intermediate y's to zero was a tossup. + // slightly different, equivalent results + //#ifdef DEBUG_HKM + // y1[j] = MAX(0.0, y1[j]); + //#endif + } + calc_ydot(m_order, y1, ydot1); + + doResidualCalc(time_curr, NSOLN_TYPE_STEADY_STATE, y1, ydot1, step1, loglevel); + + // compute the next undamped step, step1[], that would result + // if y1[] were accepted. + + doNewtonSolve(time_curr, y1, ydot1, step1, jac, loglevel); + + + // compute the weighted norm of step1 + s1 = soln_error_norm(step1); + + // write log information + if (loglevel > 3) { + print_solnDelta_norm_contrib((const double *) step0, + "DeltaSolnTrial", + (const double *) step1, + "DeltaSolnTrialTest", + "dampNewt: Important Entries for " + "Weighted Soln Updates:", + y0, y1, ff, 5); + } + if (loglevel > 1) { + printf("\t\t\tdampNewt: s0 = %g, s1 = %g, fbound = %g," + "damp = %g\n", s0, s1, fbound, damp); + } + + + // if the norm of s1 is less than the norm of s0, then + // accept this damping coefficient. Also accept it if this + // step would result in a converged solution. Otherwise, + // decrease the damping coefficient and try again. + + if (s1 < 1.0E-5 || s1 < s0) { + if (loglevel > 2) { + if (s1 > s0) { + if (s1 > 1.0) { + printf("\t\t\tdampStep: current trial step and damping" + " coefficient accepted because test step < 1\n"); + printf("\t\t\t s1 = %g, s0 = %g\n", s1, s0); + } + } + } + break; + } else { + if (loglevel > 1) { + printf("\t\t\tdampStep: current step rejected: (s1 = %g > " + "s0 = %g)", s1, s0); + if (m < (NDAMP-1)) { + printf(" Decreasing damping factor and retrying"); + } else { + printf(" Giving up!!!"); + } + printf("\n"); + } + } + num_backtracks++; + damp /= DampFactor; + } + + // If a damping coefficient was found, return 1 if the + // solution after stepping by the damped step would represent + // a converged solution, and return 0 otherwise. If no damping + // coefficient could be found, return -2. + if (m < NDAMP) { + if (s1 > 1.0) return 0; + else return 1; + } else { + if (s1 < 0.5 && (s0 < 0.5)) return 1; + if (s1 < 1.0) return 0; + return -2; + } + } + + /************************************************************************** + * + * solve_nonlinear_problem(): + * + * Find the solution to F(X) = 0 by damped Newton iteration. On + * entry, x0 contains an initial estimate of the solution. On + * successful return, x1 contains the converged solution. + * + * SolnType = TRANSIENT -> we will assume we are relaxing a transient + * equation system for now. Will make it more general later, + * if an application comes up. + * + */ + int NonlinearSolver::solve_nonlinear_problem(int SolnType, double* y_comm, + double* ydot_comm, double CJ, + double time_curr, + SquareMatrix& jac, + int &num_newt_its, + int &num_linear_solves, + int &num_backtracks, + int loglevelInput) + { + double t0 = second(); + bool m_residCurrent = false; + int m = 0; + bool forceNewJac = false; + double s1=1.e30; + + std::vector y_curr(neq_, 0.0); + std::vector ydot_curr(neq_, 0.0); + std::vector stp(neq_, 0.0); + std::vector stp1(neq_, 0.0); + + std::vector y_new(neq_, 0.0); + std::vector ydot_new(neq_, 0.0); + + mdp::mdp_copy_dbl_1(DATA_PTR(y_curr), y_comm, neq_); + // copyn((size_t)neq_, y_comm, y_curr); + mdp::mdp_copy_dbl_1(DATA_PTR(ydot_curr), ydot_comm, neq_); + + + + bool frst = true; + num_newt_its = 0; + num_linear_solves = - m_numTotalLinearSolves; + num_backtracks = 0; + int i_backtracks; + int loglevel = loglevelInput; + + while (1 > 0) { + + /* + * Increment Newton Solve counter + */ + m_numTotalNewtIts++; + num_newt_its++; + + + if (loglevel > 1) { + printf("\t\tSolve_Nonlinear_Problem: iteration %d:\n", + num_newt_its); + } + + // Check whether the Jacobian should be re-evaluated. + + forceNewJac = true; + + if (forceNewJac) { + if (loglevel > 1) { + printf("\t\t\tGetting a new Jacobian and solving system\n"); + } + beuler_jac(jac, DATA_PTR(m_resid), time_curr, CJ, DATA_PTR(y_curr), DATA_PTR(ydot_curr), + num_newt_its); + m_residCurrent = true; + } else { + if (loglevel > 1) { + printf("\t\t\tSolving system with old jacobian\n"); + } + m_residCurrent = false; + } + /* + * Go get new scales + */ + setColumnScales(); + + + doResidualCalc(time_curr, NSOLN_TYPE_STEADY_STATE, + DATA_PTR(y_curr), DATA_PTR(ydot_curr), DATA_PTR(stp), loglevel); + + // compute the undamped Newton step + doNewtonSolve(time_curr, DATA_PTR(y_curr), DATA_PTR(ydot_curr), DATA_PTR(stp), + jac, loglevel); + + // damp the Newton step + m = dampStep(time_curr, DATA_PTR(y_curr), DATA_PTR(ydot_curr), + DATA_PTR(stp), DATA_PTR(y_new), DATA_PTR(ydot_new), + DATA_PTR(stp1), s1, jac, loglevel, frst, i_backtracks); + frst = false; + num_backtracks += i_backtracks; + + /* + * Impose the minimum number of newton iterations critera + */ + if (num_newt_its < m_min_newt_its) { + if (m == 1) m = 0; + } + /* + * Impose max newton iteration + */ + if (num_newt_its > 20) { + m = -1; + if (loglevel > 1) { + printf("\t\t\tDampnewton unsuccessful (max newts exceeded) sfinal = %g\n", s1); + } + } + + if (loglevel > 1) { + if (m == 1) { + printf("\t\t\tDampNewton iteration successful, nonlin " + "converged sfinal = %g\n", s1); + } else if (m == 0) { + printf("\t\t\tDampNewton iteration successful, get new" + "direction, sfinal = %g\n", s1); + } else { + printf("\t\t\tDampnewton unsuccessful sfinal = %g\n", s1); + } + } + + // If we are converged, then let's use the best solution possible + // for an end result. We did a resolve in dampStep(). Let's update + // the solution to reflect that. + // HKM 5/16 -> Took this out, since if the last step was a + // damped step, then adding stp1[j] is undamped, and + // may lead to oscillations. It kind of defeats the + // purpose of dampStep() anyway. + // if (m == 1) { + // for (int j = 0; j < neq_; j++) { + // y_new[j] += stp1[j]; + // HKM setting intermediate y's to zero was a tossup. + // slightly different, equivalent results + // #ifdef DEBUG_HKM + // y_new[j] = MAX(0.0, y_new[j]); + // #endif + // } + // } + + bool m_filterIntermediate = false; + if (m_filterIntermediate) { + if (m == 0) { + (void) filterNewStep(time_n, DATA_PTR(y_new), DATA_PTR(ydot_new)); + } + } + // Exchange new for curr solutions + if (m == 0 || m == 1) { + mdp::mdp_copy_dbl_1(DATA_PTR(y_curr), DATA_PTR(y_new), neq_); + calc_ydot(m_order, DATA_PTR(y_curr), DATA_PTR(ydot_curr)); + } + + // convergence + if (m == 1) goto done; + + // If dampStep fails, first try a new Jacobian if an old + // one was being used. If it was a new Jacobian, then + // return -1 to signify failure. + else if (m < 0) { + goto done; + } + } + + done: + mdp::mdp_copy_dbl_1(y_comm, DATA_PTR(y_curr), neq_); + mdp::mdp_copy_dbl_1(ydot_comm, DATA_PTR(ydot_curr), neq_); + + + num_linear_solves += m_numTotalLinearSolves; + + double time_elapsed = second() - t0; + if (loglevel > 1) { + if (m == 1) { + printf("\t\tNonlinear problem solved successfully in " + "%d its, time elapsed = %g sec\n", + num_newt_its, time_elapsed); + } + } + return m; + } + + /***************************************************************8 + * + * + */ + void NonlinearSolver:: + print_solnDelta_norm_contrib(const double * const solnDelta0, + const char * const s0, + const double * const solnDelta1, + const char * const s1, + const char * const title, + const double * const y0, + const double * const y1, + double damp, + int num_entries) { + int i, j, jnum; + bool used; + double dmax0, dmax1, error, rel_norm; + printf("\t\t%s currentDamp = %g\n", title, damp); + printf("\t\t I ysoln %10s ysolnTrial " + "%10s weight relSoln0 relSoln1\n", s0, s1); + int *imax = mdp::mdp_alloc_int_1(num_entries, -1); + printf("\t\t "); print_line("-", 90); + for (jnum = 0; jnum < num_entries; jnum++) { + dmax1 = -1.0; + for (i = 0; i < neq_; i++) { + used = false; + for (j = 0; j < jnum; j++) { + if (imax[j] == i) used = true; + } + if (!used) { + error = solnDelta0[i] / m_ewt[i]; + rel_norm = sqrt(error * error); + error = solnDelta1[i] / m_ewt[i]; + rel_norm += sqrt(error * error); + if (rel_norm > dmax1) { + imax[jnum] = i; + dmax1 = rel_norm; + } + } + } + if (imax[jnum] >= 0) { + i = imax[jnum]; + error = solnDelta0[i] / m_ewt[i]; + dmax0 = sqrt(error * error); + error = solnDelta1[i] / m_ewt[i]; + dmax1 = sqrt(error * error); + printf("\t\t %4d %12.4e %12.4e %12.4e %12.4e " + "%12.4e %12.4e %12.4e\n", + i, y0[i], solnDelta0[i], y1[i], + solnDelta1[i], m_ewt[i], dmax0, dmax1); + } + } + printf("\t\t "); print_line("-", 90); + mdp::mdp_safe_free((void **) &imax); + } + + + +} + diff --git a/Cantera/src/numerics/NonlinearSolver.h b/Cantera/src/numerics/NonlinearSolver.h new file mode 100644 index 000000000..9099aa9e1 --- /dev/null +++ b/Cantera/src/numerics/NonlinearSolver.h @@ -0,0 +1,270 @@ +/** + * @file NonlinearSolve.h + * Class that calculates the solution to a nonlinear, dense, set + * of equations (see \ref numerics + * and class \link Cantera::NonlinearSolver NonlinearSolver\endlink). + */ + +/* + * $Date$ + * $Revision$ + */ +/* + * Copywrite 2004 Sandia Corporation. Under the terms of Contract + * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government + * retains certain rights in this software. + * See file License.txt for licensing information. + */ + +#ifndef CT_NONLINEARSOLVER_H +#define CT_NONLINEARSOLVER_H + +#include "ResidJacEval.h" + +namespace Cantera { + +#define NSOLN_TYPE_PSEUDO_TIME_DEPENDENT 2 +#define NSOLN_TYPE_TIME_DEPENDENT 1 +#define NSOLN_TYPE_STEADY_STATE 0 + + //! Class that calculates the solution to a nonlinear system + /*! + * + * @ingroup numerics + */ + class NonlinearSolver { + + //! Default constructor + /*! + * @param func Residual and jacobian evaluator function object + */ + NonlinearSolver(ResidJacEval *func); + + //!Copy Constructor for the %ThermoPhase object. + /*! + * @param right Item to be copied + */ + NonlinearSolver(const NonlinearSolver &right); + + //! Destructor + ~NonlinearSolver(); + + //! Assignment operator + /*! + * This is NOT a virtual function. + * + * @param right Reference to %NonlinearSolver object to be + * copied into the + * current one. + */ + NonlinearSolver& operator=(const NonlinearSolver &right); + + /** + * L2 Norm of a delta in the solution + * + * The second argument has a default of false. However, + * if true, then a table of the largest values is printed + * out to standard output. + */ + double soln_error_norm(const double * const delta_y, + bool printLargest = false); + + //! Compute the current Residual + /*! + * Compute the time dependent residual of + * the set of equations. + */ + void doTDResidualCalc(const double time_curr, const int typeCalc, + const double * const y_curr, + const double * const ydot_curr, double* const residual, + int loglevel); + + //! Compute the current Residual + /*! + * Compute the steady state residual of + * the set of equations. + */ + void doSteadyResidualCalc(const double time_curr, const int typeCalc, + const double * const y_curr, + double* const residual, int loglevel); + + void doResidualCalc(const double time_curr, const int typeCalc, + const double * const y_curr, + const double * const ydot_curr, double* const residual, + int loglevel); + + //! Compute the undamped Newton step + /*! + * + * Compute the undamped Newton step. The residual function is + * evaluated at the current time, t_n, at the current values of the + * solution vector, m_y_n, and the solution time derivative, m_ydot_n. + * The Jacobian is not recomputed. + * + * A factored jacobian is reused, if available. If a factored jacobian + * is not available, then the jacobian is factored. Before factoring, + * the jacobian is row and column-scaled. Column scaling is not + * recomputed. The row scales are recomputed here, after column + * scaling has been implemented. + * + */ + void doNewtonSolve(const double time_curr, const double * const y_curr, + const double * const ydot_curr, double* const delta_y, + SquareMatrix& jac, int loglevel); + + //! + /*! + * + * Return the factor by which the undamped Newton step 'step0' + * must be multiplied in order to keep all solution components in + * all domains between their specified lower and upper bounds. + * Other bounds may be applied here as well. + * + * Currently the bounds are hard coded into this routine: + * + * Minimum value for all variables: - 0.01 * m_ewt[i] + * Maximum value = none. + * + * Thus, this means that all solution components are expected + * to be numerical greater than zero in the limit of time step + * truncation errors going to zero. + * + * Delta bounds: The idea behind these is that the Jacobian + * couldn't possibly be representative if the + * variable is changed by a lot. (true for + * nonlinear systems, false for linear systems) + * Maximum increase in variable in any one newton iteration: + * factor of 2 + * Maximum decrease in variable in any one newton iteration: + * factor of 5 + */ + double boundStep(const double* y, + const double* step0, int loglevel); + + /** + * Internal function to calculate the predicted solution + * at a time step. + */ + void calc_y_pred(int); + + /** + * Internal function to calculate the time derivative at the + * new step + */ + void calc_ydot(int order, double * const y_curr, double * const ydot_curr); + + void beuler_jac(SquareMatrix &, double * const, + double, double, double * const, double * const, int); + + + double filterNewStep(double, double *, double *); + + //! + /*! + * On entry, step0 must contain an undamped Newton step for the + * solution x0. This method attempts to find a damping coefficient + * such that the next undamped step would have a norm smaller than + * that of step0. If successful, the new solution after taking the + * damped step is returned in y1, and the undamped step at y1 is + * returned in step1. + */ + int dampStep(double time_curr, const double* y0, + const double *ydot0, const double* step0, + double* y1, double* ydot1, double* step1, + double& s1, SquareMatrix& jac, + int& loglevel, bool writetitle, + int& num_backtracks); + + + + // Compute the weighted norm of the undamped step size step0 + + //! Find the solution to F(X) = 0 by damped Newton iteration. + /*! + * On + * entry, x0 contains an initial estimate of the solution. On + * successful return, x1 contains the converged solution. + * + * SolnType = TRANSIENT -> we will assume we are relaxing a transient + * equation system for now. Will make it more general later, + * if an application comes up. + * + */ + int solve_nonlinear_problem(int SolnType, double* y_comm, + double* ydot_comm, double CJ, + double time_curr, + SquareMatrix& jac, + int &num_newt_its, + int &num_linear_solves, + int &num_backtracks, + int loglevelInput); + + + void setColumnScales(); + + void + print_solnDelta_norm_contrib(const double * const solnDelta0, + const char * const s0, + const double * const solnDelta1, + const char * const s1, + const char * const title, + const double * const y0, + const double * const y1, + double damp, + int num_entries); + + + + //! Pointer to the residual and jacobian evaluator for the + //! function + /*! + * See ResidJacEval.h for an evaluator. + */ + ResidJacEval *m_func; + + //! Local copy of the number of equations + int neq_; + + std::vector m_ewt; + + std::vector m_y_n; + std::vector m_y_nm1; + std::vector m_colScales; + std::vector m_rowScales; + + std::vector m_resid; + + double delta_t_n; + + //! Counter for the total number of function evaluations + int m_nfe; + + //! The type of column scaled used in the solution of the problem + bool m_colScaling; + + //! int indicating whether row scaling is turned on (1) or not (0) + int m_rowScaling; + + int m_numTotalLinearSolves; + + int m_numTotalNewtIts; + + int m_min_newt_its; + + int filterNewstep; + + //! Current system time + /*! + * Note, we assume even for steady state problems that the residual + * is a function of a system time. + */ + double time_n; + + int m_matrixConditioning; + + int m_order; + }; + +} + +#endif diff --git a/Cantera/src/numerics/ResidEval.h b/Cantera/src/numerics/ResidEval.h index d2c3be01f..ec61bd2d2 100755 --- a/Cantera/src/numerics/ResidEval.h +++ b/Cantera/src/numerics/ResidEval.h @@ -2,7 +2,10 @@ * @file ResidEval.h * */ - +/* + * $Date$ + * $Revision$ + */ // Copyright 2006 California Institute of Technology #ifndef CT_RESIDEVAL_H @@ -14,86 +17,94 @@ #endif #include "ct_defs.h" +#include "ctexceptions.h" namespace Cantera { - const int c_NONE = 0; - const int c_GE_ZERO = 1; - const int c_GT_ZERO = 2; - const int c_LE_ZERO = -1; - const int c_LT_ZERO = -2; + const int c_NONE = 0; + const int c_GE_ZERO = 1; + const int c_GT_ZERO = 2; + const int c_LE_ZERO = -1; + const int c_LT_ZERO = -2; + + /** + * Virtual base class for DAE residual function evaluators. + * Classes derived from ResidEval evaluate the residual function + * \f[ + \vec{F}(t,\vec{y}, \vec{y^\prime}) + * \f] + * The DAE solver attempts to find a solution y(t) such that F = 0. + * @ingroup DAE_Group + */ + class ResidEval { + + public: + + ResidEval() {} + virtual ~ResidEval() {} /** - * Virtual base class for DAE residual function evaluators. - * Classes derived from ResidEval evaluate the residual function - * \f[ - \vec{F}(t,\vec{y}, \vec{y^\prime}) - * \f] - * The DAE solver attempts to find a solution y(t) such that F = 0. - * @ingroup DAE_Group + * Constrain solution component k. Possible values for + * 'flag' are: + * - c_NONE no constraint + * - c_GE_ZERO >= 0 + * - c_GT_ZERO > 0 + * - c_LE_ZERO <= 0 + * - c_LT_ZERO < 0 */ - class ResidEval { + virtual void constrain(const int k, const int flag) { m_constrain[k] = flag; } + int constraint(const int k) const { + std::map::const_iterator i = m_constrain.find(k); + if (i != m_constrain.end()) return i->second; + return c_NONE; + } - public: - - ResidEval() {} - virtual ~ResidEval() {} - - /** - * Constrain solution component k. Possible values for - * 'flag' are: - * - c_NONE no constraint - * - c_GE_ZERO >= 0 - * - c_GT_ZERO > 0 - * - c_LE_ZERO <= 0 - * - c_LT_ZERO < 0 - */ - virtual void constrain(int k, int flag) { m_constrain[k] = flag; } - int constraint(int k) { return m_constrain[k]; } - - /** - * Specify that solution component k is purely algebraic - - * that is, the derivative of this component does not appear - * in the residual function. - */ - virtual void setAlgebraic(int k) { m_alg[k] = 1; } - virtual bool isAlgebraic(int k) {return (m_alg[k] == 1); } + /** + * Specify that solution component k is purely algebraic - + * that is, the derivative of this component does not appear + * in the residual function. + */ + virtual void setAlgebraic(const int k) { m_alg[k] = 1; } + virtual bool isAlgebraic(const int k) {return (m_alg[k] == 1); } - /** - * Evaluate the residual function. Called by the - * integrator. - * @param t time. (input) - * @param y solution vector. (input) - * @param ydot rate of change of solution vector. (input) - * @param r residual vector (output) - */ - virtual int eval(double t, const double* y, - const double* ydot, double* r)=0; + /** + * Evaluate the residual function. Called by the + * integrator. + * @param t time. (input) + * @param y solution vector. (input) + * @param ydot rate of change of solution vector. (input) + * @param r residual vector (output) + */ + virtual int eval(const doublereal t, const doublereal * const y, + const doublereal * const ydot, + doublereal * const r) { + throw CanteraError("ResidEval::eval()", "base class called"); + } - /** - * Fill the solution and derivative vectors with the initial - * conditions at initial time t0. If these do not satisfy the - * residual equation, call one of the "corrrectInitial_xxx" - * methods before calling solve. - */ - virtual void getInitialConditions(double t0, double* y, - doublereal* ydot)=0; + /** + * Fill the solution and derivative vectors with the initial + * conditions at initial time t0. If these do not satisfy the + * residual equation, call one of the "corrrectInitial_xxx" + * methods before calling solve. + */ + virtual void getInitialConditions(const doublereal t0, doublereal * const y, + doublereal * const ydot) { + throw CanteraError("ResidEval::GetInitialConditions()", "base class called"); + } - /** - * Number of equations. - */ - virtual int nEquations()=0; + //! Return the number of equations in the equation system + virtual int nEquations() const = 0; - protected: + protected: - std::map m_alg; - std::map m_constrain; + std::map m_alg; + std::map m_constrain; - private: + private: - }; + }; } diff --git a/Cantera/src/numerics/ResidJacEval.cpp b/Cantera/src/numerics/ResidJacEval.cpp new file mode 100644 index 000000000..bb54b9423 --- /dev/null +++ b/Cantera/src/numerics/ResidJacEval.cpp @@ -0,0 +1,285 @@ +/** + * @file ResidJacEval.cpp + * + */ +/* + * $Revision$ + * $Date$ + */ +/* + * Copywrite 2004 Sandia Corporation. Under the terms of Contract + * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government + * retains certain rights in this software. + * See file License.txt for licensing information. + */ + +#include "ct_defs.h" +#include "ctlapack.h" +#include "ResidJacEval.h" + +#include +#include + +using namespace std; + +namespace Cantera { + + /************************************************************************* + * + * ResidJacEval(): + * + * Default constructor for the ResidJacEval class. + * + * atol has a default of 1.0E-13. + */ + ResidJacEval::ResidJacEval(doublereal atol) : + ResidEval(), + m_atol(atol) + { + } + + // Copy Constructor for the %ResidJacEval object + /* + */ + ResidJacEval::ResidJacEval(const ResidJacEval &right) : + ResidEval() + { + *this = operator=(right); + } + + /* + * + */ + ResidJacEval::~ResidJacEval() + { + } + + ResidJacEval& ResidJacEval::operator=(const ResidJacEval &right) { + if (this == &right) { + return *this; + } + + ResidEval::operator=(right); + + m_atol = right.m_atol; + neq_ = right.neq_; + + return *this; + } + + // Duplication routine for objects which inherit from + // %ResidJacEval + /* + * This virtual routine can be used to duplicate %ResidJacEval objects + * inherited from %ResidJacEval even if the application only has + * a pointer to %ResidJacEval to work with. + * + * These routines are basically wrappers around the derived copy + * constructor. + */ + ResidJacEval *ResidJacEval::duplMyselfAsResidJacEval() const { + ResidJacEval *ff = new ResidJacEval(*this); + return ff; + } + + int ResidJacEval::nEquations() const { + return neq_; + } + + /* + * + * setAtol(): + * + * Set the absolute tolerance value + */ + void ResidJacEval::setAtol(doublereal atol) + { + m_atol = atol; + if (m_atol <= 0.0) { + throw CanteraError("ResidJacEval::setAtol", + "atol must be greater than zero"); + } + } + + /************************************************************************** + * + * + * + * Fill the solution vector with the initial conditions + * at initial time t0. + */ + void ResidJacEval:: + getInitialConditionsDot(const doublereal t0, const size_t leny, + doublereal * const y, doublereal * const ydot) { + for (int i = 0; i < neq_; i++) { + y[i] = 0.0; + } + if (ydot) { + for (int i = 0; i < neq_; i++) { + ydot[i] = 0.0; + } + } + } + + /************************************************************************** + * + * + * + * Fill the solution vector with the initial conditions + * at initial time t0. + * + */ + void ResidJacEval:: + getInitialConditions(doublereal t0, size_t leny, + doublereal * y) { + getInitialConditionsDot(t0, leny, y, 0); + } + + /************************************************************************** + * + * user_out(): + * + * This function may be used to create output at various points in the + * execution of an application. + * + */ + void ResidJacEval:: + user_out2(const int ifunc, const doublereal t, const doublereal deltaT, + const doublereal *y, const doublereal *ydot) { + + } + + void ResidJacEval:: + user_out(const int ifunc, const doublereal t, + const doublereal *y, const doublereal *ydot) { + user_out2(ifunc, t, 0.0, y, ydot); + } + + /************************************************************************** + * + * + */ + void ResidJacEval:: + evalTimeTrackingEqns(const doublereal t, const doublereal deltaT, + const doublereal *y, + const doublereal *ydot) { + + } + + /******************************************************************** + * + * + * + * Return a vector of delta y's for calculation of the + * numerical Jacobian + */ + void ResidJacEval:: + calcDeltaSolnVariables(const doublereal t, + const doublereal * const ySoln, + const doublereal * const ySolnDot, + doublereal * const deltaYSoln, + const doublereal *const solnWeights) + { + if (!solnWeights) { + for (int i = 0; i < neq_; i++) { + deltaYSoln[i] = m_atol + fabs(1.0E-6 * ySoln[i]); + } + } else { + for (int i = 0; i < neq_; i++) { + deltaYSoln[i] = m_atol + + fmaxx(1.0E-2 * solnWeights[i], 1.0E-6 * fabs(ySoln[i])); + } + } + } + + /****************************************************************** + * + * calcSolnScales(): + * + * Returns a vector of ysolnScales[] that can be used to column scale + * Jacobians. + */ + void ResidJacEval:: + calcSolnScales(const doublereal t, + const doublereal * const ysoln, + const doublereal * const ysolnOld, + doublereal * const ysolnScales) + { + for (int i = 0; i < neq_; i++) { + ysolnScales[i] = 1.0; + } + } + + void ResidJacEval::filterSolnPrediction(doublereal t, + doublereal * const y) { + + } + + /************************************************************************** + * + * evalStoppingCriteria() + * + * If there is a stopping critera other than time set it here. + * + */ + bool ResidJacEval:: + evalStoppingCritera(doublereal &time_current, + doublereal &delta_t_n, + doublereal *y_n, + doublereal *ydot_n) + { + return false; + } + + /************************************************************************** + * + * matrixConditioning() + * + * Multiply the matrix by the inverse of a matrix which lead to a + * better conditioned system. The default, specified here, is to + * do nothing. + */ + void ResidJacEval:: + matrixConditioning(doublereal * const matrix, const int nrows, + doublereal * const rhs) + { + } + + /************************************************************************** + * + */ + void ResidJacEval:: + evalResidNJ(doublereal t, const doublereal deltaT, + const doublereal * y, + const doublereal * ydot, + doublereal * resid, + bool NJevaluation, + int id_x, + doublereal delta_x) + { + printf("Not implemented\n"); + std::exit(-1); + } + + /************************************************************************** + * + * evalJacobian() + * + * Calculate the jacobian and the residual at the current + * time and values. + * Backwards Euler is assumed. + */ + void ResidJacEval:: + evalJacobian(const doublereal t, const doublereal deltaT, + const doublereal * const y, + const doublereal * const ydot, + SquareMatrix &J, + doublereal * const resid) + { + printf("Not implemented\n"); + std::exit(-1); + } + + +} + diff --git a/Cantera/src/numerics/ResidJacEval.h b/Cantera/src/numerics/ResidJacEval.h new file mode 100644 index 000000000..eb33daabf --- /dev/null +++ b/Cantera/src/numerics/ResidJacEval.h @@ -0,0 +1,181 @@ +/** + * @file ResidJacEval.h + * + * Dense, Square (not sparse) matrices. + */ + +/* + * $Date$ + * $Revision$ + * + */ +/* + * Copywrite 2004 Sandia Corporation. Under the terms of Contract + * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government + * retains certain rights in this software. + * See file License.txt for licensing information. + */ + + +#ifndef CT_RESIDJACEVAL_H +#define CT_RESIDJACEVAL_H + +#include "ResidEval.h" +#include "SquareMatrix.h" + +namespace Cantera { + + /** + * A class for full (non-sparse) matrices with Fortran-compatible + * data storage. Adds matrix operations to class Array2D. + */ + class ResidJacEval : public ResidEval { + + public: + + /** + * Default constructor + */ + ResidJacEval(doublereal atol = 1.0e-13); + + //!Copy Constructor for the %ResidJacEval object + /*! + * @param right Item to be copied + */ + ResidJacEval(const ResidJacEval &right); + + /// Destructor. Does nothing. + virtual ~ResidJacEval(); + + //! Assignment operator + /*! + * This is NOT a virtual function. + * + * @param right Reference to %ResidJacEval object to be copied into the + * current one. + */ + ResidJacEval& ResidJacEval::operator=(const ResidJacEval &right); + + //! Duplication routine for objects which inherit from + //! residJacEval + /*! + * This virtual routine can be used to duplicate %ResidJacEval objects + * inherited from %ResidJacEval even if the application only has + * a pointer to %ResidJacEval to work with. + * + * These routines are basically wrappers around the derived copy + * constructor. + */ + virtual ResidJacEval *duplMyselfAsResidJacEval() const; + + //! Return the number of equations in the equation system + virtual int nEquations() const; + + /** + * Evaluate the residual function. + * @param t time (input, do not modify) + * @param y solution vector (input, do not modify) + * @param ydot rate of change of solution vector. (input, do + * not modify) + */ + virtual void evalResidNJ(doublereal t, const doublereal deltaT, + const doublereal * const y, + const doublereal * const ydot, + doublereal * const resid, + bool NJevaluation = false, + int id_x = 0, + doublereal delta_x = 0.0); + + /** + * Fill the solution vector with the initial conditions + * at initial time t0. + */ + virtual void getInitialConditionsDot(const doublereal t0, size_t leny, + doublereal * const y, + doublereal * const ydot); + + virtual void getInitialConditions(const doublereal t0, const size_t leny, + doublereal * const y); + + virtual void filterSolnPrediction(doublereal t, + doublereal * const y); + + + void setAtol(doublereal atol); + + virtual void evalTimeTrackingEqns(const doublereal t, const doublereal deltaT, + const doublereal * const y, + const doublereal * const ydot); + + virtual bool evalStoppingCritera(doublereal &time_current, + doublereal &delta_t_n, + doublereal *y_n, + doublereal *ydot_n); + /** + * Return a vector of delta y's for calculation of the + * numerical Jacobian + */ + virtual void + calcDeltaSolnVariables(const doublereal t, + const doublereal * const ysoln, + const doublereal * const ysolnDot, + doublereal * const deltaYsoln, + const doublereal * const solnWeights=0); + + /** + * Returns a vector of ysolnScales[] that can be used to column + * scale Jacobians. + */ + virtual void calcSolnScales(const doublereal t, + const doublereal * const ysoln, + const doublereal * const ysolnOld, + doublereal * const ysolnScales); + + /** + * This function may be used to create output at various points in the + * execution of an application. + * + */ + virtual void user_out2(const int ifunc, const doublereal t, + const doublereal deltaT, + const doublereal * const y, + const doublereal * const ydot); + + virtual void user_out(const int ifunc, const doublereal t, + const doublereal *y, + const doublereal *ydot); + + + virtual void matrixConditioning(doublereal * const matrix, const int nrows, + doublereal * const rhs); + + /********************************************************************* + * + * evalJacobian() + * + * Calculate the jacobian and the residual at the current + * time and values. + * Backwards Euler is assumed. + */ + virtual void evalJacobian(const doublereal t, const doublereal deltaT, + + const double* const y, + const double* const ydot, + SquareMatrix &J, + doublereal * const resid); + + + + protected: + + doublereal m_atol; + + int neq_; + + }; +} + +#endif + + +