Cantera  3.1.0a1
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 https://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 onedGroup
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 shared pointers to the domains to be linked together.
34  * The domain pointers must be entered in left-to-right order --- that is,
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(vector<shared_ptr<Domain1D>>& domains);
39 
40  //! @name Setting initial values
41  //!
42  //! These methods are used to set the initial values of solution components.
43  //! @{
44 
45  //! Set initial guess for one component for all domains
46  /**
47  * @param component component name
48  * @param locs A vector of relative positions, beginning with 0.0 at the
49  * left of the domain, and ending with 1.0 at the right of the domain.
50  * @param vals A vector of values corresponding to the relative position
51  * locations.
52  */
53  void setInitialGuess(const string& component, vector<double>& locs,
54  vector<double>& vals);
55 
56  /**
57  * Set a single value in the solution vector.
58  * @param dom domain number, beginning with 0 for the leftmost domain.
59  * @param comp component number
60  * @param localPoint grid point within the domain, beginning with 0 for
61  * the leftmost grid point in the domain.
62  * @param value the value.
63  */
64  void setValue(size_t dom, size_t comp, size_t localPoint, double value);
65 
66  /**
67  * Get one entry in the solution vector.
68  * @param dom domain number, beginning with 0 for the leftmost domain.
69  * @param comp component number
70  * @param localPoint grid point within the domain, beginning with 0 for
71  * the leftmost grid point in the domain.
72  */
73  double value(size_t dom, size_t comp, size_t localPoint) const;
74 
75  double workValue(size_t dom, size_t comp, size_t localPoint) const;
76 
77  /**
78  * Specify a profile for one component of one domain.
79  * @param dom domain number, beginning with 0 for the leftmost domain.
80  * @param comp component number
81  * @param pos A vector of relative positions, beginning with 0.0 at the
82  * left of the domain, and ending with 1.0 at the right of the domain.
83  * @param values A vector of values corresponding to the relative position
84  * locations.
85  *
86  * Note that the vector pos and values can have lengths different than the
87  * number of grid points, but their lengths must be equal. The values at
88  * the grid points will be linearly interpolated based on the (pos,
89  * values) specification.
90  */
91  void setProfile(size_t dom, size_t comp, const vector<double>& pos,
92  const vector<double>& values);
93 
94  //! Set component 'comp' of domain 'dom' to value 'v' at all points.
95  void setFlatProfile(size_t dom, size_t comp, double v);
96 
97  //! @}
98 
99  //! @name Logging, saving and restoring of solutions
100  //!
101  //! @{
102 
103  /**
104  * Output information on current solution for all domains to stream.
105  * @param s Output stream
106  * @since New in %Cantera 3.0.
107  */
108  void show(std::ostream& s);
109 
110  /**
111  * Show logging information on current solution for all domains.
112  * @since New in %Cantera 3.0.
113  */
114  void show();
115 
116  /**
117  * Save current simulation data to a container file or CSV format.
118  *
119  * In order to save the content of a Sim1D object, individual domains are
120  * converted to SolutionArray objects and saved using the SolutionArray::save()
121  * method. For HDF and YAML output, all domains are written to a single container
122  * file with shared header information. Simulation settings of individual domains
123  * are preserved as meta data of the corresponding SolutionArray objects.
124  * For CSV files, only state and auxiliary data of the main 1D domain are saved.
125  *
126  * The complete state of the current object can be restored from HDF and YAML
127  * container files using the restore() method, while individual domains can be
128  * loaded using SolutionArray::restore() for further analysis. While CSV do not
129  * contain complete information, they can still be used for setting initial states
130  * of individual simulation objects for some %Cantera API's.
131  *
132  * @param fname Name of output file (CSV, YAML or HDF)
133  * @param name Identifier of storage location within the container file; this
134  * node/group contains header information and multiple subgroups holding
135  * domain-specific SolutionArray data (YAML/HDF only)
136  * @param desc Custom comment describing the dataset to be stored (YAML/HDF only)
137  * @param overwrite Force overwrite if file/name exists; optional (default=false)
138  * @param compression Compression level (0-9); optional (default=0; HDF only)
139  * @param basis Output mass ("Y"/"mass") or mole ("X"/"mole") fractions;
140  * if not specified (default=""), the native basis of the underlying
141  * ThermoPhase manager is used - @see nativeState (CSV only)
142  */
143  void save(const string& fname, const string& name, const string& desc,
144  bool overwrite=false, int compression=0, const string& basis="");
145 
146  /**
147  * Save the residual of the current solution to a container file.
148  * @param fname Name of output container file
149  * @param name Identifier of solution within the container file
150  * @param desc Description of the solution
151  * @param overwrite Force overwrite if name exists; optional (default=false)
152  * @param compression Compression level (optional; HDF only)
153  */
154  void saveResidual(const string& fname, const string& name,
155  const string& desc, bool overwrite=false, int compression=0);
156 
157  /**
158  * Retrieve data and settings from a previously saved simulation.
159  *
160  * This method restores a simulation object from YAML or HDF data previously saved
161  * using the save() method.
162  *
163  * @param fname Name of container file (YAML or HDF)
164  * @param name Identifier of location within the container file; this node/group
165  * contains header information and subgroups with domain-specific SolutionArray
166  * data
167  * @return AnyMap containing header information
168  */
169  AnyMap restore(const string& fname, const string& name);
170 
171  //! @}
172 
173  void setTimeStep(double stepsize, size_t n, const int* tsteps);
174 
175  void solve(int loglevel = 0, bool refine_grid = true);
176 
177  void eval(double rdt=-1.0, int count = 1) {
178  OneDim::eval(npos, m_state->data(), m_xnew.data(), rdt, count);
179  }
180 
181  // Evaluate the governing equations and return the vector of residuals
182  void getResidual(double rdt, double* resid) {
183  OneDim::eval(npos, m_state->data(), resid, rdt, 0);
184  }
185 
186  //! Refine the grid in all domains.
187  int refine(int loglevel=0);
188 
189  //! Add node for fixed temperature point of freely propagating flame
190  int setFixedTemperature(double t);
191 
192  //! Return temperature at the point used to fix the flame location
193  double fixedTemperature();
194 
195  //! Return location of the point where temperature is fixed
196  double fixedTemperatureLocation();
197 
198  /**
199  * Set grid refinement criteria. If dom >= 0, then the settings
200  * apply only to the specified domain. If dom < 0, the settings
201  * are applied to each domain. @see Refiner::setCriteria.
202  */
203  void setRefineCriteria(int dom = -1, double ratio = 10.0,
204  double slope = 0.8, double curve = 0.8,
205  double prune = -0.1);
206 
207  /**
208  * Get the grid refinement criteria. dom must be greater than
209  * or equal to zero (that is, the domain must be specified).
210  * @see Refiner::getCriteria
211  */
212  vector<double> getRefineCriteria(int dom);
213 
214  /**
215  * Set the maximum number of grid points in the domain. If dom >= 0,
216  * then the settings apply only to the specified domain. If dom < 0,
217  * the settings are applied to each domain. @see Refiner::setMaxPoints.
218  */
219  void setMaxGridPoints(int dom, int npoints);
220 
221  /**
222  * Get the maximum number of grid points in this domain. @see Refiner::maxPoints
223  *
224  * @param dom domain number, beginning with 0 for the leftmost domain.
225  */
226  size_t maxGridPoints(size_t dom);
227 
228  //! Set the minimum grid spacing in the specified domain(s).
229  /*!
230  * @param dom Domain index. If dom == -1, the specified spacing is applied
231  * to all domains.
232  * @param gridmin The minimum allowable grid spacing [m]
233  */
234  void setGridMin(int dom, double gridmin);
235 
236  //! Set the current solution vector to the last successful time-stepping
237  //! solution. This can be used to examine the solver progress after a failed
238  //! integration.
240 
241  //! Set the current solution vector and grid to the last successful steady-
242  //! state solution. This can be used to examine the solver progress after a
243  //! failure during grid refinement.
244  void restoreSteadySolution();
245 
246  void getInitialSoln();
247 
248  double jacobian(int i, int j);
249 
250  void evalSSJacobian();
251 
252  //! Solve the equation @f$ J^T \lambda = b @f$.
253  /**
254  * Here, @f$ J = \partial f/\partial x @f$ is the Jacobian matrix of the
255  * system of equations @f$ f(x,p)=0 @f$. This can be used to efficiently
256  * solve for the sensitivities of a scalar objective function @f$ g(x,p) @f$
257  * to a vector of parameters @f$ p @f$ by solving:
258  * @f[ J^T \lambda = \left( \frac{\partial g}{\partial x} \right)^T @f]
259  * for @f$ \lambda @f$ and then computing:
260  * @f[
261  * \left.\frac{dg}{dp}\right|_{f=0} = \frac{\partial g}{\partial p}
262  * - \lambda^T \frac{\partial f}{\partial p}
263  * @f]
264  */
265  void solveAdjoint(const double* b, double* lambda);
266 
267  void resize() override;
268 
269  //! Set a function that will be called after each successful steady-state
270  //! solve, before regridding. Intended to be used for observing solver
271  //! progress for debugging purposes.
272  void setSteadyCallback(Func1* callback) {
273  m_steady_callback = callback;
274  }
275 
276 protected:
277  //! the solution vector after the last successful timestepping
278  vector<double> m_xlast_ts;
279 
280  //! the solution vector after the last successful steady-state solve (stored
281  //! before grid refinement)
282  vector<double> m_xlast_ss;
283 
284  //! the grids for each domain after the last successful steady-state solve
285  //! (stored before grid refinement)
286  vector<vector<double>> m_grid_last_ss;
287 
288  //! a work array used to hold the residual or the new solution
289  vector<double> m_xnew;
290 
291  //! timestep
292  double m_tstep;
293 
294  //! array of number of steps to take before re-attempting the steady-state
295  //! solution
296  vector<int> m_steps;
297 
298  //! User-supplied function called after a successful steady-state solve.
300 
301 private:
302  //! Calls method _finalize in each domain.
303  void finalize();
304 
305  //! Wrapper around the Newton solver
306  /*!
307  * @return 0 if successful, -1 on failure
308  */
309  int newtonSolve(int loglevel);
310 };
311 
312 }
313 #endif
A map of string keys to values whose type can vary at runtime.
Definition: AnyMap.h:427
Base class for 'functor' classes that evaluate a function of one variable.
Definition: Func1.h:75
Container class for multiple-domain 1D problems.
Definition: OneDim.h:27
void eval(size_t j, double *x, double *r, double rdt=-1.0, int count=1)
Evaluate the multi-domain residual function.
Definition: OneDim.cpp:246
double rdt() const
Reciprocal of the time step.
Definition: OneDim.h:153
std::tuple< string, size_t, 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:50
shared_ptr< vector< double > > m_state
Solution vector.
Definition: OneDim.h:332
One-dimensional simulations.
Definition: Sim1D.h:22
void restoreTimeSteppingSolution()
Set the current solution vector to the last successful time-stepping solution.
Definition: Sim1D.cpp:352
void resize() override
Call after one or more grids has changed size, for example after being refined.
Definition: Sim1D.cpp:818
void saveResidual(const string &fname, const string &name, const string &desc, bool overwrite=false, int compression=0)
Save the residual of the current solution to a container file.
Definition: Sim1D.cpp:147
vector< double > m_xnew
a work array used to hold the residual or the new solution
Definition: Sim1D.h:289
void setProfile(size_t dom, size_t comp, const vector< double > &pos, const vector< double > &values)
Specify a profile for one component of one domain.
Definition: Sim1D.cpp:79
double fixedTemperatureLocation()
Return location of the point where temperature is fixed.
Definition: Sim1D.cpp:714
vector< vector< double > > m_grid_last_ss
the grids for each domain after the last successful steady-state solve (stored before grid refinement...
Definition: Sim1D.h:286
void finalize()
Calls method _finalize in each domain.
Definition: Sim1D.cpp:381
void setValue(size_t dom, size_t comp, size_t localPoint, double value)
Set a single value in the solution vector.
Definition: Sim1D.cpp:55
void setSteadyCallback(Func1 *callback)
Set a function that will be called after each successful steady-state solve, before regridding.
Definition: Sim1D.h:272
int refine(int loglevel=0)
Refine the grid in all domains.
Definition: Sim1D.cpp:522
void show()
Show logging information on current solution for all domains.
Definition: Sim1D.cpp:341
double fixedTemperature()
Return temperature at the point used to fix the flame location.
Definition: Sim1D.cpp:701
vector< double > m_xlast_ss
the solution vector after the last successful steady-state solve (stored before grid refinement)
Definition: Sim1D.h:282
void setMaxGridPoints(int dom, int npoints)
Set the maximum number of grid points in the domain.
Definition: Sim1D.cpp:765
int setFixedTemperature(double t)
Add node for fixed temperature point of freely propagating flame.
Definition: Sim1D.cpp:611
void setInitialGuess(const string &component, vector< double > &locs, vector< double > &vals)
Set initial guess for one component for all domains.
Definition: Sim1D.cpp:41
int newtonSolve(int loglevel)
Wrapper around the Newton solver.
Definition: Sim1D.cpp:397
vector< int > m_steps
array of number of steps to take before re-attempting the steady-state solution
Definition: Sim1D.h:296
double m_tstep
timestep
Definition: Sim1D.h:292
void solveAdjoint(const double *b, double *lambda)
Solve the equation .
Definition: Sim1D.cpp:794
AnyMap restore(const string &fname, const string &name)
Retrieve data and settings from a previously saved simulation.
Definition: Sim1D.cpp:247
Func1 * m_steady_callback
User-supplied function called after a successful steady-state solve.
Definition: Sim1D.h:299
void restoreSteadySolution()
Set the current solution vector and grid to the last successful steady- state solution.
Definition: Sim1D.cpp:361
size_t maxGridPoints(size_t dom)
Get the maximum number of grid points in this domain.
Definition: Sim1D.cpp:778
void setFlatProfile(size_t dom, size_t comp, double v)
Set component 'comp' of domain 'dom' to value 'v' at all points.
Definition: Sim1D.cpp:324
double value(size_t dom, size_t comp, size_t localPoint) const
Get one entry in the solution vector.
Definition: Sim1D.cpp:63
vector< double > getRefineCriteria(int dom)
Get the grid refinement criteria.
Definition: Sim1D.cpp:741
Sim1D()
Default constructor.
Definition: Sim1D.h:29
void setGridMin(int dom, double gridmin)
Set the minimum grid spacing in the specified domain(s).
Definition: Sim1D.cpp:752
void setRefineCriteria(int dom=-1, double ratio=10.0, double slope=0.8, double curve=0.8, double prune=-0.1)
Set grid refinement criteria.
Definition: Sim1D.cpp:727
vector< double > m_xlast_ts
the solution vector after the last successful timestepping
Definition: Sim1D.h:278
void save(const string &fname, const string &name, const string &desc, bool overwrite=false, int compression=0, const string &basis="")
Save current simulation data to a container file or CSV format.
Definition: Sim1D.cpp:98
MultiJac & jacobian()
Return a reference to the Jacobian evaluator of an OneDim object.
Definition: OneDim.cpp:87
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:564
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:180