1 //! @file Domain1D.h 2 3 // This file is part of Cantera. See License.txt in the top-level directory or 4 // at https://cantera.org/license.txt for license and copyright information. 5 6 #ifndef CT_DOMAIN1D_H 7 #define CT_DOMAIN1D_H 8 9 #include "cantera/base/ctexceptions.h" 10 11 namespace Cantera 12 { 13 14 // domain types 15 const int cFlowType = 50; 16 const int cFreeFlow = 51; 17 const int cAxisymmetricStagnationFlow = 52; 18 const int cConnectorType = 100; 19 const int cSurfType = 102; 20 const int cInletType = 104; 21 const int cSymmType = 105; 22 const int cOutletType = 106; 23 const int cEmptyType = 107; 24 const int cOutletResType = 108; 25 const int cPorousType = 109; 26 27 class MultiJac; 28 class OneDim; 29 class Refiner; 30 class XML_Node; 31 32 /** 33 * Base class for one-dimensional domains. 34 * @ingroup onedim 35 */ 36 class Domain1D 37 { 38 public: 39 /** 40 * Constructor. 41 * @param nv Number of variables at each grid point. 42 * @param points Number of grid points. 43 * @param time (unused) 44 */ 45 Domain1D(size_t nv=1, size_t points=1, double time=0.0); 46 47 virtual ~Domain1D(); 48 Domain1D(const Domain1D&) = delete; 49 Domain1D& operator=(const Domain1D&) = delete; 50 51 //! Domain type flag. domainType()52 int domainType() { 53 return m_type; 54 } 55 56 //! The left-to-right location of this domain. domainIndex()57 size_t domainIndex() { 58 return m_index; 59 } 60 61 //! True if the domain is a connector domain. isConnector()62 bool isConnector() { 63 return (m_type >= cConnectorType); 64 } 65 66 //! The container holding this domain. container()67 const OneDim& container() const { 68 return *m_container; 69 } 70 71 //! Specify the container object for this domain, and the position of this 72 //! domain in the list. setContainer(OneDim * c,size_t index)73 void setContainer(OneDim* c, size_t index) { 74 m_container = c; 75 m_index = index; 76 } 77 78 //! Set the Jacobian bandwidth. See the discussion of method bandwidth(). 79 void setBandwidth(int bw = -1) { 80 m_bw = bw; 81 } 82 83 //! Set the Jacobian bandwidth for this domain. 84 /** 85 * When class OneDim computes the bandwidth of the overall multi-domain 86 * problem (in OneDim::resize()), it calls this method for the bandwidth 87 * of each domain. If setBandwidth has not been called, then a negative 88 * bandwidth is returned, in which case OneDim assumes that this domain is 89 * dense -- that is, at each point, all components depend on the value of 90 * all other components at that point. In this case, the bandwidth is bw = 91 * 2*nComponents() - 1. However, if this domain contains some components 92 * that are uncoupled from other components at the same point, then this 93 * default bandwidth may greatly overestimate the true bandwidth, with a 94 * substantial penalty in performance. For such domains, use method 95 * setBandwidth to specify the bandwidth before passing this domain to the 96 * Sim1D or OneDim constructor. 97 */ bandwidth()98 size_t bandwidth() { 99 return m_bw; 100 } 101 102 /*! 103 * Initialize. This method is called by OneDim::init() for each domain once 104 * at the beginning of a simulation. Base class method does nothing, but may 105 * be overloaded. 106 */ init()107 virtual void init() { } 108 109 virtual void setInitialState(doublereal* xlocal = 0) {} setState(size_t point,const doublereal * state,doublereal * x)110 virtual void setState(size_t point, const doublereal* state, doublereal* x) {} 111 112 /*! 113 * When called, this function should reset "bad" values in the state vector 114 * such as negative species concentrations. This function may be called 115 * after a failed solution attempt. 116 */ resetBadValues(double * xg)117 virtual void resetBadValues(double* xg) {} 118 119 /*! 120 * Resize the domain to have nv components and np grid points. This method 121 * is virtual so that subclasses can perform other actions required to 122 * resize the domain. 123 */ 124 virtual void resize(size_t nv, size_t np); 125 126 //! Return a reference to the grid refiner. refiner()127 Refiner& refiner() { 128 return *m_refiner; 129 } 130 131 //! Number of components at each grid point. nComponents()132 size_t nComponents() const { 133 return m_nv; 134 } 135 136 //! Check that the specified component index is in range. 137 //! Throws an exception if n is greater than nComponents()-1 checkComponentIndex(size_t n)138 void checkComponentIndex(size_t n) const { 139 if (n >= m_nv) { 140 throw IndexError("Domain1D::checkComponentIndex", "points", n, m_nv-1); 141 } 142 } 143 144 //! Check that an array size is at least nComponents(). 145 //! Throws an exception if nn is less than nComponents(). Used before calls 146 //! which take an array pointer. checkComponentArraySize(size_t nn)147 void checkComponentArraySize(size_t nn) const { 148 if (m_nv > nn) { 149 throw ArraySizeError("Domain1D::checkComponentArraySize", nn, m_nv); 150 } 151 } 152 153 //! Number of grid points in this domain. nPoints()154 size_t nPoints() const { 155 return m_points; 156 } 157 158 //! Check that the specified point index is in range. 159 //! Throws an exception if n is greater than nPoints()-1 checkPointIndex(size_t n)160 void checkPointIndex(size_t n) const { 161 if (n >= m_points) { 162 throw IndexError("Domain1D::checkPointIndex", "points", n, m_points-1); 163 } 164 } 165 166 //! Check that an array size is at least nPoints(). 167 //! Throws an exception if nn is less than nPoints(). Used before calls 168 //! which take an array pointer. checkPointArraySize(size_t nn)169 void checkPointArraySize(size_t nn) const { 170 if (m_points > nn) { 171 throw ArraySizeError("Domain1D::checkPointArraySize", nn, m_points); 172 } 173 } 174 175 //! Name of the nth component. May be overloaded. 176 virtual std::string componentName(size_t n) const; 177 setComponentName(size_t n,const std::string & name)178 void setComponentName(size_t n, const std::string& name) { 179 m_name[n] = name; 180 } 181 182 //! index of component with name \a name. 183 virtual size_t componentIndex(const std::string& name) const; 184 setBounds(size_t n,doublereal lower,doublereal upper)185 void setBounds(size_t n, doublereal lower, doublereal upper) { 186 m_min[n] = lower; 187 m_max[n] = upper; 188 } 189 190 //! Set tolerances for time-stepping mode 191 /*! 192 * @param rtol Relative tolerance 193 * @param atol Absolute tolerance 194 * @param n component index these tolerances apply to. If set to -1 (the 195 * default), these tolerances will be applied to all solution 196 * components. 197 */ 198 void setTransientTolerances(doublereal rtol, doublereal atol, size_t n=npos); 199 200 //! Set tolerances for steady-state mode 201 /*! 202 * @param rtol Relative tolerance 203 * @param atol Absolute tolerance 204 * @param n component index these tolerances apply to. If set to -1 (the 205 * default), these tolerances will be applied to all solution 206 * components. 207 */ 208 void setSteadyTolerances(doublereal rtol, doublereal atol, size_t n=npos); 209 210 //! Relative tolerance of the nth component. rtol(size_t n)211 doublereal rtol(size_t n) { 212 return (m_rdt == 0.0 ? m_rtol_ss[n] : m_rtol_ts[n]); 213 } 214 215 //! Absolute tolerance of the nth component. atol(size_t n)216 doublereal atol(size_t n) { 217 return (m_rdt == 0.0 ? m_atol_ss[n] : m_atol_ts[n]); 218 } 219 220 //! Steady relative tolerance of the nth component steady_rtol(size_t n)221 double steady_rtol(size_t n) { 222 return m_rtol_ss[n]; 223 } 224 225 //! Steady absolute tolerance of the nth component steady_atol(size_t n)226 double steady_atol(size_t n) { 227 return m_atol_ss[n]; 228 } 229 230 //! Transient relative tolerance of the nth component transient_rtol(size_t n)231 double transient_rtol(size_t n) { 232 return m_rtol_ts[n]; 233 } 234 235 //! Transient absolute tolerance of the nth component transient_atol(size_t n)236 double transient_atol(size_t n) { 237 return m_atol_ts[n]; 238 } 239 240 //! Upper bound on the nth component. upperBound(size_t n)241 doublereal upperBound(size_t n) const { 242 return m_max[n]; 243 } 244 245 //! Lower bound on the nth component lowerBound(size_t n)246 doublereal lowerBound(size_t n) const { 247 return m_min[n]; 248 } 249 250 //! Prepare to do time stepping with time step dt 251 /*! 252 * Copy the internally-stored solution at the last time step to array x0. 253 */ initTimeInteg(doublereal dt,const doublereal * x0)254 void initTimeInteg(doublereal dt, const doublereal* x0) { 255 std::copy(x0 + loc(), x0 + loc() + size(), m_slast.begin()); 256 m_rdt = 1.0/dt; 257 } 258 259 //! Prepare to solve the steady-state problem 260 /*! 261 * Set the internally-stored reciprocal of the time step to 0.0 262 */ setSteadyMode()263 void setSteadyMode() { 264 m_rdt = 0.0; 265 } 266 267 //! True if in steady-state mode steady()268 bool steady() { 269 return (m_rdt == 0.0); 270 } 271 272 //! True if not in steady-state mode transient()273 bool transient() { 274 return (m_rdt != 0.0); 275 } 276 277 /*! 278 * Set this if something has changed in the governing 279 * equations (e.g. the value of a constant has been changed, 280 * so that the last-computed Jacobian is no longer valid. 281 */ 282 void needJacUpdate(); 283 284 //! Evaluate the residual function at point j. If j == npos, 285 //! evaluate the residual function at all points. 286 /*! 287 * This function must be implemented in classes derived from Domain1D. 288 * 289 * @param j Grid point at which to update the residual 290 * @param[in] x State vector 291 * @param[out] r residual vector 292 * @param[out] mask Boolean mask indicating whether each solution 293 * component has a time derivative (1) or not (0). 294 * @param[in] rdt Reciprocal of the timestep (`rdt=0` implies steady- 295 * state.) 296 */ 297 virtual void eval(size_t j, doublereal* x, doublereal* r, 298 integer* mask, doublereal rdt=0.0) { 299 throw NotImplementedError("Domain1D::eval"); 300 } 301 index(size_t n,size_t j)302 size_t index(size_t n, size_t j) const { 303 return m_nv*j + n; 304 } value(const doublereal * x,size_t n,size_t j)305 doublereal value(const doublereal* x, size_t n, size_t j) const { 306 return x[index(n,j)]; 307 } 308 setJac(MultiJac * jac)309 virtual void setJac(MultiJac* jac) {} 310 311 //! Save the current solution for this domain into an XML_Node 312 /*! 313 * Base class version of the general domain1D save function. Derived classes 314 * should call the base class method in addition to saving their own data. 315 * 316 * @param o XML_Node to save the solution to. 317 * @param sol Current value of the solution vector. The object will pick 318 * out which part of the solution vector pertains to this 319 * object. 320 * @return XML_Node created to represent this domain 321 * 322 * @deprecated The XML output format is deprecated and will be removed in 323 * Cantera 3.0. 324 */ 325 virtual XML_Node& save(XML_Node& o, const doublereal* const sol); 326 327 //! Restore the solution for this domain from an XML_Node 328 /*! 329 * Base class version of the general Domain1D restore function. Derived 330 * classes should call the base class method in addition to restoring 331 * their own data. 332 * 333 * @param dom XML_Node for this domain 334 * @param soln Current value of the solution vector, local to this object. 335 * @param loglevel 0 to suppress all output; 1 to show warnings; 2 for 336 * verbose output 337 * 338 * @deprecated The XML input format is deprecated and will be removed in 339 * Cantera 3.0. 340 */ 341 virtual void restore(const XML_Node& dom, doublereal* soln, int loglevel); 342 size()343 size_t size() const { 344 return m_nv*m_points; 345 } 346 347 /** 348 * Find the index of the first grid point in this domain, and 349 * the start of its variables in the global solution vector. 350 */ 351 void locate(); 352 353 /** 354 * Location of the start of the local solution vector in the global 355 * solution vector, 356 */ 357 virtual size_t loc(size_t j = 0) const { 358 return m_iloc; 359 } 360 361 /** 362 * The index of the first (i.e., left-most) grid point belonging to this 363 * domain. 364 */ firstPoint()365 size_t firstPoint() const { 366 return m_jstart; 367 } 368 369 /** 370 * The index of the last (i.e., right-most) grid point belonging to this 371 * domain. 372 */ lastPoint()373 size_t lastPoint() const { 374 return m_jstart + m_points - 1; 375 } 376 377 /** 378 * Set the left neighbor to domain 'left.' Method 'locate' is called to 379 * update the global positions of this domain and all those to its right. 380 */ linkLeft(Domain1D * left)381 void linkLeft(Domain1D* left) { 382 m_left = left; 383 locate(); 384 } 385 386 //! Set the right neighbor to domain 'right.' linkRight(Domain1D * right)387 void linkRight(Domain1D* right) { 388 m_right = right; 389 } 390 391 //! Append domain 'right' to this one, and update all links. append(Domain1D * right)392 void append(Domain1D* right) { 393 linkRight(right); 394 right->linkLeft(this); 395 } 396 397 //! Return a pointer to the left neighbor. left()398 Domain1D* left() const { 399 return m_left; 400 } 401 402 //! Return a pointer to the right neighbor. right()403 Domain1D* right() const { 404 return m_right; 405 } 406 407 //! Value of component n at point j in the previous solution. prevSoln(size_t n,size_t j)408 double prevSoln(size_t n, size_t j) const { 409 return m_slast[m_nv*j + n]; 410 } 411 412 //! Specify an identifying tag for this domain. setID(const std::string & s)413 void setID(const std::string& s) { 414 m_id = s; 415 } 416 id()417 std::string id() const { 418 if (m_id != "") { 419 return m_id; 420 } else { 421 return fmt::format("domain {}", m_index); 422 } 423 } 424 showSolution_s(std::ostream & s,const doublereal * x)425 virtual void showSolution_s(std::ostream& s, const doublereal* x) {} 426 427 //! Print the solution. 428 virtual void showSolution(const doublereal* x); 429 z(size_t jlocal)430 doublereal z(size_t jlocal) const { 431 return m_z[jlocal]; 432 } zmin()433 doublereal zmin() const { 434 return m_z[0]; 435 } zmax()436 doublereal zmax() const { 437 return m_z[m_points - 1]; 438 } 439 440 void setProfile(const std::string& name, double* values, double* soln); 441 grid()442 vector_fp& grid() { 443 return m_z; 444 } grid()445 const vector_fp& grid() const { 446 return m_z; 447 } grid(size_t point)448 doublereal grid(size_t point) const { 449 return m_z[point]; 450 } 451 452 //! called to set up initial grid, and after grid refinement 453 virtual void setupGrid(size_t n, const doublereal* z); 454 455 /** 456 * Writes some or all initial solution values into the global solution 457 * array, beginning at the location pointed to by x. This method is called 458 * by the Sim1D constructor, and allows default values or ones that have 459 * been set locally prior to installing this domain into the container to be 460 * written to the global solution vector. 461 */ 462 virtual void _getInitialSoln(doublereal* x); 463 464 //! Initial value of solution component \a n at grid point \a j. 465 virtual doublereal initialValue(size_t n, size_t j); 466 467 /** 468 * In some cases, a domain may need to set parameters that depend on the 469 * initial solution estimate. In such cases, the parameters may be set in 470 * method _finalize. This method is called just before the Newton solver is 471 * called, and the x array is guaranteed to be the local solution vector for 472 * this domain that will be used as the initial guess. If no such parameters 473 * need to be set, then method _finalize does not need to be overloaded. 474 */ _finalize(const doublereal * x)475 virtual void _finalize(const doublereal* x) {} 476 477 /** 478 * In some cases, for computational efficiency some properties (e.g. 479 * transport coefficients) may not be updated during Jacobian evaluations. 480 * Set this to `true` to force these properties to be udpated even while 481 * calculating Jacobian elements. 482 */ forceFullUpdate(bool update)483 void forceFullUpdate(bool update) { 484 m_force_full_update = update; 485 } 486 487 protected: 488 doublereal m_rdt; 489 size_t m_nv; 490 size_t m_points; 491 vector_fp m_slast; 492 vector_fp m_max; 493 vector_fp m_min; 494 vector_fp m_rtol_ss, m_rtol_ts; 495 vector_fp m_atol_ss, m_atol_ts; 496 vector_fp m_z; 497 OneDim* m_container; 498 size_t m_index; 499 int m_type; 500 501 //! Starting location within the solution vector for unknowns that 502 //! correspond to this domain 503 /*! 504 * Remember there may be multiple domains associated with this problem 505 */ 506 size_t m_iloc; 507 508 size_t m_jstart; 509 510 Domain1D* m_left, *m_right; 511 512 //! Identity tag for the domain 513 std::string m_id; 514 std::string m_desc; 515 std::unique_ptr<Refiner> m_refiner; 516 std::vector<std::string> m_name; 517 int m_bw; 518 bool m_force_full_update; 519 }; 520 } 521 522 #endif 523