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