Cantera  3.1.0a1
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 
8 #include "cantera/oneD/Sim1D.h"
10 #include "cantera/oneD/StFlow.h"
12 #include "cantera/oneD/refine.h"
13 #include "cantera/numerics/funcs.h"
16 #include "cantera/numerics/Func1.h"
17 #include <limits>
18 #include <fstream>
19 
20 using namespace std;
21 
22 namespace Cantera
23 {
24 
25 Sim1D::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 
41 void 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 
55 void 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 
63 double 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 
71 double 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 
79 void 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 
98 void 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 
147 void 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 
160 namespace { // restrict scope of helper function to local translation unit
161 
162 //! convert data format used by Python h5py export (Cantera < 3.0)
163 AnyMap 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 
247 AnyMap 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 
324 void 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 
332 void 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  }
358  *m_state = m_xlast_ts;
359 }
360 
362 {
363  if (m_xlast_ss.empty()) {
364  throw CanteraError("Sim1D::restoreSteadySolution",
365  "No successful steady state solution");
366  }
367  *m_state = m_xlast_ss;
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 
374 void 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 
388 void 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 
397 int 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 
411 void 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
428  setSteadyMode();
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", "debug",
451  "After successful Newton solve");
452  }
453  if (loglevel > 7) {
454  saveResidual("debug_sim1d.yaml", "residual",
455  "After successful Newton solve");
456  }
457  ok = true;
458  } else {
459  debuglog(" failure. \n", loglevel);
460  if (loglevel > 6) {
461  save("debug_sim1d.yaml", "debug",
462  "After unsuccessful Newton solve");
463  }
464  if (loglevel > 7) {
465  saveResidual("debug_sim1d.yaml", "residual",
466  "After unsuccessful Newton solve");
467  }
468  if (loglevel > 0) {
469  writelog("Take {} timesteps ", nsteps);
470  }
471  dt = timeStep(nsteps, dt, m_state->data(), m_xnew.data(), loglevel-1);
472  m_xlast_ts = *m_state;
473  if (loglevel > 6) {
474  save("debug_sim1d.yaml", "debug", "After timestepping");
475  }
476  if (loglevel > 7) {
477  saveResidual("debug_sim1d.yaml", "residual",
478  "After timestepping");
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", "debug", "After regridding");
510  }
511  if (new_points && loglevel > 7) {
512  saveResidual("debug_sim1d.yaml", "residual",
513  "After regridding");
514  }
515  } else {
516  debuglog("grid refinement disabled.\n", loglevel);
517  new_points = 0;
518  }
519  }
520 }
521 
522 int Sim1D::refine(int loglevel)
523 {
524  int ianalyze, np = 0;
525  vector<double> znew, xnew;
526  vector<size_t> dsize;
527 
528  m_xlast_ss = *m_state;
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 
727 void 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 
741 vector<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 
752 void 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 
765 void 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 
778 size_t Sim1D::maxGridPoints(size_t dom)
779 {
780  Refiner& r = domain(dom).refiner();
781  return r.maxPoints();
782 }
783 
784 double Sim1D::jacobian(int i, int j)
785 {
786  return OneDim::jacobian().value(i,j);
787 }
788 
789 void Sim1D::evalSSJacobian()
790 {
791  OneDim::evalSSJacobian(m_state->data(), m_xnew.data());
792 }
793 
794 void 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 {
820  OneDim::resize();
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.
Definition: BandMatrix.cpp:253
double & value(size_t i, size_t j)
Return a changeable reference to element (i,j).
Definition: BandMatrix.cpp:148
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:66
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
virtual string componentName(size_t n) const
Name of the nth component. May be overloaded.
Definition: Domain1D.cpp:49
Refiner & refiner()
Return a reference to the grid refiner.
Definition: Domain1D.h:140
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
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
vector< double > getCriteria()
Get the grid refinement criteria.
Definition: refine.h:42
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
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.
Definition: SolutionArray.h:51
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.
Definition: ctexceptions.h:278
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.