Cantera  3.2.0a5
Loading...
Searching...
No Matches
Domain1D.h
Go to the documentation of this file.
1 //! @file Domain1D.h
2
3// This file is part of Cantera. See License.txt in the top-level directory or
4// at https://cantera.org/license.txt for license and copyright information.
5
6#ifndef CT_DOMAIN1D_H
7#define CT_DOMAIN1D_H
8
10#include "cantera/base/global.h"
11
12namespace Cantera
13{
14
15class MultiJac;
16class OneDim;
17class Refiner;
18class AnyMap;
19class Kinetics;
20class Transport;
21class Solution;
22class SolutionArray;
23
24/**
25 * Base class for one-dimensional domains.
26 * @ingroup flowGroup
27 */
29{
30public:
31 /**
32 * Constructor.
33 * @param nv Number of variables at each grid point.
34 * @param points Number of grid points.
35 * @param time (unused)
36 * @deprecated After %Cantera 3.2, this constructor will become protected
37 */
38 explicit Domain1D(size_t nv=1, size_t points=1, double time=0.0);
39
40 virtual ~Domain1D();
41 Domain1D(const Domain1D&) = delete;
42 Domain1D& operator=(const Domain1D&) = delete;
43
44 //! Domain type flag.
45 //! @since Starting in %Cantera 3.1, the return type is a `string`.
46 virtual string domainType() const { return "domain"; }
47
48 //! String indicating the domain implemented.
49 //! @since New in %Cantera 3.0.
50 //! @deprecated Transitional method. Use domainType() instead.
51 string type() const { return domainType(); }
52
53 //! The left-to-right location of this domain.
54 size_t domainIndex() {
55 return m_index;
56 }
57
58 //! True if the domain is a connector domain.
59 virtual bool isConnector() {
60 return false;
61 }
62
63 //! Set the solution manager.
64 //! @since New in %Cantera 3.0.
65 //! @deprecated To be removed after %Cantera 3.2. Solution objects must be
66 //! provided to constructors instead.
67 void setSolution(shared_ptr<Solution> sol);
68
69 //! Set the kinetics manager.
70 //! @since New in %Cantera 3.0.
71 //! @deprecated To be removed after %Cantera 3.2. A change of domain contents after
72 //! instantiation will be disabled.
73 virtual void setKinetics(shared_ptr<Kinetics> kin) {
74 warn_deprecated("Domain1D::setKinetics",
75 "After Cantera 3.2, a change of domain contents after instantiation "
76 "will be disabled.");
77 throw NotImplementedError("Domain1D::setKinetics");
78 }
79
80 //! Set transport model to existing instance
81 //! @since New in %Cantera 3.0.
82 //! @todo Convert warning message to new exception and remove virtual.
83 virtual void setTransport(shared_ptr<Transport> trans) {
84 warn_deprecated("Domain1D::setTransport",
85 "After Cantera 3.2, a change of domain contents after instantiation "
86 "will be disabled.");
87 throw NotImplementedError("Domain1D::setTransport");
88 }
89
90 //! Set transport model by name.
91 //! @param model String specifying model name.
92 //! @since New in %Cantera 3.2.
93 virtual void setTransportModel(const string& model) {
94 throw NotImplementedError("Domain1D::setTransportModel");
95 }
96
97protected:
98 //! Update transport model to existing instance
99 //! @since New in %Cantera 3.2.
100 virtual void _setKinetics(shared_ptr<Kinetics> kin) {
101 throw NotImplementedError("Domain1D::_setKinetics");
102 }
103
104 //! Update transport model to existing instance
105 //! @since New in %Cantera 3.2.
106 virtual void _setTransport(shared_ptr<Transport> trans) {
107 throw NotImplementedError("Domain1D::_setTransport");
108 }
109
110public:
111 //! The container holding this domain.
112 const OneDim& container() const {
113 return *m_container;
114 }
115
116 //! Specify the container object for this domain, and the position of this
117 //! domain in the list.
118 void setContainer(OneDim* c, size_t index) {
119 m_container = c;
120 m_index = index;
121 }
122
123 //! Set the Jacobian bandwidth. See the discussion of method bandwidth().
124 void setBandwidth(int bw = -1) {
125 m_bw = bw;
126 }
127
128 //! Set the Jacobian bandwidth for this domain.
129 /**
130 * When class OneDim computes the bandwidth of the overall multi-domain
131 * problem (in OneDim::resize()), it calls this method for the bandwidth
132 * of each domain. If setBandwidth has not been called, then a negative
133 * bandwidth is returned, in which case OneDim assumes that this domain is
134 * dense -- that is, at each point, all components depend on the value of
135 * all other components at that point. In this case, the bandwidth is bw =
136 * 2*nComponents() - 1. However, if this domain contains some components
137 * that are uncoupled from other components at the same point, then this
138 * default bandwidth may greatly overestimate the true bandwidth, with a
139 * substantial penalty in performance. For such domains, use method
140 * setBandwidth to specify the bandwidth before passing this domain to the
141 * Sim1D or OneDim constructor.
142 */
143 size_t bandwidth() {
144 return m_bw;
145 }
146
147 /**
148 * Initialize. This method is called by OneDim::init() for each domain once
149 * at the beginning of a simulation. Base class method does nothing, but may
150 * be overloaded.
151 */
152 virtual void init() { }
153
154 /**
155 * When called, this function should reset "bad" values in the state vector
156 * such as negative species concentrations. This function may be called
157 * after a failed solution attempt.
158 */
159 virtual void resetBadValues(double* xg) {}
160
161 /**
162 * Resize the domain to have nv components and np grid points. This method
163 * is virtual so that subclasses can perform other actions required to
164 * resize the domain.
165 */
166 virtual void resize(size_t nv, size_t np);
167
168 //! Return a reference to the grid refiner.
170 return *m_refiner;
171 }
172
173 //! Number of components at each grid point.
174 size_t nComponents() const {
175 return m_nv;
176 }
177
178 //! Check that the specified component index is in range.
179 //! Throws an exception if n is greater than nComponents()-1
180 //! @deprecated To be removed after %Cantera 3.2. Only used by legacy CLib.
181 void checkComponentIndex(size_t n) const {
182 warn_deprecated("Domain1D::checkComponentIndex",
183 "To be removed after Cantera 3.2. Only used by legacy CLib.");
184 if (n >= m_nv) {
185 throw IndexError("Domain1D::checkComponentIndex", "points", n, m_nv);
186 }
187 }
188
189 //! Check that an array size is at least nComponents().
190 //! Throws an exception if nn is less than nComponents(). Used before calls
191 //! which take an array pointer.
192 //! @deprecated To be removed after %Cantera 3.2. Unused.
193 void checkComponentArraySize(size_t nn) const {
194 warn_deprecated("Domain1D::checkComponentArraySize",
195 "To be removed after Cantera 3.2. Unused.");
196 if (m_nv > nn) {
197 throw ArraySizeError("Domain1D::checkComponentArraySize", nn, m_nv);
198 }
199 }
200
201 //! Number of grid points in this domain.
202 size_t nPoints() const {
203 return m_points;
204 }
205
206 //! Check that the specified point index is in range.
207 //! Throws an exception if n is greater than nPoints()-1
208 //! @deprecated To be removed after %Cantera 3.2. Only used by legacy CLib.
209 void checkPointIndex(size_t n) const {
210 warn_deprecated("Domain1D::checkPointIndex",
211 "To be removed after Cantera 3.2. Only used by legacy CLib.");
212 if (n >= m_points) {
213 throw IndexError("Domain1D::checkPointIndex", "points", n, m_points);
214 }
215 }
216
217 //! Check that an array size is at least nPoints().
218 //! Throws an exception if nn is less than nPoints(). Used before calls
219 //! which take an array pointer.
220 //! @deprecated To be removed after %Cantera 3.2. Unused.
221 void checkPointArraySize(size_t nn) const {
222 warn_deprecated("Domain1D::checkPointArraySize",
223 "To be removed after Cantera 3.2. Unused.");
224 if (m_points > nn) {
225 throw ArraySizeError("Domain1D::checkPointArraySize", nn, m_points);
226 }
227 }
228
229 //! Name of component `n`. May be overloaded.
230 virtual string componentName(size_t n) const;
231
232 //! Set the name of the component `n` to `name`.
233 void setComponentName(size_t n, const string& name) {
234 m_name[n] = name;
235 }
236
237 /**
238 * Index of component with name `name`.
239 * @param name name of component
240 * @param checkAlias if `true` (default), check alias mapping
241 */
242 virtual size_t componentIndex(const string& name, bool checkAlias=true) const;
243
244 /**
245 * Check whether the Domain contains a component.
246 * @param name name of component
247 * @param checkAlias if `true` (default), check alias mapping
248 *
249 * @since New in %Cantera 3.2.
250 */
251 virtual bool hasComponent(const string& name, bool checkAlias=true) const;
252
253 /**
254 * Update state at given location to state of associated Solution object.
255 */
256 virtual void updateState(size_t loc) {
257 throw NotImplementedError("Domain1D::updateState",
258 "Not implemented for domain type '{}'.", domainType());
259 }
260
261 /**
262 * Set the upper and lower bounds for a solution component, n.
263 *
264 * @param n solution component index
265 * @param lower lower bound on component n
266 * @param upper upper bound on component n
267 */
268 void setBounds(size_t n, double lower, double upper) {
269 m_min[n] = lower;
270 m_max[n] = upper;
271 }
272
273 //! Set tolerances for time-stepping mode
274 /*!
275 * @param rtol Relative tolerance
276 * @param atol Absolute tolerance
277 * @param n component index these tolerances apply to. If set to -1 (the
278 * default), these tolerances will be applied to all solution
279 * components.
280 */
281 void setTransientTolerances(double rtol, double atol, size_t n=npos);
282
283 //! Set tolerances for steady-state mode
284 /*!
285 * @param rtol Relative tolerance
286 * @param atol Absolute tolerance
287 * @param n component index these tolerances apply to. If set to -1 (the
288 * default), these tolerances will be applied to all solution
289 * components.
290 */
291 void setSteadyTolerances(double rtol, double atol, size_t n=npos);
292
293 //! Relative tolerance of the nth component.
294 double rtol(size_t n) {
295 return (m_rdt == 0.0 ? m_rtol_ss[n] : m_rtol_ts[n]);
296 }
297
298 //! Absolute tolerance of the nth component.
299 double atol(size_t n) {
300 return (m_rdt == 0.0 ? m_atol_ss[n] : m_atol_ts[n]);
301 }
302
303 //! Steady relative tolerance of the nth component
304 double steady_rtol(size_t n) {
305 return m_rtol_ss[n];
306 }
307
308 //! Steady absolute tolerance of the nth component
309 double steady_atol(size_t n) {
310 return m_atol_ss[n];
311 }
312
313 //! Transient relative tolerance of the nth component
314 double transient_rtol(size_t n) {
315 return m_rtol_ts[n];
316 }
317
318 //! Transient absolute tolerance of the nth component
319 double transient_atol(size_t n) {
320 return m_atol_ts[n];
321 }
322
323 //! Upper bound on the nth component.
324 double upperBound(size_t n) const {
325 return m_max[n];
326 }
327
328 //! Lower bound on the nth component
329 double lowerBound(size_t n) const {
330 return m_min[n];
331 }
332
333 /**
334 * Set grid refinement criteria. @see Refiner::setCriteria.
335 * @since New in %Cantera 3.2
336 */
337 void setRefineCriteria(double ratio = 10.0,
338 double slope = 0.8, double curve = 0.8,
339 double prune = -0.1);
340
341 /**
342 * Get the grid refinement criteria. @see Refiner::getCriteria
343 * @since New in %Cantera 3.2
344 */
345 vector<double> getRefineCriteria();
346
347 /**
348 * Performs the setup required before starting a time-stepping solution.
349 * Stores the solution provided in `x0` to the internal storage, and sets
350 * the reciprocal of the time step to `1/dt`.
351 *
352 * @param[in] dt Time step
353 * @param[in] x0 Array to store the solution at the last time step
354 */
355 void initTimeInteg(double dt, const double* x0) {
356 std::copy(x0 + loc(), x0 + loc() + size(), m_slast.begin());
357 m_rdt = 1.0/dt;
358 }
359
360 /**
361 * Set the internally-stored reciprocal of the time step to 0.0, which is
362 * used to indicate that the problem is in steady-state mode.
363 */
365 m_rdt = 0.0;
366 }
367
368 //! True if in steady-state mode
369 bool steady() {
370 return (m_rdt == 0.0);
371 }
372
373 //! True if not in steady-state mode
374 bool transient() {
375 return (m_rdt != 0.0);
376 }
377
378 /**
379 * Set this if something has changed in the governing
380 * equations (for example, the value of a constant has been changed,
381 * so that the last-computed Jacobian is no longer valid.
382 */
383 void needJacUpdate();
384
385 //! Evaluate the residual function at point j. If j == npos,
386 //! evaluate the residual function at all points.
387 /*!
388 * This function must be implemented in classes derived from Domain1D.
389 *
390 * @param[in] j Grid point at which to update the residual
391 * @param[in] x State vector
392 * @param[out] r residual vector
393 * @param[out] mask Boolean mask indicating whether each solution
394 * component has a time derivative (1) or not (0).
395 * @param[in] rdt Reciprocal of the timestep (`rdt=0` implies steady-state.)
396 */
397 virtual void eval(size_t j, double* x, double* r, integer* mask, double rdt=0.0) {
398 throw NotImplementedError("Domain1D::eval");
399 }
400
401 /**
402 * Returns the index of the solution vector, which corresponds to component
403 * n at grid point j.
404 *
405 * @param n component index
406 * @param j grid point index
407 */
408 size_t index(size_t n, size_t j) const {
409 return m_nv*j + n;
410 }
411
412 /**
413 * Returns the value of solution component n at grid point j of the solution
414 * vector x.
415 *
416 * @param x solution vector
417 * @param n component index
418 * @param j grid point index
419 *
420 * @deprecated To be removed after %Cantera 3.2. Replaceable with version accessing
421 * component by name.
422 */
423 double value(const double* x, size_t n, size_t j) const {
424 warn_deprecated("Domain1D::value",
425 "Replaceable with version accessing component by name");
426 return x[index(n,j)]; // caution: this assumes m_iloc = 0
427 }
428
429 /**
430 * Set a single component value at a boundary.
431 * @param component Name of the component.
432 *
433 * @since New in %Cantera 3.2.
434 */
435 virtual double value(const string& component) const {
436 throw NotImplementedError("Domain1D::value",
437 "Not implemented for domain type '{}'.", domainType());
438 }
439
440 /**
441 * Set a single component value in a flow domain or at a boundary.
442 * @param component Name of the component.
443 * @param value Value of the component.
444 *
445 * @since New in %Cantera 3.2.
446 */
447 virtual void setValue(const string& component, double value) {
448 throw NotImplementedError("Domain1D::setValue",
449 "Not implemented for domain type '{}'.", domainType());
450 }
451
452 /**
453 * Retrieve component values.
454 * @param component Name of the component.
455 * @returns Vector of length nPoints() containing values at grid points.
456 *
457 * @since New in %Cantera 3.2.
458 */
459 vector<double> values(const string& component) const {
460 vector<double> data(nPoints());
461 getValues(component, data);
462 return data;
463 }
464
465 /**
466 * Retrieve component values.
467 * @param component Name of the component.
468 * @param[out] values Vector of length nPoints() containing values at grid points.
469 *
470 * @since New in %Cantera 3.2.
471 */
472 virtual void getValues(const string& component, vector<double>& values) const {
473 throw NotImplementedError("Domain1D::getValues",
474 "Not implemented for domain type '{}'.", domainType());
475 }
476
477 /**
478 * Specify component values.
479 * @param component Name of the component.
480 * @param[in] values Vector of length nPoints() containing values at grid points.
481 *
482 * @since New in %Cantera 3.2.
483 */
484 virtual void setValues(const string& component, const vector<double>& values) {
485 throw NotImplementedError("Domain1D::setValues",
486 "Not implemented for domain type '{}'.", domainType());
487 }
488
489 /**
490 * Retrieve internal work array values for a component.
491 * After calling Sim1D::eval(), this array contains the values of the residual
492 * function.
493 * @param component Name of the component.
494 * @returns Vector of length nPoints() containing residuals at grid points.
495 *
496 * @since New in %Cantera 3.2.
497 */
498 vector<double> residuals(const string& component) const {
499 vector<double> data(nPoints());
500 getResiduals(component, data);
501 return data;
502 }
503
504 /**
505 * Retrieve internal work array values for a component.
506 * After calling Sim1D::eval(), this array contains the values of the residual
507 * function.
508 * @param component Name of the component.
509 * @param[out] values Vector of length nPoints() containing residuals at grid
510 * points.
511 *
512 * @since New in %Cantera 3.2.
513 */
514 virtual void getResiduals(const string& component, vector<double>& values) const {
515 throw NotImplementedError("Domain1D::getResiduals",
516 "Not applicable or not implemented for domain type '{}'.", domainType());
517 }
518
519 /**
520 * Specify a profile for a component.
521 * @param component Name of the component.
522 * @param[in] pos A vector of relative positions, beginning with 0.0 at the
523 * left of the domain, and ending with 1.0 at the right of the domain.
524 * @param[in] values A vector of values corresponding to the relative position
525 * locations.
526 *
527 * Note that the vector pos and values can have lengths different than the
528 * number of grid points, but their lengths must be equal. The values at
529 * the grid points will be linearly interpolated based on the (pos,
530 * values) specification.
531 *
532 * @since New in %Cantera 3.2.
533 */
534 virtual void setProfile(const string& component,
535 const vector<double>& pos, const vector<double>& values) {
536 throw NotImplementedError("Domain1D::setProfile",
537 "Not implemented for domain type '{}'.", domainType());
538 }
539
540 /**
541 * Specify a flat profile for a component.
542 * @param component Name of the component.
543 * @param value Constant value.
544 *
545 * @since New in %Cantera 3.2.
546 */
547 virtual void setFlatProfile(const string& component, double value) {
548 throw NotImplementedError("Domain1D::setFlatProfile",
549 "Not implemented for domain type '{}'.", domainType());
550 }
551
552 //! Save the state of this domain as a SolutionArray.
553 /*!
554 * @param soln local solution vector for this domain
555 * @todo Despite the method's name, data are copied; the intent is to access data
556 * directly in future revisions, where a non-const version will be implemented.
557 *
558 * @since New in %Cantera 3.0.
559 * @deprecated To be removed after %Cantera 3.2. Replaceable by toArray().
560 */
561 shared_ptr<SolutionArray> asArray(const double* soln) {
562 warn_deprecated("Domain1D::asArray",
563 "To be removed after Cantera 3.2. Replaceable by 'toArray'.");
564 return toArray();
565 }
566
567 //! Save the state of this domain to a SolutionArray.
568 /*!
569 * This method serves as an external interface for high-level API's; it does not
570 * provide direct access to memory.
571 * @param normalize If true, normalize concentrations (default=false)
572 *
573 * @since New in %Cantera 3.0.
574 */
575 virtual shared_ptr<SolutionArray> toArray(bool normalize=false) {
576 throw NotImplementedError("Domain1D::toArray", "Needs to be overloaded.");
577 }
578
579 //! Restore the solution for this domain from a SolutionArray
580 /*!
581 * @param[in] arr SolutionArray defining the state of this domain
582 * @param[out] soln Value of the solution vector, local to this domain
583 *
584 * @since New in %Cantera 3.0.
585 * @deprecated To be removed after %Cantera 3.2.
586 * Replaceable by version that does not require solution vector.
587 */
588 void fromArray(SolutionArray& arr, double* soln) {
589 warn_deprecated("Domain1D::fromArray",
590 "To be removed after Cantera 3.2. Replaceable by version that does not "
591 "require solution vector.");
592 // create a shared_ptr that skips deletion of the held pointer
593 shared_ptr<SolutionArray> arr_ptr(&arr, [](SolutionArray*){});
594 fromArray(arr_ptr);
595 }
596
597 //! Restore the solution for this domain from a SolutionArray.
598 /*!
599 * This method serves as an external interface for high-level API's.
600 * @param arr SolutionArray defining the state of this domain
601 * @since New in %Cantera 3.0.
602 */
603 virtual void fromArray(const shared_ptr<SolutionArray>& arr) {
604 throw NotImplementedError("Domain1D::fromArray", "Needs to be overloaded.");
605 }
606
607 /**
608 * Return a concise summary of a Domain.
609 * @see SolutionArray.info()
610 * @param keys List of components to be displayed; if empty, all components are
611 * considered.
612 * @param rows Maximum number of rendered rows.
613 * @param width Maximum width of rendered output.
614 * @since New in %Cantera 3.2
615 */
616 string info(const vector<string>& keys, int rows=10, int width=80);
617
618 /**
619 * Return a concise summary of a Domain.
620 * Skips keys input while `vector<string>` is not implemented in sourcegen.
621 * @see SolutionArray.info()
622 * @param rows Maximum number of rendered rows.
623 * @param width Maximum width of rendered output.
624 * @since New in %Cantera 3.2
625 */
626 string _info(int rows=10, int width=80) {
627 return info({}, rows, width);
628 }
629
630 //! Return thermo/kinetics/transport manager used in the domain
631 //! @since New in %Cantera 3.0.
632 //! @deprecated To be removed after %Cantera 3.2. Renamed to phase().
633 shared_ptr<Solution> solution() const {
634 warn_deprecated("Domain1D::solution",
635 "To be removed after Cantera 3.2. Renamed to phase().");
636 return m_solution;
637 }
638
639 //! Return thermo/kinetics/transport manager used in the domain
640 //! @since New in %Cantera 3.2.
641 shared_ptr<Solution> phase() const {
642 return m_solution;
643 }
644
645 //! Return the size of the solution vector (the product of #m_nv and #m_points).
646 size_t size() const {
647 return m_nv*m_points;
648 }
649
650 /**
651 * Find the index of the first grid point in this domain, and
652 * the start of its variables in the global solution vector.
653 */
654 void locate();
655
656 //! Location of the start of the local solution vector in the global solution vector
657 virtual size_t loc(size_t j = 0) const {
658 return m_iloc;
659 }
660
661 //! The index of the first (that is, left-most) grid point belonging to this domain.
662 size_t firstPoint() const {
663 return m_jstart;
664 }
665
666 //! The index of the last (that is, right-most) grid point belonging to this domain.
667 size_t lastPoint() const {
668 return m_jstart + m_points - 1;
669 }
670
671 /**
672 * Set the left neighbor to domain 'left.' Method 'locate' is called to
673 * update the global positions of this domain and all those to its right.
674 */
676 m_left = left;
677 if (!m_solution && left && left->phase()) {
679 }
680 locate();
681 }
682
683 //! Set the right neighbor to domain 'right.'
685 m_right = right;
686 if (!m_solution && right && right->phase()) {
688 }
689 }
690
691 //! Append domain 'right' to this one, and update all links.
694 right->linkLeft(this);
695 }
696
697 //! Return a pointer to the left neighbor.
698 Domain1D* left() const {
699 return m_left;
700 }
701
702 //! Return a pointer to the right neighbor.
703 Domain1D* right() const {
704 return m_right;
705 }
706
707 //! Value of component n at point j in the previous solution.
708 double prevSoln(size_t n, size_t j) const {
709 return m_slast[m_nv*j + n];
710 }
711
712 //! Specify an identifying tag for this domain.
713 void setID(const string& s) {
714 m_id = s;
715 }
716
717 //! Returns the identifying tag for this domain.
718 string id() const {
719 if (m_id != "") {
720 return m_id;
721 } else {
722 return fmt::format("domain {}", m_index);
723 }
724 }
725
726 //! Print the solution.
727 //! @param x Pointer to the local portion of the system state vector
728 virtual void show(const double* x);
729
730 //! Get the coordinate [m] of the point with local index `jlocal`
731 double z(size_t jlocal) const {
732 return m_z[jlocal];
733 }
734
735 //! Get the coordinate [m] of the first (leftmost) grid point in this domain
736 double zmin() const {
737 return m_z[0];
738 }
739
740 //! Get the coordinate [m] of the last (rightmost) grid point in this domain
741 double zmax() const {
742 return m_z[m_points - 1];
743 }
744
745 //! Set initial values for a component at each grid point
746 //! @param name Name of the component
747 //! @param values Array of length nPoints() containing the initial values
748 //! @param soln Pointer to the local portion of the system state vector
749 //! @deprecated To be removed after %Cantera 3.2.
750 //! Replaceable by version using vectors arguments.
751 void setProfile(const string& name, double* values, double* soln);
752
753 //! Access the array of grid coordinates [m]
754 vector<double>& grid() {
755 return m_z;
756 }
757
758 //! Access the array of grid coordinates [m]
759 const vector<double>& grid() const {
760 return m_z;
761 }
762
763 //! Set up initial grid.
764 //! @since New in %Cantera 3.2.
765 void setupGrid(const vector<double>& grid);
766
767 //! called to set up initial grid, and after grid refinement
768 virtual void setupGrid(size_t n, const double* z);
769
770 //! Set up uniform grid.
771 //! @param points Number of grid points
772 //! @param length Length of domain
773 //! @param start Start position of domain (default=0.)
774 //! @since New in %Cantera 3.2.
775 void setupUniformGrid(size_t points, double length, double start=0.);
776
777 /**
778 * Writes some or all initial solution values into the global solution
779 * array, beginning at the location pointed to by x. This method is called
780 * by the Sim1D constructor, and allows default values or ones that have
781 * been set locally prior to installing this domain into the container to be
782 * written to the global solution vector.
783 */
784 virtual void _getInitialSoln(double* x);
785
786 //! Initial value of solution component @e n at grid point @e j.
787 virtual double initialValue(size_t n, size_t j);
788
789 /**
790 * In some cases, a domain may need to set parameters that depend on the
791 * initial solution estimate. In such cases, the parameters may be set in
792 * method _finalize. This method is called just before the Newton solver is
793 * called, and the x array is guaranteed to be the local solution vector for
794 * this domain that will be used as the initial guess. If no such parameters
795 * need to be set, then method _finalize does not need to be overloaded.
796 */
797 virtual void _finalize(const double* x) {}
798
799 /**
800 * In some cases, for computational efficiency some properties (such as
801 * transport coefficients) may not be updated during Jacobian evaluations.
802 * Set this to `true` to force these properties to be updated even while
803 * calculating Jacobian elements.
804 */
805 void forceFullUpdate(bool update) {
806 m_force_full_update = update;
807 }
808
809 //! Set shared data pointer
810 void setData(shared_ptr<vector<double>>& data) {
811 m_state = data;
812 }
813
814protected:
815 //! Retrieve meta data
816 virtual AnyMap getMeta() const;
817
818 //! Retrieve meta data
819 virtual void setMeta(const AnyMap& meta);
820
821 shared_ptr<vector<double>> m_state; //!< data pointer shared from OneDim
822
823 double m_rdt = 0.0; //!< Reciprocal of the time step
824 size_t m_nv = 0; //!< Number of solution components
825 size_t m_points; //!< Number of grid points
826 vector<double> m_slast; //!< Solution vector at the last time step
827 vector<double> m_max; //!< Upper bounds on solution components
828 vector<double> m_min; //!< Lower bounds on solution components
829 vector<double> m_rtol_ss; //!< Relative tolerances for steady mode
830 vector<double> m_rtol_ts; //!< Relative tolerances for transient mode
831 vector<double> m_atol_ss; //!< Absolute tolerances for steady mode
832 vector<double> m_atol_ts; //!< Absolute tolerances for transient mode
833 vector<double> m_z; //!< 1D spatial grid coordinates
834
835 //! Parent OneDim simulation containing this and adjacent domains
836 OneDim* m_container = nullptr;
837
838 size_t m_index; //!< Left-to-right location of this domain
839
840 //! Starting location within the solution vector for unknowns that
841 //! correspond to this domain
842 /*!
843 * Remember there may be multiple domains associated with this problem
844 */
845 size_t m_iloc = 0;
846
847 //! Index of the first point in this domain in the global point list.
848 //! @see firstPoint(), lastPoint()
849 size_t m_jstart = 0;
850
851 Domain1D* m_left = nullptr; //!< Pointer to the domain to the left
852 Domain1D* m_right = nullptr; //!< Pointer to the domain to the right
853
854 //! Identity tag for the domain
855 string m_id;
856 unique_ptr<Refiner> m_refiner; //!< Refiner object used for placing grid points
857 vector<string> m_name; //!< Names of solution components
858 int m_bw = -1; //!< See bandwidth()
859 bool m_force_full_update = false; //!< see forceFullUpdate()
860
861 //! Composite thermo/kinetics/transport handler
862 shared_ptr<Solution> m_solution;
863};
864}
865
866#endif
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:431
Array size error.
Base class for one-dimensional domains.
Definition Domain1D.h:29
void setTransientTolerances(double rtol, double atol, size_t n=npos)
Set tolerances for time-stepping mode.
Definition Domain1D.cpp:109
virtual void resetBadValues(double *xg)
When called, this function should reset "bad" values in the state vector such as negative species con...
Definition Domain1D.h:159
size_t lastPoint() const
The index of the last (that is, right-most) grid point belonging to this domain.
Definition Domain1D.h:667
size_t m_iloc
Starting location within the solution vector for unknowns that correspond to this domain.
Definition Domain1D.h:845
void checkPointArraySize(size_t nn) const
Check that an array size is at least nPoints().
Definition Domain1D.h:221
size_t domainIndex()
The left-to-right location of this domain.
Definition Domain1D.h:54
shared_ptr< Solution > m_solution
Composite thermo/kinetics/transport handler.
Definition Domain1D.h:862
Domain1D * m_left
Pointer to the domain to the left.
Definition Domain1D.h:851
OneDim * m_container
Parent OneDim simulation containing this and adjacent domains.
Definition Domain1D.h:836
bool transient()
True if not in steady-state mode.
Definition Domain1D.h:374
void setComponentName(size_t n, const string &name)
Set the name of the component n to name.
Definition Domain1D.h:233
size_t nComponents() const
Number of components at each grid point.
Definition Domain1D.h:174
shared_ptr< SolutionArray > asArray(const double *soln)
Save the state of this domain as a SolutionArray.
Definition Domain1D.h:561
size_t bandwidth()
Set the Jacobian bandwidth for this domain.
Definition Domain1D.h:143
double rtol(size_t n)
Relative tolerance of the nth component.
Definition Domain1D.h:294
shared_ptr< Solution > solution() const
Return thermo/kinetics/transport manager used in the domain.
Definition Domain1D.h:633
vector< double > m_atol_ss
Absolute tolerances for steady mode.
Definition Domain1D.h:831
void setupUniformGrid(size_t points, double length, double start=0.)
Set up uniform grid.
Definition Domain1D.cpp:238
vector< double > m_rtol_ts
Relative tolerances for transient mode.
Definition Domain1D.h:830
size_t size() const
Return the size of the solution vector (the product of m_nv and m_points).
Definition Domain1D.h:646
vector< double > m_atol_ts
Absolute tolerances for transient mode.
Definition Domain1D.h:832
virtual void updateState(size_t loc)
Update state at given location to state of associated Solution object.
Definition Domain1D.h:256
virtual bool isConnector()
True if the domain is a connector domain.
Definition Domain1D.h:59
string _info(int rows=10, int width=80)
Return a concise summary of a Domain.
Definition Domain1D.h:626
virtual void _setKinetics(shared_ptr< Kinetics > kin)
Update transport model to existing instance.
Definition Domain1D.h:100
virtual void setTransport(shared_ptr< Transport > trans)
Set transport model to existing instance.
Definition Domain1D.h:83
virtual void setMeta(const AnyMap &meta)
Retrieve meta data.
Definition Domain1D.cpp:171
virtual void _finalize(const double *x)
In some cases, a domain may need to set parameters that depend on the initial solution estimate.
Definition Domain1D.h:797
size_t m_index
Left-to-right location of this domain.
Definition Domain1D.h:838
void setupGrid(const vector< double > &grid)
Set up initial grid.
Definition Domain1D.cpp:223
Domain1D * left() const
Return a pointer to the left neighbor.
Definition Domain1D.h:698
vector< double > residuals(const string &component) const
Retrieve internal work array values for a component.
Definition Domain1D.h:498
string id() const
Returns the identifying tag for this domain.
Definition Domain1D.h:718
vector< double > & grid()
Access the array of grid coordinates [m].
Definition Domain1D.h:754
double zmin() const
Get the coordinate [m] of the first (leftmost) grid point in this domain.
Definition Domain1D.h:736
size_t m_jstart
Index of the first point in this domain in the global point list.
Definition Domain1D.h:849
shared_ptr< Solution > phase() const
Return thermo/kinetics/transport manager used in the domain.
Definition Domain1D.h:641
size_t m_nv
Number of solution components.
Definition Domain1D.h:824
void setContainer(OneDim *c, size_t index)
Specify the container object for this domain, and the position of this domain in the list.
Definition Domain1D.h:118
size_t nPoints() const
Number of grid points in this domain.
Definition Domain1D.h:202
vector< string > m_name
Names of solution components.
Definition Domain1D.h:857
double m_rdt
Reciprocal of the time step.
Definition Domain1D.h:823
virtual string domainType() const
Domain type flag.
Definition Domain1D.h:46
bool m_force_full_update
see forceFullUpdate()
Definition Domain1D.h:859
string info(const vector< string > &keys, int rows=10, int width=80)
Return a concise summary of a Domain.
Definition Domain1D.cpp:47
double lowerBound(size_t n) const
Lower bound on the nth component.
Definition Domain1D.h:329
virtual void setTransportModel(const string &model)
Set transport model by name.
Definition Domain1D.h:93
virtual size_t componentIndex(const string &name, bool checkAlias=true) const
Index of component with name name.
Definition Domain1D.cpp:87
void checkComponentIndex(size_t n) const
Check that the specified component index is in range.
Definition Domain1D.h:181
shared_ptr< vector< double > > m_state
data pointer shared from OneDim
Definition Domain1D.h:821
void linkLeft(Domain1D *left)
Set the left neighbor to domain 'left.
Definition Domain1D.h:675
vector< double > values(const string &component) const
Retrieve component values.
Definition Domain1D.h:459
virtual void resize(size_t nv, size_t np)
Resize the domain to have nv components and np grid points.
Definition Domain1D.cpp:52
double z(size_t jlocal) const
Get the coordinate [m] of the point with local index jlocal
Definition Domain1D.h:731
virtual void setValue(const string &component, double value)
Set a single component value in a flow domain or at a boundary.
Definition Domain1D.h:447
double upperBound(size_t n) const
Upper bound on the nth component.
Definition Domain1D.h:324
Refiner & refiner()
Return a reference to the grid refiner.
Definition Domain1D.h:169
Domain1D * right() const
Return a pointer to the right neighbor.
Definition Domain1D.h:703
virtual void show(const double *x)
Print the solution.
Definition Domain1D.cpp:248
vector< double > m_rtol_ss
Relative tolerances for steady mode.
Definition Domain1D.h:829
void fromArray(SolutionArray &arr, double *soln)
Restore the solution for this domain from a SolutionArray.
Definition Domain1D.h:588
vector< double > m_slast
Solution vector at the last time step.
Definition Domain1D.h:826
virtual double value(const string &component) const
Set a single component value at a boundary.
Definition Domain1D.h:435
virtual void fromArray(const shared_ptr< SolutionArray > &arr)
Restore the solution for this domain from a SolutionArray.
Definition Domain1D.h:603
Domain1D * m_right
Pointer to the domain to the right.
Definition Domain1D.h:852
double steady_atol(size_t n)
Steady absolute tolerance of the nth component.
Definition Domain1D.h:309
void setSteadyTolerances(double rtol, double atol, size_t n=npos)
Set tolerances for steady-state mode.
Definition Domain1D.cpp:122
virtual string componentName(size_t n) const
Name of component n. May be overloaded.
Definition Domain1D.cpp:75
double transient_atol(size_t n)
Transient absolute tolerance of the nth component.
Definition Domain1D.h:319
void setSolution(shared_ptr< Solution > sol)
Set the solution manager.
Definition Domain1D.cpp:31
virtual void getValues(const string &component, vector< double > &values) const
Retrieve component values.
Definition Domain1D.h:472
virtual void _setTransport(shared_ptr< Transport > trans)
Update transport model to existing instance.
Definition Domain1D.h:106
virtual void init()
Initialize.
Definition Domain1D.h:152
double atol(size_t n)
Absolute tolerance of the nth component.
Definition Domain1D.h:299
void setID(const string &s)
Specify an identifying tag for this domain.
Definition Domain1D.h:713
vector< double > m_z
1D spatial grid coordinates
Definition Domain1D.h:833
void forceFullUpdate(bool update)
In some cases, for computational efficiency some properties (such as transport coefficients) may not ...
Definition Domain1D.h:805
const vector< double > & grid() const
Access the array of grid coordinates [m].
Definition Domain1D.h:759
const OneDim & container() const
The container holding this domain.
Definition Domain1D.h:112
virtual void setKinetics(shared_ptr< Kinetics > kin)
Set the kinetics manager.
Definition Domain1D.h:73
virtual void setValues(const string &component, const vector< double > &values)
Specify component values.
Definition Domain1D.h:484
size_t m_points
Number of grid points.
Definition Domain1D.h:825
bool steady()
True if in steady-state mode.
Definition Domain1D.h:369
void setBandwidth(int bw=-1)
Set the Jacobian bandwidth. See the discussion of method bandwidth().
Definition Domain1D.h:124
double steady_rtol(size_t n)
Steady relative tolerance of the nth component.
Definition Domain1D.h:304
string m_id
Identity tag for the domain.
Definition Domain1D.h:855
string type() const
String indicating the domain implemented.
Definition Domain1D.h:51
vector< double > m_max
Upper bounds on solution components.
Definition Domain1D.h:827
unique_ptr< Refiner > m_refiner
Refiner object used for placing grid points.
Definition Domain1D.h:856
void initTimeInteg(double dt, const double *x0)
Performs the setup required before starting a time-stepping solution.
Definition Domain1D.h:355
void setBounds(size_t n, double lower, double upper)
Set the upper and lower bounds for a solution component, n.
Definition Domain1D.h:268
vector< double > m_min
Lower bounds on solution components.
Definition Domain1D.h:828
virtual void setProfile(const string &component, const vector< double > &pos, const vector< double > &values)
Specify a profile for a component.
Definition Domain1D.h:534
void setData(shared_ptr< vector< double > > &data)
Set shared data pointer.
Definition Domain1D.h:810
virtual void eval(size_t j, double *x, double *r, integer *mask, double rdt=0.0)
Evaluate the residual function at point j.
Definition Domain1D.h:397
void append(Domain1D *right)
Append domain 'right' to this one, and update all links.
Definition Domain1D.h:692
double value(const double *x, size_t n, size_t j) const
Returns the value of solution component n at grid point j of the solution vector x.
Definition Domain1D.h:423
void setSteadyMode()
Set the internally-stored reciprocal of the time step to 0.0, which is used to indicate that the prob...
Definition Domain1D.h:364
virtual double initialValue(size_t n, size_t j)
Initial value of solution component n at grid point j.
Definition Domain1D.cpp:317
void checkPointIndex(size_t n) const
Check that the specified point index is in range.
Definition Domain1D.h:209
double prevSoln(size_t n, size_t j) const
Value of component n at point j in the previous solution.
Definition Domain1D.h:708
double zmax() const
Get the coordinate [m] of the last (rightmost) grid point in this domain.
Definition Domain1D.h:741
size_t firstPoint() const
The index of the first (that is, left-most) grid point belonging to this domain.
Definition Domain1D.h:662
virtual shared_ptr< SolutionArray > toArray(bool normalize=false)
Save the state of this domain to a SolutionArray.
Definition Domain1D.h:575
void needJacUpdate()
Set this if something has changed in the governing equations (for example, the value of a constant ha...
Definition Domain1D.cpp:135
void linkRight(Domain1D *right)
Set the right neighbor to domain 'right.'.
Definition Domain1D.h:684
virtual void _getInitialSoln(double *x)
Writes some or all initial solution values into the global solution array, beginning at the location ...
Definition Domain1D.cpp:308
size_t index(size_t n, size_t j) const
Returns the index of the solution vector, which corresponds to component n at grid point j.
Definition Domain1D.h:408
void checkComponentArraySize(size_t nn) const
Check that an array size is at least nComponents().
Definition Domain1D.h:193
vector< double > getRefineCriteria()
Get the grid refinement criteria.
Definition Domain1D.cpp:303
double transient_rtol(size_t n)
Transient relative tolerance of the nth component.
Definition Domain1D.h:314
int m_bw
See bandwidth()
Definition Domain1D.h:858
virtual size_t loc(size_t j=0) const
Location of the start of the local solution vector in the global solution vector.
Definition Domain1D.h:657
void locate()
Find the index of the first grid point in this domain, and the start of its variables in the global s...
Definition Domain1D.cpp:203
virtual AnyMap getMeta() const
Retrieve meta data.
Definition Domain1D.cpp:143
virtual bool hasComponent(const string &name, bool checkAlias=true) const
Check whether the Domain contains a component.
Definition Domain1D.cpp:98
void setRefineCriteria(double ratio=10.0, double slope=0.8, double curve=0.8, double prune=-0.1)
Set grid refinement criteria.
Definition Domain1D.cpp:298
virtual void setFlatProfile(const string &component, double value)
Specify a flat profile for a component.
Definition Domain1D.h:547
virtual void getResiduals(const string &component, vector< double > &values) const
Retrieve internal work array values for a component.
Definition Domain1D.h:514
An array index is out of range.
An error indicating that an unimplemented function has been called.
Container class for multiple-domain 1D problems.
Definition OneDim.h:27
Refine Domain1D grids so that profiles satisfy adaptation tolerances.
Definition refine.h:17
A container class holding arrays of state information.
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
Namespace for the Cantera kernel.
Definition AnyMap.cpp:595
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:180
void warn_deprecated(const string &source, const AnyBase &node, const string &message)
A deprecation warning for syntax in an input file.
Definition AnyMap.cpp:1997