Cantera  3.2.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/Flow1D.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 };
207 for (const auto& [newName, oldName] : header_pairs) {
208 if (header.hasKey(oldName)) {
209 out[newName] = header[oldName];
210 }
211 }
212
213 map<string, string> refiner_pairs = {
214 {"ratio", "ratio"},
215 {"slope", "slope"},
216 {"curve", "curve"},
217 {"prune", "prune"},
218 // {"grid-min", "???"}, // missing
219 {"max-points", "max_grid_points"},
220 };
221 for (const auto& [newName, oldName] : refiner_pairs) {
222 if (header.hasKey(oldName)) {
223 out["refine-criteria"][newName] = header[oldName];
224 }
225 }
226
227 if (header.hasKey("fixed_temperature")) {
228 double temp = header.getDouble("fixed_temperature", -1.);
229 auto profile = arr->getComponent("T").as<vector<double>>();
230 int ix = 0;
231 while (profile[ix] <= temp && ix < arr->size()) {
232 ix++;
233 }
234 if (ix != 0) {
235 auto grid = arr->getComponent("grid").as<vector<double>>();
236 out["fixed-point"]["location"] = grid[ix - 1];
237 out["fixed-point"]["temperature"] = temp;
238 }
239 }
240
241 return out;
242}
243
244} // end unnamed namespace
245
246AnyMap Sim1D::restore(const string& fname, const string& name)
247{
248 size_t dot = fname.find_last_of(".");
249 string extension = (dot != npos) ? toLowerCopy(fname.substr(dot+1)) : "";
250 if (extension == "xml") {
251 throw CanteraError("Sim1D::restore",
252 "Restoring from XML is no longer supported.");
253 }
254 AnyMap header;
255 if (extension == "h5" || extension == "hdf" || extension == "hdf5") {
256 map<string, shared_ptr<SolutionArray>> arrs;
257 header = SolutionArray::readHeader(fname, name);
258
259 for (auto dom : m_dom) {
260 auto arr = SolutionArray::create(dom->solution());
261 try {
262 arr->readEntry(fname, name, dom->id());
263 } catch (CanteraError& err) {
264 throw CanteraError("Sim1D::restore",
265 "Encountered exception when reading entry '{}' from '{}':\n{}",
266 name, fname, err.getMessage());
267 }
268 dom->resize(dom->nComponents(), arr->size());
269 if (!header.hasKey("generator")) {
270 arr->meta() = legacyH5(arr, header);
271 }
272 arrs[dom->id()] = arr;
273 }
274 resize();
275 m_xlast_ts.clear();
276 for (auto dom : m_dom) {
277 try {
278 dom->fromArray(*arrs[dom->id()], m_state->data() + dom->loc());
279 } catch (CanteraError& err) {
280 throw CanteraError("Sim1D::restore",
281 "Encountered exception when restoring domain '{}' from HDF:\n{}",
282 dom->id(), err.getMessage());
283 }
284 }
285 finalize();
286 } else if (extension == "yaml" || extension == "yml") {
287 AnyMap root = AnyMap::fromYamlFile(fname);
288 map<string, shared_ptr<SolutionArray>> arrs;
289 header = SolutionArray::readHeader(root, name);
290
291 for (auto dom : m_dom) {
292 auto arr = SolutionArray::create(dom->solution());
293 try {
294 arr->readEntry(root, name, dom->id());
295 } catch (CanteraError& err) {
296 throw CanteraError("Sim1D::restore",
297 "Encountered exception when reading entry '{}' from '{}':\n{}",
298 name, fname, err.getMessage());
299 }
300 dom->resize(dom->nComponents(), arr->size());
301 arrs[dom->id()] = arr;
302 }
303 resize();
304 m_xlast_ts.clear();
305 for (auto dom : m_dom) {
306 try {
307 dom->fromArray(*arrs[dom->id()], m_state->data() + dom->loc());
308 } catch (CanteraError& err) {
309 throw CanteraError("Sim1D::restore",
310 "Encountered exception when restoring domain '{}' from YAML:\n{}",
311 dom->id(), err.getMessage());
312 }
313 }
314 finalize();
315 } else {
316 throw CanteraError("Sim1D::restore",
317 "Unknown file extension '{}'; supported extensions include "
318 "'h5'/'hdf'/'hdf5' and 'yml'/'yaml'.", extension);
319 }
320 return header;
321}
322
323void Sim1D::setFlatProfile(size_t dom, size_t comp, double v)
324{
325 size_t np = domain(dom).nPoints();
326 for (size_t n = 0; n < np; n++) {
327 setValue(dom, comp, n, v);
328 }
329}
330
332{
333 for (size_t n = 0; n < nDomains(); n++) {
334 if (domain(n).type() != "empty") {
335 writelog("\n\n>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> "+domain(n).id()
336 +" <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n\n");
337 domain(n).show(m_state->data() + start(n));
338 }
339 }
340}
341
343{
344 if (m_xlast_ts.empty()) {
345 throw CanteraError("Sim1D::restoreTimeSteppingSolution",
346 "No successful time steps taken on this grid.");
347 }
349}
350
352{
353 if (m_xlast_ss.empty()) {
354 throw CanteraError("Sim1D::restoreSteadySolution",
355 "No successful steady state solution");
356 }
358 for (size_t n = 0; n < nDomains(); n++) {
359 vector<double>& z = m_grid_last_ss[n];
360 domain(n).setupGrid(z.size(), z.data());
361 }
362}
363
365{
366 for (size_t n = 0; n < nDomains(); n++) {
367 domain(n)._getInitialSoln(m_state->data() + start(n));
368 }
369}
370
372{
373 for (size_t n = 0; n < nDomains(); n++) {
374 domain(n)._finalize(m_state->data() + start(n));
375 }
376}
377
378void Sim1D::setTimeStep(double stepsize, size_t n, const int* tsteps)
379{
380 m_tstep = stepsize;
381 m_steps.resize(n);
382 for (size_t i = 0; i < n; i++) {
383 m_steps[i] = tsteps[i];
384 }
385}
386
387int Sim1D::newtonSolve(int loglevel)
388{
389 int m = OneDim::solve(m_state->data(), m_xnew.data(), loglevel);
390 if (m >= 0) {
391 *m_state = m_xnew;
392 return 0;
393 } else if (m > -10) {
394 return -1;
395 } else {
396 throw CanteraError("Sim1D::newtonSolve",
397 "ERROR: OneDim::solve returned m = {}", m);
398 }
399}
400
401void Sim1D::solve(int loglevel, bool refine_grid)
402{
403 int new_points = 1;
404 double dt = m_tstep;
405 m_nsteps = 0;
406 finalize();
407
408 // For debugging outputs
409 int attempt_counter = 0;
410 const int max_history = 10; // Store up to 10 previous solutions
411 if (loglevel > 6) {
413 }
414
415 while (new_points > 0) {
416 size_t istep = 0;
417 int nsteps = m_steps[istep];
418
419 bool ok = false;
420 if (loglevel > 0) {
421 writeline('.', 78, true, true);
422 }
423 while (!ok) {
424 // Keep the attempt_counter in the range of [1, max_history]
425 attempt_counter = (attempt_counter % max_history) + 1;
426
427 // Attempt to solve the steady problem
430 debuglog("\nAttempt Newton solution of steady-state problem.", loglevel);
431 int status = newtonSolve(loglevel);
432
433 if (status == 0) {
434 if (loglevel > 0) {
435 writelog("\nNewton steady-state solve succeeded.\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 writeDebugInfo("NewtonSuccess", "After successful Newton solve",
449 loglevel, attempt_counter);
450
451 ok = true;
452 } else {
453 debuglog("\nNewton steady-state solve failed.\n", loglevel);
454 writeDebugInfo("NewtonFail", "After unsuccessful Newton solve",
455 loglevel, attempt_counter);
456
457 if (loglevel > 0) {
458 writelog("\nAttempt {} timesteps.", nsteps);
459 }
460
461 dt = timeStep(nsteps, dt, m_state->data(), m_xnew.data(), loglevel-1);
463 writeDebugInfo("Timestepping", "After timestepping", loglevel,
464 attempt_counter);
465
466 // Repeat the last timestep's data for logging purposes
467 if (loglevel == 1) {
468 writelog("\nFinal timestep info: dt= {:<10.4g} log(ss)= {:<10.4g}\n", dt,
469 log10(ssnorm(m_state->data(), m_xnew.data())));
470 }
471 istep++;
472 if (istep >= m_steps.size()) {
473 nsteps = m_steps.back();
474 } else {
475 nsteps = m_steps[istep];
476 }
477 dt = std::min(dt, m_tmax);
478 }
479 }
480 if (loglevel > 0) {
481 writeline('.', 78, true, true);
482 }
483 if (loglevel > 3) {
484 show();
485 }
486
487 if (refine_grid) {
488 new_points = refine(loglevel);
489 if (new_points) {
490 // If the grid has changed, preemptively reduce the timestep
491 // to avoid multiple successive failed time steps.
492 dt = m_tstep;
493 }
494 writeDebugInfo("Regridding", "After regridding", loglevel, attempt_counter);
495 } else {
496 debuglog("grid refinement disabled.\n", loglevel);
497 new_points = 0;
498 }
499 }
500 if (new_points < 0) {
501 // If the solver finished after removing grid points, do one final evaluation
502 // of the governing equations to update internal arrays in each domain that may
503 // be used for data saved in output files.
504 for (auto dom : m_dom) {
505 dom->eval(npos, m_state->data() + dom->loc(), m_xnew.data() + dom->loc(),
506 m_mask.data());
507 }
508 }
509}
510
511int Sim1D::refine(int loglevel)
512{
513 int added = 0;
514 int discarded = 0;
515 vector<double> znew, xnew;
516 vector<size_t> dsize;
517
519 m_grid_last_ss.clear();
520
521 for (size_t n = 0; n < nDomains(); n++) {
522 Domain1D& d = domain(n);
523 Refiner& r = d.refiner();
524
525 // Save the old grid corresponding to the converged solution
526 m_grid_last_ss.push_back(d.grid());
527
528 // determine where new points are needed
529 r.analyze(d.grid().size(), d.grid().data(), m_state->data() + start(n));
530
531 if (loglevel > 0) {
532 r.show();
533 }
534
535 added += r.nNewPoints();
536 size_t comp = d.nComponents();
537
538 // loop over points in the current grid
539 size_t npnow = d.nPoints();
540 size_t nstart = znew.size();
541 for (size_t m = 0; m < npnow; m++) {
542 if (r.keepPoint(m)) {
543 // add the current grid point to the new grid
544 znew.push_back(d.z(m));
545
546 // do the same for the solution at this point
547 for (size_t i = 0; i < comp; i++) {
548 xnew.push_back(value(n, i, m));
549 }
550
551 // now check whether a new point is needed in the interval to
552 // the right of point m, and if so, add entries to znew and xnew
553 // for this new point
554 if (r.newPointNeeded(m) && m + 1 < npnow) {
555 // add new point at midpoint
556 double zmid = 0.5*(d.z(m) + d.z(m+1));
557 znew.push_back(zmid);
558 added++;
559
560 // for each component, linearly interpolate the solution to this
561 // point
562 for (size_t i = 0; i < comp; i++) {
563 double xmid = 0.5*(value(n, i, m) + value(n, i, m+1));
564 xnew.push_back(xmid);
565 }
566 }
567 } else {
568 discarded++;
569 if (loglevel > 0) {
570 writelog("refine: discarding point at {}\n", d.z(m));
571 }
572 }
573 }
574 dsize.push_back(znew.size() - nstart);
575 }
576
577 // At this point, the new grid znew and the new solution vector xnew have
578 // been constructed, but the domains themselves have not yet been modified.
579 // Now update each domain with the new grid.
580
581 size_t gridstart = 0, gridsize;
582 for (size_t n = 0; n < nDomains(); n++) {
583 Domain1D& d = domain(n);
584 gridsize = dsize[n];
585 if (gridsize != 0) {
586 d.setupGrid(gridsize, &znew[gridstart]);
587 }
588 gridstart += gridsize;
589 }
590
591 // Replace the current solution vector with the new one
592 *m_state = xnew;
593 resize();
594 finalize();
595 return added || -discarded;
596}
597
599{
600 std::filesystem::remove("debug_sim1d.yaml");
601}
602
603void Sim1D::writeDebugInfo(const string& header_suffix, const string& message,
604 int loglevel, int attempt_counter)
605{
606 string file_header;
607 if (loglevel > 6) {
608 file_header = fmt::format("solution_{}_{}", attempt_counter, header_suffix);
609 save("debug_sim1d.yaml", file_header, message, true);
610 }
611 if (loglevel > 7) {
612 file_header = fmt::format("residual_{}_{}", attempt_counter, header_suffix);
613 saveResidual("debug_sim1d.yaml", file_header, message, true);
614 }
615}
616
618{
619 int np = 0;
620 vector<double> znew, xnew;
621 double zfixed = 0.0;
622 double z1 = 0.0, z2 = 0.0;
623 vector<size_t> dsize;
624
625 for (size_t n = 0; n < nDomains(); n++) {
626 Domain1D& d = domain(n);
627 size_t comp = d.nComponents();
628 size_t mfixed = npos;
629
630 // loop over current grid to determine where new point is needed
631 Flow1D* d_free = dynamic_cast<Flow1D*>(&domain(n));
632 size_t npnow = d.nPoints();
633 size_t nstart = znew.size();
634 if (d_free && d_free->isFree()) {
635 for (size_t m = 0; m < npnow - 1; m++) {
636 bool fixedpt = false;
637 double t1 = value(n, c_offset_T, m);
638 double t2 = value(n, c_offset_T, m + 1);
639 // threshold to avoid adding new point too close to existing point
640 double thresh = min(1., 1.e-1 * (t2 - t1));
641 z1 = d.z(m);
642 z2 = d.z(m + 1);
643 if (fabs(t - t1) <= thresh) {
644 zfixed = z1;
645 fixedpt = true;
646 } else if (fabs(t2 - t) <= thresh) {
647 zfixed = z2;
648 fixedpt = true;
649 } else if ((t1 < t) && (t < t2)) {
650 mfixed = m;
651 zfixed = (z1 - z2) / (t1 - t2) * (t - t2) + z2;
652 fixedpt = true;
653 }
654
655 if (fixedpt) {
656 d_free->m_zfixed = zfixed;
657 d_free->m_tfixed = t;
658 break;
659 }
660 }
661 }
662
663 // copy solution domain and push back values
664 for (size_t m = 0; m < npnow; m++) {
665 // add the current grid point to the new grid
666 znew.push_back(d.z(m));
667
668 // do the same for the solution at this point
669 for (size_t i = 0; i < comp; i++) {
670 xnew.push_back(value(n, i, m));
671 }
672 if (m == mfixed) {
673 // add new point at zfixed (mfixed is not npos)
674 znew.push_back(zfixed);
675 np++;
676 double interp_factor = (zfixed - z2) / (z1 - z2);
677 // for each component, linearly interpolate
678 // the solution to this point
679 for (size_t i = 0; i < comp; i++) {
680 double xmid = interp_factor*(value(n, i, m) - value(n, i, m+1)) + value(n,i,m+1);
681 xnew.push_back(xmid);
682 }
683 }
684 }
685 dsize.push_back(znew.size() - nstart);
686 }
687
688 // At this point, the new grid znew and the new solution vector xnew have
689 // been constructed, but the domains themselves have not yet been modified.
690 // Now update each domain with the new grid.
691 size_t gridstart = 0;
692 for (size_t n = 0; n < nDomains(); n++) {
693 Domain1D& d = domain(n);
694 size_t gridsize = dsize[n];
695 d.setupGrid(gridsize, &znew[gridstart]);
696 gridstart += gridsize;
697 }
698
699 // Replace the current solution vector with the new one
700 *m_state = xnew;
701
702 resize();
703 finalize();
704 return np;
705}
706
708{
709 double t_fixed = std::numeric_limits<double>::quiet_NaN();
710 for (size_t n = 0; n < nDomains(); n++) {
711 Flow1D* d = dynamic_cast<Flow1D*>(&domain(n));
712 if (d && d->isFree() && d->m_tfixed > 0) {
713 t_fixed = d->m_tfixed;
714 break;
715 }
716 }
717 return t_fixed;
718}
719
721{
722 double z_fixed = std::numeric_limits<double>::quiet_NaN();
723 for (size_t n = 0; n < nDomains(); n++) {
724 Flow1D* d = dynamic_cast<Flow1D*>(&domain(n));
725 if (d && d->isFree() && d->m_tfixed > 0) {
726 z_fixed = d->m_zfixed;
727 break;
728 }
729 }
730 return z_fixed;
731}
732
733void Sim1D::setLeftControlPoint(double temperature)
734{
735 bool two_point_domain_found = false;
736 for (size_t n = 0; n < nDomains(); n++) {
737 Domain1D& d = domain(n);
738
739 // Skip if the domain type doesn't match
740 if (d.domainType() != "axisymmetric-flow") {
741 continue;
742 }
743
744 Flow1D& d_axis = dynamic_cast<Flow1D&>(domain(n));
745 size_t np = d_axis.nPoints();
746
747 // Skip if two-point control is not enabled
748 if (!d_axis.twoPointControlEnabled()) {
749 continue;
750 }
751 two_point_domain_found = true;
752
753 double current_val, next_val;
754 for (size_t m = 0; m < np-1; m++) {
755 current_val = value(n,c_offset_T,m);
756 next_val = value(n,c_offset_T,m+1);
757 if ((current_val - temperature) * (next_val - temperature) < 0.0) {
758 // Pick the coordinate of the point with the temperature closest
759 // to the desired temperature
760 size_t index = 0;
761 if (std::abs(current_val - temperature) <
762 std::abs(next_val - temperature)) {
763 index = m;
764 } else {
765 index = m+1;
766 }
767 d_axis.setLeftControlPointCoordinate(d_axis.z(index));
769 return;
770 }
771 }
772 }
773
774 if (!two_point_domain_found) {
775 throw CanteraError("Sim1D::setLeftControlPoint",
776 "No domain with two-point control enabled was found.");
777 } else {
778 throw CanteraError("Sim1D::setLeftControlPoint",
779 "No control point with temperature {} was able to be found in the"
780 "flame's temperature range.", temperature);
781 }
782}
783
784void Sim1D::setRightControlPoint(double temperature)
785{
786 bool two_point_domain_found = false;
787 for (size_t n = 0; n < nDomains(); n++) {
788 Domain1D& d = domain(n);
789
790 // Skip if the domain type doesn't match
791 if (d.domainType() != "axisymmetric-flow") {
792 continue;
793 }
794
795 Flow1D& d_axis = dynamic_cast<Flow1D&>(domain(n));
796 size_t np = d_axis.nPoints();
797
798 // Skip if two-point control is not enabled
799 if (!d_axis.twoPointControlEnabled()) {
800 continue;
801 }
802 two_point_domain_found = true;
803
804 double current_val, next_val;
805 for (size_t m = np-1; m > 0; m--) {
806 current_val = value(n,c_offset_T,m);
807 next_val = value(n,c_offset_T,m-1);
808 if ((current_val - temperature) * (next_val - temperature) < 0.0) {
809 // Pick the coordinate of the point with the temperature closest
810 // to the desired temperature
811 size_t index = 0;
812 if (std::abs(current_val - temperature) <
813 std::abs(next_val - temperature)) {
814 index = m;
815 } else {
816 index = m-1;
817 }
818 d_axis.setRightControlPointCoordinate(d_axis.z(index));
820 return;
821 }
822 }
823 }
824
825 if (!two_point_domain_found) {
826 throw CanteraError("Sim1D::setRightControlPoint",
827 "No domain with two-point control enabled was found.");
828 } else {
829 throw CanteraError("Sim1D::setRightControlPoint",
830 "No control point with temperature {} was able to be found in the"
831 "flame's temperature range.", temperature);
832 }
833
834}
835
836void Sim1D::setRefineCriteria(int dom, double ratio,
837 double slope, double curve, double prune)
838{
839 if (dom >= 0) {
840 Refiner& r = domain(dom).refiner();
841 r.setCriteria(ratio, slope, curve, prune);
842 } else {
843 for (size_t n = 0; n < nDomains(); n++) {
844 Refiner& r = domain(n).refiner();
845 r.setCriteria(ratio, slope, curve, prune);
846 }
847 }
848}
849
850vector<double> Sim1D::getRefineCriteria(int dom)
851{
852 if (dom >= 0) {
853 Refiner& r = domain(dom).refiner();
854 return r.getCriteria();
855 } else {
856 throw CanteraError("Sim1D::getRefineCriteria",
857 "Must specify domain to get criteria from");
858 }
859}
860
861void Sim1D::setGridMin(int dom, double gridmin)
862{
863 if (dom >= 0) {
864 Refiner& r = domain(dom).refiner();
865 r.setGridMin(gridmin);
866 } else {
867 for (size_t n = 0; n < nDomains(); n++) {
868 Refiner& r = domain(n).refiner();
869 r.setGridMin(gridmin);
870 }
871 }
872}
873
874void Sim1D::setMaxGridPoints(int dom, int npoints)
875{
876 if (dom >= 0) {
877 Refiner& r = domain(dom).refiner();
878 r.setMaxPoints(npoints);
879 } else {
880 for (size_t n = 0; n < nDomains(); n++) {
881 Refiner& r = domain(n).refiner();
882 r.setMaxPoints(npoints);
883 }
884 }
885}
886
887size_t Sim1D::maxGridPoints(size_t dom)
888{
889 Refiner& r = domain(dom).refiner();
890 return r.maxPoints();
891}
892
893double Sim1D::jacobian(int i, int j)
894{
895 warn_deprecated("Sim1D::jacobian", "To be removed after Cantera 3.2.");
896 return OneDim::jacobian().value(i,j);
897}
898
900{
901 OneDim::evalSSJacobian(m_state->data(), m_xnew.data());
902}
903
904void Sim1D::solveAdjoint(const double* b, double* lambda)
905{
906 for (auto& D : m_dom) {
907 D->forceFullUpdate(true);
908 }
910 for (auto& D : m_dom) {
911 D->forceFullUpdate(false);
912 }
913
914 auto multijac = dynamic_pointer_cast<MultiJac>(m_jac);
915 if (!multijac) {
916 throw CanteraError("Sim1D::solveAdjoint", "Banded (MultiJac) required");
917 }
918 // Form J^T
919 size_t bw = bandwidth();
920 BandMatrix Jt(size(), bw, bw);
921 for (size_t i = 0; i < size(); i++) {
922 size_t j1 = (i > bw) ? i - bw : 0;
923 size_t j2 = (i + bw >= size()) ? size() - 1: i + bw;
924 for (size_t j = j1; j <= j2; j++) {
925 Jt(j,i) = multijac->value(i,j);
926 }
927 }
928
929 Jt.solve(b, lambda);
930}
931
933{
935 m_xnew.resize(size(), 0.0);
936}
937
938}
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:432
double getDouble(const string &key, double default_) const
If key exists, return it as a double, otherwise return default_.
Definition AnyMap.cpp:1580
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
Definition AnyMap.cpp:1477
static void clearCachedFile(const string &filename)
Remove the specified file from the input cache if it is present.
Definition AnyMap.cpp:1816
static AnyMap fromYamlFile(const string &name, const string &parent_name="")
Create an AnyMap from a YAML file.
Definition AnyMap.cpp:1841
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:29
size_t nComponents() const
Number of components at each grid point.
Definition Domain1D.h:148
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:537
vector< double > & grid()
Access the array of grid coordinates [m].
Definition Domain1D.h:505
double zmin() const
Get the coordinate [m] of the first (leftmost) grid point in this domain.
Definition Domain1D.h:489
size_t nPoints() const
Number of grid points in this domain.
Definition Domain1D.h:170
virtual string domainType() const
Domain type flag.
Definition Domain1D.h:45
double z(size_t jlocal) const
Get the coordinate [m] of the point with local index jlocal
Definition Domain1D.h:484
Refiner & refiner()
Return a reference to the grid refiner.
Definition Domain1D.h:143
virtual void show(const double *x)
Print the solution.
Definition Domain1D.cpp:235
virtual string componentName(size_t n) const
Name of component n. May be overloaded.
Definition Domain1D.cpp:67
string type() const
String indicating the domain implemented.
Definition Domain1D.h:50
double zmax() const
Get the coordinate [m] of the last (rightmost) grid point in this domain.
Definition Domain1D.h:494
virtual void _getInitialSoln(double *x)
Writes some or all initial solution values into the global solution array, beginning at the location ...
Definition Domain1D.cpp:284
size_t index(size_t n, size_t j) const
Returns the index of the solution vector, which corresponds to component n at grid point j.
Definition Domain1D.h:335
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:410
virtual void setupGrid(size_t n, const double *z)
called to set up initial grid, and after grid refinement
Definition Domain1D.cpp:225
This class represents 1D flow domains that satisfy the one-dimensional similarity solution for chemic...
Definition Flow1D.h:46
void setLeftControlPointTemperature(double temperature)
Sets the temperature of the left control point.
Definition Flow1D.cpp:1178
void setLeftControlPointCoordinate(double z_left)
Sets the coordinate of the left control point.
Definition Flow1D.cpp:1193
void setRightControlPointCoordinate(double z_right)
Sets the coordinate of the right control point.
Definition Flow1D.cpp:1248
bool twoPointControlEnabled() const
Returns the status of the two-point control.
Definition Flow1D.h:328
void setRightControlPointTemperature(double temperature)
Sets the temperature of the right control point.
Definition Flow1D.cpp:1233
double m_tfixed
Temperature at the point used to fix the flame location.
Definition Flow1D.h:999
double m_zfixed
Location of the point where temperature is fixed.
Definition Flow1D.h:996
bool isFree() const
Retrieve flag indicating whether flow is freely propagating.
Definition Flow1D.h:360
virtual double eval(double t) const
Evaluate the function.
Definition Func1.cpp:28
void setOptions(int maxJacAge=5)
Set options.
Container class for multiple-domain 1D problems.
Definition OneDim.h:28
int solve(double *x0, double *x1, int loglevel)
Solve F(x) = 0, where F(x) is the multi-domain residual function.
Definition OneDim.cpp:230
size_t start(size_t i) const
The index of the start of domain i in the solution vector.
Definition OneDim.h:107
int m_nsteps
Number of time steps taken in the current call to solve()
Definition OneDim.h:448
virtual void resize()
Call after one or more grids has changed size, for example after being refined.
Definition OneDim.cpp:172
size_t size() const
Total solution vector length;.
Definition OneDim.h:118
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:261
double ssnorm(double *x, double *r)
Steady-state max norm (infinity norm) of the residual evaluated using solution x.
Definition OneDim.cpp:339
size_t nDomains() const
Number of domains.
Definition OneDim.h:76
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:51
size_t bandwidth() const
Jacobian bandwidth.
Definition OneDim.h:147
shared_ptr< SystemJacobian > m_jac
Jacobian evaluator.
Definition OneDim.h:401
shared_ptr< vector< double > > m_state
Solution vector.
Definition OneDim.h:399
vector< int > m_mask
Transient mask. See transientMask().
Definition OneDim.h:433
double timeStep(int nsteps, double dt, double *x, double *r, int loglevel)
Take time steps using Backward Euler.
Definition OneDim.cpp:397
void evalSSJacobian(double *x, double *rsd)
Evaluate the steady-state Jacobian, accessible via jacobian()
Definition OneDim.cpp:240
double m_tmax
maximum timestep size
Definition OneDim.h:394
int m_ss_jac_age
Maximum age of the Jacobian in steady-state mode.
Definition OneDim.h:438
void setSteadyMode()
Prepare to solve the steady-state problem.
Definition OneDim.cpp:368
Domain1D & domain(size_t i) const
Return a reference to domain i.
Definition OneDim.h:81
vector< shared_ptr< Domain1D > > m_dom
All domains comprising the system.
Definition OneDim.h:413
MultiNewton & newton()
Return a reference to the Newton iterator.
Definition OneDim.cpp:100
Refine Domain1D grids so that profiles satisfy adaptation tolerances.
Definition refine.h:17
bool newPointNeeded(size_t j)
Returns true if a new grid point is needed to the right of grid index j.
Definition refine.h:130
size_t maxPoints() const
Returns the maximum number of points allowed in the domain.
Definition refine.h:81
int nNewPoints()
Returns the number of new grid points that were needed.
Definition refine.h:109
void show()
Displays the results of the grid refinement analysis.
Definition refine.cpp:216
vector< double > getCriteria()
Get the grid refinement criteria.
Definition refine.h:60
void setMaxPoints(int npmax)
Set the maximum number of points allowed in the domain.
Definition refine.h:76
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
int analyze(size_t n, const double *z, const double *x)
Determine locations in the grid that need additional grid points and update the internal state of the...
Definition refine.cpp:43
bool keepPoint(size_t j)
Returns true if the grid point at index j should be kept.
Definition refine.h:139
void setGridMin(double gridmin)
Set the minimum allowable spacing between adjacent grid points [m].
Definition refine.h:86
void getInitialSoln()
Get the initial value of the system state from each domain in the simulation.
Definition Sim1D.cpp:364
void restoreTimeSteppingSolution()
Set the current solution vector to the last successful time-stepping solution.
Definition Sim1D.cpp:342
void resize() override
Call after one or more grids has changed size, for example after being refined.
Definition Sim1D.cpp:932
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:377
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:720
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:374
void finalize()
Calls method _finalize in each domain.
Definition Sim1D.cpp:371
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
void writeDebugInfo(const string &header_suffix, const string &message, int loglevel, int attempt_counter)
Write solver debugging information to a YAML file based on the specified log level.
Definition Sim1D.cpp:603
int refine(int loglevel=0)
Refine the grid in all domains.
Definition Sim1D.cpp:511
void show()
Show logging information on current solution for all domains.
Definition Sim1D.cpp:331
double fixedTemperature()
Return temperature at the point used to fix the flame location.
Definition Sim1D.cpp:707
vector< double > m_xlast_ss
the solution vector after the last successful steady-state solve (stored before grid refinement)
Definition Sim1D.h:370
void setMaxGridPoints(int dom, int npoints)
Set the maximum number of grid points in the domain.
Definition Sim1D.cpp:874
int setFixedTemperature(double t)
Add node for fixed temperature point of freely propagating flame.
Definition Sim1D.cpp:617
void clearDebugFile()
Deletes a debug_sim1d.yaml file if it exists.
Definition Sim1D.cpp:598
void setInitialGuess(const string &component, vector< double > &locs, vector< double > &vals)
Set initial guess for one component for all domains.
Definition Sim1D.cpp:41
void solve(int loglevel=0, bool refine_grid=true)
Performs the hybrid Newton steady/time-stepping solution.
Definition Sim1D.cpp:401
int newtonSolve(int loglevel)
Wrapper around the Newton solver.
Definition Sim1D.cpp:387
vector< int > m_steps
array of number of steps to take before re-attempting the steady-state solution
Definition Sim1D.h:384
void evalSSJacobian()
Evaluate the Jacobian in steady-state mode.
Definition Sim1D.cpp:899
double m_tstep
timestep
Definition Sim1D.h:380
void solveAdjoint(const double *b, double *lambda)
Solve the equation .
Definition Sim1D.cpp:904
AnyMap restore(const string &fname, const string &name)
Retrieve data and settings from a previously saved simulation.
Definition Sim1D.cpp:246
Func1 * m_steady_callback
User-supplied function called after a successful steady-state solve.
Definition Sim1D.h:387
void restoreSteadySolution()
Set the current solution vector and grid to the last successful steady- state solution.
Definition Sim1D.cpp:351
size_t maxGridPoints(size_t dom)
Get the maximum number of grid points in this domain.
Definition Sim1D.cpp:887
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:323
void setTimeStep(double stepsize, size_t n, const int *tsteps)
Set the number of time steps to try when the steady Newton solver is unsuccessful.
Definition Sim1D.cpp:378
void setRightControlPoint(double temperature)
Set the right control point location using the specified temperature.
Definition Sim1D.cpp:784
double value(size_t dom, size_t comp, size_t localPoint) const
Get one entry in the solution vector.
Definition Sim1D.cpp:63
double workValue(size_t dom, size_t comp, size_t localPoint) const
Get an entry in the work vector, which may contain either a new system state or the current residual ...
Definition Sim1D.cpp:71
vector< double > getRefineCriteria(int dom)
Get the grid refinement criteria.
Definition Sim1D.cpp:850
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:861
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:836
vector< double > m_xlast_ts
the solution vector after the last successful timestepping
Definition Sim1D.h:366
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
void setLeftControlPoint(double temperature)
Set the left control point location using the specified temperature.
Definition Sim1D.cpp:733
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.
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:88
#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:154
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition global.h:171
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:263
Namespace for the Cantera kernel.
Definition AnyMap.cpp:595
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:180
@ c_offset_T
temperature [kelvin]
Definition Flow1D.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:1997
Contains declarations for string manipulation functions within Cantera.