Cantera  3.0.0
Loading...
Searching...
No Matches
Sim1D.cpp
Go to the documentation of this file.
1/**
2 * @file Sim1D.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/StFlow.h"
12#include "cantera/oneD/refine.h"
17#include <limits>
18#include <fstream>
19
20using namespace std;
21
22namespace Cantera
23{
24
25Sim1D::Sim1D(vector<shared_ptr<Domain1D>>& domains) :
26 OneDim(domains),
27 m_steady_callback(0)
28{
29 // resize the internal solution vector and the work array, and perform
30 // domain-specific initialization of the solution vector.
31 resize();
32 for (size_t n = 0; n < nDomains(); n++) {
33 domain(n)._getInitialSoln(m_state->data() + start(n));
34 }
35
36 // set some defaults
37 m_tstep = 1.0e-5;
38 m_steps = { 10 };
39}
40
41Sim1D::Sim1D(vector<Domain1D*>& domains) :
42 OneDim(domains),
43 m_steady_callback(0)
44{
45 warn_deprecated("Sim1D::Sim1D(vector<Domain1D*>&)",
46 "To be removed after Cantera 3.0; superseded by "
47 "Sim1D::Sim1D(vector<shared_ptr<Domain1D>>&).");
48
49 // resize the internal solution vector and the work array, and perform
50 // domain-specific initialization of the solution vector.
51 resize();
52 for (size_t n = 0; n < nDomains(); n++) {
53 domain(n)._getInitialSoln(m_state->data() + start(n));
54 }
55
56 // set some defaults
57 m_tstep = 1.0e-5;
58 m_steps = { 10 };
59}
60
61void Sim1D::setInitialGuess(const string& component, vector<double>& locs,
62 vector<double>& vals)
63{
64 for (size_t dom=0; dom<nDomains(); dom++) {
65 Domain1D& d = domain(dom);
66 size_t ncomp = d.nComponents();
67 for (size_t comp=0; comp<ncomp; comp++) {
68 if (d.componentName(comp)==component) {
69 setProfile(dom,comp,locs,vals);
70 }
71 }
72 }
73}
74
75void Sim1D::setValue(size_t dom, size_t comp, size_t localPoint, double value)
76{
77 size_t iloc = domain(dom).loc() + domain(dom).index(comp, localPoint);
78 AssertThrowMsg(iloc < m_state->size(), "Sim1D::setValue",
79 "Index out of bounds: {} > {}", iloc, m_state->size());
80 (*m_state)[iloc] = value;
81}
82
83double Sim1D::value(size_t dom, size_t comp, size_t localPoint) const
84{
85 size_t iloc = domain(dom).loc() + domain(dom).index(comp, localPoint);
86 AssertThrowMsg(iloc < m_state->size(), "Sim1D::value",
87 "Index out of bounds: {} > {}", iloc, m_state->size());
88 return (*m_state)[iloc];
89}
90
91double Sim1D::workValue(size_t dom, size_t comp, size_t localPoint) const
92{
93 size_t iloc = domain(dom).loc() + domain(dom).index(comp, localPoint);
94 AssertThrowMsg(iloc < m_state->size(), "Sim1D::workValue",
95 "Index out of bounds: {} > {}", iloc, m_state->size());
96 return m_xnew[iloc];
97}
98
99void Sim1D::setProfile(size_t dom, size_t comp,
100 const vector<double>& pos, const vector<double>& values)
101{
102 if (pos.front() != 0.0 || pos.back() != 1.0) {
103 throw CanteraError("Sim1D::setProfile",
104 "`pos` vector must span the range [0, 1]. Got a vector spanning "
105 "[{}, {}] instead.", pos.front(), pos.back());
106 }
107 Domain1D& d = domain(dom);
108 double z0 = d.zmin();
109 double z1 = d.zmax();
110 for (size_t n = 0; n < d.nPoints(); n++) {
111 double zpt = d.z(n);
112 double frac = (zpt - z0)/(z1 - z0);
113 double v = linearInterp(frac, pos, values);
114 setValue(dom, comp, n, v);
115 }
116}
117
118void Sim1D::save(const string& fname, const string& id,
119 const string& desc, int loglevel)
120{
121 warn_deprecated("Sim1D::save",
122 "To be removed after Cantera 3.0; use version without loglevel instead.");
123 save(fname, id, desc, true); // legacy version overwrites data
124 if (loglevel > 0) {
125 writelog("Solution saved to file '{}' as entry '{}'.\n", fname, id);
126 }
127}
128
129void Sim1D::save(const string& fname, const string& name, const string& desc,
130 bool overwrite, int compression, const string& basis)
131{
132 size_t dot = fname.find_last_of(".");
133 string extension = (dot != npos) ? toLowerCopy(fname.substr(dot+1)) : "";
134 if (extension == "csv") {
135 for (auto dom : m_dom) {
136 auto arr = dom->asArray(m_state->data() + dom->loc());
137 if (dom->size() > 1) {
138 arr->writeEntry(fname, overwrite, basis);
139 break;
140 }
141 }
142 return;
143 }
144 if (basis != "") {
145 warn_user("Sim1D::save",
146 "Species basis '{}' not implemented for HDF5 or YAML output.", basis);
147 }
148 if (extension == "h5" || extension == "hdf" || extension == "hdf5") {
149 SolutionArray::writeHeader(fname, name, desc, overwrite);
150 for (auto dom : m_dom) {
151 auto arr = dom->asArray(m_state->data() + dom->loc());
152 arr->writeEntry(fname, name, dom->id(), overwrite, compression);
153 }
154 return;
155 }
156 if (extension == "yaml" || extension == "yml") {
157 // Check for an existing file and load it if present
158 AnyMap data;
159 if (std::ifstream(fname).good()) {
160 data = AnyMap::fromYamlFile(fname);
161 }
162 SolutionArray::writeHeader(data, name, desc, overwrite);
163
164 for (auto dom : m_dom) {
165 auto arr = dom->asArray(m_state->data() + dom->loc());
166 arr->writeEntry(data, name, dom->id(), overwrite);
167 }
168
169 // Write the output file and remove the now-outdated cached file
170 std::ofstream out(fname);
171 out << data.toYamlString();
173 return;
174 }
175 throw CanteraError("Sim1D::save", "Unsupported file format '{}'.", extension);
176}
177
178void Sim1D::saveResidual(const string& fname, const string& id,
179 const string& desc, int loglevel)
180{
181 warn_deprecated("Sim1D::saveResidual",
182 "To be removed after Cantera 3.0; use version without loglevel instead.");
183 saveResidual(fname, id, desc, true); // legacy version overwrites data
184 if (loglevel > 0) {
185 writelog("Solution saved to file '{}' as entry '{}'.\n", fname, id);
186 }
187}
188
189void Sim1D::saveResidual(const string& fname, const string& name,
190 const string& desc, bool overwrite, int compression)
191{
192 vector<double> res(m_state->size(), -999);
193 OneDim::eval(npos, m_state->data(), &res[0], 0.0);
194 // Temporarily put the residual into m_state, since this is the vector that the
195 // save() function reads.
196 vector<double> backup(*m_state);
197 *m_state = res;
198 save(fname, name, desc, overwrite, compression);
199 *m_state = backup;
200}
201
202namespace { // restrict scope of helper function to local translation unit
203
204//! convert data format used by Python h5py export (Cantera < 3.0)
205AnyMap legacyH5(shared_ptr<SolutionArray> arr, const AnyMap& header={})
206{
207 auto meta = arr->meta();
208 AnyMap out;
209
210 map<string, string> meta_pairs = {
211 {"type", "Domain1D_type"},
212 {"name", "name"},
213 {"emissivity-left", "emissivity_left"},
214 {"emissivity-right", "emissivity_right"},
215 };
216 for (const auto& [newName, oldName] : meta_pairs) {
217 if (meta.hasKey(oldName)) {
218 out[newName] = meta[oldName];
219 }
220 }
221
222 map<string, string> tol_pairs = {
223 {"transient-abstol", "transient_abstol"},
224 {"steady-abstol", "steady_abstol"},
225 {"transient-reltol", "transient_reltol"},
226 {"steady-reltol", "steady_reltol"},
227 };
228 for (const auto& [newName, oldName] : tol_pairs) {
229 if (meta.hasKey(oldName)) {
230 out["tolerances"][newName] = meta[oldName];
231 }
232 }
233
234 if (meta.hasKey("phase")) {
235 out["phase"]["name"] = meta["phase"]["name"];
236 out["phase"]["source"] = meta["phase"]["source"];
237 }
238
239 if (arr->size() <= 1) {
240 return out;
241 }
242
243 map<string, string> header_pairs = {
244 {"transport-model", "transport_model"},
245 {"radiation-enabled", "radiation_enabled"},
246 {"energy-enabled", "energy_enabled"},
247 {"Soret-enabled", "soret_enabled"},
248 {"species-enabled", "species_enabled"},
249 };
250 for (const auto& [newName, oldName] : header_pairs) {
251 if (header.hasKey(oldName)) {
252 out[newName] = header[oldName];
253 }
254 }
255
256 map<string, string> refiner_pairs = {
257 {"ratio", "ratio"},
258 {"slope", "slope"},
259 {"curve", "curve"},
260 {"prune", "prune"},
261 // {"grid-min", "???"}, // missing
262 {"max-points", "max_grid_points"},
263 };
264 for (const auto& [newName, oldName] : refiner_pairs) {
265 if (header.hasKey(oldName)) {
266 out["refine-criteria"][newName] = header[oldName];
267 }
268 }
269
270 if (header.hasKey("fixed_temperature")) {
271 double temp = header.getDouble("fixed_temperature", -1.);
272 auto profile = arr->getComponent("T").as<vector<double>>();
273 int ix = 0;
274 while (profile[ix] <= temp && ix < arr->size()) {
275 ix++;
276 }
277 if (ix != 0) {
278 auto grid = arr->getComponent("grid").as<vector<double>>();
279 out["fixed-point"]["location"] = grid[ix - 1];
280 out["fixed-point"]["temperature"] = temp;
281 }
282 }
283
284 return out;
285}
286
287} // end unnamed namespace
288
289AnyMap Sim1D::restore(const string& fname, const string& id, int loglevel)
290{
291 warn_deprecated("Sim1D::saveResidual",
292 "To be removed after Cantera 3.0; use version without loglevel instead.");
293 return restore(fname, id);
294}
295
296AnyMap Sim1D::restore(const string& fname, const string& name)
297{
298 size_t dot = fname.find_last_of(".");
299 string extension = (dot != npos) ? toLowerCopy(fname.substr(dot+1)) : "";
300 if (extension == "xml") {
301 throw CanteraError("Sim1D::restore",
302 "Restoring from XML is no longer supported.");
303 }
304 AnyMap header;
305 if (extension == "h5" || extension == "hdf" || extension == "hdf5") {
306 map<string, shared_ptr<SolutionArray>> arrs;
307 header = SolutionArray::readHeader(fname, name);
308
309 for (auto dom : m_dom) {
310 auto arr = SolutionArray::create(dom->solution());
311 try {
312 arr->readEntry(fname, name, dom->id());
313 } catch (CanteraError& err) {
314 throw CanteraError("Sim1D::restore",
315 "Encountered exception when reading entry '{}' from '{}':\n{}",
316 name, fname, err.getMessage());
317 }
318 dom->resize(dom->nComponents(), arr->size());
319 if (!header.hasKey("generator")) {
320 arr->meta() = legacyH5(arr, header);
321 }
322 arrs[dom->id()] = arr;
323 }
324 resize();
325 m_xlast_ts.clear();
326 for (auto dom : m_dom) {
327 try {
328 dom->fromArray(*arrs[dom->id()], m_state->data() + dom->loc());
329 } catch (CanteraError& err) {
330 throw CanteraError("Sim1D::restore",
331 "Encountered exception when restoring domain '{}' from HDF:\n{}",
332 dom->id(), err.getMessage());
333 }
334 }
335 finalize();
336 } else if (extension == "yaml" || extension == "yml") {
337 AnyMap root = AnyMap::fromYamlFile(fname);
338 map<string, shared_ptr<SolutionArray>> arrs;
339 header = SolutionArray::readHeader(root, name);
340
341 for (auto dom : m_dom) {
342 auto arr = SolutionArray::create(dom->solution());
343 try {
344 arr->readEntry(root, name, dom->id());
345 } catch (CanteraError& err) {
346 throw CanteraError("Sim1D::restore",
347 "Encountered exception when reading entry '{}' from '{}':\n{}",
348 name, fname, err.getMessage());
349 }
350 dom->resize(dom->nComponents(), arr->size());
351 arrs[dom->id()] = arr;
352 }
353 resize();
354 m_xlast_ts.clear();
355 for (auto dom : m_dom) {
356 try {
357 dom->fromArray(*arrs[dom->id()], m_state->data() + dom->loc());
358 } catch (CanteraError& err) {
359 throw CanteraError("Sim1D::restore",
360 "Encountered exception when restoring domain '{}' from YAML:\n{}",
361 dom->id(), err.getMessage());
362 }
363 }
364 finalize();
365 } else {
366 throw CanteraError("Sim1D::restore",
367 "Unknown file extension '{}'; supported extensions include "
368 "'h5'/'hdf'/'hdf5' and 'yml'/'yaml'.", extension);
369 }
370 return header;
371}
372
373void Sim1D::setFlatProfile(size_t dom, size_t comp, double v)
374{
375 size_t np = domain(dom).nPoints();
376 for (size_t n = 0; n < np; n++) {
377 setValue(dom, comp, n, v);
378 }
379}
380
381void Sim1D::showSolution(ostream& s)
382{
383 warn_deprecated("Sim1D::showSolution",
384 "To be removed after Cantera 3.0; replaced by 'show'.");
385 show(s);
386}
387
389{
390 warn_deprecated("Sim1D::showSolution",
391 "To be removed after Cantera 3.0; replaced by 'show'.");
392 show();
393}
394
395void Sim1D::show(ostream& s)
396{
397 for (size_t n = 0; n < nDomains(); n++) {
398 if (domain(n).type() != "empty") {
399 domain(n).show(s, m_state->data() + start(n));
400 }
401 }
402}
403
405{
406 for (size_t n = 0; n < nDomains(); n++) {
407 if (domain(n).type() != "empty") {
408 writelog("\n\n>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> "+domain(n).id()
409 +" <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n\n");
410 domain(n).show(m_state->data() + start(n));
411 }
412 }
413}
414
416{
417 if (m_xlast_ts.empty()) {
418 throw CanteraError("Sim1D::restoreTimeSteppingSolution",
419 "No successful time steps taken on this grid.");
420 }
422}
423
425{
426 if (m_xlast_ss.empty()) {
427 throw CanteraError("Sim1D::restoreSteadySolution",
428 "No successful steady state solution");
429 }
431 for (size_t n = 0; n < nDomains(); n++) {
432 vector<double>& z = m_grid_last_ss[n];
433 domain(n).setupGrid(z.size(), z.data());
434 }
435}
436
437void Sim1D::getInitialSoln()
438{
439 for (size_t n = 0; n < nDomains(); n++) {
440 domain(n)._getInitialSoln(m_state->data() + start(n));
441 }
442}
443
445{
446 for (size_t n = 0; n < nDomains(); n++) {
447 domain(n)._finalize(m_state->data() + start(n));
448 }
449}
450
451void Sim1D::setTimeStep(double stepsize, size_t n, const int* tsteps)
452{
453 m_tstep = stepsize;
454 m_steps.resize(n);
455 for (size_t i = 0; i < n; i++) {
456 m_steps[i] = tsteps[i];
457 }
458}
459
460int Sim1D::newtonSolve(int loglevel)
461{
462 int m = OneDim::solve(m_state->data(), m_xnew.data(), loglevel);
463 if (m >= 0) {
464 *m_state = m_xnew;
465 return 0;
466 } else if (m > -10) {
467 return -1;
468 } else {
469 throw CanteraError("Sim1D::newtonSolve",
470 "ERROR: OneDim::solve returned m = {}", m);
471 }
472}
473
474void Sim1D::solve(int loglevel, bool refine_grid)
475{
476 int new_points = 1;
477 double dt = m_tstep;
478 m_nsteps = 0;
479 finalize();
480
481 while (new_points > 0) {
482 size_t istep = 0;
483 int nsteps = m_steps[istep];
484
485 bool ok = false;
486 if (loglevel > 0) {
487 writeline('.', 78, true, true);
488 }
489 while (!ok) {
490 // Attempt to solve the steady problem
492 newton().setOptions(m_ss_jac_age);
493 debuglog("Attempt Newton solution of steady-state problem...", loglevel);
494 int status = newtonSolve(loglevel-1);
495
496 if (status == 0) {
497 if (loglevel > 0) {
498 writelog(" success.\n\n");
499 writelog("Problem solved on [");
500 for (size_t mm = 1; mm < nDomains(); mm+=2) {
501 writelog("{}", domain(mm).nPoints());
502 if (mm + 2 < nDomains()) {
503 writelog(", ");
504 }
505 }
506 writelog("] point grid(s).\n");
507 }
508 if (m_steady_callback) {
510 }
511
512 if (loglevel > 6) {
513 save("debug_sim1d.yaml", "debug",
514 "After successful Newton solve");
515 }
516 if (loglevel > 7) {
517 saveResidual("debug_sim1d.yaml", "residual",
518 "After successful Newton solve");
519 }
520 ok = true;
521 } else {
522 debuglog(" failure. \n", loglevel);
523 if (loglevel > 6) {
524 save("debug_sim1d.yaml", "debug",
525 "After unsuccessful Newton solve");
526 }
527 if (loglevel > 7) {
528 saveResidual("debug_sim1d.yaml", "residual",
529 "After unsuccessful Newton solve");
530 }
531 if (loglevel > 0) {
532 writelog("Take {} timesteps ", nsteps);
533 }
534 dt = timeStep(nsteps, dt, m_state->data(), m_xnew.data(), loglevel-1);
536 if (loglevel > 6) {
537 save("debug_sim1d.yaml", "debug", "After timestepping");
538 }
539 if (loglevel > 7) {
540 saveResidual("debug_sim1d.yaml", "residual",
541 "After timestepping");
542 }
543
544 if (loglevel == 1) {
545 writelog(" {:10.4g} {:10.4g}\n", dt,
546 log10(ssnorm(m_state->data(), m_xnew.data())));
547 }
548 istep++;
549 if (istep >= m_steps.size()) {
550 nsteps = m_steps.back();
551 } else {
552 nsteps = m_steps[istep];
553 }
554 dt = std::min(dt, m_tmax);
555 }
556 }
557 if (loglevel > 0) {
558 writeline('.', 78, true, true);
559 }
560 if (loglevel > 2) {
561 show();
562 }
563
564 if (refine_grid) {
565 new_points = refine(loglevel);
566 if (new_points) {
567 // If the grid has changed, preemptively reduce the timestep
568 // to avoid multiple successive failed time steps.
569 dt = m_tstep;
570 }
571 if (new_points && loglevel > 6) {
572 save("debug_sim1d.yaml", "debug", "After regridding");
573 }
574 if (new_points && loglevel > 7) {
575 saveResidual("debug_sim1d.yaml", "residual",
576 "After regridding");
577 }
578 } else {
579 debuglog("grid refinement disabled.\n", loglevel);
580 new_points = 0;
581 }
582 }
583}
584
585int Sim1D::refine(int loglevel)
586{
587 int ianalyze, np = 0;
588 vector<double> znew, xnew;
589 vector<size_t> dsize;
590
592 m_grid_last_ss.clear();
593
594 for (size_t n = 0; n < nDomains(); n++) {
595 Domain1D& d = domain(n);
596 Refiner& r = d.refiner();
597
598 // Save the old grid corresponding to the converged solution
599 m_grid_last_ss.push_back(d.grid());
600
601 // determine where new points are needed
602 ianalyze = r.analyze(d.grid().size(), d.grid().data(),
603 m_state->data() + start(n));
604 if (ianalyze < 0) {
605 return ianalyze;
606 }
607
608 if (loglevel > 0) {
609 r.show();
610 }
611
612 np += r.nNewPoints();
613 size_t comp = d.nComponents();
614
615 // loop over points in the current grid
616 size_t npnow = d.nPoints();
617 size_t nstart = znew.size();
618 for (size_t m = 0; m < npnow; m++) {
619 if (r.keepPoint(m)) {
620 // add the current grid point to the new grid
621 znew.push_back(d.grid(m));
622
623 // do the same for the solution at this point
624 for (size_t i = 0; i < comp; i++) {
625 xnew.push_back(value(n, i, m));
626 }
627
628 // now check whether a new point is needed in the interval to
629 // the right of point m, and if so, add entries to znew and xnew
630 // for this new point
631 if (r.newPointNeeded(m) && m + 1 < npnow) {
632 // add new point at midpoint
633 double zmid = 0.5*(d.grid(m) + d.grid(m+1));
634 znew.push_back(zmid);
635 np++;
636
637 // for each component, linearly interpolate
638 // the solution to this point
639 for (size_t i = 0; i < comp; i++) {
640 double xmid = 0.5*(value(n, i, m) + value(n, i, m+1));
641 xnew.push_back(xmid);
642 }
643 }
644 } else {
645 if (loglevel > 0) {
646 writelog("refine: discarding point at {}\n", d.grid(m));
647 }
648 }
649 }
650 dsize.push_back(znew.size() - nstart);
651 }
652
653 // At this point, the new grid znew and the new solution vector xnew have
654 // been constructed, but the domains themselves have not yet been modified.
655 // Now update each domain with the new grid.
656
657 size_t gridstart = 0, gridsize;
658 for (size_t n = 0; n < nDomains(); n++) {
659 Domain1D& d = domain(n);
660 gridsize = dsize[n];
661 if (gridsize != 0) {
662 d.setupGrid(gridsize, &znew[gridstart]);
663 }
664 gridstart += gridsize;
665 }
666
667 // Replace the current solution vector with the new one
668 *m_state = xnew;
669 resize();
670 finalize();
671 return np;
672}
673
675{
676 int np = 0;
677 vector<double> znew, xnew;
678 double zfixed = 0.0;
679 double z1 = 0.0, z2 = 0.0;
680 vector<size_t> dsize;
681
682 for (size_t n = 0; n < nDomains(); n++) {
683 Domain1D& d = domain(n);
684 size_t comp = d.nComponents();
685 size_t mfixed = npos;
686
687 // loop over current grid to determine where new point is needed
688 StFlow* d_free = dynamic_cast<StFlow*>(&domain(n));
689 size_t npnow = d.nPoints();
690 size_t nstart = znew.size();
691 if (d_free && d_free->isFree()) {
692 for (size_t m = 0; m < npnow - 1; m++) {
693 bool fixedpt = false;
694 double t1 = value(n, c_offset_T, m);
695 double t2 = value(n, c_offset_T, m + 1);
696 // threshold to avoid adding new point too close to existing point
697 double thresh = min(1., 1.e-1 * (t2 - t1));
698 z1 = d.grid(m);
699 z2 = d.grid(m + 1);
700 if (fabs(t - t1) <= thresh) {
701 zfixed = z1;
702 fixedpt = true;
703 } else if (fabs(t2 - t) <= thresh) {
704 zfixed = z2;
705 fixedpt = true;
706 } else if ((t1 < t) && (t < t2)) {
707 mfixed = m;
708 zfixed = (z1 - z2) / (t1 - t2) * (t - t2) + z2;
709 fixedpt = true;
710 }
711
712 if (fixedpt) {
713 d_free->m_zfixed = zfixed;
714 d_free->m_tfixed = t;
715 break;
716 }
717 }
718 }
719
720 // copy solution domain and push back values
721 for (size_t m = 0; m < npnow; m++) {
722 // add the current grid point to the new grid
723 znew.push_back(d.grid(m));
724
725 // do the same for the solution at this point
726 for (size_t i = 0; i < comp; i++) {
727 xnew.push_back(value(n, i, m));
728 }
729 if (m == mfixed) {
730 // add new point at zfixed (mfixed is not npos)
731 znew.push_back(zfixed);
732 np++;
733 double interp_factor = (zfixed - z2) / (z1 - z2);
734 // for each component, linearly interpolate
735 // the solution to this point
736 for (size_t i = 0; i < comp; i++) {
737 double xmid = interp_factor*(value(n, i, m) - value(n, i, m+1)) + value(n,i,m+1);
738 xnew.push_back(xmid);
739 }
740 }
741 }
742 dsize.push_back(znew.size() - nstart);
743 }
744
745 // At this point, the new grid znew and the new solution vector xnew have
746 // been constructed, but the domains themselves have not yet been modified.
747 // Now update each domain with the new grid.
748 size_t gridstart = 0;
749 for (size_t n = 0; n < nDomains(); n++) {
750 Domain1D& d = domain(n);
751 size_t gridsize = dsize[n];
752 d.setupGrid(gridsize, &znew[gridstart]);
753 gridstart += gridsize;
754 }
755
756 // Replace the current solution vector with the new one
757 *m_state = xnew;
758
759 resize();
760 finalize();
761 return np;
762}
763
765{
766 double t_fixed = std::numeric_limits<double>::quiet_NaN();
767 for (size_t n = 0; n < nDomains(); n++) {
768 StFlow* d = dynamic_cast<StFlow*>(&domain(n));
769 if (d && d->isFree() && d->m_tfixed > 0) {
770 t_fixed = d->m_tfixed;
771 break;
772 }
773 }
774 return t_fixed;
775}
776
778{
779 double z_fixed = std::numeric_limits<double>::quiet_NaN();
780 for (size_t n = 0; n < nDomains(); n++) {
781 StFlow* d = dynamic_cast<StFlow*>(&domain(n));
782 if (d && d->isFree() && d->m_tfixed > 0) {
783 z_fixed = d->m_zfixed;
784 break;
785 }
786 }
787 return z_fixed;
788}
789
790void Sim1D::setRefineCriteria(int dom, double ratio,
791 double slope, double curve, double prune)
792{
793 if (dom >= 0) {
794 Refiner& r = domain(dom).refiner();
795 r.setCriteria(ratio, slope, curve, prune);
796 } else {
797 for (size_t n = 0; n < nDomains(); n++) {
798 Refiner& r = domain(n).refiner();
799 r.setCriteria(ratio, slope, curve, prune);
800 }
801 }
802}
803
804vector<double> Sim1D::getRefineCriteria(int dom)
805{
806 if (dom >= 0) {
807 Refiner& r = domain(dom).refiner();
808 return r.getCriteria();
809 } else {
810 throw CanteraError("Sim1D::getRefineCriteria",
811 "Must specify domain to get criteria from");
812 }
813}
814
815void Sim1D::setGridMin(int dom, double gridmin)
816{
817 if (dom >= 0) {
818 Refiner& r = domain(dom).refiner();
819 r.setGridMin(gridmin);
820 } else {
821 for (size_t n = 0; n < nDomains(); n++) {
822 Refiner& r = domain(n).refiner();
823 r.setGridMin(gridmin);
824 }
825 }
826}
827
828void Sim1D::setMaxGridPoints(int dom, int npoints)
829{
830 if (dom >= 0) {
831 Refiner& r = domain(dom).refiner();
832 r.setMaxPoints(npoints);
833 } else {
834 for (size_t n = 0; n < nDomains(); n++) {
835 Refiner& r = domain(n).refiner();
836 r.setMaxPoints(npoints);
837 }
838 }
839}
840
841size_t Sim1D::maxGridPoints(size_t dom)
842{
843 Refiner& r = domain(dom).refiner();
844 return r.maxPoints();
845}
846
847double Sim1D::jacobian(int i, int j)
848{
849 return OneDim::jacobian().value(i,j);
850}
851
852void Sim1D::evalSSJacobian()
853{
854 OneDim::evalSSJacobian(m_state->data(), m_xnew.data());
855}
856
857void Sim1D::solveAdjoint(const double* b, double* lambda)
858{
859 for (auto& D : m_dom) {
860 D->forceFullUpdate(true);
861 }
862 evalSSJacobian();
863 for (auto& D : m_dom) {
864 D->forceFullUpdate(false);
865 }
866
867 // Form J^T
868 size_t bw = bandwidth();
869 BandMatrix Jt(size(), bw, bw);
870 for (size_t i = 0; i < size(); i++) {
871 size_t j1 = (i > bw) ? i - bw : 0;
872 size_t j2 = (i + bw >= size()) ? size() - 1: i + bw;
873 for (size_t j = j1; j <= j2; j++) {
874 Jt(j,i) = m_jac->value(i,j);
875 }
876 }
877
878 Jt.solve(b, lambda);
879}
880
882{
884 m_xnew.resize(size(), 0.0);
885}
886
887}
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:427
double getDouble(const string &key, double default_) const
If key exists, return it as a double, otherwise return default_.
Definition AnyMap.cpp:1520
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
Definition AnyMap.cpp:1423
static void clearCachedFile(const string &filename)
Remove the specified file from the input cache if it is present.
Definition AnyMap.cpp:1747
static AnyMap fromYamlFile(const string &name, const string &parent_name="")
Create an AnyMap from a YAML file.
Definition AnyMap.cpp:1771
A class for banded matrices, involving matrix inversion processes.
Definition BandMatrix.h:37
int solve(const double *const b, double *const x)
Solve the matrix problem Ax = b.
double & value(size_t i, size_t j)
Return a changeable reference to element (i,j).
Base class for exceptions thrown by Cantera classes.
virtual string getMessage() const
Method overridden by derived classes to format the error message.
Base class for one-dimensional domains.
Definition Domain1D.h:41
size_t nComponents() const
Number of components at each grid point.
Definition Domain1D.h:160
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
Refiner & refiner()
Return a reference to the grid refiner.
Definition Domain1D.h:155
virtual string componentName(size_t n) const
Name of the nth component. May be overloaded.
Definition Domain1D.cpp:56
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
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
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
virtual double eval(double t) const
Evaluate the function.
Definition Func1.cpp:61
void setOptions(int maxJacAge=5)
Set options.
Definition MultiNewton.h:68
Container class for multiple-domain 1D problems.
Definition OneDim.h:27
int solve(double *x0, double *x1, int loglevel)
Solve F(x) = 0, where F(x) is the multi-domain residual function.
Definition OneDim.cpp:258
size_t start(size_t i) const
The index of the start of domain i in the solution vector.
Definition OneDim.h:95
int m_nsteps
Number of time steps taken in the current call to solve()
Definition OneDim.h:377
virtual void resize()
Call after one or more grids has changed size, for example after being refined.
Definition OneDim.cpp:200
size_t size() const
Total solution vector length;.
Definition OneDim.h:106
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:292
double ssnorm(double *x, double *r)
Steady-state max norm (infinity norm) of the residual evaluated using solution x.
Definition OneDim.cpp:324
size_t nDomains() const
Number of domains.
Definition OneDim.h:64
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:66
size_t bandwidth() const
Jacobian bandwidth.
Definition OneDim.h:136
shared_ptr< vector< double > > m_state
Solution vector.
Definition OneDim.h:342
unique_ptr< MultiJac > m_jac
Jacobian evaluator.
Definition OneDim.h:344
vector< Domain1D * > m_dom
Definition OneDim.h:356
double timeStep(int nsteps, double dt, double *x, double *r, int loglevel)
Take time steps using Backward Euler.
Definition OneDim.cpp:382
double m_tmax
maximum timestep size
Definition OneDim.h:337
void setSteadyMode()
Prepare to solve the steady-state problem.
Definition OneDim.cpp:353
Domain1D & domain(size_t i) const
Return a reference to domain i.
Definition OneDim.h:69
MultiNewton & newton()
Return a reference to the Newton iterator.
Definition OneDim.cpp:137
Refine Domain1D grids so that profiles satisfy adaptation tolerances.
Definition refine.h:17
size_t maxPoints() const
Returns the maximum number of points allowed in the domain.
Definition refine.h:57
vector< double > getCriteria()
Get the grid refinement criteria.
Definition refine.h:42
void setMaxPoints(int npmax)
Set the maximum number of points allowed in the domain.
Definition refine.h:52
void setCriteria(double ratio=10.0, double slope=0.8, double curve=0.8, double prune=-0.1)
Set grid refinement criteria.
Definition refine.cpp:21
void setGridMin(double gridmin)
Set the minimum allowable spacing between adjacent grid points [m].
Definition refine.h:62
void restoreTimeSteppingSolution()
Set the current solution vector to the last successful time-stepping solution.
Definition Sim1D.cpp:415
void resize() override
Call after one or more grids has changed size, for example after being refined.
Definition Sim1D.cpp:881
vector< double > m_xnew
a work array used to hold the residual or the new solution
Definition Sim1D.h:358
void saveResidual(const string &fname, const string &id, const string &desc, int loglevel)
Save the residual of the current solution to a container file.
Definition Sim1D.cpp:178
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:99
double fixedTemperatureLocation()
Return location of the point where temperature is fixed.
Definition Sim1D.cpp:777
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:355
void finalize()
Calls method _finalize in each domain.
Definition Sim1D.cpp:444
void setValue(size_t dom, size_t comp, size_t localPoint, double value)
Set a single value in the solution vector.
Definition Sim1D.cpp:75
int refine(int loglevel=0)
Refine the grid in all domains.
Definition Sim1D.cpp:585
void show()
Show logging information on current solution for all domains.
Definition Sim1D.cpp:404
double fixedTemperature()
Return temperature at the point used to fix the flame location.
Definition Sim1D.cpp:764
vector< double > m_xlast_ss
the solution vector after the last successful steady-state solve (stored before grid refinement)
Definition Sim1D.h:351
void save(const string &fname, const string &id, const string &desc, int loglevel)
Save the current solution to a container file.
Definition Sim1D.cpp:118
void setMaxGridPoints(int dom, int npoints)
Set the maximum number of grid points in the domain.
Definition Sim1D.cpp:828
int setFixedTemperature(double t)
Add node for fixed temperature point of freely propagating flame.
Definition Sim1D.cpp:674
AnyMap restore(const string &fname, const string &id, int loglevel)
Initialize the solution with a previously-saved solution.
Definition Sim1D.cpp:289
void setInitialGuess(const string &component, vector< double > &locs, vector< double > &vals)
Set initial guess for one component for all domains.
Definition Sim1D.cpp:61
void showSolution()
Show logging information on current solution for all domains.
Definition Sim1D.cpp:388
int newtonSolve(int loglevel)
Wrapper around the Newton solver.
Definition Sim1D.cpp:460
vector< int > m_steps
array of number of steps to take before re-attempting the steady-state solution
Definition Sim1D.h:365
double m_tstep
timestep
Definition Sim1D.h:361
void solveAdjoint(const double *b, double *lambda)
Solve the equation .
Definition Sim1D.cpp:857
Func1 * m_steady_callback
User-supplied function called after a successful steady-state solve.
Definition Sim1D.h:368
void restoreSteadySolution()
Set the current solution vector and grid to the last successful steady- state solution.
Definition Sim1D.cpp:424
size_t maxGridPoints(size_t dom)
Get the maximum number of grid points in this domain.
Definition Sim1D.cpp:841
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:373
double value(size_t dom, size_t comp, size_t localPoint) const
Get one entry in the solution vector.
Definition Sim1D.cpp:83
vector< double > getRefineCriteria(int dom)
Get the grid refinement criteria.
Definition Sim1D.cpp:804
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:815
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:790
vector< double > m_xlast_ts
the solution vector after the last successful timestepping
Definition Sim1D.h:347
static AnyMap readHeader(const string &fname, const string &name)
Read header information from a HDF container file.
static void writeHeader(const string &fname, const string &name, const string &desc, bool overwrite=false)
Write header data to a HDF container file.
static shared_ptr< SolutionArray > create(const shared_ptr< Solution > &sol, int size=0, const AnyMap &meta={})
Instantiate a new SolutionArray reference.
This class represents 1D flow domains that satisfy the one-dimensional similarity solution for chemic...
Definition StFlow.h:45
double m_tfixed
Temperature at the point used to fix the flame location.
Definition StFlow.h:566
double m_zfixed
Location of the point where temperature is fixed.
Definition StFlow.h:563
bool isFree() const
Retrieve flag indicating whether flow is freely propagating.
Definition StFlow.h:304
Header for a file containing miscellaneous numerical functions.
string toLowerCopy(const string &input)
Convert to lower case.
MultiJac & jacobian()
Return a reference to the Jacobian evaluator of an OneDim object.
Definition OneDim.cpp:133
#define AssertThrowMsg(expr, procedure,...)
Assertion must be true or an error is thrown.
void debuglog(const string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
Definition global.h:158
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition global.h:175
double dot(InputIter x_begin, InputIter x_end, InputIter2 y_begin)
Function that calculates a templated inner product.
Definition utilities.h:82
double linearInterp(double x, const vector< double > &xpts, const vector< double > &fpts)
Linearly interpolate a function defined on a discrete grid.
Definition funcs.cpp:13
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
@ c_offset_T
temperature
Definition StFlow.h:27
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
Contains declarations for string manipulation functions within Cantera.