1 /** 2 * @file ctonedim.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 #define CANTERA_USE_INTERNAL 9 10 #include "cantera/clib/ctonedim.h" 11 12 // Cantera includes 13 #include "cantera/oneD/Sim1D.h" 14 #include "cantera/oneD/Boundary1D.h" 15 #include "cantera/transport/TransportBase.h" 16 #include "Cabinet.h" 17 #include "cantera/base/stringUtils.h" 18 19 #include <fstream> 20 21 using namespace std; 22 using namespace Cantera; 23 24 typedef Cabinet<Sim1D> SimCabinet; 25 typedef Cabinet<Domain1D> DomainCabinet; 26 template<> SimCabinet* SimCabinet::s_storage = 0; 27 template<> DomainCabinet* DomainCabinet::s_storage = 0; 28 29 typedef Cabinet<ThermoPhase> ThermoCabinet; 30 typedef Cabinet<Kinetics> KineticsCabinet; 31 typedef Cabinet<Transport> TransportCabinet; 32 33 extern "C" { 34 ct_clearOneDim()35 int ct_clearOneDim() 36 { 37 try { 38 DomainCabinet::clear(); 39 SimCabinet::clear(); 40 return 0; 41 } catch (...) { 42 return handleAllExceptions(-1, ERR); 43 } 44 } 45 domain_del(int i)46 int domain_del(int i) 47 { 48 try { 49 DomainCabinet::del(i); 50 return 0; 51 } catch (...) { 52 return handleAllExceptions(-1 , ERR); 53 } 54 } 55 domain_type(int i)56 int domain_type(int i) 57 { 58 try { 59 return DomainCabinet::item(i).domainType(); 60 } catch (...) { 61 return handleAllExceptions(-1, ERR); 62 } 63 } 64 domain_index(int i)65 size_t domain_index(int i) 66 { 67 try { 68 return DomainCabinet::item(i).domainIndex(); 69 } catch (...) { 70 return handleAllExceptions(npos, npos); 71 } 72 } 73 domain_nComponents(int i)74 size_t domain_nComponents(int i) 75 { 76 try { 77 return DomainCabinet::item(i).nComponents(); 78 } catch (...) { 79 return handleAllExceptions(npos, npos); 80 } 81 } 82 domain_nPoints(int i)83 size_t domain_nPoints(int i) 84 { 85 try { 86 return DomainCabinet::item(i).nPoints(); 87 } catch (...) { 88 return handleAllExceptions(npos, npos); 89 } 90 } 91 domain_componentName(int i,int n,int sz,char * buf)92 int domain_componentName(int i, int n, int sz, char* buf) 93 { 94 try { 95 Domain1D& dom = DomainCabinet::item(i); 96 dom.checkComponentIndex(n); 97 return static_cast<int>(copyString(dom.componentName(n), buf, sz)); 98 } catch (...) { 99 return handleAllExceptions(-1, ERR); 100 } 101 } 102 domain_componentIndex(int i,const char * name)103 size_t domain_componentIndex(int i, const char* name) 104 { 105 try { 106 return DomainCabinet::item(i).componentIndex(name); 107 } catch (...) { 108 return handleAllExceptions(-1, ERR); 109 } 110 } 111 domain_grid(int i,int n)112 double domain_grid(int i, int n) 113 { 114 try { 115 Domain1D& dom = DomainCabinet::item(i); 116 dom.checkPointIndex(n); 117 return dom.grid(n); 118 } catch (...) { 119 return handleAllExceptions(DERR, DERR); 120 } 121 } 122 domain_setBounds(int i,int n,double lower,double upper)123 int domain_setBounds(int i, int n, double lower, double upper) 124 { 125 try { 126 Domain1D& dom = DomainCabinet::item(i); 127 dom.checkComponentIndex(n); 128 dom.setBounds(n, lower, upper); 129 return 0; 130 } catch (...) { 131 return handleAllExceptions(-1, ERR); 132 } 133 } 134 domain_upperBound(int i,int n)135 double domain_upperBound(int i, int n) 136 { 137 try { 138 Domain1D& dom = DomainCabinet::item(i); 139 dom.checkComponentIndex(n); 140 return dom.upperBound(n); 141 } catch (...) { 142 return handleAllExceptions(DERR, DERR); 143 } 144 } 145 domain_lowerBound(int i,int n)146 double domain_lowerBound(int i, int n) 147 { 148 try { 149 Domain1D& dom = DomainCabinet::item(i); 150 dom.checkComponentIndex(n); 151 return dom.lowerBound(n); 152 } catch (...) { 153 return handleAllExceptions(DERR, DERR); 154 } 155 } 156 domain_setSteadyTolerances(int i,int n,double rtol,double atol)157 int domain_setSteadyTolerances(int i, int n, double rtol, 158 double atol) 159 { 160 try { 161 Domain1D& dom = DomainCabinet::item(i); 162 dom.checkComponentIndex(n); 163 dom.setSteadyTolerances(rtol, atol, n); 164 return 0; 165 } catch (...) { 166 return handleAllExceptions(-1, ERR); 167 } 168 } 169 domain_setTransientTolerances(int i,int n,double rtol,double atol)170 int domain_setTransientTolerances(int i, int n, double rtol, 171 double atol) 172 { 173 try { 174 Domain1D& dom = DomainCabinet::item(i); 175 dom.checkComponentIndex(n); 176 dom.setTransientTolerances(rtol, atol, n); 177 return 0; 178 } catch (...) { 179 return handleAllExceptions(-1, ERR); 180 } 181 } 182 domain_rtol(int i,int n)183 double domain_rtol(int i, int n) 184 { 185 try { 186 Domain1D& dom = DomainCabinet::item(i); 187 dom.checkComponentIndex(n); 188 return dom.rtol(n); 189 } catch (...) { 190 return handleAllExceptions(DERR, DERR); 191 } 192 } 193 domain_atol(int i,int n)194 double domain_atol(int i, int n) 195 { 196 try { 197 Domain1D& dom = DomainCabinet::item(i); 198 dom.checkComponentIndex(n); 199 return dom.atol(n); 200 } catch (...) { 201 return handleAllExceptions(DERR, DERR); 202 } 203 } 204 domain_setupGrid(int i,size_t npts,const double * grid)205 int domain_setupGrid(int i, size_t npts, const double* grid) 206 { 207 try { 208 DomainCabinet::item(i).setupGrid(npts, grid); 209 return 0; 210 } catch (...) { 211 return handleAllExceptions(-1, ERR); 212 } 213 } 214 domain_setID(int i,const char * id)215 int domain_setID(int i, const char* id) 216 { 217 try { 218 DomainCabinet::item(i).setID(id); 219 return 0; 220 } catch (...) { 221 return handleAllExceptions(-1, ERR); 222 } 223 } 224 inlet_new()225 int inlet_new() 226 { 227 try { 228 return DomainCabinet::add(new Inlet1D()); 229 } catch (...) { 230 return handleAllExceptions(-1, ERR); 231 } 232 } 233 surf_new()234 int surf_new() 235 { 236 try { 237 return DomainCabinet::add(new Surf1D()); 238 } catch (...) { 239 return handleAllExceptions(-1, ERR); 240 } 241 } 242 reactingsurf_new()243 int reactingsurf_new() 244 { 245 try { 246 return DomainCabinet::add(new ReactingSurf1D()); 247 } catch (...) { 248 return handleAllExceptions(-1, ERR); 249 } 250 } 251 symm_new()252 int symm_new() 253 { 254 try { 255 return DomainCabinet::add(new Symm1D()); 256 } catch (...) { 257 return handleAllExceptions(-1, ERR); 258 } 259 } 260 outlet_new()261 int outlet_new() 262 { 263 try { 264 return DomainCabinet::add(new Outlet1D()); 265 } catch (...) { 266 return handleAllExceptions(-1, ERR); 267 } 268 } 269 outletres_new()270 int outletres_new() 271 { 272 try { 273 return DomainCabinet::add(new OutletRes1D()); 274 } catch (...) { 275 return handleAllExceptions(-1, ERR); 276 } 277 } 278 bdry_setMdot(int i,double mdot)279 int bdry_setMdot(int i, double mdot) 280 { 281 try { 282 DomainCabinet::get<Boundary1D>(i).setMdot(mdot); 283 return 0; 284 } catch (...) { 285 return handleAllExceptions(-1, ERR); 286 } 287 } 288 bdry_setTemperature(int i,double t)289 int bdry_setTemperature(int i, double t) 290 { 291 try { 292 DomainCabinet::get<Boundary1D>(i).setTemperature(t); 293 return 0; 294 } catch (...) { 295 return handleAllExceptions(-1, ERR); 296 } 297 } 298 bdry_setMoleFractions(int i,const char * x)299 int bdry_setMoleFractions(int i, const char* x) 300 { 301 try { 302 DomainCabinet::get<Boundary1D>(i).setMoleFractions(x); 303 return 0; 304 } catch (...) { 305 return handleAllExceptions(-1, ERR); 306 } 307 } 308 bdry_temperature(int i)309 double bdry_temperature(int i) 310 { 311 try { 312 return DomainCabinet::get<Boundary1D>(i).temperature(); 313 } catch (...) { 314 return handleAllExceptions(DERR, DERR); 315 } 316 } 317 bdry_massFraction(int i,int k)318 double bdry_massFraction(int i, int k) 319 { 320 try { 321 return DomainCabinet::get<Boundary1D>(i).massFraction(k); 322 } catch (...) { 323 return handleAllExceptions(DERR, DERR); 324 } 325 } 326 bdry_mdot(int i)327 double bdry_mdot(int i) 328 { 329 try { 330 return DomainCabinet::get<Boundary1D>(i).mdot(); 331 } catch (...) { 332 return handleAllExceptions(DERR, DERR); 333 } 334 } 335 reactingsurf_setkineticsmgr(int i,int j)336 int reactingsurf_setkineticsmgr(int i, int j) 337 { 338 try { 339 InterfaceKinetics& k = Cabinet<Kinetics>::get<InterfaceKinetics>(j); 340 DomainCabinet::get<ReactingSurf1D>(i).setKineticsMgr(&k); 341 return 0; 342 } catch (...) { 343 return handleAllExceptions(-1, ERR); 344 } 345 } 346 reactingsurf_enableCoverageEqs(int i,int onoff)347 int reactingsurf_enableCoverageEqs(int i, int onoff) 348 { 349 try { 350 DomainCabinet::get<ReactingSurf1D>(i).enableCoverageEquations(onoff != 0); 351 return 0; 352 } catch (...) { 353 return handleAllExceptions(-1, ERR); 354 } 355 } 356 inlet_setSpreadRate(int i,double v)357 int inlet_setSpreadRate(int i, double v) 358 { 359 try { 360 DomainCabinet::get<Inlet1D>(i).setSpreadRate(v); 361 return 0; 362 } catch (...) { 363 return handleAllExceptions(-1, ERR); 364 } 365 } 366 367 //------------------ stagnation flow domains -------------------- 368 stflow_new(int iph,int ikin,int itr,int itype)369 int stflow_new(int iph, int ikin, int itr, int itype) 370 { 371 try { 372 IdealGasPhase& ph = ThermoCabinet::get<IdealGasPhase>(iph); 373 StFlow* x = new StFlow(&ph, ph.nSpecies(), 2); 374 if (itype == 1) { 375 x->setAxisymmetricFlow(); 376 } else if (itype == 2) { 377 x->setFreeFlow(); 378 } else { 379 delete x; 380 return -2; 381 } 382 x->setKinetics(KineticsCabinet::item(ikin)); 383 x->setTransport(TransportCabinet::item(itr)); 384 return DomainCabinet::add(x); 385 } catch (...) { 386 return handleAllExceptions(-1, ERR); 387 } 388 } 389 390 stflow_setTransport(int i,int itr)391 int stflow_setTransport(int i, int itr) 392 { 393 try { 394 DomainCabinet::get<StFlow>(i).setTransport(TransportCabinet::item(itr)); 395 return 0; 396 } catch (...) { 397 return handleAllExceptions(-1, ERR); 398 } 399 } 400 stflow_enableSoret(int i,int iSoret)401 int stflow_enableSoret(int i, int iSoret) 402 { 403 try { 404 bool withSoret = (iSoret > 0); 405 DomainCabinet::get<StFlow>(i).enableSoret(withSoret); 406 return 0; 407 } catch (...) { 408 return handleAllExceptions(-1, ERR); 409 } 410 } 411 stflow_setPressure(int i,double p)412 int stflow_setPressure(int i, double p) 413 { 414 try { 415 DomainCabinet::get<StFlow>(i).setPressure(p); 416 return 0; 417 } catch (...) { 418 return handleAllExceptions(-1, ERR); 419 } 420 } 421 stflow_pressure(int i)422 double stflow_pressure(int i) 423 { 424 try { 425 return DomainCabinet::get<StFlow>(i).pressure(); 426 } catch (...) { 427 return handleAllExceptions(DERR, DERR); 428 } 429 } 430 stflow_setFixedTempProfile(int i,size_t n,const double * pos,size_t m,const double * temp)431 int stflow_setFixedTempProfile(int i, size_t n, const double* pos, 432 size_t m, const double* temp) 433 { 434 try { 435 vector_fp vpos(n), vtemp(n); 436 for (size_t j = 0; j < n; j++) { 437 vpos[j] = pos[j]; 438 vtemp[j] = temp[j]; 439 } 440 DomainCabinet::get<StFlow>(i).setFixedTempProfile(vpos, vtemp); 441 return 0; 442 } catch (...) { 443 return handleAllExceptions(-1, ERR); 444 } 445 } 446 stflow_solveEnergyEqn(int i,int flag)447 int stflow_solveEnergyEqn(int i, int flag) 448 { 449 try { 450 if (flag > 0) { 451 DomainCabinet::get<StFlow>(i).solveEnergyEqn(npos); 452 } else { 453 DomainCabinet::get<StFlow>(i).fixTemperature(npos); 454 } 455 return 0; 456 } catch (...) { 457 return handleAllExceptions(-1, ERR); 458 } 459 } 460 461 //------------------- Sim1D -------------------------------------- 462 sim1D_new(size_t nd,const int * domains)463 int sim1D_new(size_t nd, const int* domains) 464 { 465 try { 466 vector<Domain1D*> d; 467 for (size_t n = 0; n < nd; n++) { 468 d.push_back(&DomainCabinet::item(domains[n])); 469 } 470 return SimCabinet::add(new Sim1D(d)); 471 } catch (...) { 472 return handleAllExceptions(-1, ERR); 473 } 474 } 475 sim1D_del(int i)476 int sim1D_del(int i) 477 { 478 try { 479 SimCabinet::del(i); 480 return 0; 481 } catch (...) { 482 return handleAllExceptions(-1, ERR); 483 } 484 } 485 sim1D_setValue(int i,int dom,int comp,int localPoint,double value)486 int sim1D_setValue(int i, int dom, int comp, int localPoint, double value) 487 { 488 try { 489 SimCabinet::item(i).setValue(dom, comp, localPoint, value); 490 return 0; 491 } catch (...) { 492 return handleAllExceptions(-1, ERR); 493 } 494 } 495 sim1D_setProfile(int i,int dom,int comp,size_t np,const double * pos,size_t nv,const double * v)496 int sim1D_setProfile(int i, int dom, int comp, size_t np, const double* pos, 497 size_t nv, const double* v) 498 { 499 try { 500 Sim1D& sim = SimCabinet::item(i); 501 sim.checkDomainIndex(dom); 502 sim.domain(dom).checkComponentIndex(comp); 503 vector_fp vv, pv; 504 for (size_t n = 0; n < np; n++) { 505 vv.push_back(v[n]); 506 pv.push_back(pos[n]); 507 } 508 sim.setProfile(dom, comp, pv, vv); 509 return 0; 510 } catch (...) { 511 return handleAllExceptions(-1, ERR); 512 } 513 } 514 sim1D_setFlatProfile(int i,int dom,int comp,double v)515 int sim1D_setFlatProfile(int i, int dom, int comp, double v) 516 { 517 try { 518 Sim1D& sim = SimCabinet::item(i); 519 sim.checkDomainIndex(dom); 520 sim.domain(dom).checkComponentIndex(comp); 521 sim.setFlatProfile(dom, comp, v); 522 return 0; 523 } catch (...) { 524 return handleAllExceptions(-1, ERR); 525 } 526 } 527 sim1D_showSolution(int i,const char * fname)528 int sim1D_showSolution(int i, const char* fname) 529 { 530 try { 531 string fn = string(fname); 532 if (fn == "-") { 533 SimCabinet::item(i).showSolution(); 534 } else { 535 ofstream fout(fname); 536 SimCabinet::item(i).showSolution(fout); 537 } 538 return 0; 539 } catch (...) { 540 return handleAllExceptions(-1, ERR); 541 } 542 } 543 sim1D_setTimeStep(int i,double stepsize,size_t ns,const int * nsteps)544 int sim1D_setTimeStep(int i, double stepsize, size_t ns, const int* nsteps) 545 { 546 try { 547 SimCabinet::item(i).setTimeStep(stepsize, ns, nsteps); 548 return 0; 549 } catch (...) { 550 return handleAllExceptions(-1, ERR); 551 } 552 } 553 sim1D_getInitialSoln(int i)554 int sim1D_getInitialSoln(int i) 555 { 556 try { 557 SimCabinet::item(i).getInitialSoln(); 558 return 0; 559 } catch (...) { 560 return handleAllExceptions(-1, ERR); 561 } 562 } 563 sim1D_solve(int i,int loglevel,int refine_grid)564 int sim1D_solve(int i, int loglevel, int refine_grid) 565 { 566 try { 567 bool r = (refine_grid == 0 ? false : true); 568 SimCabinet::item(i).solve(loglevel, r); 569 return 0; 570 } catch (...) { 571 return handleAllExceptions(-1, ERR); 572 } 573 } 574 sim1D_refine(int i,int loglevel)575 int sim1D_refine(int i, int loglevel) 576 { 577 try { 578 SimCabinet::item(i).refine(loglevel); 579 return 0; 580 } catch (...) { 581 return handleAllExceptions(-1, ERR); 582 } 583 } 584 sim1D_setRefineCriteria(int i,int dom,double ratio,double slope,double curve,double prune)585 int sim1D_setRefineCriteria(int i, int dom, double ratio, 586 double slope, double curve, double prune) 587 { 588 try { 589 SimCabinet::item(i).setRefineCriteria(dom, ratio, slope, curve, prune); 590 return 0; 591 } catch (...) { 592 return handleAllExceptions(-1, ERR); 593 } 594 } 595 sim1D_setGridMin(int i,int dom,double gridmin)596 int sim1D_setGridMin(int i, int dom, double gridmin) 597 { 598 try { 599 SimCabinet::item(i).setGridMin(dom, gridmin); 600 return 0; 601 } catch (...) { 602 return handleAllExceptions(-1, ERR); 603 } 604 } 605 sim1D_save(int i,const char * fname,const char * id,const char * desc)606 int sim1D_save(int i, const char* fname, const char* id, const char* desc) 607 { 608 try { 609 SimCabinet::item(i).save(fname, id, desc); 610 return 0; 611 } catch (...) { 612 return handleAllExceptions(-1, ERR); 613 } 614 } 615 sim1D_restore(int i,const char * fname,const char * id)616 int sim1D_restore(int i, const char* fname, const char* id) 617 { 618 try { 619 SimCabinet::item(i).restore(fname, id); 620 return 0; 621 } catch (...) { 622 return handleAllExceptions(-1, ERR); 623 } 624 } 625 sim1D_writeStats(int i,int printTime)626 int sim1D_writeStats(int i, int printTime) 627 { 628 try { 629 SimCabinet::item(i).writeStats(printTime); 630 return 0; 631 } catch (...) { 632 return handleAllExceptions(-1, ERR); 633 } 634 } 635 sim1D_domainIndex(int i,const char * name)636 int sim1D_domainIndex(int i, const char* name) 637 { 638 try { 639 return (int) SimCabinet::item(i).domainIndex(name); 640 } catch (...) { 641 return handleAllExceptions(-1, ERR); 642 } 643 } 644 sim1D_value(int i,int idom,int icomp,int localPoint)645 double sim1D_value(int i, int idom, int icomp, int localPoint) 646 { 647 try { 648 Sim1D& sim = SimCabinet::item(i); 649 sim.checkDomainIndex(idom); 650 sim.domain(idom).checkComponentIndex(icomp); 651 return sim.value(idom, icomp, localPoint); 652 } catch (...) { 653 return handleAllExceptions(DERR, DERR); 654 } 655 } 656 sim1D_workValue(int i,int idom,int icomp,int localPoint)657 double sim1D_workValue(int i, int idom, int icomp, int localPoint) 658 { 659 try { 660 Sim1D& sim = SimCabinet::item(i); 661 sim.checkDomainIndex(idom); 662 sim.domain(idom).checkComponentIndex(icomp); 663 return sim.workValue(idom, icomp, localPoint); 664 } catch (...) { 665 return handleAllExceptions(DERR, DERR); 666 } 667 } 668 sim1D_eval(int i,double rdt,int count)669 int sim1D_eval(int i, double rdt, int count) 670 { 671 try { 672 SimCabinet::item(i).eval(rdt, count); 673 return 0; 674 } catch (...) { 675 return handleAllExceptions(-1, ERR); 676 } 677 } 678 sim1D_setMaxJacAge(int i,int ss_age,int ts_age)679 int sim1D_setMaxJacAge(int i, int ss_age, int ts_age) 680 { 681 try { 682 SimCabinet::item(i).setJacAge(ss_age, ts_age); 683 return 0; 684 } catch (...) { 685 return handleAllExceptions(-1, ERR); 686 } 687 } 688 sim1D_setFixedTemperature(int i,double temp)689 int sim1D_setFixedTemperature(int i, double temp) 690 { 691 try { 692 SimCabinet::item(i).setFixedTemperature(temp); 693 return 0; 694 } catch (...) { 695 return handleAllExceptions(-1, ERR); 696 } 697 } 698 } 699