Cantera  2.4.0
Domain1D.cpp
Go to the documentation of this file.
1 /**
2  * @file Domain1D.cpp
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 
10 #include "cantera/base/ctml.h"
11 
12 using namespace std;
13 
14 namespace Cantera
15 {
16 
17 Domain1D::Domain1D(size_t nv, size_t points, double time) :
18  m_rdt(0.0),
19  m_nv(0),
20  m_container(0),
21  m_index(npos),
22  m_type(0),
23  m_iloc(0),
24  m_jstart(0),
25  m_left(0),
26  m_right(0),
27  m_bw(-1),
28  m_force_full_update(false)
29 {
30  resize(nv, points);
31 }
32 
33 void Domain1D::resize(size_t nv, size_t np)
34 {
35  // if the number of components is being changed, then a
36  // new grid refiner is required.
37  if (nv != m_nv || !m_refiner) {
38  m_nv = nv;
39  m_refiner.reset(new Refiner(*this));
40  }
41  m_nv = nv;
42  m_name.resize(m_nv,"");
43  m_max.resize(m_nv, 0.0);
44  m_min.resize(m_nv, 0.0);
45  // Default error tolerances for all domains
46  m_rtol_ss.resize(m_nv, 1.0e-4);
47  m_atol_ss.resize(m_nv, 1.0e-9);
48  m_rtol_ts.resize(m_nv, 1.0e-4);
49  m_atol_ts.resize(m_nv, 1.0e-11);
50  m_points = np;
51  m_z.resize(np, 0.0);
52  m_slast.resize(m_nv * m_points, 0.0);
53  locate();
54 }
55 
56 std::string Domain1D::componentName(size_t n) const
57 {
58  if (m_name[n] != "") {
59  return m_name[n];
60  } else {
61  return fmt::format("component {}", n);
62  }
63 }
64 
65 size_t Domain1D::componentIndex(const std::string& name) const
66 {
67  for (size_t n = 0; n < nComponents(); n++) {
68  if (name == componentName(n)) {
69  return n;
70  }
71  }
72  throw CanteraError("Domain1D::componentIndex",
73  "no component named "+name);
74 }
75 
76 void Domain1D::setTransientTolerances(doublereal rtol, doublereal atol, size_t n)
77 {
78  if (n == npos) {
79  for (n = 0; n < m_nv; n++) {
80  m_rtol_ts[n] = rtol;
81  m_atol_ts[n] = atol;
82  }
83  } else {
84  m_rtol_ts[n] = rtol;
85  m_atol_ts[n] = atol;
86  }
87 }
88 
89 void Domain1D::setSteadyTolerances(doublereal rtol, doublereal atol, size_t n)
90 {
91  if (n == npos) {
92  for (n = 0; n < m_nv; n++) {
93  m_rtol_ss[n] = rtol;
94  m_atol_ss[n] = atol;
95  }
96  } else {
97  m_rtol_ss[n] = rtol;
98  m_atol_ss[n] = atol;
99  }
100 }
101 
103 {
104  if (m_container) {
105  m_container->jacobian().setAge(10000);
106  m_container->saveStats();
107  }
108 }
109 
110 XML_Node& Domain1D::save(XML_Node& o, const doublereal* const sol)
111 {
112  XML_Node& d = o.addChild("domain");
113  d.addAttribute("points", nPoints());
114  d.addAttribute("components", nComponents());
115  d.addAttribute("id", id());
116  addFloatArray(d, "abstol_transient", nComponents(), &m_atol_ts[0]);
117  addFloatArray(d, "reltol_transient", nComponents(), &m_rtol_ts[0]);
118  addFloatArray(d, "abstol_steady", nComponents(), &m_atol_ss[0]);
119  addFloatArray(d, "reltol_steady", nComponents(), &m_rtol_ss[0]);
120  return d;
121 }
122 
123 void Domain1D::restore(const XML_Node& dom, doublereal* soln, int loglevel)
124 {
125  vector_fp values;
126  vector<XML_Node*> nodes = dom.getChildren("floatArray");
127  for (size_t i = 0; i < nodes.size(); i++) {
128  string title = nodes[i]->attrib("title");
129  getFloatArray(*nodes[i], values, false);
130  if (values.size() != nComponents()) {
131  if (loglevel > 0) {
132  writelog("Warning: Domain1D::restore: Got an array of length {}"
133  " when one of length {} was expected. "
134  "Tolerances for individual species may not be preserved.\n",
135  values.size(), nComponents());
136  }
137  // The number of components will differ when restoring from a
138  // mechanism with a different number of species. Assuming that
139  // tolerances are the same for all species, we can just copy the
140  // tolerance from the last species.
141  if (!values.empty()) {
142  values.resize(nComponents(), values[values.size()-1]);
143  } else {
144  // If the tolerance vector is empty, just leave the defaults
145  // in place.
146  continue;
147  }
148  }
149  if (title == "abstol_transient") {
150  m_atol_ts = values;
151  } else if (title == "reltol_transient") {
152  m_rtol_ts = values;
153  } else if (title == "abstol_steady") {
154  m_atol_ss = values;
155  } else if (title == "reltol_steady") {
156  m_rtol_ss = values;
157  } else {
158  throw CanteraError("Domain1D::restore",
159  "Got an unexpected array, '" + title + "'");
160  }
161  }
162 }
163 
165 {
166  if (m_left) {
167  // there is a domain on the left, so the first grid point in this domain
168  // is one more than the last one on the left
169  m_jstart = m_left->lastPoint() + 1;
170 
171  // the starting location in the solution vector
172  m_iloc = m_left->loc() + m_left->size();
173  } else {
174  // this is the left-most domain
175  m_jstart = 0;
176  m_iloc = 0;
177  }
178  // if there is a domain to the right of this one, then repeat this for it
179  if (m_right) {
180  m_right->locate();
181  }
182 }
183 
184 void Domain1D::setupGrid(size_t n, const doublereal* z)
185 {
186  if (n > 1) {
187  resize(m_nv, n);
188  for (size_t j = 0; j < m_points; j++) {
189  m_z[j] = z[j];
190  }
191  }
192 }
193 
194 void Domain1D::showSolution(const doublereal* x)
195 {
196  size_t nn = m_nv/5;
197  for (size_t i = 0; i < nn; i++) {
198  writeline('-', 79, false, true);
199  writelog("\n z ");
200  for (size_t n = 0; n < 5; n++) {
201  writelog(" {:>10s} ", componentName(i*5 + n));
202  }
203  writeline('-', 79, false, true);
204  for (size_t j = 0; j < m_points; j++) {
205  writelog("\n {:10.4g} ", m_z[j]);
206  for (size_t n = 0; n < 5; n++) {
207  double v = value(x, i*5+n, j);
208  writelog(" {:10.4g} ", v);
209  }
210  }
211  writelog("\n");
212  }
213  size_t nrem = m_nv - 5*nn;
214  writeline('-', 79, false, true);
215  writelog("\n z ");
216  for (size_t n = 0; n < nrem; n++) {
217  writelog(" {:>10s} ", componentName(nn*5 + n));
218  }
219  writeline('-', 79, false, true);
220  for (size_t j = 0; j < m_points; j++) {
221  writelog("\n {:10.4g} ", m_z[j]);
222  for (size_t n = 0; n < nrem; n++) {
223  double v = value(x, nn*5+n, j);
224  writelog(" {:10.4g} ", v);
225  }
226  }
227  writelog("\n");
228 }
229 
230 void Domain1D::setProfile(const std::string& name, double* values, double* soln)
231 {
232  for (size_t n = 0; n < m_nv; n++) {
233  if (name == componentName(n)) {
234  for (size_t j = 0; j < m_points; j++) {
235  soln[index(n, j) + m_iloc] = values[j];
236  }
237  return;
238  }
239  }
240  throw CanteraError("Domain1D::setProfile", "unknown component: "+name);
241 }
242 
243 void Domain1D::_getInitialSoln(doublereal* x)
244 {
245  for (size_t j = 0; j < m_points; j++) {
246  for (size_t n = 0; n < m_nv; n++) {
247  x[index(n,j)] = initialValue(n,j);
248  }
249  }
250 }
251 
252 doublereal Domain1D::initialValue(size_t n, size_t j)
253 {
254  throw CanteraError("Domain1D::initialValue",
255  "base class method called!");
256  return 0.0;
257 }
258 
259 } // namespace
void addFloatArray(XML_Node &node, const std::string &title, const size_t n, const doublereal *const vals, const std::string &units, const std::string &type, const doublereal minval, const doublereal maxval)
This function adds a child node with the name, "floatArray", with a value consisting of a comma separ...
Definition: ctml.cpp:40
std::vector< XML_Node * > getChildren(const std::string &name) const
Get a vector of pointers to XML_Node containing all of the children of the current node which match t...
Definition: xml.cpp:864
virtual void _getInitialSoln(doublereal *x)
Writes some or all initial solution values into the global solution array, beginning at the location ...
Definition: Domain1D.cpp:243
CTML ("Cantera Markup Language") is the variant of XML that Cantera uses to store data...
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:164
size_t getFloatArray(const XML_Node &node, vector_fp &v, const bool convert, const std::string &unitsString, const std::string &nodeName)
This function reads the current node or a child node of the current node with the default name...
Definition: ctml.cpp:256
MultiJac & jacobian()
Return a reference to the Jacobian evaluator.
Definition: OneDim.cpp:101
size_t m_iloc
Starting location within the solution vector for unknowns that correspond to this domain...
Definition: Domain1D.h:502
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:153
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:97
STL namespace.
virtual void restore(const XML_Node &dom, doublereal *soln, int loglevel)
Restore the solution for this domain from an XML_Node.
Definition: Domain1D.cpp:123
virtual XML_Node & save(XML_Node &o, const doublereal *const sol)
Save the current solution for this domain into an XML_Node.
Definition: Domain1D.cpp:110
void setAge(int age)
Set the Jacobian age.
Definition: MultiJac.h:58
size_t lastPoint() const
The index of the last (i.e., right-most) grid point belonging to this domain.
Definition: Domain1D.h:369
virtual void resize(size_t nv, size_t np)
Definition: Domain1D.cpp:33
virtual void showSolution(const doublereal *x)
Print the solution.
Definition: Domain1D.cpp:194
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
void saveStats()
Save statistics on function and Jacobian evaluation, and reset the counters.
Definition: OneDim.cpp:137
virtual std::string componentName(size_t n) const
Name of the nth component. May be overloaded.
Definition: Domain1D.cpp:56
void addAttribute(const std::string &attrib, const std::string &value)
Add or modify an attribute of the current node.
Definition: xml.cpp:474
size_t componentIndex(const std::string &name) const
index of component with name name.
Definition: Domain1D.cpp:65
size_t nComponents() const
Number of components at each grid point.
Definition: Domain1D.h:134
void setSteadyTolerances(doublereal rtol, doublereal atol, size_t n=npos)
Set tolerances for steady-state mode.
Definition: Domain1D.cpp:89
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
doublereal atol(size_t n)
Absolute tolerance of the nth component.
Definition: Domain1D.h:218
Refine Domain1D grids so that profiles satisfy adaptation tolerances.
Definition: refine.h:16
doublereal rtol(size_t n)
Relative tolerance of the nth component.
Definition: Domain1D.h:213
void needJacUpdate()
Definition: Domain1D.cpp:102
size_t nPoints() const
Number of grid points in this domain.
Definition: Domain1D.h:156
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
void setTransientTolerances(doublereal rtol, doublereal atol, size_t n=npos)
Set tolerances for time-stepping mode.
Definition: Domain1D.cpp:76
virtual void setupGrid(size_t n, const doublereal *z)
called to set up initial grid, and after grid refinement
Definition: Domain1D.cpp:184
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:353
virtual doublereal initialValue(size_t n, size_t j)
Initial value of solution component n at grid point j.
Definition: Domain1D.cpp:252