Cantera  3.0.0
Loading...
Searching...
No Matches
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 https://cantera.org/license.txt for license and copyright information.
7
10#include "cantera/oneD/refine.h"
11#include "cantera/base/AnyMap.h"
13
14namespace Cantera
15{
16
17Domain1D::Domain1D(size_t nv, size_t points, double time)
18{
19 resize(nv, points);
20}
21
22Domain1D::~Domain1D()
23{
24}
25
27{
28 warn_deprecated("Domain1D::domainType",
29 "To be changed after Cantera 3.0; for new behavior, see 'type'.");
30 return m_type;
31}
32
33void 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 = make_unique<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
56string 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
65size_t Domain1D::componentIndex(const 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
76void Domain1D::setTransientTolerances(double rtol, double 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
89void Domain1D::setSteadyTolerances(double rtol, double 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
111{
112 auto wrap_tols = [this](const vector<double>& tols) {
113 // If all tolerances are the same, just store the scalar value.
114 // Otherwise, store them by component name
115 set<double> unique_tols(tols.begin(), tols.end());
116 if (unique_tols.size() == 1) {
117 return AnyValue(tols[0]);
118 } else {
119 AnyMap out;
120 for (size_t i = 0; i < tols.size(); i++) {
121 out[componentName(i)] = tols[i];
122 }
123 return AnyValue(out);
124 }
125 };
126 AnyMap state;
127 state["type"] = type();
128 state["points"] = static_cast<long int>(nPoints());
129 if (nComponents() && nPoints()) {
130 state["tolerances"]["transient-abstol"] = wrap_tols(m_atol_ts);
131 state["tolerances"]["steady-abstol"] = wrap_tols(m_atol_ss);
132 state["tolerances"]["transient-reltol"] = wrap_tols(m_rtol_ts);
133 state["tolerances"]["steady-reltol"] = wrap_tols(m_rtol_ss);
134 }
135 return state;
136}
137
138AnyMap Domain1D::serialize(const double* soln) const
139{
140 warn_deprecated("Domain1D::serialize",
141 "To be removed after Cantera 3.0; superseded by asArray.");
142 AnyMap out;
143 auto arr = asArray(soln);
144 arr->writeEntry(out, "", "");
145 return out;
146}
147
148shared_ptr<SolutionArray> Domain1D::toArray(bool normalize) const {
149 if (!m_state) {
150 throw CanteraError("Domain1D::toArray",
151 "Domain needs to be installed in a container before calling asArray.");
152 }
153 auto ret = asArray(m_state->data() + m_iloc);
154 if (normalize) {
155 ret->normalize();
156 }
157 return ret;
158}
159
160void Domain1D::fromArray(const shared_ptr<SolutionArray>& arr)
161{
162 if (!m_state) {
163 throw CanteraError("Domain1D::fromArray",
164 "Domain needs to be installed in a container before calling fromArray.");
165 }
166 resize(nComponents(), arr->size());
167 m_container->resize();
168 fromArray(*arr, m_state->data() + m_iloc);
169 _finalize(m_state->data() + m_iloc);
170}
171
172void Domain1D::setMeta(const AnyMap& meta)
173{
174 auto set_tols = [&](const AnyValue& tols, const string& which, vector<double>& out)
175 {
176 if (!tols.hasKey(which)) {
177 return;
178 }
179 const auto& tol = tols[which];
180 if (tol.isScalar()) {
181 out.assign(nComponents(), tol.asDouble());
182 } else {
183 for (size_t i = 0; i < nComponents(); i++) {
184 string name = componentName(i);
185 if (tol.hasKey(name)) {
186 out[i] = tol[name].asDouble();
187 } else {
188 warn_user("Domain1D::setMeta", "No {} found for component '{}'",
189 which, name);
190 }
191 }
192 }
193 };
194
195 if (meta.hasKey("tolerances")) {
196 const auto& tols = meta["tolerances"];
197 set_tols(tols, "transient-abstol", m_atol_ts);
198 set_tols(tols, "transient-reltol", m_rtol_ts);
199 set_tols(tols, "steady-abstol", m_atol_ss);
200 set_tols(tols, "steady-reltol", m_rtol_ss);
201 }
202}
203
204void Domain1D::restore(const AnyMap& state, double* soln, int loglevel)
205{
206 warn_deprecated("Domain1D::restore",
207 "To be removed after Cantera 3.0; restore from SolutionArray instead.");
208 auto arr = SolutionArray::create(solution());
209 arr->readEntry(state, "", "");
210 fromArray(*arr, soln);
211}
212
214{
215 if (m_left) {
216 // there is a domain on the left, so the first grid point in this domain
217 // is one more than the last one on the left
218 m_jstart = m_left->lastPoint() + 1;
219
220 // the starting location in the solution vector
221 m_iloc = m_left->loc() + m_left->size();
222 } else {
223 // this is the left-most domain
224 m_jstart = 0;
225 m_iloc = 0;
226 }
227 // if there is a domain to the right of this one, then repeat this for it
228 if (m_right) {
229 m_right->locate();
230 }
231}
232
233void Domain1D::setupGrid(size_t n, const double* z)
234{
235 if (n > 1) {
236 resize(m_nv, n);
237 for (size_t j = 0; j < m_points; j++) {
238 m_z[j] = z[j];
239 }
240 }
241}
242
243void Domain1D::showSolution_s(std::ostream& s, const double* x)
244{
245 warn_deprecated("Domain1D::showSolution_s",
246 "To be removed after Cantera 3.0; replaced by 'show'.");
247 show(s, x);
248}
249
250void Domain1D::showSolution(const double* x)
251{
252 warn_deprecated("Domain1D::showSolution",
253 "To be removed after Cantera 3.0; replaced by 'show'.");
254 show(x);
255}
256
257void Domain1D::show(const double* x)
258{
259 size_t nn = m_nv/5;
260 for (size_t i = 0; i < nn; i++) {
261 writeline('-', 79, false, true);
262 writelog("\n z ");
263 for (size_t n = 0; n < 5; n++) {
264 writelog(" {:>10s} ", componentName(i*5 + n));
265 }
266 writeline('-', 79, false, true);
267 for (size_t j = 0; j < m_points; j++) {
268 writelog("\n {:10.4g} ", m_z[j]);
269 for (size_t n = 0; n < 5; n++) {
270 double v = value(x, i*5+n, j);
271 writelog(" {:10.4g} ", v);
272 }
273 }
274 writelog("\n");
275 }
276 size_t nrem = m_nv - 5*nn;
277 writeline('-', 79, false, true);
278 writelog("\n z ");
279 for (size_t n = 0; n < nrem; n++) {
280 writelog(" {:>10s} ", componentName(nn*5 + n));
281 }
282 writeline('-', 79, false, true);
283 for (size_t j = 0; j < m_points; j++) {
284 writelog("\n {:10.4g} ", m_z[j]);
285 for (size_t n = 0; n < nrem; n++) {
286 double v = value(x, nn*5+n, j);
287 writelog(" {:10.4g} ", v);
288 }
289 }
290 writelog("\n");
291}
292
293void Domain1D::setProfile(const string& name, double* values, double* soln)
294{
295 for (size_t n = 0; n < m_nv; n++) {
296 if (name == componentName(n)) {
297 for (size_t j = 0; j < m_points; j++) {
298 soln[index(n, j) + m_iloc] = values[j];
299 }
300 return;
301 }
302 }
303 throw CanteraError("Domain1D::setProfile", "unknown component: "+name);
304}
305
307{
308 for (size_t j = 0; j < m_points; j++) {
309 for (size_t n = 0; n < m_nv; n++) {
310 x[index(n,j)] = initialValue(n,j);
311 }
312 }
313}
314
315double Domain1D::initialValue(size_t n, size_t j)
316{
317 throw NotImplementedError("Domain1D::initialValue");
318}
319
320} // namespace
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:427
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
Definition AnyMap.cpp:1423
A wrapper for a variable whose type is determined at runtime.
Definition AnyMap.h:86
bool hasKey(const string &key) const
Returns true if this AnyValue is an AnyMap and that map contains a key with the given name.
Definition AnyMap.cpp:617
Base class for exceptions thrown by Cantera classes.
void setTransientTolerances(double rtol, double atol, size_t n=npos)
Set tolerances for time-stepping mode.
Definition Domain1D.cpp:76
size_t lastPoint() const
The index of the last (that is, right-most) grid point belonging to this domain.
Definition Domain1D.h:434
size_t m_iloc
Starting location within the solution vector for unknowns that correspond to this domain.
Definition Domain1D.h:594
size_t nComponents() const
Number of components at each grid point.
Definition Domain1D.h:160
double rtol(size_t n)
Relative tolerance of the nth component.
Definition Domain1D.h:239
AnyMap serialize(const double *soln) const
Save the state of this domain as an AnyMap.
Definition Domain1D.cpp:138
shared_ptr< Solution > solution() const
Return thermo/kinetics/transport manager used in the domain.
Definition Domain1D.h:400
virtual void setMeta(const AnyMap &meta)
Retrieve meta data.
Definition Domain1D.cpp:172
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:550
virtual string type() const
String indicating the domain implemented.
Definition Domain1D.h:62
size_t nPoints() const
Number of grid points in this domain.
Definition Domain1D.h:182
shared_ptr< vector< double > > m_state
data pointer shared from OneDim
Definition Domain1D.h:574
virtual void resize(size_t nv, size_t np)
Resize the domain to have nv components and np grid points.
Definition Domain1D.cpp:33
void setSteadyTolerances(double rtol, double atol, size_t n=npos)
Set tolerances for steady-state mode.
Definition Domain1D.cpp:89
virtual shared_ptr< SolutionArray > asArray(const double *soln) const
Save the state of this domain as a SolutionArray.
Definition Domain1D.h:354
virtual string componentName(size_t n) const
Name of the nth component. May be overloaded.
Definition Domain1D.cpp:56
double atol(size_t n)
Absolute tolerance of the nth component.
Definition Domain1D.h:244
void restore(const AnyMap &state, double *soln, int loglevel)
Restore the solution for this domain from an AnyMap.
Definition Domain1D.cpp:204
shared_ptr< SolutionArray > toArray(bool normalize=false) const
Save the state of this domain to a SolutionArray.
Definition Domain1D.cpp:148
size_t m_points
Number of grid points.
Definition Domain1D.h:578
virtual void showSolution(const double *x)
Print the solution.
Definition Domain1D.cpp:250
int domainType()
Domain type flag.
Definition Domain1D.cpp:26
virtual double initialValue(size_t n, size_t j)
Initial value of solution component n at grid point j.
Definition Domain1D.cpp:315
virtual size_t componentIndex(const string &name) const
index of component with name name.
Definition Domain1D.cpp:65
virtual void fromArray(SolutionArray &arr, double *soln)
Restore the solution for this domain from a SolutionArray.
Definition Domain1D.h:386
virtual void showSolution_s(std::ostream &s, const double *x)
Definition Domain1D.cpp:243
void needJacUpdate()
Set this if something has changed in the governing equations (for example, the value of a constant ha...
Definition Domain1D.cpp:102
virtual void _getInitialSoln(double *x)
Writes some or all initial solution values into the global solution array, beginning at the location ...
Definition Domain1D.cpp:306
Domain1D(size_t nv=1, size_t points=1, double time=0.0)
Constructor.
Definition Domain1D.cpp:17
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:418
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:213
virtual AnyMap getMeta() const
Retrieve meta data.
Definition Domain1D.cpp:110
virtual void show(std::ostream &s, const double *x)
Print the solution.
Definition Domain1D.h:500
virtual void setupGrid(size_t n, const double *z)
called to set up initial grid, and after grid refinement
Definition Domain1D.cpp:233
void setAge(int age)
Set the Jacobian age.
Definition MultiJac.h:59
An error indicating that an unimplemented function has been called.
virtual void resize()
Call after one or more grids has changed size, for example after being refined.
Definition OneDim.cpp:200
void saveStats()
Save statistics on function and Jacobian evaluation, and reset the counters.
Definition OneDim.cpp:169
static shared_ptr< SolutionArray > create(const shared_ptr< Solution > &sol, int size=0, const AnyMap &meta={})
Instantiate a new SolutionArray reference.
MultiJac & jacobian()
Return a reference to the Jacobian evaluator of an OneDim object.
Definition OneDim.cpp:133
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition global.h:175
void warn_user(const string &method, const string &msg, const Args &... args)
Print a user warning raised from method as CanteraWarning.
Definition global.h:267
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:195
void warn_deprecated(const string &source, const AnyBase &node, const string &message)
A deprecation warning for syntax in an input file.
Definition AnyMap.cpp:1926