Cantera  2.3.0
Sim1D.h
Go to the documentation of this file.
1 /**
2  * @file Sim1D.h
3  */
4 
5 // This file is part of Cantera. See License.txt in the top-level directory or
6 // at http://www.cantera.org/license.txt for license and copyright information.
7 
8 #ifndef CT_SIM1D_H
9 #define CT_SIM1D_H
10 
11 #include "OneDim.h"
12 
13 namespace Cantera
14 {
15 
16 /**
17  * One-dimensional simulations. Class Sim1D extends class OneDim by storing
18  * the solution vector, and by adding a hybrid Newton/time-stepping solver.
19  * @ingroup onedim
20  */
21 class Sim1D : public OneDim
22 {
23 public:
24  //! Default constructor.
25  /*!
26  * This constructor is provided to make the class default-constructible, but
27  * is not meant to be used in most applications. Use the next constructor
28  */
29  Sim1D() {}
30 
31  /**
32  * Standard constructor.
33  * @param domains A vector of pointers to the domains to be linked together.
34  * The domain pointers must be entered in left-to-right order --- i.e.,
35  * the pointer to the leftmost domain is domain[0], the pointer to the
36  * domain to its right is domain[1], etc.
37  */
38  Sim1D(std::vector<Domain1D*>& domains);
39 
40  /**
41  * @name Setting initial values
42  *
43  * These methods are used to set the initial values of solution components.
44  */
45  //@{
46 
47  /// Set initial guess for one component for all domains
48  /**
49  * @param component component name
50  * @param locs A vector of relative positions, beginning with 0.0 at the
51  * left of the domain, and ending with 1.0 at the right of the domain.
52  * @param vals A vector of values corresponding to the relative position
53  * locations.
54  */
55  void setInitialGuess(const std::string& component, vector_fp& locs,
56  vector_fp& vals);
57 
58  /**
59  * Set a single value in the solution vector.
60  * @param dom domain number, beginning with 0 for the leftmost domain.
61  * @param comp component number
62  * @param localPoint grid point within the domain, beginning with 0 for
63  * the leftmost grid point in the domain.
64  * @param value the value.
65  */
66  void setValue(size_t dom, size_t comp, size_t localPoint, doublereal value);
67 
68  /**
69  * Get one entry in the solution vector.
70  * @param dom domain number, beginning with 0 for the leftmost domain.
71  * @param comp component number
72  * @param localPoint grid point within the domain, beginning with 0 for
73  * the leftmost grid point in the domain.
74  */
75  doublereal value(size_t dom, size_t comp, size_t localPoint) const;
76 
77  doublereal workValue(size_t dom, size_t comp, size_t localPoint) const;
78 
79  /**
80  * Specify a profile for one component of one domain.
81  * @param dom domain number, beginning with 0 for the leftmost domain.
82  * @param comp component number
83  * @param pos A vector of relative positions, beginning with 0.0 at the
84  * left of the domain, and ending with 1.0 at the right of the domain.
85  * @param values A vector of values corresponding to the relative position
86  * locations.
87  *
88  * Note that the vector pos and values can have lengths different than the
89  * number of grid points, but their lengths must be equal. The values at
90  * the grid points will be linearly interpolated based on the (pos,
91  * values) specification.
92  */
93  void setProfile(size_t dom, size_t comp, const vector_fp& pos,
94  const vector_fp& values);
95 
96  /// Set component 'comp' of domain 'dom' to value 'v' at all points.
97  void setFlatProfile(size_t dom, size_t comp, doublereal v);
98 
99  //@}
100 
101  void save(const std::string& fname, const std::string& id,
102  const std::string& desc, int loglevel=1);
103 
104  void saveResidual(const std::string& fname, const std::string& id,
105  const std::string& desc, int loglevel=1);
106 
107  /// Print to stream s the current solution for all domains.
108  void showSolution(std::ostream& s);
109  void showSolution();
110 
111  const doublereal* solution() {
112  return m_x.data();
113  }
114 
115  void setTimeStep(double stepsize, size_t n, const int* tsteps);
116 
117  void solve(int loglevel = 0, bool refine_grid = true);
118 
119  void eval(doublereal rdt=-1.0, int count = 1) {
120  OneDim::eval(npos, m_x.data(), m_xnew.data(), rdt, count);
121  }
122 
123  // Evaluate the governing equations and return the vector of residuals
124  void getResidual(double rdt, double* resid) {
125  OneDim::eval(npos, m_x.data(), resid, rdt, 0);
126  }
127 
128  /// Refine the grid in all domains.
129  int refine(int loglevel=0);
130 
131  //! Add node for fixed temperature point of freely propagating flame
132  int setFixedTemperature(doublereal t);
133 
134  /**
135  * Set grid refinement criteria. If dom >= 0, then the settings
136  * apply only to the specified domain. If dom < 0, the settings
137  * are applied to each domain. @see Refiner::setCriteria.
138  */
139  void setRefineCriteria(int dom = -1, doublereal ratio = 10.0,
140  doublereal slope = 0.8, doublereal curve = 0.8, doublereal prune = -0.1);
141 
142  /**
143  * Set the maximum number of grid points in the domain. If dom >= 0,
144  * then the settings apply only to the specified domain. If dom < 0,
145  * the settings are applied to each domain. @see Refiner::setMaxPoints.
146  */
147  void setMaxGridPoints(int dom, int npoints);
148 
149  /**
150  * Get the maximum number of grid points in this domain. @see Refiner::maxPoints
151  *
152  * @param dom domain number, beginning with 0 for the leftmost domain.
153  */
154  size_t maxGridPoints(size_t dom);
155 
156  //! Set the minimum grid spacing in the specified domain(s).
157  /*!
158  * @param dom Domain index. If dom == -1, the specified spacing is applied
159  * to all domains.
160  * @param gridmin The minimum allowable grid spacing [m]
161  */
162  void setGridMin(int dom, double gridmin);
163 
164  //! Initialize the solution with a previously-saved solution.
165  void restore(const std::string& fname, const std::string& id, int loglevel=2);
166 
167  //! Set the current solution vector to the last successful time-stepping
168  //! solution. This can be used to examine the solver progress after a failed
169  //! integration.
171 
172  //! Set the current solution vector and grid to the last successful steady-
173  //! state solution. This can be used to examine the solver progress after a
174  //! failure during grid refinement.
175  void restoreSteadySolution();
176 
177  void getInitialSoln();
178 
179  void setSolution(const doublereal* soln) {
180  std::copy(soln, soln + m_x.size(), m_x.data());
181  }
182 
183  const doublereal* solution() const {
184  return m_x.data();
185  }
186 
187  doublereal jacobian(int i, int j);
188 
189  void evalSSJacobian();
190 
191  //! Solve the equation \f$ J^T \lambda = b \f$.
192  /**
193  * Here, \f$ J = \partial f/\partial x \f$ is the Jacobian matrix of the
194  * system of equations \f$ f(x,p)=0 \f$. This can be used to efficiently
195  * solve for the sensitivities of a scalar objective function \f$ g(x,p) \f$
196  * to a vector of parameters \f$ p \f$ by solving:
197  * \f[ J^T \lambda = \left( \frac{\partial g}{\partial x} \right)^T \f]
198  * for \f$ \lambda \f$ and then computing:
199  * \f[
200  * \left.\frac{dg}{dp}\right|_{f=0} = \frac{\partial g}{\partial p}
201  * - \lambda^T \frac{\partial f}{\partial p}
202  * \f]
203  */
204  void solveAdjoint(const double* b, double* lambda);
205 
206  virtual void resize();
207 
208  //! Set a function that will be called after each successful steady-state
209  //! solve, before regridding. Intended to be used for observing solver
210  //! progress for debugging purposes.
211  void setSteadyCallback(Func1* callback) {
212  m_steady_callback = callback;
213  }
214 
215 protected:
216  //! the solution vector
218 
219  //! the solution vector after the last successful timestepping
221 
222  //! the solution vector after the last successful steady-state solve (stored
223  //! before grid refinement)
225 
226  //! the grids for each domain after the last successful steady-state solve
227  //! (stored before grid refinement)
228  std::vector<vector_fp> m_grid_last_ss;
229 
230  //! a work array used to hold the residual or the new solution
232 
233  //! timestep
234  doublereal m_tstep;
235 
236  //! array of number of steps to take before re-attempting the steady-state
237  //! solution
239 
240  //! User-supplied function called after a successful steady-state solve.
242 
243 private:
244  /// Calls method _finalize in each domain.
245  void finalize();
246 
247  //! Wrapper around the Newton solver
248  /*!
249  * @return 0 if successful, -1 on failure
250  */
251  int newtonSolve(int loglevel);
252 };
253 
254 }
255 #endif
Container class for multiple-domain 1D problems.
Definition: OneDim.h:25
std::vector< vector_fp > m_grid_last_ss
the grids for each domain after the last successful steady-state solve (stored before grid refinement...
Definition: Sim1D.h:228
size_t maxGridPoints(size_t dom)
Get the maximum number of grid points in this domain.
Definition: Sim1D.cpp:554
void restore(const std::string &fname, const std::string &id, int loglevel=2)
Initialize the solution with a previously-saved solution.
Definition: Sim1D.cpp:107
vector_fp m_x
the solution vector
Definition: Sim1D.h:217
vector_fp m_xlast_ts
the solution vector after the last successful timestepping
Definition: Sim1D.h:220
void eval(size_t j, double *x, double *r, doublereal rdt=-1.0, int count=1)
Evaluate the multi-domain residual function.
Definition: OneDim.cpp:258
vector_fp m_xlast_ss
the solution vector after the last successful steady-state solve (stored before grid refinement) ...
Definition: Sim1D.h:224
MultiJac & jacobian()
Return a reference to the Jacobian evaluator.
Definition: OneDim.cpp:101
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
void setValue(size_t dom, size_t comp, size_t localPoint, doublereal value)
Set a single value in the solution vector.
Definition: Sim1D.cpp:50
void setRefineCriteria(int dom=-1, doublereal ratio=10.0, doublereal slope=0.8, doublereal curve=0.8, doublereal prune=-0.1)
Set grid refinement criteria.
Definition: Sim1D.cpp:514
void setInitialGuess(const std::string &component, vector_fp &locs, vector_fp &vals)
Set initial guess for one component for all domains.
Definition: Sim1D.cpp:37
void solveAdjoint(const double *b, double *lambda)
Solve the equation .
Definition: Sim1D.cpp:570
vector_int m_steps
array of number of steps to take before re-attempting the steady-state solution
Definition: Sim1D.h:238
int newtonSolve(int loglevel)
Wrapper around the Newton solver.
Definition: Sim1D.cpp:213
std::vector< int > vector_int
Vector of ints.
Definition: ct_defs.h:159
int setFixedTemperature(doublereal t)
Add node for fixed temperature point of freely propagating flame.
Definition: Sim1D.cpp:427
One-dimensional simulations.
Definition: Sim1D.h:21
void finalize()
Calls method _finalize in each domain.
Definition: Sim1D.cpp:197
void restoreTimeSteppingSolution()
Set the current solution vector to the last successful time-stepping solution.
Definition: Sim1D.cpp:168
void setSteadyCallback(Func1 *callback)
Set a function that will be called after each successful steady-state solve, before regridding...
Definition: Sim1D.h:211
virtual void resize()
Call after one or more grids has changed size, e.g. after being refined.
Definition: Sim1D.cpp:594
std::tuple< std::string, size_t, std::string > component(size_t i)
Return the domain, local point index, and component name for the i-th component of the global solutio...
Definition: OneDim.cpp:65
Func1 * m_steady_callback
User-supplied function called after a successful steady-state solve.
Definition: Sim1D.h:241
Base class for &#39;functor&#39; classes that evaluate a function of one variable.
Definition: Func1.h:41
doublereal m_tstep
timestep
Definition: Sim1D.h:234
void setFlatProfile(size_t dom, size_t comp, doublereal v)
Set component &#39;comp&#39; of domain &#39;dom&#39; to value &#39;v&#39; at all points.
Definition: Sim1D.cpp:140
Sim1D()
Default constructor.
Definition: Sim1D.h:29
vector_fp m_xnew
a work array used to hold the residual or the new solution
Definition: Sim1D.h:231
doublereal value(size_t dom, size_t comp, size_t localPoint) const
Get one entry in the solution vector.
Definition: Sim1D.cpp:58
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:157
int refine(int loglevel=0)
Refine the grid in all domains.
Definition: Sim1D.cpp:341
void restoreSteadySolution()
Set the current solution vector and grid to the last successful steady- state solution.
Definition: Sim1D.cpp:177
void setProfile(size_t dom, size_t comp, const vector_fp &pos, const vector_fp &values)
Specify a profile for one component of one domain.
Definition: Sim1D.cpp:74
doublereal rdt() const
Reciprocal of the time step.
Definition: OneDim.h:140
void setMaxGridPoints(int dom, int npoints)
Set the maximum number of grid points in the domain.
Definition: Sim1D.cpp:541
Namespace for the Cantera kernel.
Definition: application.cpp:29
void setGridMin(int dom, double gridmin)
Set the minimum grid spacing in the specified domain(s).
Definition: Sim1D.cpp:528