Skip to content

Commit

Permalink
[Doc] Document remaining members of MultiJac, OneDim, and Sim1D
Browse files Browse the repository at this point in the history
  • Loading branch information
@speth @ischoegl
speth authored and ischoegl committed Oct 28, 2024
1 parent 6e5f848 commit 3ba21c4
Show file tree
Hide file tree
Showing 5 changed files with 136 additions and 58 deletions.
16 changes: 11 additions & 5 deletions include/cantera/oneD/MultiJac.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);

/**
Expand All @@ -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);

Expand Down Expand Up @@ -88,14 +90,18 @@ class MultiJac : public BandMatrix

vector<double> 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<double>::epsilon());

double m_elapsed = 0.0; //!< Elapsed CPU time taken to compute the Jacobian
vector<double> m_ssdiag; //!< Diagonal of the steady-state Jacobian

//! Transient mask for transient terms, 1 if transient, 0 if steady-state
vector<int> 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)
};
}

Expand Down
34 changes: 23 additions & 11 deletions include/cantera/oneD/MultiNewton.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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;
}
Expand All @@ -137,8 +143,14 @@ class MultiNewton
void resize(size_t points);

protected:
//! Work arrays of size #m_n used in solve().
vector<double> m_x, m_stp, m_stp1;
//! Work array holding the system state after the last successful step. Size #m_n.
vector<double> m_x;

//! Work array holding the undamped Newton step or the system residual. Size #m_n.
vector<double> m_stp;

//! Work array holding the damped Newton step. Size #m_n.
vector<double> m_stp1;

//! Maximum allowable Jacobian age before it is recomputed.
int m_maxAge = 5;
Expand Down
85 changes: 65 additions & 20 deletions include/cantera/oneD/OneDim.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ class AnyMap;
class OneDim
{
public:
//! Default constructor
OneDim();

//! Construct a OneDim container for the domains in the list *domains*.
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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];
}
Expand Down Expand Up @@ -204,22 +205,26 @@ 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<int>& transientMask() {
return m_mask;
}

/**
* 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
Expand All @@ -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;
}
Expand All @@ -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;
Expand All @@ -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);

/**
Expand Down Expand Up @@ -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
Expand All @@ -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<shared_ptr<Domain1D>> m_dom;

//! All connector and boundary domains
vector<shared_ptr<Domain1D>> m_connect;

//! All bulk/flow domains
vector<shared_ptr<Domain1D>> 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<size_t> m_nvars;

//! Location in the state vector of the first component of each point, across all
//! domains. Accessed with loc().
vector<size_t> m_loc;

//! Transient mask. See transientMask().
vector<int> 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;
Expand All @@ -379,18 +413,29 @@ class OneDim
int m_nsteps_max = 500;

private:
// statistics
int m_nevals = 0;
double m_evaltime = 0;
vector<size_t> m_gridpts;
vector<int> m_jacEvals;
vector<double> 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<size_t> m_gridpts; //!< Number of grid points in this grid
vector<int> m_jacEvals; //!< Number of Jacobian evaluations on this grid
vector<double> 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<int> m_funcEvals;

//! Time [s] spent on residual function evaluations on this grid (not counting
//! evaluations used to construct Jacobians).
vector<double> m_funcElapsed;

//! Number of time steps taken in each call to solve() (for example, for each
//! successive grid refinement)
vector<int> m_timeSteps;

//! @}
};

}
Expand Down
Loading

0 comments on commit 3ba21c4

Please sign in to comment.