diff --git a/include/cantera/oneD/MultiJac.h b/include/cantera/oneD/MultiJac.h index ca29f217ee..08a848f52f 100644 --- a/include/cantera/oneD/MultiJac.h +++ b/include/cantera/oneD/MultiJac.h @@ -23,6 +23,8 @@ namespace Cantera class MultiJac : public BandMatrix { public: + //! Constructor. + //! @param r The nonlinear system for which to compute the Jacobian. MultiJac(OneDim& r); /** @@ -36,9 +38,9 @@ class MultiJac : public BandMatrix * perturbing each variable, evaluating the residual function, and then * estimating the partial derivatives numerically using finite differences. * - * @param x0 Solution vector at which to evaluate the Jacobian - * @param resid0 Residual vector at x0 - * @param rdt Reciprocal of the time step + * @param x0 Solution vector at which to evaluate the Jacobian + * @param resid0 Residual vector at x0 + * @param rdt Reciprocal of the time step */ void eval(double* x0, double* resid0, double rdt); @@ -88,14 +90,18 @@ class MultiJac : public BandMatrix vector m_r1; //!< Perturbed residual vector double m_rtol = 1e-5; //!< Relative tolerance for perturbing solution components + + //! Absolute tolerance for perturbing solution components double m_atol = sqrt(std::numeric_limits::epsilon()); + double m_elapsed = 0.0; //!< Elapsed CPU time taken to compute the Jacobian vector m_ssdiag; //!< Diagonal of the steady-state Jacobian //! Transient mask for transient terms, 1 if transient, 0 if steady-state vector m_mask; - int m_nevals = 0; - int m_age = 100000; //!< Age of the Jacobian (times 'incrementAge' has been called) + + int m_nevals = 0; //!< Number of Jacobian evaluations. + int m_age = 100000; //!< Age of the Jacobian (times incrementAge() has been called) }; } diff --git a/include/cantera/oneD/MultiNewton.h b/include/cantera/oneD/MultiNewton.h index da2a3604e3..b90f622c9a 100644 --- a/include/cantera/oneD/MultiNewton.h +++ b/include/cantera/oneD/MultiNewton.h @@ -23,11 +23,14 @@ namespace Cantera class MultiNewton { public: + //! Constructor + //! @param sz Number of variables in the system MultiNewton(int sz); virtual ~MultiNewton() {}; MultiNewton(const MultiNewton&) = delete; MultiNewton& operator=(const MultiNewton&) = delete; + //! Get the number of variables in the system. size_t size() { return m_n; } @@ -89,15 +92,16 @@ class MultiNewton * where @f$ \Delta x_k = J^{-1}(x_k) F(x_k) @f$, and * @f$ \Delta x_{k+1} = J^{-1}(x_{k}) F(x_{k+1}) @f$. * - * @param[in] x0 initial solution about which a Newton step will be taken - * @param[in] step0 initial undamped Newton step - * @param[out] x1 solution after taking the damped Newton step - * @param[out] step1 Newton step after taking the damped Newton step - * @param[out] s1 norm of the subsequent Newton step after taking the damped Newton step - * @param[in] r domain object, used for evaluating residuals over all domains - * @param[in] jac Jacobian evaluator - * @param[in] loglevel controls amount of printed diagnostics - * @param[in] writetitle controls if logging title is printed + * @param[in] x0 initial solution about which a Newton step will be taken + * @param[in] step0 initial undamped Newton step + * @param[out] x1 solution after taking the damped Newton step + * @param[out] step1 Newton step after taking the damped Newton step + * @param[out] s1 norm of the subsequent Newton step after taking the damped Newton + * step + * @param[in] r domain object, used for evaluating residuals over all domains + * @param[in] jac Jacobian evaluator + * @param[in] loglevel controls amount of printed diagnostics + * @param[in] writetitle controls if logging title is printed * * @returns * - `1` a damping coefficient was found and the solution converges. @@ -129,6 +133,8 @@ class MultiNewton int solve(double* x0, double* x1, OneDim& r, MultiJac& jac, int loglevel); //! Set options. + //! @param maxJacAge Maximum number of steps that can be taken before requiring + //! a Jacobian update void setOptions(int maxJacAge = 5) { m_maxAge = maxJacAge; } @@ -137,8 +143,14 @@ class MultiNewton void resize(size_t points); protected: - //! Work arrays of size #m_n used in solve(). - vector m_x, m_stp, m_stp1; + //! Work array holding the system state after the last successful step. Size #m_n. + vector m_x; + + //! Work array holding the undamped Newton step or the system residual. Size #m_n. + vector m_stp; + + //! Work array holding the damped Newton step. Size #m_n. + vector m_stp1; //! Maximum allowable Jacobian age before it is recomputed. int m_maxAge = 5; diff --git a/include/cantera/oneD/OneDim.h b/include/cantera/oneD/OneDim.h index c1dd6b5827..1a73ab74aa 100644 --- a/include/cantera/oneD/OneDim.h +++ b/include/cantera/oneD/OneDim.h @@ -26,6 +26,7 @@ class AnyMap; class OneDim { public: + //! Default constructor OneDim(); //! Construct a OneDim container for the domains in the list *domains*. @@ -69,6 +70,7 @@ class OneDim return *m_dom[i]; } + //! Get the index of the domain named `name`. size_t domainIndex(const string& name); //! Check that the specified domain index is in range. @@ -116,13 +118,12 @@ class OneDim return m_dom.back().get(); } - //! Number of solution components at global point jg. + //! Number of solution components at global point `jg`. size_t nVars(size_t jg) { return m_nvars[jg]; } - //! Location in the solution vector of the first component of global point - //! jg. + //! Location in the solution vector of the first component of global point `jg`. size_t loc(size_t jg) { return m_loc[jg]; } @@ -204,6 +205,8 @@ class OneDim //! Call after one or more grids has changed size, for example after being refined. virtual void resize(); + //! Access the vector indicating which equations contain a transient term. + //! Elements are 1 for equations with a transient terms and 0 otherwise. vector& transientMask() { return m_mask; } @@ -211,15 +214,17 @@ class OneDim /** * Take time steps using Backward Euler. * - * @param nsteps number of steps - * @param dt initial step size - * @param x current solution vector - * @param r solution vector after time stepping - * @param loglevel controls amount of printed diagnostics + * @param nsteps number of steps + * @param dt initial step size + * @param x current solution vector + * @param r solution vector after time stepping + * @param loglevel controls amount of printed diagnostics * @returns size of last timestep taken */ double timeStep(int nsteps, double dt, double* x, double* r, int loglevel); + //! Reset values such as negative species concentrations in each domain. + //! @see Domain1D::resetBadValues void resetBadValues(double* x); //! Write statistics about the number of iterations and Jacobians at each @@ -231,10 +236,15 @@ class OneDim */ void writeStats(int printTime = 1); - // options + //! @name Options + //! @{ + + //! Set the minimum time step allowed during time stepping void setMinTimeStep(double tmin) { m_tmin = tmin; } + + //! Set the maximum time step allowed during time stepping void setMaxTimeStep(double tmax) { m_tmax = tmax; } @@ -243,7 +253,7 @@ class OneDim * Sets a factor by which the time step is reduced if the time stepping * fails. The default value is 0.5. * - * @param tfactor factor time step is multiplied by if time stepping fails + * @param tfactor factor time step is multiplied by if time stepping fails */ void setTimeStepFactor(double tfactor) { m_tfactor = tfactor; @@ -260,7 +270,13 @@ class OneDim int maxTimeStepCount() const { return m_nsteps_max; } + //! @} + //! Set the maximum number of steps that can be taken using the same Jacobian + //! before it must be re-evaluated. + //! @param ss_age Age limit during steady-state mode + //! @param ts_age Age limit during time stepping mode. If not specified, the + //! steady-state age is also used during time stepping. void setJacAge(int ss_age, int ts_age=-1); /** @@ -334,7 +350,10 @@ class OneDim } protected: - void evalSSJacobian(double* x, double* xnew); + //! Evaluate the steady-state Jacobian, accessible via jacobian() + //! @param[in] x Current state vector, length size() + //! @param[out] rsd Storage for the residual, length size() + void evalSSJacobian(double* x, double* rsd); double m_tmin = 1e-16; //!< minimum timestep size double m_tmax = 1e+08; //!< maximum timestep size @@ -352,19 +371,34 @@ class OneDim size_t m_bw = 0; //!< Jacobian bandwidth size_t m_size = 0; //!< solution vector size + //! All domains comprising the system vector> m_dom; + + //! All connector and boundary domains vector> m_connect; + + //! All bulk/flow domains vector> m_bulk; + //! Indicates whether one-time initialization for each domain has been completed. bool m_init = false; + + //! Number of variables at each point, across all domains. Length points(). + //! Accessed with nVars(). vector m_nvars; + + //! Location in the state vector of the first component of each point, across all + //! domains. Accessed with loc(). vector m_loc; + + //! Transient mask. See transientMask(). vector m_mask; + + //! Total number of points. size_t m_pts = 0; - // options - int m_ss_jac_age = 20; - int m_ts_jac_age = 20; + int m_ss_jac_age = 20; //!< Maximum age of the Jacobian in steady-state mode. + int m_ts_jac_age = 20; //!< Maximum age of the Jacobian in time-stepping mode. //! Function called at the start of every call to #eval. Func1* m_interrupt = nullptr; @@ -379,18 +413,29 @@ class OneDim int m_nsteps_max = 500; private: - // statistics - int m_nevals = 0; - double m_evaltime = 0; - vector m_gridpts; - vector m_jacEvals; - vector m_jacElapsed; + //! @name Statistics + //! Solver stats are collected after successfully solving on a particular grid. + //! @{ + int m_nevals = 0; //!< Number of calls to eval() + double m_evaltime = 0; //!< Total time [s] spent in eval() + + vector m_gridpts; //!< Number of grid points in this grid + vector m_jacEvals; //!< Number of Jacobian evaluations on this grid + vector m_jacElapsed; //!< Time [s] spent evaluating Jacobians on this grid + + //! Number of residual function evaluations on this grid (not counting evaluations + //! used to construct Jacobians). vector m_funcEvals; + + //! Time [s] spent on residual function evaluations on this grid (not counting + //! evaluations used to construct Jacobians). vector m_funcElapsed; //! Number of time steps taken in each call to solve() (for example, for each //! successive grid refinement) vector m_timeSteps; + + //! @} }; } diff --git a/include/cantera/oneD/Sim1D.h b/include/cantera/oneD/Sim1D.h index 8a9a8d1169..2448bd5ada 100644 --- a/include/cantera/oneD/Sim1D.h +++ b/include/cantera/oneD/Sim1D.h @@ -30,7 +30,7 @@ class Sim1D : public OneDim /** * Standard constructor. - * @param domains A vector of shared pointers to the domains to be linked together. + * @param domains A vector of shared pointers to the domains to be linked together. * The domain pointers must be entered in left-to-right order --- that is, * the pointer to the leftmost domain is domain[0], the pointer to the * domain to its right is domain[1], etc. @@ -44,10 +44,10 @@ class Sim1D : public OneDim //! Set initial guess for one component for all domains /** - * @param component component name - * @param locs A vector of relative positions, beginning with 0.0 at the + * @param component component name + * @param locs A vector of relative positions, beginning with 0.0 at the * left of the domain, and ending with 1.0 at the right of the domain. - * @param vals A vector of values corresponding to the relative position + * @param vals A vector of values corresponding to the relative position * locations. */ void setInitialGuess(const string& component, vector& locs, @@ -55,32 +55,37 @@ class Sim1D : public OneDim /** * Set a single value in the solution vector. - * @param dom domain number, beginning with 0 for the leftmost domain. - * @param comp component number - * @param localPoint grid point within the domain, beginning with 0 for + * @param dom domain number, beginning with 0 for the leftmost domain. + * @param comp component number + * @param localPoint grid point within the domain, beginning with 0 for * the leftmost grid point in the domain. - * @param value the value. + * @param value the value. */ void setValue(size_t dom, size_t comp, size_t localPoint, double value); /** * Get one entry in the solution vector. - * @param dom domain number, beginning with 0 for the leftmost domain. - * @param comp component number - * @param localPoint grid point within the domain, beginning with 0 for + * @param dom domain number, beginning with 0 for the leftmost domain. + * @param comp component number + * @param localPoint grid point within the domain, beginning with 0 for * the leftmost grid point in the domain. */ double value(size_t dom, size_t comp, size_t localPoint) const; + //! Get an entry in the work vector, which may contain either a new system state + //! or the current residual of the system. + //! @param dom domain index + //! @param comp component index + //! @param localPoint grid point within the domain double workValue(size_t dom, size_t comp, size_t localPoint) const; /** * Specify a profile for one component of one domain. - * @param dom domain number, beginning with 0 for the leftmost domain. - * @param comp component number - * @param pos A vector of relative positions, beginning with 0.0 at the + * @param dom domain number, beginning with 0 for the leftmost domain. + * @param comp component number + * @param pos A vector of relative positions, beginning with 0.0 at the * left of the domain, and ending with 1.0 at the right of the domain. - * @param values A vector of values corresponding to the relative position + * @param values A vector of values corresponding to the relative position * locations. * * Note that the vector pos and values can have lengths different than the @@ -202,6 +207,13 @@ class Sim1D : public OneDim //! @} + //! Set the number of time steps to try when the steady Newton solver is + //! unsuccessful. + //! @param stepsize Initial time step size [s] + //! @param n Length of `tsteps` array + //! @param tsteps A sequence of time step counts to take after subsequent failures + //! of the steady-state solver. The last value in `tsteps` will be used again + //! after further unsuccessful solution attempts. void setTimeStep(double stepsize, size_t n, const int* tsteps); /** @@ -301,15 +313,15 @@ class Sim1D : public OneDim /** * Get the maximum number of grid points in this domain. @see Refiner::maxPoints * - * @param dom domain number, beginning with 0 for the leftmost domain. + * @param dom domain number, beginning with 0 for the leftmost domain. */ size_t maxGridPoints(size_t dom); //! Set the minimum grid spacing in the specified domain(s). /*! - * @param dom Domain index. If dom == -1, the specified spacing is applied - * to all domains. - * @param gridmin The minimum allowable grid spacing [m] + * @param dom Domain index. If dom == -1, the specified spacing is applied + * to all domains. + * @param gridmin The minimum allowable grid spacing [m] */ void setGridMin(int dom, double gridmin); @@ -323,10 +335,13 @@ class Sim1D : public OneDim //! failure during grid refinement. void restoreSteadySolution(); + //! Get the initial value of the system state from each domain in the simulation. void getInitialSoln(); + //! Get the Jacobian element @f$ J_{ij} = \partial f_i / \partial x_j \f$ double jacobian(int i, int j); + //! Evaluate the Jacobian in steady-state mode. void evalSSJacobian(); //! Solve the equation @f$ J^T \lambda = b @f$. diff --git a/src/oneD/OneDim.cpp b/src/oneD/OneDim.cpp index a7c9b0071d..6e8290e393 100644 --- a/src/oneD/OneDim.cpp +++ b/src/oneD/OneDim.cpp @@ -216,13 +216,13 @@ int OneDim::solve(double* x, double* xnew, int loglevel) return m_newt->solve(x, xnew, *this, *m_jac, loglevel); } -void OneDim::evalSSJacobian(double* x, double* xnew) +void OneDim::evalSSJacobian(double* x, double* rsd) { double rdt_save = m_rdt; m_jac_ok = false; setSteadyMode(); - eval(npos, x, xnew, 0.0, 0); - m_jac->eval(x, xnew, 0.0); + eval(npos, x, rsd, 0.0, 0); + m_jac->eval(x, rsd, 0.0); m_rdt = rdt_save; }