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