diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 2d07f896a9..7b46b929c1 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -690,7 +690,7 @@ jobs: echo "figure.max_open_warning: 100" > matplotlibrc export LD_LIBRARY_PATH=build/lib find samples/python -type f -iname "*.py" \ - -exec sh -c 'for n; do echo "$n" | tee -a results.txt && python3 "$n" >> results.txt || exit 1; done' sh {} + + -exec sh -c 'for n; do echo "$n" | tee -a results.txt && python3 -u "$n" >> results.txt || exit 1; done' sh {} + env: # The pyparsing ignore setting is due to a new warning introduced in Matplotlib==3.6.0 # Ignore NasaPoly2 warnings from n-hexane-NUIG-2015.yaml diff --git a/doc/sphinx/develop/reactor-integration.md b/doc/sphinx/develop/reactor-integration.md index e5f97fdc15..a74c47dc48 100644 --- a/doc/sphinx/develop/reactor-integration.md +++ b/doc/sphinx/develop/reactor-integration.md @@ -24,10 +24,11 @@ reac = ct.IdealGasReactor(gas, clone=False) sim = ct.ReactorNet([reac]) ``` -The `__init__` method of the Python {py:class}`ReactorNet` class calls -{ct}`ReactorNet::addReactor` for each {py:class}`Reactor` object provided in the list -supplied. When the first {ct}`Reactor` is added to the network, the {ct}`ReactorNet` -creates a new {ct}`Integrator` used to integrate the governing equations. +The constructor for the {py:class}`ReactorNet` class builds the network from the +provided list of {py:class}`Reactor` and {py:class}`ReactorSurface` objects and creates +a new {ct}`Integrator` used to integrate the governing equations. If there are +connections such as walls and mass flow controllers between different reactors, the +adjacent reactors are automatically added to the reactor network. The {ct}`Integrator` class is Cantera's interface for ODE/DAE system integrators. {ct}`Integrator` is a [polymorphic base class](http://www.cplusplus.com/doc/tutorial/polymorphism/); diff --git a/include/cantera/kinetics/ImplicitSurfChem.h b/include/cantera/kinetics/ImplicitSurfChem.h deleted file mode 100644 index 109884a1f7..0000000000 --- a/include/cantera/kinetics/ImplicitSurfChem.h +++ /dev/null @@ -1,304 +0,0 @@ -/** - * @file ImplicitSurfChem.h - * Declarations for the implicit integration of surface site density equations - * (see @ref kineticsmgr and class - * @link Cantera::ImplicitSurfChem ImplicitSurfChem@endlink). - */ - -// This file is part of Cantera. See License.txt in the top-level directory or -// at https://cantera.org/license.txt for license and copyright information. - -#ifndef CT_IMPSURFCHEM_H -#define CT_IMPSURFCHEM_H - -#include "cantera/numerics/Integrator.h" -#include "cantera/kinetics/InterfaceKinetics.h" -#include "cantera/kinetics/solveSP.h" - -namespace Cantera -{ - -//! @defgroup surfSolverGroup Surface Problem Solver - -//! Advances the surface coverages of the associated set of SurfacePhase -//! objects in time -/*! - * This function advances a set of SurfPhase objects, each associated with one - * InterfaceKinetics object, in time. The following equation is used for each - * surface phase, *i*. - * - * @f[ - * \dot \theta_k = \dot s_k (\sigma_k / s_0) - * @f] - * - * In this equation, - * - @f$ \theta_k @f$ is the site coverage for the kth species. - * - @f$ \dot s_k @f$ is the source term for the kth species - * - @f$ \sigma_k @f$ is the number of surface sites covered by each species k. - * - @f$ s_0 @f$ is the total site density of the interfacial phase. - * - * Additionally, the 0'th equation in the set is discarded. Instead the - * alternate equation is solved for - * - * @f[ - * \sum_{k=0}^{N-1} \dot \theta_k = 0 - * @f] - * - * This last equation serves to ensure that sum of the @f$ \theta_k @f$ values - * stays constant. - * - * The object uses the CVODE software to advance the surface equations. - * - * The solution vector used by this object is as follows: For each surface - * phase with @f$ N_s @f$ surface sites, it consists of the surface coverages - * @f$ \theta_k @f$ for @f$ k = 0, N_s - 1 @f$ - * - * @ingroup surfSolverGroup - */ -class ImplicitSurfChem : public FuncEval -{ -public: - //! Constructor for multiple surfaces. - /*! - * @param k Vector of pointers to InterfaceKinetics objects Each object - * consists of a surface or an edge containing internal degrees of - * freedom representing the concentration of surface adsorbates. - * @param rtol The relative tolerance for the integrator - * @param atol The absolute tolerance for the integrator - * @param maxStepSize The maximum step-size the integrator is allowed to take. - * If zero, this option is disabled. - * @param maxSteps The maximum number of time-steps the integrator can take. - * If not supplied, uses the default value in the CVodesIntegrator (20000). - * @param maxErrTestFails the maximum permissible number of error test failures - * If not supplied, uses the default value in CVODES (7). - */ - ImplicitSurfChem(vector k, - double rtol=1.e-7, double atol=1.e-14, - double maxStepSize=0, size_t maxSteps=20000, - size_t maxErrTestFails=7); - - /** - * Must be called before calling method 'advance' - */ - void initialize(double t0=0.0); - - /** - * Set the maximum integration step-size. Note, setting this value to zero - * disables this option - */ - void setMaxStepSize(double maxstep=0.0); - - /** - * Set the relative and absolute integration tolerances. - */ - void setTolerances(double rtol=1.e-7, double atol=1.e-14); - - /** - * Set the maximum number of CVODES integration steps. - */ - void setMaxSteps(size_t maxsteps=20000); - - /** - * Set the maximum number of CVODES error test failures - */ - void setMaxErrTestFails(size_t maxErrTestFails=7); - - //! Integrate from t0 to t1. The integrator is reinitialized first. - /*! - * This routine does a time accurate solve from t = t0 to t = t1. - * of the surface problem. - * - * @param t0 Initial Time -> this is an input - * @param t1 Final Time -> This is an input - */ - void integrate(double t0, double t1); - - //! Integrate from t0 to t1 without reinitializing the integrator. - /*! - * Use when the coverages have not changed from their values on return - * from the last call to integrate or integrate0. - * - * @param t0 Initial Time -> this is an input - * @param t1 Final Time -> This is an input - */ - void integrate0(double t0, double t1); - - //! Solve for the pseudo steady-state of the surface problem - /*! - * Solve for the steady state of the surface problem. - * This is the same thing as the advanceCoverages() function, - * but at infinite times. - * - * Note, a direct solve is carried out under the hood here, - * to reduce the computational time. - * - * @param ifuncOverride One of the values defined in @ref solvesp_methods. - * The default is -1, which means that the program will decide. - * @param timeScaleOverride When a pseudo transient is selected this value - * can be used to override the default time scale for - * integration which is one. When SFLUX_TRANSIENT is used, this - * is equal to the time over which the equations are integrated. - * When SFLUX_INITIALIZE is used, this is equal to the time used - * in the initial transient algorithm, before the equation - * system is solved directly. - */ - void solvePseudoSteadyStateProblem(int ifuncOverride = -1, - double timeScaleOverride = 1.0); - - // overloaded methods of class FuncEval - - //! Return the number of equations - size_t neq() const override { - return m_nv; - } - - //! Evaluate the value of ydot[k] at the current conditions - /*! - * @param t Time (seconds) - * @param y Vector containing the current solution vector - * @param ydot Output vector containing the value of the - * derivative of the surface coverages. - * @param p Unused parameter pass-through parameter vector - */ - void eval(double t, double* y, double* ydot, double* p) override; - - //! Get the current state of the solution vector - /*! - * @param y Value of the solution vector to be used. - * On output, this contains the initial value - * of the solution. - */ - void getState(double* y) override; - - /** - * Get the specifications for the problem from the values - * in the ThermoPhase objects for all phases. - * - * 1. concentrations of all species in all phases, #m_concSpecies - * 2. Temperature and pressure - * - * @param vecConcSpecies Vector of concentrations. The phase concentration - * vectors are contiguous within the object, in the same - * order as the unknown vector. - */ - void getConcSpecies(double* const vecConcSpecies) const; - - //! Sets the concentrations within phases that are unknowns in - //! the surface problem - /*! - * Fills the local concentration vector for all of the species in all of - * the phases that are unknowns in the surface problem. - * - * @param vecConcSpecies Vector of concentrations. The phase concentration - * vectors are contiguous within the object, in the same - * order as the unknown vector. - */ - void setConcSpecies(const double* const vecConcSpecies); - - //! Sets the state variable in all thermodynamic phases (surface and - //! surrounding bulk phases) to the input temperature and pressure - /*! - * @param TKelvin input temperature (kelvin) - * @param PresPa input pressure in pascal. - */ - void setCommonState_TP(double TKelvin, double PresPa); - - //! Returns a reference to the vector of pointers to the - //! InterfaceKinetics objects - /*! - * This should probably go away in the future, as it opens up the class. - */ - vector & getObjects() { - return m_vecKinPtrs; - } - - int checkMatch(vector m_vec, ThermoPhase* thPtr); - - void setIOFlag(int ioFlag) { - m_ioFlag = ioFlag; - } - -protected: - //! Set the mixture to a state consistent with solution - //! vector y. - /*! - * This function will set the surface site factions in the underlying - * SurfPhase objects to the current value of the solution vector. - * - * @param y Current value of the solution vector. The length is equal to - * the sum of the number of surface sites in all the surface phases. - */ - void updateState(double* y); - - //! vector of pointers to surface phases. - vector m_surf; - - //! Vector of pointers to bulk phases - vector m_bulkPhases; - - //! vector of pointers to InterfaceKinetics objects - vector m_vecKinPtrs; - - //! Vector of number of species in each Surface Phase - vector m_nsp; - - vector m_specStartIndex; - - //! Total number of surface species in all surface phases - /*! - * This is the total number of unknowns in m_mode 0 problem - */ - size_t m_nv = 0; - - size_t m_numTotalBulkSpecies = 0; - size_t m_numTotalSpecies = 0; - - vector> pLocVec; - //! Pointer to the CVODE integrator - unique_ptr m_integ; - double m_atol, m_rtol; // tolerances - double m_maxstep; //!< max step size - size_t m_nmax; //!< maximum number of steps allowed - size_t m_maxErrTestFails; //!< maximum number of error test failures allowed - vector m_work; - - /** - * Temporary vector - length num species in the Kinetics object. This is - * the sum of the number of species in each phase included in the kinetics - * object. - */ - vector m_concSpecies; - vector m_concSpeciesSave; - - /** - * Index into the species vector of the kinetics manager, - * pointing to the first species from the surrounding medium. - */ - int m_mediumSpeciesStart = -1; - /** - * Index into the species vector of the kinetics manager, pointing to the - * first species from the condensed phase of the particles. - */ - int m_bulkSpeciesStart = -1; - /** - * Index into the species vector of the kinetics manager, pointing to the - * first species from the surface of the particles - */ - int m_surfSpeciesStart = -1; - /** - * Pointer to the helper method, Placid, which solves the surface problem. - */ - unique_ptr m_surfSolver; - - //! If true, a common temperature and pressure for all surface and bulk - //! phases associated with the surface problem is imposed - bool m_commonTempPressForPhases = true; - -private: - //! Controls the amount of printing from this routine - //! and underlying routines. - int m_ioFlag = 0; -}; -} - -#endif diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h index 029c76f153..874857e5ad 100644 --- a/include/cantera/kinetics/InterfaceKinetics.h +++ b/include/cantera/kinetics/InterfaceKinetics.h @@ -15,7 +15,7 @@ namespace Cantera // forward declarations class SurfPhase; -class ImplicitSurfChem; +class ReactorNet; //! A kinetics manager for heterogeneous reaction mechanisms. The reactions are //! assumed to occur at a 2D interface between two 3D phases. @@ -56,7 +56,6 @@ class InterfaceKinetics : public Kinetics public: //! Constructor InterfaceKinetics() = default; - ~InterfaceKinetics() override; void resizeReactions() override; @@ -162,6 +161,10 @@ class InterfaceKinetics : public Kinetics * \dot {\theta}_k = \dot s_k (\sigma_k / s_0) * @f] * + * For the purpose of this integration, the temperature and pressure of the adjacent + * bulk phases set equal to those of the surface phase, and the composition of the + * bulk phases is held fixed. + * * @param tstep Time value to advance the surface coverages * @param rtol The relative tolerance for the integrator * @param atol The absolute tolerance for the integrator @@ -179,25 +182,13 @@ class InterfaceKinetics : public Kinetics //! Solve for the pseudo steady-state of the surface problem /*! * This is the same thing as the advanceCoverages() function, - * but at infinite times. - * - * Note, a direct solve is carried out under the hood here, - * to reduce the computational time. + * but at infinite times. Carries out a direct solution of the nonlinear + * algebraic system using SteadyReactorSolver. * - * @param ifuncOverride One of the values defined in @ref solvesp_methods. - * The default is -1, which means that the program will decide. - * @param timeScaleOverride When a pseudo transient is selected this value - * can be used to override the default time scale for - * integration which is one. When SFLUX_TRANSIENT is used, this - * is equal to the time over which the equations are integrated. - * When SFLUX_INITIALIZE is used, this is equal to the time used - * in the initial transient algorithm, before the equation - * system is solved directly. + * @param loglevel Log level argument passed to SteadyReactorSolver + * @since In %Cantera 4.0, the previous arguments were replaced with `loglevel`. */ - void solvePseudoSteadyStateProblem(int ifuncOverride = -1, - double timeScaleOverride = 1.0); - - void setIOFlag(int ioFlag); + void solvePseudoSteadyStateProblem(int loglevel=0); //! Update the standard state chemical potentials and species equilibrium //! constant entries @@ -311,6 +302,10 @@ class InterfaceKinetics : public Kinetics //! included void assertDerivativesValid(const string& name); + //! Create a ReactorNet where only the surface species are active, for use with + //! advanceCoverages() and solvePseudoSteadyStateProblem(). + void buildNetwork(); + //! @} //! Temporary work vector of length m_kk @@ -394,7 +389,7 @@ class InterfaceKinetics : public Kinetics * be used to solve this single InterfaceKinetics object's surface problem * uncoupled from other surface phases. */ - ImplicitSurfChem* m_integrator = nullptr; + ReactorNet* m_integrator = nullptr; bool m_ROP_ok = false; @@ -449,8 +444,6 @@ class InterfaceKinetics : public Kinetics */ vector> m_rxnPhaseIsProduct; - int m_ioFlag = 0; - //! Number of dimensions of reacting phase (2 for InterfaceKinetics, 1 for //! EdgeKinetics) size_t m_nDim = 2; diff --git a/include/cantera/kinetics/solveSP.h b/include/cantera/kinetics/solveSP.h deleted file mode 100644 index b4536c8261..0000000000 --- a/include/cantera/kinetics/solveSP.h +++ /dev/null @@ -1,489 +0,0 @@ -/** - * @file solveSP.h Header file for implicit surface problem solver (see @ref - * chemkinetics and class @link Cantera::solveSP solveSP@endlink). - */ - -// This file is part of Cantera. See License.txt in the top-level directory or -// at https://cantera.org/license.txt for license and copyright information. - -#ifndef SOLVESP_H -#define SOLVESP_H - -#include "cantera/kinetics/InterfaceKinetics.h" -#include "cantera/numerics/DenseMatrix.h" - -//! @defgroup solvesp_methods Surface Problem Solver Methods -//! @ingroup surfSolverGroup -//! @{ - -//! This assumes that the initial guess supplied to the routine is far from -//! the correct one. Substantial work plus transient time-stepping is to be -//! expected to find a solution. -const int SFLUX_INITIALIZE = 1; - -//! Need to solve the surface problem in order to calculate the surface fluxes -//! of gas-phase species. (Can expect a moderate change in the solution -//! vector; try to solve the system by direct methods with no damping first, -//! then try time-stepping if the first method fails). A "time_scale" supplied -//! here is used in the algorithm to determine when to shut off time-stepping. -const int SFLUX_RESIDUAL = 2; - -//! Calculation of the surface problem is due to the need for a numerical -//! Jacobian for the gas-problem. The solution is expected to be very close to -//! the initial guess, and accuracy is needed because solution variables have -//! been perturbed from nominal values to create Jacobian entries. -const int SFLUX_JACOBIAN = 3; - -//! The transient calculation is performed here for an amount of time -//! specified by "time_scale". It is not guaranteed to be time-accurate - -//! just stable and fairly fast. The solution after del_t time is returned, -//! whether it's converged to a steady state or not. This is a poor man's time -//! stepping algorithm. -const int SFLUX_TRANSIENT = 4; -//! @} End of solvesp_method - -//! @defgroup solvesp_bulkFunc Surface Problem Bulk Phase Mode -//! Functionality expected from the bulk phase. This changes the equations -//! that will be used to solve for the bulk mole fractions. -//! @ingroup surfSolverGroup -//! @{ - -//! Deposition of a bulk phase is to be expected. Bulk mole fractions are -//! determined from ratios of growth rates of bulk species. -const int BULK_DEPOSITION = 1; - -//! Etching of a bulk phase is to be expected. Bulk mole fractions are assumed -//! constant, and given by the initial conditions. This is also used whenever -//! the condensed phase is part of the larger solution. -const int BULK_ETCH = 2; -//! @} End of solvesp_bulkFunc - -namespace Cantera -{ - -//! @addtogroup surfSolverGroup -//! @{ - -//! Method to solve a pseudo steady state surface problem -/*! - * The following class handles solving the surface problem. The calculation - * uses Newton's method to obtain the surface fractions of the surface and - * bulk species by requiring that the surface species production rate = 0 and - * that the either the bulk fractions are proportional to their production - * rates or they are constants. - * - * Currently, the bulk mole fractions are treated as constants. Implementation - * of their being added to the unknown solution vector is delayed. - * - * Lets introduce the unknown vector for the "surface problem". The surface - * problem is defined as the evaluation of the surface site fractions for - * multiple surface phases. The unknown vector will consist of the vector of - * surface concentrations for each species in each surface vector. Species - * are grouped first by their surface phases - * - * - C_i_j = Concentration of the ith species in the jth surface phase - * - Nj = number of surface species in the jth surface phase - * - * The unknown solution vector is defined as follows: - * - * C_i_j | kindexSP - * --------- | ---------- - * C_0_0 | 0 - * C_1_0 | 1 - * C_2_0 | 2 - * . . . | ... - * C_N0-1_0 | N0-1 - * C_0_1 | N0 - * C_1_1 | N0+1 - * C_2_1 | N0+2 - * . . . | ... - * C_N1-1_1 | NO+N1-1 - * - * Note there are a couple of different types of species indices floating - * around in the formulation of this object. - * - * kindexSP: This is the species index in the contiguous vector of unknowns - * for the surface problem. - * - * Note, in the future, BULK_DEPOSITION systems will be added, and the - * solveSP unknown vector will get more complicated. It will include the mole - * fraction and growth rates of specified bulk phases - * - * Indices which relate to individual kinetics objects use the suffix KSI - * (kinetics species index). - * - * ## Solution Method - * - * This routine is typically used within a residual calculation in a large code. - * It's typically invoked millions of times for large calculations, and it must - * work every time. Therefore, requirements demand that it be robust but also - * efficient. - * - * The solution methodology is largely determined by the `ifunc` parameter, - * that is input to the solution object. This parameter may have one of the - * values defined in @ref solvesp_methods. - * - * ### Pseudo time stepping algorithm: - * The time step is determined from sdot[], so that the time step - * doesn't ever change the value of a variable by more than 100%. - * - * This algorithm does use a damped Newton's method to relax the equations. - * Damping is based on a "delta damping" technique. The solution unknowns - * are not allowed to vary too much between iterations. - * - * `EXTRA_ACCURACY`: A constant that is the ratio of the required update norm - * in this Newton iteration compared to that in the nonlinear solver. A value - * of 0.1 is used so surface species are safely overconverged. - * - * Functions called: - * - `ct_dgetrf` -- First half of LAPACK direct solve of a full Matrix - * - `ct_dgetrs` -- Second half of LAPACK direct solve of a full matrix. - * Returns solution vector in the right-hand-side vector, resid. - */ -class solveSP -{ -public: - //! Constructor for the object - /*! - * @param surfChemPtr Pointer to the ImplicitSurfChem object that - * defines the surface problem to be solved. - * @param bulkFunc Integer representing how the bulk phases should be - * handled. See @ref solvesp_bulkFunc. Currently, - * only the default value of BULK_ETCH is supported. - */ - solveSP(ImplicitSurfChem* surfChemPtr, int bulkFunc = BULK_ETCH); - - //! Destructor - ~solveSP() = default; - - solveSP(const solveSP&) = delete; - solveSP& operator=(const solveSP&) = delete; - -public: - //! Main routine that actually calculates the pseudo steady state - //! of the surface problem - /*! - * The actual converged solution is returned as part of the internal state - * of the InterfaceKinetics objects. - * - * Uses Newton's method to get the surface fractions of the surface and - * bulk species by requiring that the surface species production rate = 0 - * and that the bulk fractions are proportional to their production rates. - * - * @param ifunc Determines the type of solution algorithm to be used. See - * @ref solvesp_methods for possible values. - * @param time_scale Time over which to integrate the surface equations, - * where applicable - * @param TKelvin Temperature (kelvin) - * @param PGas Pressure (pascals) - * @param reltol Relative tolerance to use - * @param abstol absolute tolerance. - * @return 1 if the surface problem is successfully solved. - * -1 if the surface problem wasn't solved successfully. - * Note the actual converged solution is returned as part of the - * internal state of the InterfaceKinetics objects. - */ - int solveSurfProb(int ifunc, double time_scale, double TKelvin, - double PGas, double reltol, double abstol); - -private: - //! Printing routine that optionally gets called at the start of every - //! invocation - void print_header(int ioflag, int ifunc, double time_scale, - int damping, double reltol, double abstol); - - //! Printing routine that gets called after every iteration - void printIteration(int ioflag, double damp, int label_d, int label_t, - double inv_t, double t_real, size_t iter, - double update_norm, double resid_norm, - bool do_time, bool final=false); - - //! Calculate a conservative delta T to use in a pseudo-steady state - //! algorithm - /*! - * This routine calculates a pretty conservative 1/del_t based on - * MAX_i(sdot_i/(X_i*SDen0)). This probably guarantees diagonal dominance. - * - * Small surface fractions are allowed to intervene in the del_t - * determination, no matter how small. This may be changed. - * Now minimum changed to 1.0e-12, - * - * Maximum time step set to time_scale. - * - * @param netProdRateSolnSP Output variable. Net production rate of all - * of the species in the solution vector. - * @param XMolSolnSP output variable. Mole fraction of all of the species - * in the solution vector - * @param label Output variable. Pointer to the value of the species - * index (kindexSP) that is controlling the time step - * @param label_old Output variable. Pointer to the value of the species - * index (kindexSP) that controlled the time step at the previous - * iteration - * @param label_factor Output variable. Pointer to the current factor - * that is used to indicate the same species is controlling the time - * step. - * @param ioflag Level of the output requested. - * @returns the 1. / delta T to be used on the next step - */ - double calc_t(double netProdRateSolnSP[], double XMolSolnSP[], - int* label, int* label_old, - double* label_factor, int ioflag); - - //! Calculate the solution and residual weights - /*! - * @param wtSpecies Weights to use for the soln unknowns. These are in - * concentration units - * @param wtResid Weights to sue for the residual unknowns. - * @param Jac Jacobian. Row sum scaling is used for the Jacobian - * @param CSolnSP Solution vector for the surface problem - * @param abstol Absolute error tolerance - * @param reltol Relative error tolerance - */ - void calcWeights(double wtSpecies[], double wtResid[], - const Array2D& Jac, const double CSolnSP[], - const double abstol, const double reltol); - - /** - * Update the surface states of the surface phases. - */ - void updateState(const double* cSurfSpec); - - //! Update mole fraction vector consisting of unknowns in surface problem - /*! - * @param XMolSolnSP Vector of mole fractions for the unknowns in the - * surface problem. - */ - void updateMFSolnSP(double* XMolSolnSP); - - //! Update the mole fraction vector for a specific kinetic species vector - //! corresponding to one InterfaceKinetics object - /*! - * @param XMolKinSp Mole fraction vector corresponding to a particular - * kinetic species for a single InterfaceKinetics Object - * This is a vector over all the species in all of the - * phases in the InterfaceKinetics object - * @param isp ID of the InterfaceKinetics Object. - */ - void updateMFKinSpecies(double* XMolKinSp, int isp); - - //! Update the vector that keeps track of the largest species in each - //! surface phase. - /*! - * @param CSolnSP Vector of the current values of the surface concentrations - * in all of the surface species. - */ - void evalSurfLarge(const double* CSolnSP); - - //! Main Function evaluation - /*! - * @param resid output Vector of residuals, length = m_neq - * @param CSolnSP Vector of species concentrations, unknowns in the - * problem, length = m_neq - * @param CSolnOldSP Old Vector of species concentrations, unknowns in the - * problem, length = m_neq - * @param do_time Calculate a time dependent residual - * @param deltaT Delta time for time dependent problem. - */ - void fun_eval(double* resid, const double* CSolnSP, - const double* CSolnOldSP, const bool do_time, const double deltaT); - - //! Main routine that calculates the current residual and Jacobian - /*! - * @param jac Jacobian to be evaluated. - * @param resid output Vector of residuals, length = m_neq - * @param CSolnSP Vector of species concentrations, unknowns in the - * problem, length = m_neq. These are tweaked in order - * to derive the columns of the Jacobian. - * @param CSolnSPOld Old Vector of species concentrations, unknowns in the - * problem, length = m_neq - * @param do_time Calculate a time dependent residual - * @param deltaT Delta time for time dependent problem. - */ - void resjac_eval(DenseMatrix& jac, double* resid, - double* CSolnSP, - const double* CSolnSPOld, const bool do_time, - const double deltaT); - - //! Vector of interface kinetics objects - /*! - * Each of these is associated with one and only one surface phase. - */ - vector &m_objects; - - //! Total number of equations to solve in the implicit problem. - /*! - * Note, this can be zero, and frequently is - */ - size_t m_neq = 0; - - //! This variable determines how the bulk phases are to be handled - /*! - * Possible values are given in @ref solvesp_bulkFunc. - */ - int m_bulkFunc; - - //! Number of surface phases in the surface problem - /*! - * This number is equal to the number of InterfaceKinetics objects - * in the problem. (until further noted) - */ - size_t m_numSurfPhases = 0; - - //! Total number of surface species in all surface phases. - /*! - * This is also the number of equations to solve for m_mode=0 system - * It's equal to the sum of the number of species in each of the - * m_numSurfPhases. - */ - size_t m_numTotSurfSpecies = 0; - - //! Mapping between the surface phases and the InterfaceKinetics objects - /*! - * Currently this is defined to be a 1-1 mapping (and probably assumed - * in some places) - * m_surfKinObjID[i] = i - */ - vector m_indexKinObjSurfPhase; - - //! Vector of length number of surface phases containing - //! the number of surface species in each phase - /*! - * Length is equal to the number of surface phases, m_numSurfPhases - */ - vector m_nSpeciesSurfPhase; - - //! Vector of surface phase pointers - /*! - * This is created during the constructor - * Length is equal to the number of surface phases, m_numSurfPhases - */ - vector m_ptrsSurfPhase; - - //! Index of the start of the unknowns for each solution phase - /*! - * i_eqn = m_eqnIndexStartPhase[isp] - * - * isp is the phase id in the list of phases solved by the - * surface problem. - * - * i_eqn is the equation number of the first unknown in the - * solution vector corresponding to isp'th phase. - */ - vector m_eqnIndexStartSolnPhase; - - //! Total number of volumetric condensed phases included in the steady state - //! problem handled by this routine. - /*! - * This is equal to or less than the total number of volumetric phases in - * all of the InterfaceKinetics objects. We usually do not include bulk - * phases. Bulk phases are only included in the calculation when their - * domain isn't included in the underlying continuum model conservation - * equation system. - * - * This is equal to 0, for the time being - */ - size_t m_numBulkPhasesSS = 0; - - //! Vector of number of species in the m_numBulkPhases phases. - /*! - * Length is number of bulk phases - */ - vector m_numBulkSpecies; - - //! Total number of species in all bulk phases. - /*! - * This is also the number of bulk equations to solve when bulk equation - * solving is turned on. - */ - size_t m_numTotBulkSpeciesSS = 0; - - //! Vector of bulk phase pointers, length is equal to m_numBulkPhases. - vector m_bulkPhasePtrs; - - //! Index between the equation index and the position in the kinetic - //! species array for the appropriate kinetics operator - /*! - * Length = m_neq. - * - * ksp = m_kinSpecIndex[ieq] - * ksp is the kinetic species index for the ieq'th equation. - */ - vector m_kinSpecIndex; - - //! Index between the equation index and the index of the - //! InterfaceKinetics object - /*! - * Length m_neq - */ - vector m_kinObjIndex; - - //! Vector containing the indices of the largest species - //! in each surface phase - /*! - * `k = m_spSurfLarge[i]` where `k` is the local species index, that is, it - * varies from 0 to (num species in phase - 1) and `i` is the surface - * phase index in the problem. Length is equal to #m_numSurfPhases. - */ - vector m_spSurfLarge; - - //! Maximum number of species in any single kinetics operator - //! -> also maxed wrt the total # of solution species - size_t m_maxTotSpecies = 0; - - //! Temporary vector with length equal to max m_maxTotSpecies - vector m_netProductionRatesSave; - - //! Temporary vector with length equal to max m_maxTotSpecies - vector m_numEqn1; - - //! Temporary vector with length equal to max m_maxTotSpecies - vector m_numEqn2; - - //! Temporary vector with length equal to max m_maxTotSpecies - vector m_CSolnSave; - - //! Solution vector. length MAX(1, m_neq) - vector m_CSolnSP; - - //! Saved initial solution vector. length MAX(1, m_neq) - vector m_CSolnSPInit; - - //! Saved solution vector at the old time step. length MAX(1, m_neq) - vector m_CSolnSPOld; - - //! Weights for the residual norm calculation. length MAX(1, m_neq) - vector m_wtResid; - - //! Weights for the species concentrations norm calculation - /*! - * length MAX(1, m_neq) - */ - vector m_wtSpecies; - - //! Residual for the surface problem - /*! - * The residual vector of length "dim" that, that has the value of "sdot" - * for surface species. The residuals for the bulk species are a function - * of the sdots for all species in the bulk phase. The last residual of - * each phase enforces {Sum(fractions) = 1}. After linear solve (dgetrf_ & - * dgetrs_), resid holds the update vector. - * - * length MAX(1, m_neq) - */ - vector m_resid; - - //! Vector of mole fractions. length m_maxTotSpecies - vector m_XMolKinSpecies; - - //! Jacobian. m_neq by m_neq computed Jacobian matrix for the local - //! Newton's method. - DenseMatrix m_Jac; - -public: - int m_ioflag = 0; -}; - -//! @} End of surfSolverGroup - -} -#endif diff --git a/include/cantera/numerics/CVodesIntegrator.h b/include/cantera/numerics/CVodesIntegrator.h index df6d902f42..a690542902 100644 --- a/include/cantera/numerics/CVodesIntegrator.h +++ b/include/cantera/numerics/CVodesIntegrator.h @@ -20,8 +20,7 @@ namespace Cantera /** * Wrapper class for 'cvodes' integrator from LLNL. * - * See FuncEval.h. Classes that use CVodesIntegrator: - * ImplicitSurfChem, ReactorNet + * See FuncEval.h. Classes that use CVodesIntegrator: ReactorNet */ class CVodesIntegrator : public Integrator { diff --git a/include/cantera/zeroD/ConstPressureMoleReactor.h b/include/cantera/zeroD/ConstPressureMoleReactor.h index 0ea9c0ce20..ea589ce0db 100644 --- a/include/cantera/zeroD/ConstPressureMoleReactor.h +++ b/include/cantera/zeroD/ConstPressureMoleReactor.h @@ -30,12 +30,8 @@ class ConstPressureMoleReactor : public MoleReactor void getState(double* y) override; void eval(double t, double* LHS, double* RHS) override; - - vector steadyConstraints() const override { - throw CanteraError("ConstPressureMoleReactor::steadyConstraints", - "ConstPressureMoleReactor is not currently compatible with the steady-state" - " solver.\nSee https://github.com/Cantera/enhancements/issues/234"); - } + void evalSteady(double t, double* LHS, double* RHS) override; + vector initializeSteady() override; void updateState(double* y) override; @@ -47,6 +43,7 @@ class ConstPressureMoleReactor : public MoleReactor protected: const size_t m_sidx = 1; + double m_initialMass; //!< Initial mass [kg]; used for steady-state calculations }; } diff --git a/include/cantera/zeroD/ConstPressureReactor.h b/include/cantera/zeroD/ConstPressureReactor.h index 9566289b74..75875070a6 100644 --- a/include/cantera/zeroD/ConstPressureReactor.h +++ b/include/cantera/zeroD/ConstPressureReactor.h @@ -33,7 +33,8 @@ class ConstPressureReactor : public Reactor void getState(double* y) override; void eval(double t, double* LHS, double* RHS) override; - vector steadyConstraints() const override; + void evalSteady(double t, double* LHS, double* RHS) override; + vector initializeSteady() override; void updateState(double* y) override; @@ -46,6 +47,9 @@ class ConstPressureReactor : public Reactor double upperBound(size_t k) const override; double lowerBound(size_t k) const override; void resetBadValues(double* y) override; + +protected: + double m_initialMass; //!< Initial mass [kg]; used for steady-state calculations }; } diff --git a/include/cantera/zeroD/FlowReactor.h b/include/cantera/zeroD/FlowReactor.h index ea2f4a481a..5c8a662bc2 100644 --- a/include/cantera/zeroD/FlowReactor.h +++ b/include/cantera/zeroD/FlowReactor.h @@ -47,8 +47,8 @@ class FlowReactor : public IdealGasReactor void evalDae(double t, double* y, double* ydot, double* residual) override; void getConstraints(double* constraints) override; - vector steadyConstraints() const override { - throw CanteraError("FlowReactor::steadyConstraints", + vector initializeSteady() override { + throw CanteraError("FlowReactor::initializeSteady", "FlowReactor is not compatible with time-dependent steady-state solver."); } diff --git a/include/cantera/zeroD/IdealGasConstPressureReactor.h b/include/cantera/zeroD/IdealGasConstPressureReactor.h index 7203a3fcd1..19f136845e 100644 --- a/include/cantera/zeroD/IdealGasConstPressureReactor.h +++ b/include/cantera/zeroD/IdealGasConstPressureReactor.h @@ -32,8 +32,9 @@ class IdealGasConstPressureReactor : public ConstPressureReactor void initialize(double t0=0.0) override; void eval(double t, double* LHS, double* RHS) override; + void evalSteady(double t, double* LHS, double* RHS) override; void updateState(double* y) override; - vector steadyConstraints() const override; + vector initializeSteady() override; //! Return the index in the solution vector for this reactor of the //! component named *nm*. Possible values for *nm* are "mass", @@ -46,6 +47,12 @@ class IdealGasConstPressureReactor : public ConstPressureReactor protected: vector m_hk; //!< Species molar enthalpies + + //! Initial mass [kg]; used for steady-state calculations + double m_initialMass; + + //! Initial temperature [K]; used for steady-state calculations + double m_initialTemperature; }; } diff --git a/include/cantera/zeroD/IdealGasMoleReactor.h b/include/cantera/zeroD/IdealGasMoleReactor.h index 3347eaa6f8..de2a684e1c 100644 --- a/include/cantera/zeroD/IdealGasMoleReactor.h +++ b/include/cantera/zeroD/IdealGasMoleReactor.h @@ -30,13 +30,14 @@ class IdealGasMoleReactor : public MoleReactor string componentName(size_t k) override; double upperBound(size_t k) const override; double lowerBound(size_t k) const override; - vector steadyConstraints() const override; + vector initializeSteady() override; void getState(double* y) override; void initialize(double t0=0.0) override; void eval(double t, double* LHS, double* RHS) override; + void evalSteady(double t, double* LHS, double* RHS) override; void updateState(double* y) override; @@ -53,6 +54,12 @@ class IdealGasMoleReactor : public MoleReactor protected: vector m_uk; //!< Species molar internal energies double m_TotalCv; //!< Total heat capacity (@f$ m c_v @f$) [J/K] + + //! Initial volume [m³]; used for steady-state calculations + double m_initialVolume; + + //! Initial temperature [K]; used for steady-state calculations + double m_initialTemperature; }; } diff --git a/include/cantera/zeroD/IdealGasReactor.h b/include/cantera/zeroD/IdealGasReactor.h index c8f74a5db9..eeb104e161 100644 --- a/include/cantera/zeroD/IdealGasReactor.h +++ b/include/cantera/zeroD/IdealGasReactor.h @@ -31,7 +31,8 @@ class IdealGasReactor : public Reactor void initialize(double t0=0.0) override; void eval(double t, double* LHS, double* RHS) override; - vector steadyConstraints() const override; + void evalSteady(double t, double* LHS, double* RHS) override; + vector initializeSteady() override; void updateState(double* y) override; //! Return the index in the solution vector for this reactor of the @@ -45,6 +46,12 @@ class IdealGasReactor : public Reactor protected: vector m_uk; //!< Species molar internal energies + + //! Initial volume [m³]; used for steady-state calculations + double m_initialVolume; + + //! Initial temperature [K]; used for steady-state calculations + double m_initialTemperature; }; } diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h index 11004fac09..f0db61d25e 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -83,7 +83,8 @@ class Reactor : public ReactorBase void getState(double* y) override; void initialize(double t0=0.0) override; void eval(double t, double* LHS, double* RHS) override; - vector steadyConstraints() const override; + void evalSteady(double t, double* LHS, double* RHS) override; + vector initializeSteady() override; void updateState(double* y) override; void addSensitivityReaction(size_t rxn) override; void addSensitivitySpeciesEnthalpy(size_t k) override; @@ -157,6 +158,9 @@ class Reactor : public ReactorBase bool m_energy = true; vector m_advancelimits; //!< Advance step limit + + //! Initial volume [m³]; used for steady-state calculations + double m_initialVolume; }; } diff --git a/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h index fb434b416a..e209a1e3a0 100644 --- a/include/cantera/zeroD/ReactorBase.h +++ b/include/cantera/zeroD/ReactorBase.h @@ -243,7 +243,8 @@ class ReactorBase : public std::enable_shared_from_this * @param[out] y state vector representing the initial state of the reactor */ virtual void getState(double* y) { - throw NotImplementedError("ReactorBase::getState"); + throw NotImplementedError("ReactorBase::getState", + "Not implemented for reactor type '{}'.", type()); } //! Get the current state and derivative vector of the reactor for a DAE solver @@ -277,20 +278,36 @@ class ReactorBase : public std::enable_shared_from_this throw NotImplementedError("ReactorBase::evalDae"); } + //! Evaluate the governing equations with modifications for the steady-state solver. + //! + //! This method calls the standard eval() method then modifies elements of `RHS` + //! that correspond to algebraic constraints. + //! + //! @since New in %Cantera 4.0. + virtual void evalSteady(double t, double* LHS, double* RHS) { + throw NotImplementedError("ReactorBase::evalSteady", + "Not implemented for reactor type '{}'.", type()); + } + //! Given a vector of length neq(), mark which variables should be //! considered algebraic constraints virtual void getConstraints(double* constraints) { throw NotImplementedError("ReactorBase::getConstraints"); } - //! Get the indices of equations that are algebraic constraints when solving the + //! Initialize the reactor before solving a steady-state problem. + //! + //! This method is responsible for storing the initial value for any algebraic + //! constraints and returning the indices of those constraints. + //! + //! @return Indices of equations that are algebraic constraints when solving the //! steady-state problem. //! //! @warning This method is an experimental part of the %Cantera API and may be //! changed or removed without notice. - //! @since New in %Cantera 3.2. - virtual vector steadyConstraints() const { - throw NotImplementedError("ReactorBase::steadyConstraints"); + //! @since New in %Cantera 3.2. Renamed from `steadyConstraints` in %Cantera 4.0. + virtual vector initializeSteady() { + throw NotImplementedError("ReactorBase::initializeSteady"); } //! Set the state of the reactor to correspond to the state vector *y*. diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index b4736e0668..ab01d7cc9d 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/ReactorNet.h @@ -29,16 +29,15 @@ class SystemJacobian; class ReactorNet : public FuncEval { public: - ReactorNet(); //! Create reactor network containing single reactor. //! @since New in %Cantera 3.2. ReactorNet(shared_ptr reactor); //! Create reactor network from multiple reactors. //! @since New in %Cantera 3.2. - ReactorNet(vector>& reactors); - ~ReactorNet() override; + ReactorNet(span> reactors); ReactorNet(const ReactorNet&) = delete; ReactorNet& operator=(const ReactorNet&) = delete; + ~ReactorNet(); //! @name Methods to set up a simulation //! @{ @@ -153,11 +152,8 @@ class ReactorNet : public FuncEval //! - The mass of constant pressure reactor types (such as ConstPressureReactor and //! IdealGasConstPressureReactor) must be constant; if flow devices are used, //! inlet and outlet flows must be balanced. - //! - The solver is currently not compatible with the ConstPressureMoleReactor or - //! IdealGasConstPressureMoleReactor classes. //! - Only ideal gas reactor types can be used for when the energy equation is //! disabled (fixed temperature simulations). - //! - Reacting surfaces are not yet supported. //! //! @param loglevel Print information about solver progress to aid in understanding //! cases where the solver fails to converge. Higher levels are more verbose. @@ -173,7 +169,7 @@ class ReactorNet : public FuncEval //! - 7: Print current residual vector after different solver stages //! //! @see SteadyStateSystem, MultiNewton - //! @since New in %Cantera 3.2. + //! @since New in %Cantera 3.2. Support for reacting surfaces added in %Cantera 4.0. void solveSteady(int loglevel=0); //! Get the Jacobian used by the steady-state solver. @@ -186,8 +182,8 @@ class ReactorNet : public FuncEval //! Return a reference to the *n*-th reactor in this network. The reactor //! indices are determined by the order in which the reactors were added //! to the reactor network. - Reactor& reactor(int n) { - return *m_bulkReactors[n]; + ReactorBase& reactor(int n) { + return *m_reactors[n]; } //! Returns `true` if verbose logging output is enabled. @@ -260,6 +256,15 @@ class ReactorNet : public FuncEval void eval(double t, double* y, double* ydot, double* p) override; + //! Evaluate the governing equations adapted for the steady-state solver. + //! + //! @param[in] y Current state vector, length neq() + //! @param[out] residual For differential variables, this is the time derivative; + //! for algebraic variables, this is the residual of the governing equation. + //! + //! @since New in %Cantera 4.0. + void evalSteady(double* y, double* residual); + //! eval coupling for IDA / DAEs void evalDae(double t, double* y, double* ydot, double* p, double* residual) override; @@ -329,7 +334,7 @@ class ReactorNet : public FuncEval //! Called to trigger integrator reinitialization before further //! integration. void setNeedsReinit() { - m_integrator_init = false; + m_needIntegratorInit = true; } //! Set the maximum number of internal integration steps the @@ -370,17 +375,13 @@ class ReactorNet : public FuncEval void evalRootFunctions(double t, const double* y, double* gout) override; protected: - //! Add the reactor *r* to this reactor network. - //! @since Changed in %Cantera 3.2. Previous version used a reference. - void addReactor(shared_ptr reactor); - //! Check that preconditioning is supported by all reactors in the network virtual void checkPreconditionerSupported() const; void updatePreconditioner(double gamma) override; //! Create reproducible names for reactors and walls/connectors. - void updateNames(Reactor& r); + void updateNames(ReactorBase& r); //! Estimate a future state based on current derivatives. //! The function is intended for internal use by ReactorNet::advance @@ -393,9 +394,9 @@ class ReactorNet : public FuncEval virtual int lastOrder() const; vector> m_reactors; - vector m_bulkReactors; - vector m_surfaces; - set m_reservoirs; + vector> m_bulkReactors; + vector> m_surfaces; + vector> m_reservoirs; set m_flowDevices; set m_walls; map m_counts; //!< Map used for default name generation @@ -408,8 +409,8 @@ class ReactorNet : public FuncEval //! The initial value of the independent variable in the system. double m_initial_time = 0.0; - bool m_init = false; - bool m_integrator_init = false; //!< True if integrator initialization is current + bool m_integratorInitialized = false; //!< True if the integrator has been initialized at least once + bool m_needIntegratorInit = true; //!< True if integrator needs to be (re)initialized size_t m_nv = 0; vector m_atol; @@ -473,10 +474,6 @@ class SteadyReactorSolver : public SteadyStateSystem //! Initial value of each state variable vector m_initialState; - - //! Indices of variables that are held constant in the time-stepping mode of the - //! steady-state solver. - vector m_algebraic; }; diff --git a/include/cantera/zeroD/ReactorSurface.h b/include/cantera/zeroD/ReactorSurface.h index 85362fc111..62ac42af8a 100644 --- a/include/cantera/zeroD/ReactorSurface.h +++ b/include/cantera/zeroD/ReactorSurface.h @@ -77,6 +77,18 @@ class ReactorSurface : public ReactorBase throw NotImplementedError("ReactorSurface::addSurface"); } + //! Get the number of Reactor and Reservoir objects adjacent to this surface + //! @since New in %Cantera 4.0. + size_t nAdjacent() const { + return m_reactors.size(); + } + + //! Access an adjacent Reactor or Reservoir + //! @since New in %Cantera 4.0. + shared_ptr adjacent(size_t n) { + return m_reactors.at(n); + } + //! Set the surface coverages. Array `cov` has length equal to the number of //! surface species. void setCoverages(const double* cov); @@ -93,8 +105,10 @@ class ReactorSurface : public ReactorBase void getState(double* y) override; void initialize(double t0=0.0) override; + vector initializeSteady() override; void updateState(double* y) override; void eval(double t, double* LHS, double* RHS) override; + void evalSteady(double t, double* LHS, double* RHS) override; void addSensitivityReaction(size_t rxn) override; // void addSensitivitySpeciesEnthalpy(size_t k) override; @@ -104,10 +118,9 @@ class ReactorSurface : public ReactorBase size_t componentIndex(const string& nm) const override; string componentName(size_t k) override; - // vector steadyConstraints() const override; - // double upperBound(size_t k) const override; - // double lowerBound(size_t k) const override; - // void resetBadValues(double* y) override; + double upperBound(size_t k) const override; + double lowerBound(size_t k) const override; + void resetBadValues(double* y) override; // void setDerivativeSettings(AnyMap& settings) override; protected: @@ -115,7 +128,7 @@ class ReactorSurface : public ReactorBase shared_ptr m_surf; shared_ptr m_kinetics; - vector m_reactors; + vector> m_reactors; }; //! A surface where the state variables are the total number of moles of each species. @@ -135,6 +148,10 @@ class MoleReactorSurface : public ReactorSurface void getState(double* y) override; void updateState(double* y) override; void eval(double t, double* LHS, double* RHS) override; + void evalSteady(double t, double* LHS, double* RHS) override; + double upperBound(size_t k) const override; + double lowerBound(size_t k) const override; + void resetBadValues(double* y) override; void getJacobianElements(vector>& trips) override; protected: @@ -174,6 +191,8 @@ class FlowReactorSurface : public ReactorSurface return "FlowReactorSurface"; } + bool timeIsIndependent() const override { return false; } + void evalDae(double t, double* y, double* ydot, double* residual) override; void getStateDae(double* y, double* ydot) override; void getConstraints(double* constraints) override; diff --git a/interfaces/cython/cantera/kinetics.pxd b/interfaces/cython/cantera/kinetics.pxd index a5d38f415e..300d745521 100644 --- a/interfaces/cython/cantera/kinetics.pxd +++ b/interfaces/cython/cantera/kinetics.pxd @@ -75,7 +75,7 @@ cdef extern from "cantera/kinetics/Kinetics.h" namespace "Cantera": cdef extern from "cantera/kinetics/InterfaceKinetics.h": cdef cppclass CxxInterfaceKinetics "Cantera::InterfaceKinetics": void advanceCoverages(double, double, double, double, size_t, size_t) except +translate_exception - void solvePseudoSteadyStateProblem() except +translate_exception + void solvePseudoSteadyStateProblem(int) except +translate_exception double interfaceCurrent(size_t) except +translate_exception diff --git a/interfaces/cython/cantera/kinetics.pyi b/interfaces/cython/cantera/kinetics.pyi index 2bb77d6800..6eea11575d 100644 --- a/interfaces/cython/cantera/kinetics.pyi +++ b/interfaces/cython/cantera/kinetics.pyi @@ -173,7 +173,7 @@ class InterfaceKinetics(Kinetics): max_steps: int = 20000, max_error_test_failures: int = 7, ) -> None: ... - def advance_coverages_to_steady_state(self) -> None: ... + def advance_coverages_to_steady_state(self, loglevel: int = 0) -> None: ... def phase_index(self, phase: ThermoPhase | str | int) -> int: ... def _phase_slice(self, phase: ThermoPhase | str | int) -> slice[int, int, None]: ... def get_creation_rates(self, phase: ThermoPhase | str | int) -> Array: ... diff --git a/interfaces/cython/cantera/kinetics.pyx b/interfaces/cython/cantera/kinetics.pyx index ed160988c8..c70a1ceb96 100644 --- a/interfaces/cython/cantera/kinetics.pyx +++ b/interfaces/cython/cantera/kinetics.pyx @@ -943,11 +943,11 @@ cdef class InterfaceKinetics(Kinetics): (self.kinetics).advanceCoverages( dt, rtol, atol, max_step_size, max_steps, max_error_test_failures) - def advance_coverages_to_steady_state(self): + def advance_coverages_to_steady_state(self, loglevel=0): """ This method advances the surface coverages to steady state. """ - (self.kinetics).solvePseudoSteadyStateProblem() + (self.kinetics).solvePseudoSteadyStateProblem(loglevel) def phase_index(self, phase): """ diff --git a/interfaces/cython/cantera/reactor.pxd b/interfaces/cython/cantera/reactor.pxd index db8c9ea3fd..e6ea70d40e 100644 --- a/interfaces/cython/cantera/reactor.pxd +++ b/interfaces/cython/cantera/reactor.pxd @@ -46,7 +46,7 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera": size_t componentIndex(string&) except +translate_exception string componentName(size_t) except +translate_exception void getState(double*) except +translate_exception - void setInitialVolume(double) + void setInitialVolume(double) except +translate_exception void addSensitivityReaction(size_t) except +translate_exception cdef cppclass CxxReactor "Cantera::Reactor" (CxxReactorBase): diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index 45e7b5e635..d3b3003507 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/interfaces/cython/cantera/reactor.pyx @@ -1586,7 +1586,7 @@ cdef class ReactorNet: def __init__(self, reactors=()): self._reactors = [] # prevents premature garbage collection cdef vector[shared_ptr[CxxReactorBase]] cxx_reactors - cdef Reactor r + cdef ReactorBase r for r in reactors: self._reactors.append(r) cxx_reactors.push_back(r._rbase) diff --git a/samples/python/kinetics/diamond_cvd.py b/samples/python/kinetics/diamond_cvd.py index b366fdfca4..7900216232 100644 --- a/samples/python/kinetics/diamond_cvd.py +++ b/samples/python/kinetics/diamond_cvd.py @@ -50,7 +50,8 @@ for n in range(20): x[ih] /= 1.4 g.TPX = t, p, x - d.advance_coverages(10.0) # integrate the coverages to steady state + d.TP = t, p + d.advance_coverages_to_steady_state() carbon_dot = d.net_production_rates[iC] mdot = mw * carbon_dot rate = mdot / dbulk.density diff --git a/samples/python/kinetics/sofc.py b/samples/python/kinetics/sofc.py index aaca13713f..3a25c80fa3 100644 --- a/samples/python/kinetics/sofc.py +++ b/samples/python/kinetics/sofc.py @@ -41,10 +41,6 @@ anode_gas_X = 'H2:0.97, H2O:0.03' cathode_gas_X = 'O2:1.0, H2O:0.001' -# time to integrate coverage eqs. to steady state in -# 'advance_coverages'. This should be more than enough time. -tss = 50.0 - sigma = 2.0 # electrolyte conductivity [Siemens / m] ethick = 5.0e-5 # electrolyte thickness [m] TPB_length_per_area = 1.0e7 # TPB length per unit area [1/m] @@ -173,7 +169,7 @@ def cathode_curr(E): # complex model is required. But as long as the thermal chemistry is fast # relative to charge transfer, this should be an OK approximation. for s in [anode_surf, oxide_surf_a, cathode_surf, oxide_surf_c]: - s.advance_coverages(tss) + s.advance_coverages_to_steady_state() show_coverages(s) # %% diff --git a/samples/python/onedim/catalytic_combustion.py b/samples/python/onedim/catalytic_combustion.py index a6ed9aa771..b2f59b51c9 100644 --- a/samples/python/onedim/catalytic_combustion.py +++ b/samples/python/onedim/catalytic_combustion.py @@ -63,10 +63,10 @@ gas = surf_phase.adjacent["gas"] gas.TPX = tinlet, p, comp1 -# integrate the coverage equations in time for 1 s, holding the gas -# composition fixed to generate a good starting estimate for the coverages. +# integrate the coverage equations holding the gas composition fixed to generate a good +# starting estimate for the coverages. surf_phase.coverages = {'PT(S)': 0.5, 'O(S)':0.5} -surf_phase.advance_coverages(1.0) +surf_phase.advance_coverages_to_steady_state() # create the object that simulates the stagnation flow, and specify an initial # grid diff --git a/samples/python/reactors/1D_packed_bed.py b/samples/python/reactors/1D_packed_bed.py index a11df8b197..482d22259f 100644 --- a/samples/python/reactors/1D_packed_bed.py +++ b/samples/python/reactors/1D_packed_bed.py @@ -212,11 +212,9 @@ rhou0 = gas.density * v_in # Initial surface coverages -# advancing coverages over a long period of time to get the steady state. -surf.advance_coverages(1e10) +surf.advance_coverages_to_steady_state() Zk_0 = surf.coverages - # %% # Define residual function required for IDA solver # ------------------------------------------------ diff --git a/samples/python/reactors/combustor.py b/samples/python/reactors/combustor.py index f485621a4a..49750845cc 100644 --- a/samples/python/reactors/combustor.py +++ b/samples/python/reactors/combustor.py @@ -57,6 +57,8 @@ def mdot(t): return combustor.mass / residence_time +residence_time = 0.1 # starting residence time + inlet_mfc = ct.MassFlowController(inlet, combustor, mdot=mdot) # %% @@ -75,7 +77,6 @@ def mdot(t): states = ct.SolutionArray(gas, extra=['tres']) -residence_time = 0.1 # starting residence time while combustor.T > 500: sim.initial_time = 0.0 # reset the integrator sim.solve_steady() diff --git a/samples/python/reactors/continuous_reactor.py b/samples/python/reactors/continuous_reactor.py index 3a03499900..a44b51a7e4 100644 --- a/samples/python/reactors/continuous_reactor.py +++ b/samples/python/reactors/continuous_reactor.py @@ -117,7 +117,7 @@ upstream=stirred_reactor, downstream=exhaust, primary=mass_flow_controller, - K=1e-3, + K=1e-6, ) reactor_network = ct.ReactorNet([stirred_reactor]) @@ -219,28 +219,14 @@ reactor_X = inlet_X for reactor_temperature in T: - gas.TPX = reactor_temperature, reactor_pressure, inlet_X - fuel_air_mixture_tank = ct.Reservoir(gas) - # Use composition from the previous iteration to speed up convergence - gas.TPX = reactor_temperature, reactor_pressure, reactor_X - stirred_reactor = ct.IdealGasReactor(gas, energy="off", volume=reactor_volume) - fuel_air_mixture_tank = ct.Reservoir(gas) - stirred_reactor = ct.IdealGasReactor(gas, energy="off", volume=reactor_volume) - mass_flow_controller = ct.MassFlowController( - upstream=fuel_air_mixture_tank, - downstream=stirred_reactor, - mdot=lambda t: stirred_reactor.mass / residence_time, - ) - pressure_regulator = ct.PressureController( - upstream=stirred_reactor, downstream=exhaust, primary=mass_flow_controller, - K=1e-3, - ) - reactor_network = ct.ReactorNet([stirred_reactor]) + fuel_air_mixture_tank.phase.TPX = reactor_temperature, reactor_pressure, inlet_X + stirred_reactor.phase.TPX = reactor_temperature, reactor_pressure, reactor_X # Re-run the isothermal simulations tic = time.time() counter = 0 + reactor_network.initial_time = 0.0 while reactor_network.time < max_simulation_time: reactor_network.step() counter += 1 diff --git a/samples/python/reactors/pfr.py b/samples/python/reactors/pfr.py index 92a7009dbf..d15d48f18c 100644 --- a/samples/python/reactors/pfr.py +++ b/samples/python/reactors/pfr.py @@ -112,7 +112,7 @@ # We need an outlet to the downstream reservoir. This will determine the # pressure in the reactor. The value of K will only affect the transient # pressure difference. -v = ct.PressureController(r2, downstream, primary=m, K=1e-5) +v = ct.PressureController(r2, downstream, primary=m, K=1e-12) sim2 = ct.ReactorNet([r2]) sim2.max_time_step = 1e4 @@ -129,7 +129,7 @@ upstream.phase.TDY = r2.phase.TDY # integrate the reactor forward in time until steady state is reached sim2.reinitialize() - sim2.advance_to_steady_state() + sim2.solve_steady() # compute velocity and transform into time u2[n] = mass_flow_rate2 / area / r2.phase.density t_r2[n] = r2.mass / mass_flow_rate2 # residence time in this reactor diff --git a/samples/python/reactors/surf_pfr_chain.py b/samples/python/reactors/surf_pfr_chain.py index 33be1765a9..ebe655b0b2 100644 --- a/samples/python/reactors/surf_pfr_chain.py +++ b/samples/python/reactors/surf_pfr_chain.py @@ -98,7 +98,7 @@ # We need an outlet to the downstream reservoir. This will determine the # pressure in the reactor. The value of K will only affect the transient # pressure difference. -v = ct.PressureController(r, downstream, primary=m, K=1e-5) +v = ct.PressureController(r, downstream, primary=m, K=1e-6) sim = ct.ReactorNet([r]) @@ -108,7 +108,7 @@ # Set the state of the reservoir to match that of the previous reactor upstream.phase.TDY = r.phase.TDY sim.reinitialize() - sim.advance_to_steady_state() + sim.solve_steady() dist = n * rlen * 1.0e3 # distance in mm if n % 10 == 0: diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp index 604a45901f..19859599f2 100644 --- a/src/equil/vcs_solve_TP.cpp +++ b/src/equil/vcs_solve_TP.cpp @@ -2688,6 +2688,7 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc, writelogendl(); } + m_TmpPhase.assign(m_numPhases, 0.0); double* tlogMoles = &m_TmpPhase[0]; // Might as well recalculate the phase mole vector and compare to the stored diff --git a/src/kinetics/ImplicitSurfChem.cpp b/src/kinetics/ImplicitSurfChem.cpp deleted file mode 100644 index 9630738995..0000000000 --- a/src/kinetics/ImplicitSurfChem.cpp +++ /dev/null @@ -1,296 +0,0 @@ -/** - * @file ImplicitSurfChem.cpp - * Definitions for the implicit integration of surface site density equations - * (see @ref kineticsmgr and class - * @link Cantera::ImplicitSurfChem ImplicitSurfChem@endlink). - */ - -// This file is part of Cantera. See License.txt in the top-level directory or -// at https://cantera.org/license.txt for license and copyright information. - -#include "cantera/kinetics/ImplicitSurfChem.h" -#include "cantera/kinetics/solveSP.h" -#include "cantera/thermo/SurfPhase.h" - -namespace Cantera -{ - -ImplicitSurfChem::ImplicitSurfChem( - vector k, double rtol, double atol, - double maxStepSize, size_t maxSteps, - size_t maxErrTestFails) : - m_atol(atol), - m_rtol(rtol), - m_maxstep(maxStepSize), - m_nmax(maxSteps), - m_maxErrTestFails(maxErrTestFails) -{ - size_t ntmax = 0; - size_t kinSpIndex = 0; - // Loop over the number of surface kinetics objects - for (size_t n = 0; n < k.size(); n++) { - InterfaceKinetics* kinPtr = k[n]; - m_vecKinPtrs.push_back(kinPtr); - SurfPhase* surf = dynamic_cast(&k[n]->thermo(0)); - if (surf == nullptr) { - throw CanteraError("ImplicitSurfChem::ImplicitSurfChem", - "kinetics manager contains no surface phase"); - } - m_surf.push_back(surf); - size_t nsp = m_surf.back()->nSpecies(); - m_nsp.push_back(nsp); - m_nv += m_nsp.back(); - size_t nt = k[n]->nTotalSpecies(); - ntmax = std::max(nt, ntmax); - m_specStartIndex.push_back(kinSpIndex); - kinSpIndex += nsp; - size_t nPhases = kinPtr->nPhases(); - vector pLocTmp(nPhases); - pLocTmp[0] = -int(n); - size_t imatch = npos; - for (size_t ip = 1; ip < nPhases; ip++) { - ThermoPhase* thPtr = & kinPtr->thermo(ip); - if ((imatch = checkMatch(m_bulkPhases, thPtr)) == npos) { - m_bulkPhases.push_back(thPtr); - nsp = thPtr->nSpecies(); - m_numTotalBulkSpecies += nsp; - imatch = m_bulkPhases.size() - 1; - } - pLocTmp[ip] = int(imatch); - } - pLocVec.push_back(pLocTmp); - } - m_numTotalSpecies = m_nv + m_numTotalBulkSpecies; - m_concSpecies.resize(m_numTotalSpecies, 0.0); - m_concSpeciesSave.resize(m_numTotalSpecies, 0.0); - - m_integ.reset(newIntegrator("CVODE")); - - // use backward differencing, with a full Jacobian computed - // numerically, and use a Newton linear iterator - m_integ->setMethod(BDF_Method); - m_integ->setLinearSolverType("DENSE"); - m_work.resize(ntmax); -} - -int ImplicitSurfChem::checkMatch(vector m_vec, ThermoPhase* thPtr) -{ - int retn = -1; - for (int i = 0; i < (int) m_vec.size(); i++) { - ThermoPhase* th = m_vec[i]; - if (th == thPtr) { - return i; - } - } - return retn; -} - -void ImplicitSurfChem::getState(double* c) -{ - size_t loc = 0; - for (size_t n = 0; n < m_surf.size(); n++) { - m_surf[n]->getCoverages(c + loc); - loc += m_nsp[n]; - } -} - -void ImplicitSurfChem::setMaxStepSize(double maxstep) -{ - m_maxstep = maxstep; - if (m_maxstep > 0) { - m_integ->setMaxStepSize(m_maxstep); - } -} - -void ImplicitSurfChem::setTolerances(double rtol, double atol) -{ - m_rtol = rtol; - m_atol = atol; - m_integ->setTolerances(m_rtol, m_atol); -} - -void ImplicitSurfChem::setMaxSteps(size_t maxsteps) -{ - m_nmax = maxsteps; - m_integ->setMaxSteps(static_cast(m_nmax)); -} - -void ImplicitSurfChem::setMaxErrTestFails(size_t maxErrTestFails) -{ - m_maxErrTestFails = maxErrTestFails; - m_integ->setMaxErrTestFails(static_cast(m_maxErrTestFails)); -} - -void ImplicitSurfChem::initialize(double t0) -{ - this->setTolerances(m_rtol, m_atol); - this->setMaxStepSize(m_maxstep); - this->setMaxSteps(m_nmax); - this->setMaxErrTestFails(m_maxErrTestFails); - m_integ->initialize(t0, *this); -} - -void ImplicitSurfChem::integrate(double t0, double t1) -{ - this->initialize(t0); - if (fabs(t1 - t0) < m_maxstep || m_maxstep == 0) { - // limit max step size on this run to t1 - t0 - m_integ->setMaxStepSize(t1 - t0); - } - m_integ->integrate(t1); - updateState(m_integ->solution()); -} - -void ImplicitSurfChem::integrate0(double t0, double t1) -{ - m_integ->integrate(t1); - updateState(m_integ->solution()); -} - -void ImplicitSurfChem::updateState(double* c) -{ - size_t loc = 0; - for (size_t n = 0; n < m_surf.size(); n++) { - m_surf[n]->setCoverages(c + loc); - loc += m_nsp[n]; - } -} - -void ImplicitSurfChem::eval(double time, double* y, double* ydot, double* p) -{ - updateState(y); // synchronize the surface state(s) with y - size_t loc = 0; - for (size_t n = 0; n < m_surf.size(); n++) { - double rs0 = 1.0/m_surf[n]->siteDensity(); - m_vecKinPtrs[n]->getNetProductionRates(m_work.data()); - double sum = 0.0; - for (size_t k = 1; k < m_nsp[n]; k++) { - ydot[k + loc] = m_work[k] * rs0 * m_surf[n]->size(k); - sum -= ydot[k]; - } - ydot[loc] = sum; - loc += m_nsp[n]; - } -} - -void ImplicitSurfChem::solvePseudoSteadyStateProblem(int ifuncOverride, - double timeScaleOverride) -{ - int ifunc; - // set bulkFunc. We assume that the bulk concentrations are constant. - int bulkFunc = BULK_ETCH; - // time scale - time over which to integrate equations - double time_scale = timeScaleOverride; - if (!m_surfSolver) { - m_surfSolver = make_unique(this, bulkFunc); - // set ifunc, which sets the algorithm. - ifunc = SFLUX_INITIALIZE; - } else { - ifunc = SFLUX_RESIDUAL; - } - - // Possibly override the ifunc value - if (ifuncOverride >= 0) { - ifunc = ifuncOverride; - } - - // Get the specifications for the problem from the values - // in the ThermoPhase objects for all phases. - // - // 1) concentrations of all species in all phases, m_concSpecies[] - // 2) Temperature and pressure - getConcSpecies(m_concSpecies.data()); - InterfaceKinetics* ik = m_vecKinPtrs[0]; - ThermoPhase& tp = ik->thermo(0); - double TKelvin = tp.temperature(); - double PGas = tp.pressure(); - - // Make sure that there is a common temperature and pressure for all - // ThermoPhase objects belonging to the interfacial kinetics object, if it - // is required by the problem statement. - if (m_commonTempPressForPhases) { - setCommonState_TP(TKelvin, PGas); - } - - double reltol = 1.0E-6; - double atol = 1.0E-20; - - // Install a filter for negative concentrations. One of the few ways solveSP - // can fail is if concentrations on input are below zero. - bool rset = false; - for (size_t k = 0; k < m_nv; k++) { - if (m_concSpecies[k] < 0.0) { - rset = true; - m_concSpecies[k] = 0.0; - } - } - if (rset) { - setConcSpecies(m_concSpecies.data()); - } - - m_surfSolver->m_ioflag = m_ioFlag; - - // Save the current solution - m_concSpeciesSave = m_concSpecies; - - int retn = m_surfSolver->solveSurfProb(ifunc, time_scale, TKelvin, PGas, - reltol, atol); - if (retn != 1) { - // reset the concentrations - m_concSpecies = m_concSpeciesSave; - setConcSpecies(m_concSpeciesSave.data()); - ifunc = SFLUX_INITIALIZE; - retn = m_surfSolver->solveSurfProb(ifunc, time_scale, TKelvin, PGas, - reltol, atol); - if (retn != 1) { - throw CanteraError("ImplicitSurfChem::solvePseudoSteadyStateProblem", - "solveSP return an error condition!"); - } - } -} - -void ImplicitSurfChem::getConcSpecies(double* const vecConcSpecies) const -{ - size_t kstart; - for (size_t ip = 0; ip < m_surf.size(); ip++) { - ThermoPhase* TP_ptr = m_surf[ip]; - kstart = m_specStartIndex[ip]; - TP_ptr->getConcentrations(vecConcSpecies + kstart); - } - kstart = m_nv; - for (size_t ip = 0; ip < m_bulkPhases.size(); ip++) { - ThermoPhase* TP_ptr = m_bulkPhases[ip]; - TP_ptr->getConcentrations(vecConcSpecies + kstart); - kstart += TP_ptr->nSpecies(); - } -} - -void ImplicitSurfChem::setConcSpecies(const double* const vecConcSpecies) -{ - size_t kstart; - for (size_t ip = 0; ip < m_surf.size(); ip++) { - ThermoPhase* TP_ptr = m_surf[ip]; - kstart = m_specStartIndex[ip]; - TP_ptr->setConcentrations(vecConcSpecies + kstart); - } - kstart = m_nv; - for (size_t ip = 0; ip < m_bulkPhases.size(); ip++) { - ThermoPhase* TP_ptr = m_bulkPhases[ip]; - TP_ptr->setConcentrations(vecConcSpecies + kstart); - kstart += TP_ptr->nSpecies(); - } -} - -void ImplicitSurfChem::setCommonState_TP(double TKelvin, double PresPa) -{ - for (size_t ip = 0; ip < m_surf.size(); ip++) { - ThermoPhase* TP_ptr = m_surf[ip]; - TP_ptr->setState_TP(TKelvin, PresPa); - } - for (size_t ip = 0; ip < m_bulkPhases.size(); ip++) { - ThermoPhase* TP_ptr = m_bulkPhases[ip]; - TP_ptr->setState_TP(TKelvin, PresPa); - } -} - -} diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index 276bc288f5..7fa2d916a1 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -6,7 +6,10 @@ // at https://cantera.org/license.txt for license and copyright information. #include "cantera/kinetics/InterfaceKinetics.h" -#include "cantera/kinetics/ImplicitSurfChem.h" +#include "cantera/zeroD/ReactorNet.h" +#include "cantera/zeroD/ReactorFactory.h" +#include "cantera/zeroD/ReactorSurface.h" +#include "cantera/zeroD/Reservoir.h" #include "cantera/kinetics/Reaction.h" #include "cantera/thermo/SurfPhase.h" #include "cantera/base/utilities.h" @@ -14,6 +17,8 @@ namespace Cantera { +// Constructor / destructor definitions required due to forward-declared unique_ptr +// members. InterfaceKinetics::~InterfaceKinetics() { delete m_integrator; @@ -471,14 +476,6 @@ void InterfaceKinetics::setMultiplier(size_t i, double f) m_ROP_ok = false; } -void InterfaceKinetics::setIOFlag(int ioFlag) -{ - m_ioFlag = ioFlag; - if (m_integrator) { - m_integrator->setIOFlag(ioFlag); - } -} - void InterfaceKinetics::addThermo(shared_ptr thermo) { Kinetics::addThermo(thermo); @@ -510,34 +507,72 @@ void InterfaceKinetics::resizeSpecies() m_phi.resize(nPhases(), 0.0); } +void InterfaceKinetics::buildNetwork() +{ + for (auto& phase : m_thermo) { + if (!phase->root()) { + throw CanteraError("InterfaceKinetics::buildNetwork", + "Phase '{}' is not attached to a Solution.", phase->name()); + } + } + vector> reservoirs; + for (size_t i = 1; i < nPhases(); i++) { + auto r = newReservoir(thermo(i).root(), false); + reservoirs.push_back(r); + } + auto rsurf = newReactorSurface(thermo(0).root(), reservoirs, false); + m_integrator = new ReactorNet(rsurf); +} + void InterfaceKinetics::advanceCoverages(double tstep, double rtol, double atol, double maxStepSize, size_t maxSteps, size_t maxErrTestFails) { - if (m_integrator == 0) { - vector k{this}; - m_integrator = new ImplicitSurfChem(k); + // Stash the state of adjacent phases, and set their T and P to match the surface + vector> savedStates(nPhases()); + for (size_t i = 1; i < nPhases(); i++) { + savedStates[i].resize(thermo(i).partialStateSize()); + thermo(i).savePartialState(savedStates[i].size(), savedStates[i].data()); + thermo(i).setState_TP(thermo(0).temperature(), thermo(0).pressure()); } + + if (!m_integrator) { + buildNetwork(); + } + m_integrator->setTolerances(rtol, atol); - m_integrator->setMaxStepSize(maxStepSize); + m_integrator->setMaxTimeStep(maxStepSize); m_integrator->setMaxSteps(maxSteps); m_integrator->setMaxErrTestFails(maxErrTestFails); - m_integrator->integrate(0.0, tstep); - delete m_integrator; - m_integrator = 0; + m_integrator->setInitialTime(0.0); + m_integrator->advance(tstep); + + // Restore adjacent phases to their original states + for (size_t i = 1; i < nPhases(); i++) { + thermo(i).restorePartialState(savedStates[i].size(), savedStates[i].data()); + } } -void InterfaceKinetics::solvePseudoSteadyStateProblem( - int ifuncOverride, double timeScaleOverride) +void InterfaceKinetics::solvePseudoSteadyStateProblem(int loglevel) { - // create our own solver object - if (m_integrator == 0) { - vector k{this}; - m_integrator = new ImplicitSurfChem(k); - m_integrator->initialize(); + // Stash the state of adjacent phases, and set their T and P to match the surface + vector> savedStates(nPhases()); + for (size_t i = 1; i < nPhases(); i++) { + savedStates[i].resize(thermo(i).partialStateSize()); + thermo(i).savePartialState(savedStates[i].size(), savedStates[i].data()); + thermo(i).setState_TP(thermo(0).temperature(), thermo(0).pressure()); + } + + if (!m_integrator) { + buildNetwork(); + } + + m_integrator->setVerbose(loglevel != 0); + m_integrator->solveSteady(loglevel); + + // Restore adjacent phases to their original states + for (size_t i = 1; i < nPhases(); i++) { + thermo(i).restorePartialState(savedStates[i].size(), savedStates[i].data()); } - m_integrator->setIOFlag(m_ioFlag); - // New direct method to go here - m_integrator->solvePseudoSteadyStateProblem(ifuncOverride, timeScaleOverride); } void InterfaceKinetics::setPhaseExistence(const size_t iphase, const int exists) diff --git a/src/kinetics/solveSP.cpp b/src/kinetics/solveSP.cpp deleted file mode 100644 index 7499a61f6d..0000000000 --- a/src/kinetics/solveSP.cpp +++ /dev/null @@ -1,692 +0,0 @@ -/** - * @file: solveSP.cpp Implicit surface site concentration solver - */ - -// This file is part of Cantera. See License.txt in the top-level directory or -// at https://cantera.org/license.txt for license and copyright information. - -#include "cantera/kinetics/solveSP.h" -#include "cantera/thermo/SurfPhase.h" -#include "cantera/kinetics/ImplicitSurfChem.h" - -namespace Cantera -{ - -// STATIC ROUTINES DEFINED IN THIS FILE - -static double calc_damping(double* x, double* dx, size_t dim, int*); -static double calcWeightedNorm(const double [], const double dx[], size_t); - -// solveSP Class Definitions - -solveSP::solveSP(ImplicitSurfChem* surfChemPtr, int bulkFunc) : - m_objects(surfChemPtr->getObjects()), - m_bulkFunc(bulkFunc) -{ - for (size_t n = 0; n < m_objects.size(); n++) { - InterfaceKinetics* kin = m_objects[n]; - SurfPhase* sp = dynamic_cast(&kin->thermo(0)); - if (sp == nullptr) { - throw CanteraError("solveSP::solveSP", - "InterfaceKinetics object has no surface phase"); - } - - m_numSurfPhases++; - m_indexKinObjSurfPhase.push_back(n); - - m_ptrsSurfPhase.push_back(sp); - size_t nsp = sp->nSpecies(); - m_nSpeciesSurfPhase.push_back(nsp); - m_numTotSurfSpecies += nsp; - } - - if (bulkFunc == BULK_DEPOSITION) { - m_neq = m_numTotSurfSpecies + m_numTotBulkSpeciesSS; - } else { - m_neq = m_numTotSurfSpecies; - } - - for (size_t n = 0; n < m_numSurfPhases; n++) { - size_t tsp = m_objects[n]->nTotalSpecies(); - m_maxTotSpecies = std::max(m_maxTotSpecies, tsp); - } - m_maxTotSpecies = std::max(m_maxTotSpecies, m_neq); - - m_netProductionRatesSave.resize(m_maxTotSpecies, 0.0); - m_numEqn1.resize(m_maxTotSpecies, 0.0); - m_numEqn2.resize(m_maxTotSpecies, 0.0); - m_XMolKinSpecies.resize(m_maxTotSpecies, 0.0); - m_CSolnSave.resize(m_neq, 0.0); - m_spSurfLarge.resize(m_numSurfPhases, 0); - m_kinSpecIndex.resize(m_numTotSurfSpecies + m_numTotBulkSpeciesSS, 0); - m_kinObjIndex.resize(m_numTotSurfSpecies + m_numTotBulkSpeciesSS, 0); - m_eqnIndexStartSolnPhase.resize(m_numSurfPhases + m_numBulkPhasesSS, 0); - - size_t kindexSP = 0; - for (size_t isp = 0; isp < m_numSurfPhases; isp++) { - size_t iKinObject = m_indexKinObjSurfPhase[isp]; - InterfaceKinetics* kin = m_objects[iKinObject]; - size_t kstart = kin->kineticsSpeciesIndex(0, 0); - size_t nsp = m_nSpeciesSurfPhase[isp]; - m_eqnIndexStartSolnPhase[isp] = kindexSP; - for (size_t k = 0; k < nsp; k++, kindexSP++) { - m_kinSpecIndex[kindexSP] = kstart + k; - m_kinObjIndex[kindexSP] = isp; - } - } - - // Dimension solution vector - size_t dim1 = std::max(1, m_neq); - m_CSolnSP.resize(dim1, 0.0); - m_CSolnSPInit.resize(dim1, 0.0); - m_CSolnSPOld.resize(dim1, 0.0); - m_wtResid.resize(dim1, 0.0); - m_wtSpecies.resize(dim1, 0.0); - m_resid.resize(dim1, 0.0); - m_Jac.resize(dim1, dim1, 0.0); -} - -int solveSP::solveSurfProb(int ifunc, double time_scale, double TKelvin, - double PGas, double reltol, double abstol) -{ - double EXTRA_ACCURACY = 0.001; - if (ifunc == SFLUX_JACOBIAN) { - EXTRA_ACCURACY *= 0.001; - } - int label_t=-1; // Species IDs for time control - int label_d = -1; // Species IDs for damping control - int label_t_old=-1; - double label_factor = 1.0; - int iter=0; // iteration number on nonlinear solver - int iter_max=1000; // maximum number of nonlinear iterations - double deltaT = 1.0E-10; // Delta time step - double damp=1.0; - double inv_t = 0.0; - double t_real = 0.0, update_norm = 1.0E6; - bool do_time = false, not_converged = true; - m_ioflag = std::min(m_ioflag, 1); - - // Set the initial value of the do_time parameter - if (ifunc == SFLUX_INITIALIZE || ifunc == SFLUX_TRANSIENT) { - do_time = true; - } - - // Store the initial guess for the surface problem in the soln vector, - // CSoln, and in an separate vector CSolnInit. - size_t loc = 0; - for (size_t n = 0; n < m_numSurfPhases; n++) { - m_ptrsSurfPhase[n]->getConcentrations(m_numEqn1.data()); - for (size_t k = 0; k < m_nSpeciesSurfPhase[n]; k++) { - m_CSolnSP[loc] = m_numEqn1[k]; - loc++; - } - } - - m_CSolnSPInit = m_CSolnSP; - - // Calculate the largest species in each phase - evalSurfLarge(m_CSolnSP.data()); - - if (m_ioflag) { - print_header(m_ioflag, ifunc, time_scale, true, reltol, abstol); - } - - // Quick return when there isn't a surface problem to solve - if (m_neq == 0) { - not_converged = false; - update_norm = 0.0; - } - - // Start of Newton's method - while (not_converged && iter < iter_max) { - iter++; - // Store previous iteration's solution in the old solution vector - m_CSolnSPOld = m_CSolnSP; - - // Evaluate the largest surface species for each surface phase every - // 5 iterations. - if (iter%5 == 4) { - evalSurfLarge(m_CSolnSP.data()); - } - - // Calculate the value of the time step - // - heuristics to stop large oscillations in deltaT - if (do_time) { - // don't hurry increase in time step at the same time as damping - if (damp < 1.0) { - label_factor = 1.0; - } - double tmp = calc_t(m_netProductionRatesSave.data(), - m_XMolKinSpecies.data(), - &label_t, &label_t_old, &label_factor, m_ioflag); - if (iter < 10) { - inv_t = tmp; - } else if (tmp > 2.0*inv_t) { - inv_t = 2.0*inv_t; - } else { - inv_t = tmp; - } - - // Check end condition - if (ifunc == SFLUX_TRANSIENT) { - tmp = t_real + 1.0/inv_t; - if (tmp > time_scale) { - inv_t = 1.0/(time_scale - t_real); - } - } - } else { - // make steady state calc a step of 1 million seconds to prevent - // singular Jacobians for some pathological cases - inv_t = 1.0e-6; - } - deltaT = 1.0/inv_t; - - // Call the routine to numerically evaluation the Jacobian and residual - // for the current iteration. - resjac_eval(m_Jac, m_resid.data(), m_CSolnSP.data(), - m_CSolnSPOld.data(), do_time, deltaT); - - // Calculate the weights. Make sure the calculation is carried out on - // the first iteration. - if (iter%4 == 1) { - calcWeights(m_wtSpecies.data(), m_wtResid.data(), - m_Jac, m_CSolnSP.data(), abstol, reltol); - } - - // Find the weighted norm of the residual - double resid_norm = calcWeightedNorm(m_wtResid.data(), m_resid.data(), m_neq); - - // Solve Linear system. The solution is in m_resid - solve(m_Jac, m_resid.data()); - - // Calculate the Damping factor needed to keep all unknowns between 0 - // and 1, and not allow too large a change (factor of 2) in any unknown. - damp = calc_damping(m_CSolnSP.data(), m_resid.data(), m_neq, &label_d); - - // Calculate the weighted norm of the update vector Here, resid is the - // delta of the solution, in concentration units. - update_norm = calcWeightedNorm(m_wtSpecies.data(), - m_resid.data(), m_neq); - - // Update the solution vector and real time Crop the concentrations to - // zero. - for (size_t irow = 0; irow < m_neq; irow++) { - m_CSolnSP[irow] -= damp * m_resid[irow]; - } - for (size_t irow = 0; irow < m_neq; irow++) { - m_CSolnSP[irow] = std::max(0.0, m_CSolnSP[irow]); - } - updateState(m_CSolnSP.data()); - - if (do_time) { - t_real += damp/inv_t; - } - - if (m_ioflag) { - printIteration(m_ioflag, damp, label_d, label_t, inv_t, t_real, iter, - update_norm, resid_norm, do_time); - } - - if (ifunc == SFLUX_TRANSIENT) { - not_converged = (t_real < time_scale); - } else { - if (do_time) { - if (t_real > time_scale || - (resid_norm < 1.0e-7 && - update_norm*time_scale/t_real < EXTRA_ACCURACY)) { - do_time = false; - } - } else { - not_converged = ((update_norm > EXTRA_ACCURACY) || - (resid_norm > EXTRA_ACCURACY)); - } - } - } // End of Newton's Method while statement - - // End Newton's method. If not converged, print error message and - // recalculate sdot's at equal site fractions. - if (not_converged && m_ioflag) { - writelog("#$#$#$# Error in solveSP $#$#$#$ \n"); - writelogf("Newton iter on surface species did not converge, " - "update_norm = %e \n", update_norm); - writelog("Continuing anyway\n"); - } - - // Decide on what to return in the solution vector. Right now, will always - // return the last solution no matter how bad - if (m_ioflag) { - fun_eval(m_resid.data(), m_CSolnSP.data(), m_CSolnSPOld.data(), - false, deltaT); - double resid_norm = calcWeightedNorm(m_wtResid.data(), m_resid.data(), m_neq); - printIteration(m_ioflag, damp, label_d, label_t, inv_t, t_real, iter, - update_norm, resid_norm, do_time, true); - } - - // Return with the appropriate flag - if (update_norm > 1.0) { - return -1; - } - return 1; -} - -void solveSP::updateState(const double* CSolnSP) -{ - vector X; - size_t loc = 0; - for (size_t n = 0; n < m_numSurfPhases; n++) { - X.resize(m_nSpeciesSurfPhase[n]); - for (size_t k = 0; k < X.size(); k++) { - X[k] = CSolnSP[loc + k] / m_ptrsSurfPhase[n]->siteDensity(); - } - m_ptrsSurfPhase[n]->setMoleFractions_NoNorm(X.data()); - loc += m_nSpeciesSurfPhase[n]; - } -} - -void solveSP::updateMFSolnSP(double* XMolSolnSP) -{ - for (size_t isp = 0; isp < m_numSurfPhases; isp++) { - size_t keqnStart = m_eqnIndexStartSolnPhase[isp]; - m_ptrsSurfPhase[isp]->getMoleFractions(XMolSolnSP + keqnStart); - } -} - -void solveSP::updateMFKinSpecies(double* XMolKinSpecies, int isp) -{ - InterfaceKinetics* kin = m_objects[isp]; - for (size_t iph = 0; iph < kin->nPhases(); iph++) { - size_t ksi = kin->kineticsSpeciesIndex(0, iph); - kin->thermo(iph).getMoleFractions(XMolKinSpecies + ksi); - } -} - -void solveSP::evalSurfLarge(const double* CSolnSP) -{ - size_t kindexSP = 0; - for (size_t isp = 0; isp < m_numSurfPhases; isp++) { - double Clarge = CSolnSP[kindexSP]; - m_spSurfLarge[isp] = 0; - kindexSP++; - for (size_t k = 1; k < m_nSpeciesSurfPhase[isp]; k++, kindexSP++) { - if (CSolnSP[kindexSP] > Clarge) { - Clarge = CSolnSP[kindexSP]; - m_spSurfLarge[isp] = k; - } - } - } -} - -void solveSP::fun_eval(double* resid, const double* CSoln, const double* CSolnOld, - const bool do_time, const double deltaT) -{ - size_t k; - double lenScale = 1.0E-9; - if (m_numSurfPhases > 0) { - // update the surface concentrations with the input surface - // concentration vector - updateState(CSoln); - - // Get the net production rates of all of the species in the - // surface kinetics mechanism - // - // HKM Should do it here for all kinetics objects so that - // bulk will eventually work. - if (do_time) { - size_t kindexSP = 0; - for (size_t isp = 0; isp < m_numSurfPhases; isp++) { - size_t nsp = m_nSpeciesSurfPhase[isp]; - InterfaceKinetics* kinPtr = m_objects[isp]; - size_t kins = kindexSP; - kinPtr->getNetProductionRates(m_netProductionRatesSave.data()); - for (k = 0; k < nsp; k++, kindexSP++) { - resid[kindexSP] = - (CSoln[kindexSP] - CSolnOld[kindexSP]) / deltaT - - m_netProductionRatesSave[k]; - } - - size_t kspecial = kins + m_spSurfLarge[isp]; - double sd = m_ptrsSurfPhase[isp]->siteDensity(); - resid[kspecial] = sd; - for (k = 0; k < nsp; k++) { - resid[kspecial] -= CSoln[kins + k]; - } - } - } else { - size_t kindexSP = 0; - for (size_t isp = 0; isp < m_numSurfPhases; isp++) { - size_t nsp = m_nSpeciesSurfPhase[isp]; - InterfaceKinetics* kinPtr = m_objects[isp]; - size_t kins = kindexSP; - kinPtr->getNetProductionRates(m_netProductionRatesSave.data()); - for (k = 0; k < nsp; k++, kindexSP++) { - resid[kindexSP] = - m_netProductionRatesSave[k]; - } - size_t kspecial = kins + m_spSurfLarge[isp]; - double sd = m_ptrsSurfPhase[isp]->siteDensity(); - resid[kspecial] = sd; - for (k = 0; k < nsp; k++) { - resid[kspecial] -= CSoln[kins + k]; - } - } - } - - if (m_bulkFunc == BULK_DEPOSITION) { - size_t kindexSP = m_numTotSurfSpecies; - for (size_t isp = 0; isp < m_numBulkPhasesSS; isp++) { - double* XBlk = m_numEqn1.data(); - size_t nsp = m_nSpeciesSurfPhase[isp]; - double grRate = 0.0; - for (k = 0; k < nsp; k++) { - if (m_netProductionRatesSave[k] > 0.0) { - grRate += m_netProductionRatesSave[k]; - } - } - resid[kindexSP] = m_bulkPhasePtrs[isp]->molarDensity(); - for (k = 0; k < nsp; k++) { - resid[kindexSP] -= CSoln[kindexSP + k]; - } - if (grRate > 0.0) { - for (k = 1; k < nsp; k++) { - if (m_netProductionRatesSave[k] > 0.0) { - resid[kindexSP + k] = XBlk[k] * grRate - - m_netProductionRatesSave[k]; - } else { - resid[kindexSP + k] = XBlk[k] * grRate; - } - } - } else { - grRate = 1.0E-6; - //! @todo the appearance of k in this formula is suspicious - grRate += fabs(m_netProductionRatesSave[k]); - for (k = 1; k < nsp; k++) { - resid[kindexSP + k] = grRate * (XBlk[k] - 1.0/nsp); - } - } - if (do_time) { - for (k = 1; k < nsp; k++) { - resid[kindexSP + k] += - lenScale / deltaT * - (CSoln[kindexSP + k]- CSolnOld[kindexSP + k]); - } - } - kindexSP += nsp; - } - } - } -} - -void solveSP::resjac_eval(DenseMatrix& jac, double resid[], double CSoln[], - const double CSolnOld[], const bool do_time, - const double deltaT) -{ - size_t kColIndex = 0; - // Calculate the residual - fun_eval(resid, CSoln, CSolnOld, do_time, deltaT); - // Now we will look over the columns perturbing each unknown. - for (size_t jsp = 0; jsp < m_numSurfPhases; jsp++) { - size_t nsp = m_nSpeciesSurfPhase[jsp]; - double sd = m_ptrsSurfPhase[jsp]->siteDensity(); - for (size_t kCol = 0; kCol < nsp; kCol++) { - double cSave = CSoln[kColIndex]; - double dc = std::max(1.0E-10 * sd, fabs(cSave) * 1.0E-7); - CSoln[kColIndex] += dc; - fun_eval(m_numEqn2.data(), CSoln, CSolnOld, do_time, deltaT); - for (size_t i = 0; i < m_neq; i++) { - jac(i, kColIndex) = (m_numEqn2[i] - resid[i])/dc; - } - CSoln[kColIndex] = cSave; - kColIndex++; - } - } - - if (m_bulkFunc == BULK_DEPOSITION) { - for (size_t jsp = 0; jsp < m_numBulkPhasesSS; jsp++) { - size_t nsp = m_numBulkSpecies[jsp]; - double sd = m_bulkPhasePtrs[jsp]->molarDensity(); - for (size_t kCol = 0; kCol < nsp; kCol++) { - double cSave = CSoln[kColIndex]; - double dc = std::max(1.0E-10 * sd, fabs(cSave) * 1.0E-7); - CSoln[kColIndex] += dc; - fun_eval(m_numEqn2.data(), CSoln, CSolnOld, do_time, deltaT); - for (size_t i = 0; i < m_neq; i++) { - jac(i, kColIndex) = (m_numEqn2[i] - resid[i])/dc; - } - CSoln[kColIndex] = cSave; - kColIndex++; - } - } - } -} - -/** - * This function calculates a damping factor for the Newton iteration update - * vector, dxneg, to insure that all site and bulk fractions, x, remain - * bounded between zero and one. - * - * dxneg[] = negative of the update vector. - * - * The constant "APPROACH" sets the fraction of the distance to the boundary - * that the step can take. If the full step would not force any fraction - * outside of 0-1, then Newton's method is allowed to operate normally. - */ -static double calc_damping(double x[], double dxneg[], size_t dim, int* label) -{ - const double APPROACH = 0.80; - double damp = 1.0; - static double damp_old = 1.0; //! @todo this variable breaks thread safety - *label = -1; - - for (size_t i = 0; i < dim; i++) { - // Calculate the new suggested new value of x[i] - double xnew = x[i] - damp * dxneg[i]; - - // Calculate the allowed maximum and minimum values of x[i] - // - Only going to allow x[i] to converge to zero by a - // single order of magnitude at a time - double xtop = 1.0 - 0.1*fabs(1.0-x[i]); - double xbot = fabs(x[i]*0.1) - 1.0e-16; - if (xnew > xtop) { - damp = - APPROACH * (1.0 - x[i]) / dxneg[i]; - *label = int(i); - } else if (xnew < xbot) { - damp = APPROACH * x[i] / dxneg[i]; - *label = int(i); - } else if (xnew > 3.0*std::max(x[i], 1.0E-10)) { - damp = - 2.0 * std::max(x[i], 1.0E-10) / dxneg[i]; - *label = int(i); - } - } - damp = std::max(damp, 1e-2); - - // Only allow the damping parameter to increase by a factor of three each - // iteration. Heuristic to avoid oscillations in the value of damp - if (damp > damp_old*3) { - damp = damp_old*3; - *label = -1; - } - - // Save old value of the damping parameter for use in subsequent calls. - damp_old = damp; - return damp; - -} /* calc_damping */ - -/** - * This function calculates the norm of an update, dx[], based on the - * weighted values of x. - */ -static double calcWeightedNorm(const double wtX[], const double dx[], size_t dim) -{ - double norm = 0.0; - if (dim == 0) { - return 0.0; - } - for (size_t i = 0; i < dim; i++) { - norm += pow(dx[i] / wtX[i], 2); - } - return sqrt(norm/dim); -} - -void solveSP::calcWeights(double wtSpecies[], double wtResid[], - const Array2D& Jac, const double CSoln[], - const double abstol, const double reltol) -{ - // First calculate the weighting factor for the concentrations of the - // surface species and bulk species. - size_t kindex = 0; - for (size_t isp = 0; isp < m_numSurfPhases; isp++) { - double sd = m_ptrsSurfPhase[isp]->siteDensity(); - for (size_t k = 0; k < m_nSpeciesSurfPhase[isp]; k++, kindex++) { - wtSpecies[kindex] = abstol * sd + reltol * fabs(CSoln[kindex]); - } - } - if (m_bulkFunc == BULK_DEPOSITION) { - for (size_t isp = 0; isp < m_numBulkPhasesSS; isp++) { - double sd = m_bulkPhasePtrs[isp]->molarDensity(); - for (size_t k = 0; k < m_numBulkSpecies[isp]; k++, kindex++) { - wtSpecies[kindex] = abstol * sd + reltol * fabs(CSoln[kindex]); - } - } - } - - // Now do the residual Weights. Since we have the Jacobian, we will use it - // to generate a number based on the what a significant change in a solution - // variable does to each residual. This is a row sum scale operation. - for (size_t k = 0; k < m_neq; k++) { - wtResid[k] = 0.0; - for (size_t jcol = 0; jcol < m_neq; jcol++) { - wtResid[k] += fabs(Jac(k,jcol) * wtSpecies[jcol]); - } - } -} - -double solveSP::calc_t(double netProdRateSolnSP[], double XMolSolnSP[], int* label, - int* label_old, double* label_factor, int ioflag) -{ - double inv_timeScale = 1.0E-10; - size_t kindexSP = 0; - *label = 0; - updateMFSolnSP(XMolSolnSP); - for (size_t isp = 0; isp < m_numSurfPhases; isp++) { - // Get the interface kinetics associated with this surface - InterfaceKinetics* kin = m_objects[isp]; - - kin->getNetProductionRates(m_numEqn1.data()); - double sden = kin->thermo(0).molarDensity(); - for (size_t k = 0; k < m_nSpeciesSurfPhase[isp]; k++, kindexSP++) { - netProdRateSolnSP[kindexSP] = m_numEqn1[k]; - double tmp = std::max(XMolSolnSP[kindexSP], 1.0e-10); - tmp *= sden; - tmp = fabs(netProdRateSolnSP[kindexSP]/ tmp); - if (netProdRateSolnSP[kindexSP]> 0.0) { - tmp /= 100.; - } - if (tmp > inv_timeScale) { - inv_timeScale = tmp; - *label = int(kindexSP); - } - } - } - - // Increase time step exponentially as same species repeatedly controls time - // step - if (*label == *label_old) { - *label_factor *= 1.5; - } else { - *label_old = *label; - *label_factor = 1.0; - } - return inv_timeScale / *label_factor; -} // calc_t - -void solveSP::print_header(int ioflag, int ifunc, double time_scale, - int damping, double reltol, double abstol) -{ - if (ioflag) { - writelog("\n================================ SOLVESP CALL SETUP " - "========================================\n"); - if (ifunc == SFLUX_INITIALIZE) { - writelog("\n SOLVESP Called with Initialization turned on\n"); - writelogf(" Time scale input = %9.3e\n", time_scale); - } else if (ifunc == SFLUX_RESIDUAL) { - writelog("\n SOLVESP Called to calculate steady state residual\n"); - writelog(" from a good initial guess\n"); - } else if (ifunc == SFLUX_JACOBIAN) { - writelog("\n SOLVESP Called to calculate steady state Jacobian\n"); - writelog(" from a good initial guess\n"); - } else if (ifunc == SFLUX_TRANSIENT) { - writelog("\n SOLVESP Called to integrate surface in time\n"); - writelogf(" for a total of %9.3e sec\n", time_scale); - } else { - throw CanteraError("solveSP::print_header", - "Unknown ifunc flag = {}", ifunc); - } - - if (m_bulkFunc == BULK_DEPOSITION) { - writelog(" The composition of the Bulk Phases will be calculated\n"); - } else if (m_bulkFunc == BULK_ETCH) { - writelog(" Bulk Phases have fixed compositions\n"); - } else { - throw CanteraError("solveSP::print_header", - "Unknown bulkFunc flag = {}", m_bulkFunc); - } - - if (damping) { - writelog(" Damping is ON \n"); - } else { - writelog(" Damping is OFF \n"); - } - - writelogf(" Reltol = %9.3e, Abstol = %9.3e\n", reltol, abstol); - } - - if (ioflag == 1) { - writelog("\n\n\t Iter Time Del_t Damp DelX " - " Resid Name-Time Name-Damp\n"); - writelog("\t -----------------------------------------------" - "------------------------------------\n"); - } -} - -void solveSP::printIteration(int ioflag, double damp, int label_d, - int label_t, double inv_t, double t_real, - size_t iter, double update_norm, - double resid_norm, bool do_time, bool final) -{ - if (ioflag == 1) { - if (final) { - writelogf("\tFIN%3d ", iter); - } else { - writelogf("\t%6d ", iter); - } - if (do_time) { - writelogf("%9.4e %9.4e ", t_real, 1.0/inv_t); - } else { - writeline(' ', 22, false); - } - if (damp < 1.0) { - writelogf("%9.4e ", damp); - } else { - writeline(' ', 11, false); - } - writelogf("%9.4e %9.4e", update_norm, resid_norm); - if (do_time) { - size_t k = m_kinSpecIndex[label_t]; - size_t isp = m_kinObjIndex[label_t]; - writelog(" %-16s", m_objects[isp]->kineticsSpeciesName(k)); - } else { - writeline(' ', 16, false); - } - if (label_d >= 0) { - size_t k = m_kinSpecIndex[label_d]; - size_t isp = m_kinObjIndex[label_d]; - writelogf(" %-16s", m_objects[isp]->kineticsSpeciesName(k)); - } - if (final) { - writelog(" -- success"); - } - writelog("\n"); - } -} // printIteration - -} diff --git a/src/numerics/CVodesIntegrator.cpp b/src/numerics/CVodesIntegrator.cpp index f176576bff..481d7c86f3 100644 --- a/src/numerics/CVodesIntegrator.cpp +++ b/src/numerics/CVodesIntegrator.cpp @@ -160,6 +160,10 @@ void CVodesIntegrator::setTolerances(double reltol, size_t n, double* abstol) NV_Ith_S(m_abstol, i) = abstol[i]; } m_reltol = reltol; + if (m_cvode_mem) { + int flag = CVodeSVtolerances(m_cvode_mem, m_reltol, m_abstol); + checkError(flag, "setTolerances", "CVodeSVtolerances"); + } } void CVodesIntegrator::setTolerances(double reltol, double abstol) @@ -167,12 +171,26 @@ void CVodesIntegrator::setTolerances(double reltol, double abstol) m_itol = CV_SS; m_reltol = reltol; m_abstols = abstol; + if (m_cvode_mem) { + int flag = CVodeSStolerances(m_cvode_mem, m_reltol, m_abstols); + checkError(flag, "setTolerances", "CVodeSStolerances"); + } } void CVodesIntegrator::setSensitivityTolerances(double reltol, double abstol) { m_reltolsens = reltol; m_abstolsens = abstol; + if (m_cvode_mem && m_yS) { + vector atol(m_np); + for (size_t n = 0; n < m_np; n++) { + // This scaling factor is tuned so that reaction and species enthalpy + // sensitivities can be computed simultaneously with the same abstol. + atol[n] = m_abstolsens / m_func->m_paramScales[n]; + } + int flag = CVodeSensSStolerances(m_cvode_mem, m_reltolsens, atol.data()); + checkError(flag, "setSensitivityTolerances", "CVodeSensSStolerances"); + } } void CVodesIntegrator::setMethod(MethodType t) @@ -255,15 +273,7 @@ void CVodesIntegrator::sensInit(double t0, FuncEval& func) int flag = CVodeSensInit(m_cvode_mem, static_cast(m_np), CV_STAGGERED, CVSensRhsFn(0), m_yS); checkError(flag, "sensInit", "CVodeSensInit"); - - vector atol(m_np); - for (size_t n = 0; n < m_np; n++) { - // This scaling factor is tuned so that reaction and species enthalpy - // sensitivities can be computed simultaneously with the same abstol. - atol[n] = m_abstolsens / func.m_paramScales[n]; - } - flag = CVodeSensSStolerances(m_cvode_mem, m_reltolsens, atol.data()); - checkError(flag, "sensInit", "CVodeSensSStolerances"); + setSensitivityTolerances(m_reltolsens, m_abstolsens); } void CVodesIntegrator::initialize(double t0, FuncEval& func) diff --git a/src/zeroD/ConstPressureMoleReactor.cpp b/src/zeroD/ConstPressureMoleReactor.cpp index e6801cff6a..87347d1b28 100644 --- a/src/zeroD/ConstPressureMoleReactor.cpp +++ b/src/zeroD/ConstPressureMoleReactor.cpp @@ -106,6 +106,18 @@ void ConstPressureMoleReactor::eval(double time, double* LHS, double* RHS) } } +void ConstPressureMoleReactor::evalSteady(double t, double* LHS, double* RHS) +{ + eval(t, LHS, RHS); + RHS[1] = m_mass - m_initialMass; +} + +vector ConstPressureMoleReactor::initializeSteady() +{ + m_initialMass = m_mass; + return {1}; +} + size_t ConstPressureMoleReactor::componentIndex(const string& nm) const { if (nm == "enthalpy") { diff --git a/src/zeroD/ConstPressureReactor.cpp b/src/zeroD/ConstPressureReactor.cpp index 1d43a640b4..2000219add 100644 --- a/src/zeroD/ConstPressureReactor.cpp +++ b/src/zeroD/ConstPressureReactor.cpp @@ -117,13 +117,15 @@ void ConstPressureReactor::eval(double time, double* LHS, double* RHS) } } -vector ConstPressureReactor::steadyConstraints() const { - if (nSurfs() != 0) { - throw CanteraError("ConstPressureReactor::steadyConstraints", - "Steady state solver cannot currently be used with ConstPressureReactor" - " when reactor surfaces are present.\n" - "See https://github.com/Cantera/enhancements/issues/234"); - } +void ConstPressureReactor::evalSteady(double t, double* LHS, double* RHS) +{ + eval(t, LHS, RHS); + RHS[0] = m_mass - m_initialMass; +} + +vector ConstPressureReactor::initializeSteady() +{ + m_initialMass = m_mass; return {0}; // mass } diff --git a/src/zeroD/IdealGasConstPressureReactor.cpp b/src/zeroD/IdealGasConstPressureReactor.cpp index 6ddde615fe..a6ff09925d 100644 --- a/src/zeroD/IdealGasConstPressureReactor.cpp +++ b/src/zeroD/IdealGasConstPressureReactor.cpp @@ -114,14 +114,19 @@ void IdealGasConstPressureReactor::eval(double time, double* LHS, double* RHS) } } -vector IdealGasConstPressureReactor::steadyConstraints() const +void IdealGasConstPressureReactor::evalSteady(double t, double* LHS, double* RHS) { - if (nSurfs() != 0) { - throw CanteraError("IdealGasConstPressureReactor::steadyConstraints", - "Steady state solver cannot currently be used with " - " IdealGasConstPressureReactor when reactor surfaces are present.\n" - "See https://github.com/Cantera/enhancements/issues/234"); + eval(0.0, LHS, RHS); + RHS[0] = m_mass - m_initialMass; + if (!energyEnabled()) { + RHS[1] = m_thermo->temperature() - m_initialTemperature; } +} + +vector IdealGasConstPressureReactor::initializeSteady() +{ + m_initialMass = m_mass; + m_initialTemperature = m_thermo->temperature(); if (energyEnabled()) { return {0}; // mass } else { diff --git a/src/zeroD/IdealGasMoleReactor.cpp b/src/zeroD/IdealGasMoleReactor.cpp index b6dcfdbbdd..ddd6fe03d6 100644 --- a/src/zeroD/IdealGasMoleReactor.cpp +++ b/src/zeroD/IdealGasMoleReactor.cpp @@ -85,14 +85,10 @@ double IdealGasMoleReactor::lowerBound(size_t k) const { } } -vector IdealGasMoleReactor::steadyConstraints() const +vector IdealGasMoleReactor::initializeSteady() { - if (nSurfs() != 0) { - throw CanteraError("IdealGasMoleReactor::steadyConstraints", - "Steady state solver cannot currently be used with IdealGasMoleReactor" - " when reactor surfaces are present.\n" - "See https://github.com/Cantera/enhancements/issues/234"); - } + m_initialVolume = m_vol; + m_initialTemperature = m_thermo->temperature(); if (energyEnabled()) { return {1}; // volume } else { @@ -173,6 +169,15 @@ void IdealGasMoleReactor::eval(double time, double* LHS, double* RHS) } } +void IdealGasMoleReactor::evalSteady(double t, double* LHS, double* RHS) +{ + eval(t, LHS, RHS); + if (!energyEnabled()) { + RHS[0] = m_thermo->temperature() - m_initialTemperature; + } + RHS[1] = m_vol - m_initialVolume; +} + void IdealGasMoleReactor::getJacobianElements(vector>& trips) { // dnk_dnj represents d(dot(n_k)) / d (n_j) but is first assigned as diff --git a/src/zeroD/IdealGasReactor.cpp b/src/zeroD/IdealGasReactor.cpp index d3b178e909..02031f21e5 100644 --- a/src/zeroD/IdealGasReactor.cpp +++ b/src/zeroD/IdealGasReactor.cpp @@ -122,14 +122,19 @@ void IdealGasReactor::eval(double time, double* LHS, double* RHS) } } -vector IdealGasReactor::steadyConstraints() const +void IdealGasReactor::evalSteady(double t, double* LHS, double* RHS) { - if (nSurfs() != 0) { - throw CanteraError("IdealGasReactor::steadyConstraints", - "Steady state solver cannot currently be used with IdealGasReactor" - " when reactor surfaces are present.\n" - "See https://github.com/Cantera/enhancements/issues/234"); + eval(t, LHS, RHS); + RHS[1] = m_vol - m_initialVolume; + if (!energyEnabled()) { + RHS[2] = m_thermo->temperature() - m_initialTemperature; } +} + +vector IdealGasReactor::initializeSteady() +{ + m_initialTemperature = m_thermo->temperature(); + m_initialVolume = m_vol; if (energyEnabled()) { return {1}; // volume } else { diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index aefd6e4c4e..83ffcbd7ea 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -184,6 +184,12 @@ void Reactor::eval(double time, double* LHS, double* RHS) } } +void Reactor::evalSteady(double time, double* LHS, double* RHS) +{ + eval(time, LHS, RHS); + RHS[1] = m_vol - m_initialVolume; +} + void Reactor::evalWalls(double t) { // time is currently unused @@ -196,19 +202,15 @@ void Reactor::evalWalls(double t) } } -vector Reactor::steadyConstraints() const +vector Reactor::initializeSteady() { if (!energyEnabled()) { - throw CanteraError("Reactor::steadyConstraints", "Steady state solver cannot" + throw CanteraError("Reactor::initializeSteady", "Steady state solver cannot" " be used with {0} when energy equation is disabled." "\nConsider using IdealGas{0} instead.\n" "See https://github.com/Cantera/enhancements/issues/234", type()); } - if (nSurfs() != 0) { - throw CanteraError("Reactor::steadyConstraints", "Steady state solver cannot" - " currently be used when reactor surfaces are present.\n" - "See https://github.com/Cantera/enhancements/issues/234."); - } + m_initialVolume = m_vol; return {1}; // volume } diff --git a/src/zeroD/ReactorBase.cpp b/src/zeroD/ReactorBase.cpp index 1ff243b6d2..526d2f6935 100644 --- a/src/zeroD/ReactorBase.cpp +++ b/src/zeroD/ReactorBase.cpp @@ -66,16 +66,31 @@ bool ReactorBase::setDefaultName(map& counts) void ReactorBase::addInlet(FlowDevice& inlet) { + if (m_net) { + throw CanteraError("ReactorBase::addInlet", + "Cannot add an inlet to reactor '{}' after it has been" + " added to a ReactorNet.", m_name); + } m_inlet.push_back(&inlet); } void ReactorBase::addOutlet(FlowDevice& outlet) { + if (m_net) { + throw CanteraError("ReactorBase::addOutlet", + "Cannot add an outlet to reactor '{}' after it has been" + " added to a ReactorNet.", m_name); + } m_outlet.push_back(&outlet); } void ReactorBase::addWall(WallBase& w, int lr) { + if (m_net) { + throw CanteraError("ReactorBase::addWall", + "Cannot add a wall to reactor '{}' after it has been" + " added to a ReactorNet.", m_name); + } m_wall.push_back(&w); if (lr == 0) { m_lr.push_back(0); @@ -91,6 +106,11 @@ WallBase& ReactorBase::wall(size_t n) void ReactorBase::addSurface(ReactorSurface* surf) { + if (m_net) { + throw CanteraError("ReactorBase::addSurface", + "Cannot add a surface to reactor '{}' after it has been" + " added to a ReactorNet.", m_name); + } if (find(m_surfaces.begin(), m_surfaces.end(), surf) == m_surfaces.end()) { m_surfaces.push_back(surf); } @@ -160,6 +180,10 @@ ReactorNet& ReactorBase::network() void ReactorBase::setNetwork(ReactorNet* net) { + if (m_net) { + throw CanteraError("ReactorBase::setNetwork", + "Reactor {} is already part of a ReactorNet", m_name); + } m_net = net; } diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index b7fc6f4f98..e0b6de41d3 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -7,6 +7,7 @@ #include "cantera/zeroD/FlowDevice.h" #include "cantera/zeroD/ReactorSurface.h" #include "cantera/zeroD/Wall.h" +#include "cantera/zeroD/Reservoir.h" #include "cantera/base/utilities.h" #include "cantera/base/Array.h" #include "cantera/base/Solution.h" @@ -18,38 +19,160 @@ #include "cantera/oneD/MultiNewton.h" #include +#include namespace Cantera { -ReactorNet::ReactorNet() +ReactorNet::~ReactorNet() { - suppressErrors(true); } ReactorNet::ReactorNet(shared_ptr reactor) + : ReactorNet(span>{&reactor, 1}) { - suppressErrors(true); - addReactor(reactor); } -ReactorNet::ReactorNet(vector>& reactors) +ReactorNet::ReactorNet(span> reactors) { + if (reactors.empty()) { + throw CanteraError("ReactorNet::ReactorNet", "No reactors in network!"); + } + suppressErrors(true); - for (auto& reactor : reactors) { - addReactor(reactor); + bool isOde = true; + + // Names of Reactors and ReactorSurfaces using each Solution; should be only one + // reactor per Solution. + map> solutions; + + // All ReactorBase objects connected to this network, starting with the ones + // explicitly included. + std::deque> reactorQueue(reactors.begin(), reactors.end()); + std::set> visited; + + while (!reactorQueue.empty()) { + auto R = reactorQueue.front(); + reactorQueue.pop_front(); + + if (visited.find(R) != visited.end()) { + continue; + } + visited.insert(R); + + if (auto bulk = std::dynamic_pointer_cast(R)) { + m_bulkReactors.push_back(bulk); + m_reactors.push_back(R); + } else if (auto surf = std::dynamic_pointer_cast(R)) { + m_surfaces.push_back(surf); + m_reactors.push_back(R); + for (size_t i = 0; i < surf->nAdjacent(); i++) { + reactorQueue.push_back(surf->adjacent(i)); + } + } else if (auto res = std::dynamic_pointer_cast(R)) { + m_reservoirs.push_back(res); + } + + // Discover reservoirs, flow devices, and walls connected to this reactor + for (size_t i = 0; i < R->nInlets(); i++) { + auto& inlet = R->inlet(i); + m_flowDevices.insert(&inlet); + reactorQueue.push_back(inlet.in().shared_from_this()); + } + for (size_t i = 0; i < R->nOutlets(); i++) { + auto& outlet = R->outlet(i); + m_flowDevices.insert(&outlet); + reactorQueue.push_back(outlet.out().shared_from_this()); + } + for (size_t i = 0; i < R->nWalls(); i++) { + auto& wall = R->wall(i); + m_walls.insert(&wall); + reactorQueue.push_back(wall.left().shared_from_this()); + reactorQueue.push_back(wall.right().shared_from_this()); + } + + auto phase = R->phase(); + for (size_t i = 0; i < R->nSurfs(); i++) { + if (R->surface(i)->phase()->adjacent(phase->name()) != phase) { + throw CanteraError("ReactorNet::initialize", + "Bulk phase '{}' used by interface '{}' must be the same object\n" + "as the contents of the adjacent reactor '{}'.", + phase->name(), R->surface(i)->name(), R->name()); + } + reactorQueue.push_back(R->surface(i)->shared_from_this()); + } + R->setNetwork(this); + updateNames(*R); + solutions[phase.get()].insert(R->name()); } -} -ReactorNet::~ReactorNet() -{ + for (auto& R : m_bulkReactors) { + if (R->type() == "FlowReactor" && m_bulkReactors.size() > 1) { + throw CanteraError("ReactorNet::initialize", + "FlowReactors must be used alone."); + } + m_timeIsIndependent = R->timeIsIndependent(); + R->initialize(m_time); + isOde = R->isOde(); + R->setOffset(m_nv); + size_t nv = R->neq(); + m_nv += nv; + + for (auto current : m_bulkReactors) { + if (current->isOde() != R->isOde()) { + throw CanteraError("ReactorNet::ReactorNet", + "Cannot mix Reactor types using both ODEs and DAEs ({} and {})", + current->type(), R->type()); + } + if (current->timeIsIndependent() != R->timeIsIndependent()) { + throw CanteraError("ReactorNet::ReactorNet", + "Cannot mix Reactor types using time and space as independent variables" + "\n({} and {})", current->type(), R->type()); + } + } + } + + for (auto surf : m_surfaces) { + surf->setOffset(m_nv); + surf->initialize(m_time); + m_nv += surf->neq(); + solutions[surf->phase().get()].insert(surf->name()); + } + + for (auto& [soln, reactors] : solutions) { + if (reactors.size() > 1) { + string shared; + for (auto r : reactors) { + if (r != *reactors.begin()) { + shared += ", "; + } + shared += fmt::format("'{}'", r); + } + throw CanteraError("ReactorNet::initialize", "The following reactors /" + " reactor surfaces are using the same Solution object: {}. Use" + " independent Solution objects or set the 'clone' argument to 'true'" + " when creating the Reactor or ReactorSurface objects.", shared); + } + } + + m_ydot.resize(m_nv, 0.0); + m_yest.resize(m_nv, 0.0); + m_advancelimits.resize(m_nv, -1.0); + m_atol.resize(m_nv); + + m_integ.reset(newIntegrator(isOde ? "CVODE" : "IDA")); + // use backward differencing, with a full Jacobian computed + // numerically, and use a Newton linear iterator + m_integ->setMethod(BDF_Method); + m_integ->setLinearSolverType("DENSE"); + setTolerances(-1.0, -1.0); } void ReactorNet::setInitialTime(double time) { m_time = time; m_initial_time = time; - m_integrator_init = false; + m_needIntegratorInit = true; } void ReactorNet::setMaxTimeStep(double maxstep) @@ -71,7 +194,8 @@ void ReactorNet::setTolerances(double rtol, double atol) if (atol >= 0.0) { m_atols = atol; } - m_init = false; + fill(m_atol.begin(), m_atol.end(), m_atols); + m_integ->setTolerances(m_rtol, neq(), m_atol.data()); } void ReactorNet::setSensitivityTolerances(double rtol, double atol) @@ -82,7 +206,7 @@ void ReactorNet::setSensitivityTolerances(double rtol, double atol) if (atol >= 0.0) { m_atolsens = atol; } - m_init = false; + m_integ->setSensitivityTolerances(m_rtolsens, m_atolsens); } double ReactorNet::time() { @@ -105,110 +229,21 @@ double ReactorNet::distance() { void ReactorNet::initialize() { - m_nv = 0; - debuglog("Initializing reactor network.\n", m_verbose); - if (m_reactors.empty()) { - throw CanteraError("ReactorNet::initialize", - "no reactors in network!"); - } - // Names of Reactors and ReactorSurfaces using each Solution; should be only one - map> solutions; - // Unique ReactorSurface objects. Can be attached to multiple Reactor objects - set surfaces; - m_flowDevices.clear(); - m_walls.clear(); - m_reservoirs.clear(); - m_reactors.resize(m_bulkReactors.size()); // Surfaces will be re-added - for (size_t n = 0; n < m_bulkReactors.size(); n++) { - Reactor& r = *m_bulkReactors[n]; - shared_ptr bulk = r.phase(); - r.setOffset(m_nv); - r.initialize(m_time); - size_t nv = r.neq(); - m_nv += nv; - - if (m_verbose) { - writelog("Reactor {:d}: {:d} variables.\n", n, nv); - writelog(" {:d} sensitivity params.\n", r.nSensParams()); - } - if (r.type() == "FlowReactor" && m_bulkReactors.size() > 1) { - throw CanteraError("ReactorNet::initialize", - "FlowReactors must be used alone."); - } - solutions[bulk.get()].insert(r.name()); - for (size_t i = 0; i < r.nSurfs(); i++) { - if (r.surface(i)->phase()->adjacent(bulk->name()) != bulk) { - throw CanteraError("ReactorNet::initialize", - "Bulk phase '{}' used by interface '{}' must be the same object\n" - "as the contents of the adjacent reactor '{}'.", - bulk->name(), r.surface(i)->name(), r.name()); - } - surfaces.insert(r.surface(i)); - } - // Discover reservoirs, flow devices, and walls connected to this reactor - for (size_t i = 0; i < r.nInlets(); i++) { - auto& inlet = r.inlet(i); - m_flowDevices.insert(&inlet); - if (inlet.in().type() == "Reservoir") { - m_reservoirs.insert(&inlet.in()); - inlet.in().setNetwork(this); - } - } - for (size_t i = 0; i < r.nOutlets(); i++) { - auto& outlet = r.outlet(i); - m_flowDevices.insert(&outlet); - if (outlet.out().type() == "Reservoir") { - m_reservoirs.insert(&outlet.out()); - outlet.out().setNetwork(this); - } - } - for (size_t i = 0; i < r.nWalls(); i++) { - auto& wall = r.wall(i); - m_walls.insert(&wall); - if (wall.left().type() == "Reservoir") { - m_reservoirs.insert(&wall.left()); - } else if (wall.right().type() == "Reservoir") { - m_reservoirs.insert(&wall.right()); - } + if (m_integratorInitialized) { + if (m_needIntegratorInit) { + reinitialize(); } + return; } - for (auto surf : surfaces) { - solutions[surf->phase().get()].insert(surf->name()); - m_reactors.push_back(surf->shared_from_this()); - surf->setOffset(m_nv); - surf->initialize(m_time); - m_nv += surf->neq(); - solutions[surf->phase().get()].insert(surf->name()); - } - for (auto& [soln, reactors] : solutions) { - if (reactors.size() > 1) { - string shared; - for (auto r : reactors) { - if (r != *reactors.begin()) { - shared += ", "; - } - shared += fmt::format("'{}'", r); - } - throw CanteraError("ReactorNet::initialize", "The following reactors /" - " reactor surfaces are using the same Solution object: {}. Use" - " independent Solution objects or set the 'clone' argument to 'true'" - " when creating the Reactor or ReactorSurface objects.", shared); - } - } - - m_ydot.resize(m_nv,0.0); - m_yest.resize(m_nv,0.0); - m_advancelimits.resize(m_nv,-1.0); - m_atol.resize(neq()); - fill(m_atol.begin(), m_atol.end(), m_atols); - m_integ->setTolerances(m_rtol, neq(), m_atol.data()); - m_integ->setSensitivityTolerances(m_rtolsens, m_atolsens); if (!m_linearSolverType.empty()) { m_integ->setLinearSolverType(m_linearSolverType); } if (m_precon) { m_integ->setPreconditioner(m_precon); } + for (auto& reactor : m_reactors) { + reactor->updateConnected(true); + } m_integ->initialize(m_time, *this); if (m_verbose) { writelog("Number of equations: {:d}\n", neq()); @@ -217,22 +252,22 @@ void ReactorNet::initialize() if (m_integ->preconditionerSide() != PreconditionerSide::NO_PRECONDITION) { checkPreconditionerSupported(); } - m_integrator_init = true; - m_init = true; + m_needIntegratorInit = false; + m_integratorInitialized = true; } void ReactorNet::reinitialize() { - if (m_init) { + if (m_integratorInitialized) { debuglog("Re-initializing reactor network.\n", m_verbose); - for (auto& reservoir : m_reservoirs) { - reservoir->updateConnected(true); + for (auto& reactor : m_reactors) { + reactor->updateConnected(true); } m_integ->reinitialize(m_time, *this); if (m_integ->preconditionerSide() != PreconditionerSide::NO_PRECONDITION) { checkPreconditionerSupported(); } - m_integrator_init = true; + m_needIntegratorInit = false; } else { initialize(); } @@ -241,13 +276,13 @@ void ReactorNet::reinitialize() void ReactorNet::setLinearSolverType(const string& linSolverType) { m_linearSolverType = linSolverType; - m_integrator_init = false; + m_needIntegratorInit = true; } void ReactorNet::setPreconditioner(shared_ptr preconditioner) { m_precon = preconditioner; - m_integrator_init = false; + m_needIntegratorInit = true; } void ReactorNet::setMaxSteps(int nmax) @@ -262,11 +297,7 @@ int ReactorNet::maxSteps() void ReactorNet::advance(double time) { - if (!m_init) { - initialize(); - } else if (!m_integrator_init) { - reinitialize(); - } + initialize(); m_integ->integrate(time); m_time = m_integ->currentTime(); updateState(m_integ->solution()); @@ -274,12 +305,7 @@ void ReactorNet::advance(double time) double ReactorNet::advance(double time, bool applylimit) { - if (!m_init) { - initialize(); - } else if (!m_integrator_init) { - reinitialize(); - } - + initialize(); if (!applylimit || !hasAdvanceLimits()) { // No limit enforcement requested; integrate to requested time advance(time); @@ -355,11 +381,7 @@ double ReactorNet::advance(double time, bool applylimit) double ReactorNet::step() { - if (!m_init) { - initialize(); - } else if (!m_integrator_init) { - reinitialize(); - } + initialize(); m_time = m_integ->step(m_time + 1.0); updateState(m_integ->solution()); return m_time; @@ -367,11 +389,7 @@ double ReactorNet::step() void ReactorNet::solveSteady(int loglevel) { - if (!m_init) { - initialize(); - } else if (!m_integrator_init) { - reinitialize(); - } + initialize(); vector y(neq()); getState(y.data()); SteadyReactorSolver solver(this, y.data()); @@ -383,11 +401,7 @@ void ReactorNet::solveSteady(int loglevel) Eigen::SparseMatrix ReactorNet::steadyJacobian(double rdt) { - if (!m_init) { - initialize(); - } else if (!m_integrator_init) { - reinitialize(); - } + initialize(); vector y0(neq()); vector y1(neq()); getState(y0.data()); @@ -401,10 +415,7 @@ Eigen::SparseMatrix ReactorNet::steadyJacobian(double rdt) void ReactorNet::getEstimate(double time, int k, double* yest) { - if (!m_init) { - initialize(); - } - // initialize + initialize(); double* cvode_dky = m_integ->solution(); for (size_t j = 0; j < m_nv; j++) { yest[j] = cvode_dky[j]; @@ -465,42 +476,7 @@ void ReactorNet::evalRootFunctions(double t, const double* y, double* gout) gout[0] = g; } -void ReactorNet::addReactor(shared_ptr reactor) -{ - auto r = std::dynamic_pointer_cast(reactor); - if (!r) { - throw CanteraError("ReactorNet::addReactor", - "Reactor with type '{}' cannot be added to network.", - reactor->type()); - } - - for (auto current : m_bulkReactors) { - if (current->isOde() != r->isOde()) { - throw CanteraError("ReactorNet::addReactor", - "Cannot mix Reactor types using both ODEs and DAEs ({} and {})", - current->type(), r->type()); - } - if (current->timeIsIndependent() != r->timeIsIndependent()) { - throw CanteraError("ReactorNet::addReactor", - "Cannot mix Reactor types using time and space as independent variables" - "\n({} and {})", current->type(), r->type()); - } - } - m_timeIsIndependent = r->timeIsIndependent(); - r->setNetwork(this); - m_bulkReactors.push_back(r.get()); - m_reactors.push_back(r); - if (!m_integ) { - m_integ.reset(newIntegrator(r->isOde() ? "CVODE" : "IDA")); - // use backward differencing, with a full Jacobian computed - // numerically, and use a Newton linear iterator - m_integ->setMethod(BDF_Method); - m_integ->setLinearSolverType("DENSE"); - } - updateNames(*r); -} - -void ReactorNet::updateNames(Reactor& r) +void ReactorNet::updateNames(ReactorBase& r) { // ensure that reactors and components have reproducible names r.setDefaultName(m_counts); @@ -563,6 +539,21 @@ void ReactorNet::eval(double t, double* y, double* ydot, double* p) checkFinite("ydot", ydot, m_nv); } +void ReactorNet::evalSteady(double* y, double* residual) +{ + updateState(y); + m_LHS.assign(m_nv, 1); + m_RHS.assign(m_nv, 0); + for (auto& R : m_reactors) { + size_t offset = R->offset(); + R->evalSteady(m_time, m_LHS.data() + offset, m_RHS.data() + offset); + for (size_t i = offset; i < offset + R->neq(); i++) { + residual[i] = m_RHS[i] / m_LHS[i]; + } + } + checkFinite("residual", residual, m_nv); +} + void ReactorNet::evalDae(double t, double* y, double* ydot, double* p, double* residual) { m_time = t; @@ -585,9 +576,7 @@ void ReactorNet::getConstraints(double* constraints) double ReactorNet::sensitivity(size_t k, size_t p) { - if (!m_init) { - initialize(); - } + initialize(); if (p >= m_sens_params.size()) { throw IndexError("ReactorNet::sensitivity", "m_sens_params", p, m_sens_params.size()); @@ -631,9 +620,7 @@ void ReactorNet::updateState(double* y) void ReactorNet::getDerivative(int k, double* dky) { - if (!m_init) { - initialize(); - } + initialize(); double* cvode_dky = m_integ->derivative(m_time, k); for (size_t j = 0; j < m_nv; j++) { dky[j] = cvode_dky[j]; @@ -642,9 +629,7 @@ void ReactorNet::getDerivative(int k, double* dky) void ReactorNet::setAdvanceLimits(const double *limits) { - if (!m_init) { - initialize(); - } + initialize(); for (auto& R : m_bulkReactors) { R->setAdvanceLimits(limits + R->offset()); } @@ -689,9 +674,7 @@ void ReactorNet::getStateDae(double* y, double* ydot) size_t ReactorNet::globalComponentIndex(const string& component, size_t reactor) { - if (!m_init) { - initialize(); - } + initialize(); return m_reactors[reactor]->offset() + m_reactors[reactor]->componentIndex(component); } @@ -715,7 +698,7 @@ double ReactorNet::upperBound(size_t i) const { size_t iTot = 0; size_t i0 = i; - for (auto r : m_bulkReactors) { + for (auto r : m_reactors) { if (i < r->neq()) { return r->upperBound(i); } else { @@ -730,7 +713,7 @@ double ReactorNet::lowerBound(size_t i) const { size_t iTot = 0; size_t i0 = i; - for (auto r : m_bulkReactors) { + for (auto r : m_reactors) { if (i < r->neq()) { return r->lowerBound(i); } else { @@ -750,7 +733,7 @@ void ReactorNet::resetBadValues(double* y) { size_t ReactorNet::registerSensitivityParameter( const string& name, double value, double scale) { - if (m_integrator_init) { + if (m_integratorInitialized) { throw CanteraError("ReactorNet::registerSensitivityParameter", "Sensitivity parameters cannot be added after the" "integrator has been initialized."); @@ -856,18 +839,19 @@ SteadyReactorSolver::SteadyReactorSolver(ReactorNet* net, double* x0) SteadyStateSystem::resize(); m_initialState.assign(x0, x0 + m_size); setInitialGuess(x0); + setJacobianPerturbation(m_jacobianRelPerturb, 1000 * m_net->atol(), + m_jacobianThreshold); m_mask.assign(m_size, 1); size_t start = 0; for (size_t i = 0; i < net->nReactors(); i++) { auto& R = net->reactor(i); - for (auto& m : R.steadyConstraints()) { - m_algebraic.push_back(start + m); + R.updateState(x0 + start); + auto algebraic = R.initializeSteady(); + for (auto& m : algebraic) { + m_mask[start + m] = 0; } start += R.neq(); } - for (auto& n : m_algebraic) { - m_mask[n] = 0; - } } void SteadyReactorSolver::eval(double* x, double* r, double rdt, int count) @@ -876,13 +860,9 @@ void SteadyReactorSolver::eval(double* x, double* r, double rdt, int count) rdt = m_rdt; } vector xv(x, x + size()); - m_net->eval(0.0, x, r, nullptr); + m_net->evalSteady(x, r); for (size_t i = 0; i < size(); i++) { - r[i] -= (x[i] - m_initialState[i]) * rdt; - } - // Hold algebraic constraints fixed - for (auto& n : m_algebraic) { - r[n] = x[n] - m_initialState[n]; + r[i] -= (x[i] - m_initialState[i]) * rdt * m_mask[i]; } } @@ -934,7 +914,7 @@ double SteadyReactorSolver::weightedNorm(const double* step) const double sum = 0.0; const double* x = m_state->data(); for (size_t i = 0; i < size(); i++) { - double ewt = m_net->rtol()*x[i] + m_net->atol(); + double ewt = m_net->rtol()*fabs(x[i]) + m_net->atol(); double f = step[i] / ewt; sum += f*f; } diff --git a/src/zeroD/ReactorSurface.cpp b/src/zeroD/ReactorSurface.cpp index 92a02b6df8..b7e76c7279 100644 --- a/src/zeroD/ReactorSurface.cpp +++ b/src/zeroD/ReactorSurface.cpp @@ -26,7 +26,7 @@ ReactorSurface::ReactorSurface(shared_ptr soln, vector> adjacent; for (auto R : reactors) { adjacent.push_back(R->phase()); - m_reactors.push_back(R.get()); + m_reactors.push_back(R); R->addSurface(this); } if (clone) { @@ -95,6 +95,11 @@ void ReactorSurface::initialize(double t0) m_thermo->setState_TP(m_reactors[0]->temperature(), m_reactors[0]->pressure()); } +vector ReactorSurface::initializeSteady() +{ + return {0}; // sum of coverages constraint +} + void ReactorSurface::updateState(double* y) { m_surf->setCoveragesNoNorm(y); @@ -114,6 +119,18 @@ void ReactorSurface::eval(double t, double* LHS, double* RHS) RHS[0] = sum; } +void ReactorSurface::evalSteady(double t, double* LHS, double* RHS) +{ + eval(t, LHS, RHS); + vector cov(m_nsp); + m_surf->getCoverages(cov.data()); + double sum = 0.0; + for (size_t k = 0; k < m_nsp; k++) { + sum += cov[k]; + } + RHS[0] = 1.0 - sum; +} + void ReactorSurface::applySensitivity(double* params) { if (!params) { @@ -169,6 +186,25 @@ string ReactorSurface::componentName(size_t k) return m_surf->speciesName(k); } +double ReactorSurface::upperBound(size_t k) const +{ + return 1.0; +} + +double ReactorSurface::lowerBound(size_t k) const +{ + return -Tiny; +} + +void ReactorSurface::resetBadValues(double* y) +{ + for (size_t k = 0; k < m_nsp; k++) { + y[k] = std::max(y[k], 0.0); + } + m_surf->setCoverages(y); + m_surf->getCoverages(y); +} + // ------ MoleReactorSurface methods ------ void MoleReactorSurface::initialize(double t0) @@ -189,7 +225,7 @@ void MoleReactorSurface::initialize(double t0) int offset = static_cast(R->offset() + R->speciesOffset()); for (size_t k = 0; k < nsp; k++) { m_kin2net[k + k0] = offset + k; - m_kin2reactor[k + k0] = R; + m_kin2reactor[k + k0] = R.get(); } } } @@ -222,6 +258,33 @@ void MoleReactorSurface::eval(double t, double* LHS, double* RHS) } } +void MoleReactorSurface::evalSteady(double t, double* LHS, double* RHS) +{ + eval(t, LHS, RHS); + double cov_sum = 0.0; + for (size_t k = 0; k < m_nsp; k++) { + cov_sum += m_cov_tmp[k]; + } + RHS[0] = 1.0 - cov_sum; +} + +double MoleReactorSurface::upperBound(size_t k) const +{ + return BigNumber; +} + +double MoleReactorSurface::lowerBound(size_t k) const +{ + return -Tiny; +} + +void MoleReactorSurface::resetBadValues(double* y) +{ + for (size_t k = 0; k < m_nsp; k++) { + y[k] = std::max(y[k], 0.0); + } +} + void MoleReactorSurface::getJacobianElements(vector>& trips) { auto sdot_ddC = m_kinetics->netProductionRates_ddCi(); diff --git a/test/data/haca2.yaml b/test/data/haca2.yaml index a9a668406f..3103a3a206 100644 --- a/test/data/haca2.yaml +++ b/test/data/haca2.yaml @@ -22,7 +22,7 @@ phases: species: - gri30.yaml/species: [H, N2, CH3, CH4, C2H2, H2, OH, H2O, CO, O2] state: - T: 1100.0 + T: 1400.0 P: 1.01325e+05 X: {H: 0.01, N2: 0.8899, H2: 0.04, CH4: 0.01, C2H2: 0.01, OH: 1.0e-04, H2O: 0.04, O2: 1.0e-03} @@ -30,6 +30,9 @@ phases: thermo: fixed-stoichiometry elements: [C] species: [CB-CB3] + state: + T: 1400.0 + P: 1.01325e+05 - name: soot_interface thermo: ideal-surface adjacent-phases: [gas, soot] diff --git a/test/matlab/ctTestConstPressureReactor.m b/test/matlab/ctTestConstPressureReactor.m index 2a36563158..4981735926 100644 --- a/test/matlab/ctTestConstPressureReactor.m +++ b/test/matlab/ctTestConstPressureReactor.m @@ -12,7 +12,8 @@ r2 solid resGas - env + env1 + env2 w1 w2 mfc1 @@ -60,23 +61,24 @@ function makeReactors(self, arg) self.r2.V = 0.2; self.resGas.TP = {T0 - 300, P0}; - self.env = ct.zeroD.Reservoir(self.resGas); + self.env1 = ct.zeroD.Reservoir(self.resGas); + self.env2 = ct.zeroD.Reservoir(self.resGas); U = 300 * arg.addQ; - self.w1 = ct.zeroD.Wall(self.r1, self.env); + self.w1 = ct.zeroD.Wall(self.r1, self.env1); self.w1.area = 0.1; self.w1.expansionRateCoeff = 1e3; self.w1.heatTransferCoeff = U; - self.w2 = ct.zeroD.Wall(self.r2, self.env); + self.w2 = ct.zeroD.Wall(self.r2, self.env2); self.w2.area = 0.1; self.w2.heatTransferCoeff = U; if arg.addMdot - self.mfc1 = ct.zeroD.MassFlowController(self.env, self.r1); + self.mfc1 = ct.zeroD.MassFlowController(self.env1, self.r1); self.mfc1.massFlowRate = 0.05; - self.mfc2 = ct.zeroD.MassFlowController(self.env, self.r2); + self.mfc2 = ct.zeroD.MassFlowController(self.env2, self.r2); self.mfc2.massFlowRate = 0.05; end diff --git a/test/matlab/ctTestFlowDevice.m b/test/matlab/ctTestFlowDevice.m index 3fb10d36ae..8af7e1e442 100644 --- a/test/matlab/ctTestFlowDevice.m +++ b/test/matlab/ctTestFlowDevice.m @@ -5,7 +5,6 @@ gas2 r1 r2 - net end properties (SetAccess = protected) @@ -38,16 +37,11 @@ function makeReactors(self, arg) self.gas2 = self.gas1; end - if arg.nr == 1 - self.net = ct.zeroD.ReactorNet(self.r1); - elseif arg.nr >= 2 + if arg.nr >= 2 self.gas2.TPX = {arg.T2, arg.P2, arg.X2}; self.r2 = ct.zeroD.Reactor(self.gas2); self.r2.energyEnabled = true; - self.net = ct.zeroD.ReactorNet({self.r1, self.r2}); end - - self.verifyEqual(self.net.time, 0, 'AbsTol', self.atol); end end @@ -61,7 +55,7 @@ function testMFC(self) function testMFCType(self) self.makeReactors('nr', 2); mfc = ct.zeroD.MassFlowController(self.r1, self.r2); - n1 = ct.zeroD.ReactorNet({self.r1, self.r2}); + net = ct.zeroD.ReactorNet({self.r1, self.r2}); self.verifyEqual(mfc.type, 'MassFlowController'); self.verifyTrue(startsWith(mfc.name, 'MassFlowController_')); @@ -75,11 +69,12 @@ function testMFCErrors(self) function testValve1(self) self.makeReactors('P1', ct.OneAtm, 'X1', 'AR:1.0', 'X2', 'O2:1.0'); - self.net.rtol = 1e-12; valve = ct.zeroD.Valve(self.r1, self.r2); k = 2e-5; valve.valveCoeff = k; - self.net.time = 0; + net = ct.zeroD.ReactorNet({self.r1, self.r2}); + net.rtol = 1e-12; + net.time = 0; % Skipped until Reactor.inlets and Reactor.outlets are implemented % self.verifyEqual(self.r1.outlets, self.r2.inlets); @@ -97,7 +92,7 @@ function testValve1(self) Y1a = self.r1.phase.Y; Y2a = self.r2.phase.Y; - self.net.advance(0.1); + net.advance(0.1); m1b = self.r1.D * self.r1.V; m2b = self.r2.D * self.r2.V; @@ -114,13 +109,14 @@ function testValve1(self) function testValve2(self) self.makeReactors('P1', ct.OneAtm); - self.net.rtol = 1e-11; + net.rtol = 1e-11; self.r1.energyEnabled = false; self.r2.energyEnabled = false; valve = ct.zeroD.Valve(self.r1, self.r2); k = 2e-5; valve.valveCoeff = k; - self.net.time = 0; + net = ct.zeroD.ReactorNet({self.r1, self.r2}); + net.time = 0; % Skipped until Valve.valveCoeff getter is implemented % self.verifyEqual(valve.ValveCoeff, k); @@ -137,7 +133,7 @@ function testValve2(self) B = k * (P1a / m1a + P2a / m2a); for t = linspace(1e-5, 0.5) - self.net.advance(t); + net.advance(t); m1 = self.r1.D * self.r1.V; m2 = self.r2.D * self.r2.V; @@ -160,7 +156,7 @@ function testValveType1(self) self.makeReactors(); res = ct.zeroD.Reservoir(self.gas1); v = ct.zeroD.Valve(self.r1, res); - n1 = ct.zeroD.ReactorNet(self.r1); + net = ct.zeroD.ReactorNet(self.r1); self.verifyTrue(startsWith(self.r1.name, sprintf('%s_', self.r1.type))); self.verifyTrue(startsWith(res.name, sprintf('%s_', res.type))); @@ -174,7 +170,7 @@ function testValveType2(self) self.makeReactors(); res = ct.zeroD.Reservoir(self.gas1); v = ct.zeroD.Valve(res, self.r1); - n1 = ct.zeroD.ReactorNet(self.r1); + net = ct.zeroD.ReactorNet(self.r1); self.verifyTrue(startsWith(self.r1.name, sprintf('%s_', self.r1.type))); self.verifyTrue(startsWith(res.name, sprintf('%s_', res.type))); diff --git a/test/matlab/ctTestFlowReactor2.m b/test/matlab/ctTestFlowReactor2.m index 339a51e801..2839084386 100644 --- a/test/matlab/ctTestFlowReactor2.m +++ b/test/matlab/ctTestFlowReactor2.m @@ -71,7 +71,7 @@ function testMixedReactorTypes(self) self.net = ct.zeroD.ReactorNet({r1, r2}); catch ME self.verifySubstring(ME.identifier, 'Cantera:ctError'); - self.verifySubstring(ME.message, 'Cannot mix Reactor types'); + self.verifySubstring(ME.message, 'FlowReactors must be used alone'); end end diff --git a/test/matlab/ctTestReactor.m b/test/matlab/ctTestReactor.m index 5756373a18..508950ca75 100644 --- a/test/matlab/ctTestReactor.m +++ b/test/matlab/ctTestReactor.m @@ -5,7 +5,6 @@ gas2 r1 r2 - net w end @@ -32,16 +31,12 @@ function makeReactors(self, arg) self.gas1.TPX = {arg.T1, arg.P1, arg.X1}; self.r1 = ct.zeroD.Reactor(self.gas1); - if arg.nr == 1 - self.net = ct.zeroD.ReactorNet(self.r1); - elseif arg.nr >= 2 + if arg.nr == 2 self.gas2 = ct.Solution('h2o2.yaml', '', 'none'); self.gas2.TPX = {arg.T2, arg.P2, arg.X2}; self.r2 = ct.zeroD.Reactor(self.gas2); self.r2.energyEnabled = true; - self.net = ct.zeroD.ReactorNet({self.r1, self.r2}); end - self.verifyEqual(self.net.time, 0, 'AbsTol', self.atol); end function addWall(self, arg) @@ -78,12 +73,14 @@ function testTypes(self) function testIndependentVariable(self) self.makeReactors('nr', 1); - self.verifyEqual(self.net.time, 0.0, 'AbsTol', self.atol); + net = ct.zeroD.ReactorNet(self.r1); + self.verifyEqual(net.time, 0.0, 'AbsTol', self.atol); end function testDisjoint(self) self.makeReactors('T1', 300, 'P1', 101325, 'T2', 500, 'P2', 300000); - self.net.advance(1.0); + net = ct.zeroD.ReactorNet({self.r1, self.r2}); + net.advance(1.0); self.verifyEqual(self.gas1.T, 300, 'AbsTol', self.atol); self.verifyEqual(self.gas2.T, 500, 'AbsTol', self.atol); @@ -93,21 +90,22 @@ function testDisjoint(self) function testTimeStepping(self) self.makeReactors(); + net = ct.zeroD.ReactorNet({self.r1, self.r2}); tStart = 0.3; tEnd = 10.0; dtMax = 0.07; t = tStart; - self.net.maxTimeStep = dtMax; - self.net.time = tStart; - self.verifyEqual(self.net.time, tStart); + net.maxTimeStep = dtMax; + net.time = tStart; + self.verifyEqual(net.time, tStart); while t < tEnd tPrev = t; - t = self.net.step; + t = net.step; self.verifyLessThanOrEqual(t - tPrev, 1.0001 * dtMax); - self.verifyEqual(t, self.net.time, 'AbsTol', self.atol); + self.verifyEqual(t, net.time, 'AbsTol', self.atol); end end @@ -138,9 +136,9 @@ function testWall3(self) function testEqualizePressure(self) self.makeReactors('P1', 101325, 'P2', 300000); self.addWall('K', 0.1, 'A', 1.0); - - self.net.advance(1.0); - self.verifyEqual(self.net.time, 1.0, 'AbsTol', self.atol); + net = ct.zeroD.ReactorNet({self.r1, self.r2}); + net.advance(1.0); + self.verifyEqual(net.time, 1.0, 'AbsTol', self.atol); self.verifyEqual(self.r1.phase.P, self.r2.phase.P, 'RelTol', self.rtol); self.verifyNotEqual(self.r1.T, self.r2.T); end @@ -148,9 +146,10 @@ function testEqualizePressure(self) function testHeatTransfer1(self) self.makeReactors('T1', 300, 'T2', 1000); self.addWall('U', 500, 'A', 1.0); + net = ct.zeroD.ReactorNet({self.r1, self.r2}); - self.net.advance(10.0); - self.verifyEqual(self.net.time, 10.0, 'RelTol', self.atol); + net.advance(10.0); + self.verifyEqual(net.time, 10.0, 'RelTol', self.atol); self.verifyEqual(self.r1.T, self.r2.T, 'RelTol', 5e-7); self.verifyNotEqual(self.r1.P, self.r2.P); end @@ -159,12 +158,14 @@ function testHeatTransfer2(self) self.assumeFail('Skipped until heatTransferCoeff getter is implemented'); self.makeReactors('T1', 300, 'T2', 1000); self.addWall('U', 200, 'A', 1.0); + net = ct.zeroD.ReactorNet({self.r1, self.r2}); - self.net.advance(1.0); + net.advance(1.0); T1a = self.r1.T; T2a = self.r2.T; self.makeReactors('T1', 300, 'T2', 1000); + net = ct.zeroD.ReactorNet({self.r1, self.r2}); self.r1.V = 0.25; self.r2.V = 0.25; self.addWall('U', 100, 'A', 0.5); @@ -172,7 +173,7 @@ function testHeatTransfer2(self) Qdot1 = self.w.heatTransferCoeff * self.w.area * (self.r1.T - self.r2.T); self.verifyEqual(Qdot1, self.w.heatRate, 'RelTol', self.rtol); - self.net.advance(1.0); + net.advance(1.0); Qdot2 = self.w.heatTransferCoeff * self.w.area * (self.r1.T - self.r2.T); self.verifyEqual(Qdot2, self.w.heatRate, 'RelTol', self.rtol); T1b = self.r1.T; @@ -188,18 +189,18 @@ function testTolerances(self) T0 = 1100; X0 = 'H2:1.0, O2:0.5, AR:8.0'; self.makeReactors('nr', 1, 'T1', T0, 'P1', P0, 'X1', X0); - self.net.rtol = rtol; - self.net.atol = atol; - - self.verifyEqual(self.net.rtol, rtol); - self.verifyEqual(self.net.atol, atol); + net = ct.zeroD.ReactorNet({self.r1}); + net.rtol = rtol; + net.atol = atol; + self.verifyEqual(net.rtol, rtol); + self.verifyEqual(net.atol, atol); tEnd = 1.0; nSteps = 0; t = 0; while t < tEnd - t = self.net.step(); + t = net.step(); nSteps = nSteps + 1; end @@ -216,10 +217,11 @@ function testTolerances(self) function testAdvanceReverse(self) self.makeReactors('nr', 1); - self.net.advance(0.1); + net = ct.zeroD.ReactorNet(self.r1); + net.advance(0.1); try - self.net.advance(0.09); + net.advance(0.09); catch ME self.verifySubstring(ME.identifier, 'Cantera:ctError'); self.verifySubstring(ME.message, 'backwards in time'); @@ -231,8 +233,9 @@ function testEquilibriumUV(self) T0 = 1100; X0 = 'H2:1.0, O2:0.5, AR:8.0'; self.makeReactors('nr', 1, 'T1', T0, 'P1', P0, 'X1', X0); + net = ct.zeroD.ReactorNet({self.r1}); - self.net.advance(1.0); + net.advance(1.0); gas = ct.Solution('h2o2.yaml', '', 'none'); gas.TPX = {T0, P0, X0}; @@ -254,9 +257,9 @@ function testEquilibriumHP(self) self.gas1.TPX = {T0, P0, X0}; self.r1 = ct.zeroD.IdealGasConstPressureReactor(self.gas1); - self.net = ct.zeroD.ReactorNet(self.r1); - self.net.time = 0.0; - self.net.advance(1.0); + net = ct.zeroD.ReactorNet(self.r1); + net.time = 0.0; + net.advance(1.0); self.gas2 = ct.Solution('h2o2.yaml', '', 'none'); self.gas2.TPX = {T0, P0, X0}; @@ -283,12 +286,13 @@ function testWallVelocity(self) v = ct.Func1('tabulated-linear', [0.0, 1.0, 2.0], [0.0, 1.0, 0.0]); self.w.velocity = v; - self.net.advance(1.0); + net = ct.zeroD.ReactorNet({self.r1, self.r2}); + net.advance(1.0); self.verifyEqual(self.w.velocity, v(1.0), 'AbsTol', self.atol); self.verifyEqual(self.w.expansionRate, 1.0 * A, 'AbsTol', self.atol); - self.net.advance(2.0); + net.advance(2.0); self.verifyEqual(selfw.expansionRate, 0.0, 'AbsTol', self.atol); self.verifyEqual(self.r1.V, V1 + 1.0 * A, 'RelTol', self.rtol); @@ -299,7 +303,8 @@ function testDisableEnergy(self) self.makeReactors('T1', 500); self.r1.energyEnabled = false; self.addWall('A', 1.0, 'U', 2500); - self.net.advance(11.0); + net = ct.zeroD.ReactorNet({self.r1, self.r2}); + net.advance(11.0); self.verifyEqual(self.r1.T, 500, 'RelTol', self.rtol); self.verifyEqual(self.r2.T, 500, 'RelTol', self.rtol); @@ -308,7 +313,8 @@ function testDisableEnergy(self) function testDisableChemistry(self) self.makeReactors('T1', 1000, 'nr', 1, 'X1', 'H2:2.0,O2:1.0'); self.r1.chemistryEnabled = false; - self.net.advance(11.0); + net = ct.zeroD.ReactorNet({self.r1}); + net.advance(11.0); x1 = self.r1.phase.X(self.r1.phase.speciesIndex('H2')); x2 = self.r1.phase.X(self.r1.phase.speciesIndex('O2')); @@ -332,7 +338,8 @@ function testHeatFluxFunc(self) self.w.heatFlux = f; Q = 0.3 * 60000; - self.net.advance(1.1); + net = ct.zeroD.ReactorNet({self.r1, self.r2}); + net.advance(1.1); self.verifyEqual(self.w.heatFlux, f(1.1), 'RelTol', self.rtol); U1b = self.r1.V * self.r1.D * self.r1.phase.U; U2b = self.r2.V * self.r2.D * self.r2.phase.U; diff --git a/test/python/test_reactor.py b/test/python/test_reactor.py index f1627f6d46..23d4d240ac 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -26,21 +26,17 @@ class TestReactor: reactorClass = ct.Reactor def make_reactors(self, clone=True, n_reactors=2, - T1=300, P1=101325, X1='O2:1.0', - T2=300, P2=101325, X2='O2:1.0'): + T1=300, P1=ct.one_atm, X1='O2:1.0', + T2=300, P2=ct.one_atm, X2='O2:1.0'): self.gas1 = ct.Solution('h2o2.yaml', transport_model=None) self.gas1.TPX = T1, P1, X1 self.r1 = self.reactorClass(self.gas1, clone=clone) - if n_reactors == 1: - self.net = ct.ReactorNet([self.r1]) - else: + if n_reactors == 2: self.gas1.TPX = T2, P2, X2 # Second reactor must use cloned Solution self.r2 = self.reactorClass(self.gas1, clone=True) - self.net = ct.ReactorNet([self.r1, self.r2]) - assert self.net.initial_time == 0. def add_wall(self, **kwargs): self.w = ct.Wall(self.r1, self.r2, **kwargs) @@ -48,9 +44,10 @@ def add_wall(self, **kwargs): def test_verbose(self): self.make_reactors(clone=False, n_reactors=1) - assert not self.net.verbose - self.net.verbose = True - assert self.net.verbose + net = ct.ReactorNet([self.r1]) + assert not net.verbose + net.verbose = True + assert net.verbose def test_volume(self): g = ct.Solution('h2o2.yaml', transport_model=None) @@ -62,6 +59,7 @@ def test_volume(self): def test_names(self): self.make_reactors() + net = ct.ReactorNet([self.r1, self.r2]) pattern = re.compile(r'(\d+)') digits1 = pattern.search(self.r1.name).group(0) @@ -78,37 +76,41 @@ def test_types(self): def test_component_index(self): self.make_reactors(n_reactors=1) - self.net.step() + net = ct.ReactorNet([self.r1]) + net.step() - N0 = self.net.n_vars - self.gas1.n_species + N0 = net.n_vars - self.gas1.n_species for i, name in enumerate(self.gas1.species_names): assert i + N0 == self.r1.component_index(name) def test_component_names(self): self.make_reactors(n_reactors=2) - self.net.initialize() - N = self.net.n_vars // 2 + net = ct.ReactorNet([self.r1, self.r2]) + net.initialize() + N = net.n_vars // 2 for i in range(N): assert self.r1.component_index(self.r1.component_name(i)) == i - assert (self.net.component_name(i) + assert (net.component_name(i) == '{}: {}'.format(self.r1.name, self.r1.component_name(i))) - assert (self.net.component_name(N+i) + assert (net.component_name(N+i) == '{}: {}'.format(self.r2.name, self.r2.component_name(i))) def test_independent_variable(self): self.make_reactors(clone=False, n_reactors=1) + net = ct.ReactorNet([self.r1]) with pytest.raises(ct.CanteraError, match="independent variable"): - self.net.distance + net.distance - assert self.net.time == 0.0 + assert net.time == 0.0 def test_disjoint(self): T1, P1 = 300, 101325 T2, P2 = 500, 300000 self.make_reactors(T1=T1, T2=T2, P1=P1, P2=P2) - self.net.advance(1.0) + net = ct.ReactorNet([self.r1, self.r2]) + net.advance(1.0) # Nothing should change from the initial condition assert T1 == approx(self.r1.phase.T) @@ -120,23 +122,25 @@ def test_derivative(self): T1, P1 = 300, 101325 self.make_reactors(n_reactors=1, T1=T1, P1=P1) - self.net.advance(1.0) + net = ct.ReactorNet([self.r1]) + net.advance(1.0) # compare cvode derivative to numerical derivative - dydt = self.net.get_derivative(1) - dt = -self.net.time - dy = -self.net.get_state() - self.net.step() - dt += self.net.time - dy += self.net.get_state() - for i in range(self.net.n_vars): + dydt = net.get_derivative(1) + dt = -net.time + dy = -net.get_state() + net.step() + dt += net.time + dy += net.get_state() + for i in range(net.n_vars): assert dydt[i] == approx(dy[i]/dt) def test_finite_difference_jacobian(self): self.make_reactors(n_reactors=1, T1=900, P1=101325, X1="H2:0.4, O2:0.4, N2:0.2") + net = ct.ReactorNet([self.r1]) kH2 = self.gas1.species_index("H2") while self.r1.phase.X[kH2] > 0.3: - self.net.step() + net.step() J = self.r1.finite_difference_jacobian assert J.shape == (self.r1.n_vars, self.r1.n_vars) @@ -172,39 +176,40 @@ def test_finite_difference_jacobian(self): def test_timestepping(self): self.make_reactors() + net = ct.ReactorNet([self.r1, self.r2]) tStart = 0.3 tEnd = 10.0 dt_max = 0.07 t = tStart - self.net.max_time_step = dt_max - assert self.net.max_time_step == dt_max - self.net.initial_time = tStart - assert self.net.initial_time == tStart - assert self.net.time == approx(tStart) - + net.max_time_step = dt_max + assert net.max_time_step == dt_max + net.initial_time = tStart + assert net.initial_time == tStart + assert net.time == approx(tStart) while t < tEnd: tPrev = t - t = self.net.step() + t = net.step() assert t - tPrev <= 1.0001 * dt_max - assert t == approx(self.net.time) + assert t == approx(net.time) def test_maxsteps(self): self.make_reactors() + net = ct.ReactorNet([self.r1, self.r2]) # set the up a case where we can't take # enough time-steps to reach the endtime max_steps = 10 max_step_size = 1e-07 - self.net.initial_time = 0. - self.net.max_time_step = max_step_size - self.net.max_steps = max_steps + net.initial_time = 0. + net.max_time_step = max_step_size + net.max_steps = max_steps with pytest.raises( ct.CanteraError, match='Maximum number of timesteps'): - self.net.advance(1e-04) - assert self.net.time <= max_steps * max_step_size - assert self.net.max_steps == max_steps + net.advance(1e-04) + assert net.time <= max_steps * max_step_size + assert net.max_steps == max_steps def test_wall_type1(self): self.make_reactors(P1=101325, P2=300000) @@ -244,9 +249,10 @@ def test_equalize_pressure(self): assert self.r1.walls[0] == self.w assert self.r2.walls[0] == self.w - self.net.advance(1.0) + net = ct.ReactorNet([self.r1, self.r2]) + net.advance(1.0) - assert self.net.time == approx(1.0) + assert net.time == approx(1.0) assert self.r1.phase.P == approx(self.r2.phase.P) assert self.r1.T != approx(self.r2.T) @@ -256,18 +262,19 @@ def integrate(atol, rtol): T0 = 1100 X0 = 'H2:1.0, O2:0.5, AR:8.0' self.make_reactors(n_reactors=1, T1=T0, P1=P0, X1=X0) - self.net.rtol = rtol - self.net.atol = atol + net = ct.ReactorNet([self.r1]) + net.rtol = rtol + net.atol = atol - assert self.net.rtol == rtol - assert self.net.atol == atol + assert net.rtol == rtol + assert net.atol == atol tEnd = 1.0 nSteps = 0 t = 0 while t < tEnd: - t = self.net.step() + t = net.step() nSteps += 1 return nSteps @@ -283,22 +290,22 @@ def test_advance_limits(self): T0 = 1100 X0 = 'H2:1.0, O2:0.5, AR:8.0' self.make_reactors(n_reactors=1, T1=T0, P1=P0, X1=X0) + net = ct.ReactorNet([self.r1]) limit_H2 = .01 - ix = self.net.global_component_index('H2', 0) + ix = net.global_component_index('H2', 0) self.r1.set_advance_limit('H2', limit_H2) - assert self.net.advance_limits[ix] == limit_H2 + assert net.advance_limits[ix] == limit_H2 self.r1.set_advance_limit('H2', None) - assert self.net.advance_limits[ix] == -1. + assert net.advance_limits[ix] == -1. self.r1.set_advance_limit('H2', limit_H2) - self.net.advance_limits = None - assert self.net.advance_limits[ix] == -1. - + net.advance_limits = None + assert net.advance_limits[ix] == -1. self.r1.set_advance_limit('H2', limit_H2) - self.net.advance_limits = 0 * self.net.advance_limits - 1. - assert self.net.advance_limits[ix] == -1. + net.advance_limits = 0 * net.advance_limits - 1. + assert net.advance_limits[ix] == -1. def test_advance_with_limits(self): def integrate(limit_H2 = None, apply=True): @@ -306,11 +313,11 @@ def integrate(limit_H2 = None, apply=True): T0 = 1100 X0 = 'H2:1.0, O2:0.5, AR:8.0' self.make_reactors(n_reactors=1, T1=T0, P1=P0, X1=X0) + net = ct.ReactorNet([self.r1]) if limit_H2 is not None: self.r1.set_advance_limit('H2', limit_H2) - ix = self.net.global_component_index('H2', 0) - assert self.net.advance_limits[ix] == limit_H2 - + ix = net.global_component_index('H2', 0) + assert net.advance_limits[ix] == limit_H2 tEnd = 0.1 tStep = 7e-4 nSteps = 0 @@ -318,7 +325,7 @@ def integrate(limit_H2 = None, apply=True): t = tStep while t < tEnd: - t_curr = self.net.advance(t, apply_limit=apply) + t_curr = net.advance(t, apply_limit=apply) nSteps += 1 if t_curr < t: nEvents += 1 @@ -350,20 +357,21 @@ def test_advance_limit_triggers(self): T0 = 1100 X0 = 'H2:1.0, O2:0.5, AR:8.0' self.make_reactors(n_reactors=1, T1=T0, P1=P0, X1=X0) + net = ct.ReactorNet([self.r1]) comp = 'H2' limit = 1e-3 target = 5e-2 - self.net.initialize() - baseline = np.copy(self.net.get_state()) + net.initialize() + baseline = np.copy(net.get_state()) self.r1.set_advance_limit(comp, limit) - ix = self.net.global_component_index(comp, 0) + ix = net.global_component_index(comp, 0) - reached = self.net.advance(target, apply_limit=True) + reached = net.advance(target, apply_limit=True) - assert reached == approx(self.net.time) - assert self.net.time < target - delta = abs(self.net.get_state()[ix] - baseline[ix]) + assert reached == approx(net.time) + assert net.time < target + delta = abs(net.get_state()[ix] - baseline[ix]) assert delta == approx(limit, rel=0.1, abs=1e-10) def test_advance_limit_logging(self, capsys): @@ -378,10 +386,11 @@ def test_advance_limit_logging(self, capsys): target = 5e-2 self.make_reactors(n_reactors=1, T1=T0, P1=P0, X1=X0) - self.net.verbose = True + net = ct.ReactorNet([self.r1]) + net.verbose = True self.r1.set_advance_limit('H2', limit) capsys.readouterr() - self.net.advance(target, apply_limit=True) + net.advance(target, apply_limit=True) out = capsys.readouterr().out assert "Advance limit triggered" in out assert "y_start" in out @@ -390,10 +399,11 @@ def test_advance_limit_logging(self, capsys): assert "limit =" in out self.make_reactors(n_reactors=1, T1=T0, P1=P0, X1=X0) - self.net.verbose = True + net = ct.ReactorNet([self.r1]) + net.verbose = True self.r1.set_advance_limit('H2', limit) capsys.readouterr() - self.net.advance(target, apply_limit=False) + net.advance(target, apply_limit=False) out = capsys.readouterr().out assert "Advance limit triggered" not in out @@ -404,13 +414,16 @@ def test_advance_limit_cleanup_after_failure(self): X0 = 'H2:1.0, O2:0.5, AR:8.0' self.make_reactors(n_reactors=1, T1=T0, P1=P0, X1=X0) self.r1.set_advance_limit('H2', 1e-3) - self.net.max_time_step = 1e-6 - self.net.max_steps = 1 + net = ct.ReactorNet([self.r1]) + net.max_time_step = 1e-6 + net.max_steps = 1 with pytest.raises(ct.CanteraError, match="Maximum number of timesteps"): - self.net.advance(1e-2, apply_limit=True) + net.advance(1e-2, apply_limit=True) + self.make_reactors(n_reactors=1, T1=T0, P1=P0, X1=X0) self.r1.set_advance_limit('H2', 1e-3) - self.net.advance(1e-2, apply_limit=True) + net = ct.ReactorNet([self.r1]) + net.advance(1e-2, apply_limit=True) def test_multicomponent_advance_limits(self): """Most restrictive component determines which limit triggers first""" @@ -418,33 +431,35 @@ def test_multicomponent_advance_limits(self): T0 = 1100 X0 = 'H2:1.0, O2:0.5, AR:8.0' self.make_reactors(n_reactors=1, T1=T0, P1=P0, X1=X0) + net = ct.ReactorNet([self.r1]) limits = {'H2': 3e-5, 'O2': 2e-3} - self.net.initialize() - base = np.copy(self.net.get_state()) + net.initialize() + base = np.copy(net.get_state()) for comp, lim in limits.items(): self.r1.set_advance_limit(comp, lim) - h2_ix = self.net.global_component_index('H2', 0) - o2_ix = self.net.global_component_index('O2', 0) - t_hit_h2 = self.net.advance(5e-2, apply_limit=True) - assert t_hit_h2 == approx(self.net.time) - delta_h2 = abs(self.net.get_state()[h2_ix] - base[h2_ix]) - delta_o2 = abs(self.net.get_state()[o2_ix] - base[o2_ix]) + h2_ix = net.global_component_index('H2', 0) + o2_ix = net.global_component_index('O2', 0) + t_hit_h2 = net.advance(5e-2, apply_limit=True) + assert t_hit_h2 == approx(net.time) + delta_h2 = abs(net.get_state()[h2_ix] - base[h2_ix]) + delta_o2 = abs(net.get_state()[o2_ix] - base[o2_ix]) assert delta_h2 >= 0.8 * limits['H2'] assert delta_o2 < 0.2 * limits['O2'] - base = np.copy(self.net.get_state()) + base = np.copy(net.get_state()) self.r1.set_advance_limit('H2', None) - t_hit_o2 = self.net.advance(t_hit_h2 + 5e-2, apply_limit=True) - assert t_hit_o2 == approx(self.net.time) - delta_o2 = abs(self.net.get_state()[o2_ix] - base[o2_ix]) + t_hit_o2 = net.advance(t_hit_h2 + 5e-2, apply_limit=True) + assert t_hit_o2 == approx(net.time) + delta_o2 = abs(net.get_state()[o2_ix] - base[o2_ix]) assert delta_o2 >= 0.8 * limits['O2'] def test_heat_transfer1(self): # Connected reactors reach thermal equilibrium after some time self.make_reactors(T1=300, T2=1000) self.add_wall(U=500, A=1.0) + net = ct.ReactorNet([self.r1, self.r2]) - self.net.advance(10.0) - assert self.net.time == approx(10.0) + net.advance(10.0) + assert net.time == approx(10.0) assert self.r1.T == approx(self.r2.T, rel=5e-7) assert self.r1.phase.P != approx(self.r2.phase.P) @@ -456,16 +471,18 @@ def test_advance_limits_invalid(self): def test_advance_reverse(self): self.make_reactors(n_reactors=1) - self.net.advance(0.1) + net = ct.ReactorNet([self.r1]) + net.advance(0.1) with pytest.raises(ct.CanteraError, match="backwards in time"): - self.net.advance(0.09) + net.advance(0.09) def test_heat_transfer2(self): # Result should be the same if (m * cp) / (U * A) is held constant self.make_reactors(T1=300, T2=1000) self.add_wall(U=200, A=1.0) + net = ct.ReactorNet([self.r1, self.r2]) - self.net.advance(1.0) + net.advance(1.0) T1a = self.r1.T T2a = self.r2.T @@ -473,10 +490,11 @@ def test_heat_transfer2(self): self.r1.volume = 0.25 self.r2.volume = 0.25 w = self.add_wall(U=100, A=0.5) + net = ct.ReactorNet([self.r1, self.r2]) assert (w.heat_transfer_coeff * w.area * (self.r1.T - self.r2.T) == approx(w.heat_rate)) - self.net.advance(1.0) + net.advance(1.0) assert (w.heat_transfer_coeff * w.area * (self.r1.T - self.r2.T) == approx(w.heat_rate)) T1b = self.r1.T @@ -493,8 +511,9 @@ def test_equilibrium_UV(self): T0 = 1100 X0 = 'H2:1.0, O2:0.5, AR:8.0' self.make_reactors(n_reactors=1, T1=T0, P1=P0, X1=X0) + net = ct.ReactorNet([self.r1]) - self.net.advance(1.0) + net.advance(1.0) gas = ct.Solution('h2o2.yaml', transport_model=None) gas.TPX = T0, P0, X0 @@ -539,14 +558,14 @@ def test_wall_velocity(self): self.r2.volume = V2 self.add_wall(A=A) - v = ct.Tabulated1([0.0, 1.0, 2.0], [0.0, 1.0, 0.0]) - self.w.velocity = v - self.net.advance(1.0) + + net = ct.ReactorNet([self.r1, self.r2]) + net.advance(1.0) assert self.w.velocity == approx(v(1.0)) assert self.w.expansion_rate == approx(1.0 * A, rel=1e-7) - self.net.advance(2.0) + net.advance(2.0) assert self.w.expansion_rate == approx(0.0, rel=1e-7) assert self.r1.volume == approx(V1 + 1.0 * A, rel=1e-7) @@ -556,7 +575,8 @@ def test_disable_energy(self): self.make_reactors(T1=500) self.r1.energy_enabled = False self.add_wall(A=1.0, U=2500) - self.net.advance(11.0) + net = ct.ReactorNet([self.r1, self.r2]) + net.advance(11.0) assert self.r1.T == approx(500) assert self.r2.T == approx(500) @@ -565,7 +585,8 @@ def test_disable_chemistry(self): self.make_reactors(T1=1000, n_reactors=1, X1='H2:2.0,O2:1.0') self.r1.chemistry_enabled = False - self.net.advance(11.0) + net = ct.ReactorNet([self.r1]) + net.advance(11.0) assert self.r1.T == approx(1000) assert self.r1.phase.X[self.r1.phase.species_index('H2')] == approx(2.0/3.0) @@ -586,7 +607,8 @@ def test_heat_flux_func(self): self.w.heat_flux = hfunc Q = 0.3 * 60000 - self.net.advance(1.1) + net = ct.ReactorNet([self.r1, self.r2]) + net.advance(1.1) assert self.w.heat_flux == hfunc(1.1) U1b = self.r1.volume * self.r1.density * self.r1.phase.u U2b = self.r2.volume * self.r2.density * self.r2.phase.u @@ -630,19 +652,20 @@ def mdot(t): ma = self.r1.volume * self.r1.density Ya = self.r1.Y - self.net.rtol = 1e-11 - self.net.max_time_step = 0.05 + net = ct.ReactorNet([self.r1]) + net.rtol = 1e-11 + net.max_time_step = 0.05 - self.net.advance(0.1) + net.advance(0.1) assert mfc.mass_flow_rate == approx(0.) - self.net.advance(0.3) + net.advance(0.3) assert mfc.mass_flow_rate == approx(0.04) - self.net.advance(1.0) + net.advance(1.0) assert mfc.mass_flow_rate == approx(0.08) - self.net.advance(1.2) + net.advance(1.2) assert mfc.mass_flow_rate == approx(0.) - self.net.advance(2.5) + net.advance(2.5) mb = self.r1.volume * self.r1.density Yb = self.r1.Y @@ -665,15 +688,15 @@ def test_mass_flow_controller_errors(self): mfc = ct.MassFlowController(self.r1, self.r2) mfc.mass_flow_rate = lambda t: eggs + with pytest.raises(Exception, match='eggs'): - self.net.step() + net = ct.ReactorNet([self.r1, self.r2]) with pytest.raises(NotImplementedError): mfc.pressure_function = lambda p: p**2 def test_valve1(self): self.make_reactors(P1=10*ct.one_atm, X1='AR:1.0', X2='O2:1.0') - self.net.rtol = 1e-12 valve = ct.Valve(self.r1, self.r2) k = 2e-5 valve.valve_coeff = k @@ -682,7 +705,10 @@ def test_valve1(self): assert valve.valve_coeff == k assert self.r1.energy_enabled assert self.r2.energy_enabled - self.net.initialize() + + net = ct.ReactorNet([self.r1, self.r2]) + net.rtol = 1e-12 + net.initialize() assert (self.r1.phase.P - self.r2.phase.P) * k == approx( valve.mass_flow_rate) @@ -691,7 +717,7 @@ def test_valve1(self): Y1a = self.r1.phase.Y Y2a = self.r2.phase.Y - self.net.advance(0.1) + net.advance(0.1) m1b = self.r1.phase.density * self.r1.volume m2b = self.r2.phase.density * self.r2.volume @@ -709,7 +735,6 @@ def test_valve2(self): # (constant T) we can compare with an analytical solution for # the mass of each reactor as a function of time self.make_reactors(P1=10*ct.one_atm) - self.net.rtol = 1e-11 self.r1.energy_enabled = False self.r2.energy_enabled = False valve = ct.Valve(self.r1, self.r2) @@ -729,8 +754,11 @@ def test_valve2(self): A = k * P1a * (1 + m2a/m1a) B = k * (P1a/m1a + P2a/m2a) + net = ct.ReactorNet([self.r1, self.r2]) + net.rtol = 1e-11 + for t in np.linspace(1e-5, 0.5): - self.net.advance(t) + net.advance(t) m1 = self.r1.phase.density * self.r1.volume m2 = self.r2.phase.density * self.r2.volume assert m2 == approx((m2a - A/B) * np.exp(-B * t) + A/B) @@ -742,8 +770,6 @@ def test_valve3(self): # and flow rate. self.make_reactors(P1=10*ct.one_atm, X1='AR:0.5, O2:0.5', X2='O2:1.0') - self.net.rtol = 1e-12 - self.net.atol = 1e-20 valve = ct.Valve(self.r1, self.r2) mdot = lambda dP: 5e-3 * np.sqrt(dP) if dP > 0 else 0.0 valve.pressure_function = mdot @@ -758,9 +784,13 @@ def speciesMass(k): mO2 = speciesMass(kO2) mAr = speciesMass(kAr) + net = ct.ReactorNet([self.r1, self.r2]) + net.rtol = 1e-12 + net.atol = 1e-20 + t = 0 while t < 1.0: - t = self.net.step() + t = net.step() p1 = self.r1.phase.P p2 = self.r2.phase.P assert mdot(p1-p2) == approx(valve.mass_flow_rate) @@ -771,7 +801,6 @@ def speciesMass(k): def test_valve_timing(self): # test timed valve self.make_reactors(P1=10*ct.one_atm, X1='AR:1.0', X2='O2:1.0') - self.net.rtol = 1e-12 valve = ct.Valve(self.r1, self.r2) k = 2e-5 valve.valve_coeff = k @@ -779,19 +808,23 @@ def test_valve_timing(self): delta_p = lambda: self.r1.phase.P - self.r2.phase.P mdot = lambda: valve.valve_coeff * (self.r1.phase.P - self.r2.phase.P) - self.net.initialize() + + net = ct.ReactorNet([self.r1, self.r2]) + net.rtol = 1e-12 + net.initialize() + assert valve.time_function == 0.0 assert valve.pressure_function == approx(delta_p()) assert valve.mass_flow_rate == 0.0 - self.net.advance(0.01) + net.advance(0.01) assert valve.time_function == 0.0 assert valve.pressure_function == approx(delta_p()) assert valve.mass_flow_rate == 0.0 - self.net.advance(0.01 + 1e-9) + net.advance(0.01 + 1e-9) assert valve.time_function == 1.0 assert valve.pressure_function == approx(delta_p()) assert valve.mass_flow_rate == approx(mdot()) - self.net.advance(0.02) + net.advance(0.02) assert valve.time_function == 1.0 assert valve.pressure_function == approx(delta_p()) assert valve.mass_flow_rate == approx(mdot()) @@ -800,7 +833,7 @@ def test_valve_type1(self): self.make_reactors() res = ct.Reservoir(self.gas1) v = ct.Valve(self.r1, res) - ct.ReactorNet([self.r1]) # assigns default names + net = ct.ReactorNet([self.r1]) # assigns default names assert self.r1.name.startswith(f"{self.r1.type}_") # default name assert res.name.startswith(f"{res.type}_") # default name assert v.type == "Valve" @@ -812,7 +845,7 @@ def test_valve_type2(self): self.make_reactors() res = ct.Reservoir(self.gas1) ct.Valve(res, self.r1) - ct.ReactorNet([self.r1]) # assigns default names + net = ct.ReactorNet([self.r1]) # assigns default names assert self.r1.name.startswith(f"{self.r1.type}_") # default name assert res.name.startswith(f"{res.type}_") # default name @@ -837,9 +870,10 @@ def test_pressure_controller1(self): pc.device_coefficient = 1e-5 assert pc.pressure_coeff == 1e-5 + net = ct.ReactorNet([self.r1]) t = 0 while t < 1.0: - t = self.net.step() + t = net.step() assert mdot(t) == approx(mfc.mass_flow_rate) dP = self.r1.phase.P - outlet_reservoir.phase.P assert mdot(t) + 1e-5 * dP == approx(pc.mass_flow_rate) @@ -863,9 +897,10 @@ def test_pressure_controller2(self): pc.pressure_function = pfunc assert pc.pressure_coeff == 1. + net = ct.ReactorNet([self.r1]) t = 0 while t < 1.0: - t = self.net.step() + t = net.step() assert mdot(t) == approx(mfc.mass_flow_rate) dP = self.r1.phase.P - outlet_reservoir.phase.P assert mdot(t) + pfunc(dP) == approx(pc.mass_flow_rate) @@ -902,30 +937,32 @@ def test_pressure_controller_errors(self): def test_set_initial_time(self): self.make_reactors(P1=10*ct.one_atm, X1='AR:1.0', X2='O2:1.0') - self.net.rtol = 1e-12 valve = ct.Valve(self.r1, self.r2) pfunc_a = lambda dP: 5e-3 * np.sqrt(dP) if dP > 0 else 0.0 valve.pressure_function = pfunc_a + net = ct.ReactorNet([self.r1, self.r2]) + net.rtol = 1e-12 t0 = 0.0 tf = t0 + 0.5 - self.net.advance(tf) - assert self.net.time == approx(tf) + net.advance(tf) + assert net.time == approx(tf) p1a = self.r1.phase.P p2a = self.r2.phase.P assert valve.pressure_function == approx(pfunc_a(p1a - p2a)) self.make_reactors(P1=10*ct.one_atm, X1='AR:1.0', X2='O2:1.0') - self.net.rtol = 1e-12 valve = ct.Valve(self.r1, self.r2) pfunc_b = lambda dP: 5e-3 * np.sqrt(dP) if dP > 0 else 0.0 valve.pressure_function = pfunc_b t0 = 0.2 - self.net.initial_time = t0 + net = ct.ReactorNet([self.r1, self.r2]) + net.rtol = 1e-12 + net.initial_time = t0 tf = t0 + 0.5 - self.net.advance(tf) - assert self.net.time == approx(tf) + net.advance(tf) + assert net.time == approx(tf) p1b = self.r1.phase.P p2b = self.r2.phase.P assert valve.pressure_function == approx(pfunc_b(p1b - p2b)) @@ -936,7 +973,8 @@ def test_set_initial_time(self): def test_reinitialize(self, allow_deprecated): self.make_reactors(T1=300, T2=1000) self.add_wall(U=200, A=1.0) - self.net.advance(1.0) + net = ct.ReactorNet([self.r1, self.r2]) + net.advance(1.0) T1a = self.r1.T T2a = self.r2.T @@ -949,7 +987,7 @@ def test_reinitialize(self, allow_deprecated): assert self.r1.T == approx(300) assert self.r2.T == approx(1000) - self.net.advance(2.0) + net.advance(2.0) T1b = self.r1.T T2b = self.r2.T @@ -959,7 +997,8 @@ def test_reinitialize(self, allow_deprecated): # TODO: Remove after Cantera 4.0 when syncState is removed def test_syncState_deprecated(self): self.make_reactors(n_reactors=1) - self.net.advance(0.1) + net = ct.ReactorNet([self.r1]) + net.advance(0.1) with pytest.raises(ct.CanteraError): self.r1.syncState() @@ -968,32 +1007,59 @@ def test_reservoir_sync(self): self.gas1.TP = 900, ct.one_atm reservoir = ct.Reservoir(self.gas1) wall = ct.Wall(self.r1, reservoir, U=500) - self.net.advance(1.0) + net = ct.ReactorNet([self.r1]) + net.advance(1.0) assert self.r1.T == approx(872.099, rel=1e-3) reservoir.phase.TP = 700, ct.one_atm - self.net.reinitialize() - self.net.advance(2.0) + net.reinitialize() + net.advance(2.0) assert self.r1.T == approx(747.27, rel=1e-3) + def test_invalid_modification(self): + self.make_reactors() + net = ct.ReactorNet([self.r1, self.r2]) + with pytest.raises(ct.CanteraError, match="Cannot add a wall"): + ct.Wall(self.r1, self.r2, A=2.0) + + res = ct.Reservoir(self.gas1) + with pytest.raises(ct.CanteraError, match="Cannot add an inlet"): + ct.MassFlowController(res, self.r2, mdot=0.1) + + with pytest.raises(ct.CanteraError, match="Cannot add an outlet"): + ct.MassFlowController(self.r2, res, mdot=0.1) + + def test_invalid_multiple_networks(self): + self.make_reactors() + net = ct.ReactorNet([self.r1, self.r2]) + with pytest.raises(ct.CanteraError, match="already part of a ReactorNet"): + net2 = ct.ReactorNet([self.r2]) + + def test_invalid_empty_network(self): + with pytest.raises(ct.CanteraError, match="No reactors in network"): + net = ct.ReactorNet([]) + def test_unpicklable(self): self.make_reactors() + net = ct.ReactorNet([self.r1, self.r2]) import pickle with pytest.raises(NotImplementedError): pickle.dumps(self.r1) with pytest.raises(NotImplementedError): - pickle.dumps(self.net) + pickle.dumps(net) def test_uncopyable(self): self.make_reactors() + net = ct.ReactorNet([self.r1, self.r2]) import copy with pytest.raises(NotImplementedError): copy.copy(self.r1) with pytest.raises(NotImplementedError): - copy.copy(self.net) + copy.copy(net) def test_invalid_property(self): self.make_reactors() - for x in (self.r1, self.net): + net = ct.ReactorNet([self.r1, self.r2]) + for x in (self.r1, net): with pytest.raises(AttributeError): x.foobar = 300 with pytest.raises(AttributeError): @@ -1007,11 +1073,12 @@ def test_bad_kwarg(self): def test_preconditioner_unsupported(self): self.make_reactors() - self.net.preconditioner = ct.AdaptivePreconditioner() + net = ct.ReactorNet([self.r1, self.r2]) + net.preconditioner = ct.AdaptivePreconditioner() # initialize should throw an error because the mass fraction # reactors do not support preconditioning with pytest.raises(ct.CanteraError): - self.net.initialize() + net.initialize() @pytest.mark.skipif(_graphviz is None, reason="graphviz is not installed") def test_draw_reactor(self): @@ -1069,8 +1136,9 @@ def test_draw_reactors_same_name(self): self.make_reactors() self.r1.name = 'Reactor' self.r2.name = 'Reactor' + net = ct.ReactorNet([self.r1, self.r2]) with pytest.raises(AssertionError, match="unique names"): - self.net.draw() + net.draw() @pytest.mark.skipif(_graphviz is None, reason="graphviz is not installed") def test_draw_grouped_reactors(self): @@ -1079,7 +1147,8 @@ def test_draw_grouped_reactors(self): self.r2.name = "Reactor 2" self.r1.group_name = "Group 1" self.r2.group_name = "Group 2" - graph = self.net.draw() + net = ct.ReactorNet([self.r1, self.r2]) + graph = net.draw() expected = ['\tsubgraph "cluster_Group 1" {\n', '\t\t"Reactor 1"\n', '\t\tlabel="Group 1"\n', @@ -1142,7 +1211,8 @@ def test_draw_flow_controller(self): edge_attr={'xlabel': 'MFC'}, name="MFC") ct.PressureController(self.r1, outlet_reservoir, primary=mfc, name="PC") mfc.edge_attr.update({'color': 'purple'}) - self.net.advance_to_steady_state() + net = ct.ReactorNet([self.r1]) + net.advance_to_steady_state() graph = mfc.draw(node_attr={'style': 'filled'}, edge_attr={'style': 'dotted', 'color': 'blue'}) expected = [('\tInlet -> Reactor [label="MFC\\nṁ = 2 kg/s" color=blue ' @@ -1167,10 +1237,11 @@ def test_draw_network(self): ct.PressureController(self.r1, outlet, primary=mfc_hot1, name='pc_h1') ct.PressureController(self.r1, outlet, primary=mfc_hot2, name='pc_h2') ct.PressureController(self.r2, outlet, primary=mfc_cold, name='pc_c') - self.net.advance_to_steady_state() - graph = self.net.draw(mass_flow_attr={'color': 'green'}, - heat_flow_attr={'color': 'orange'}, - print_state=True) + net = ct.ReactorNet([self.r1, self.r2]) + net.advance_to_steady_state() + graph = net.draw(mass_flow_attr={'color': 'green'}, + heat_flow_attr={'color': 'orange'}, + print_state=True) expected = { '\tRC [label="{RC|{T (K)\\n202.18|P (bar)\\n1.013}}" shape=Mrecord]\n', '\tOut [label="{Out|{T (K)\\n200.00|P (bar)\\n1.013}}" shape=Mrecord]\n', @@ -1419,16 +1490,17 @@ def create_reactors(self, add_Q=False, add_mdot=False, add_surf=False): self.r2.volume = 0.2 resGas.TP = T0 - 300, P0 - env = ct.Reservoir(resGas) + env1 = ct.Reservoir(resGas) + env2 = ct.Reservoir(resGas) U = 300 if add_Q else 0 - self.w1 = ct.Wall(self.r1, env, K=1e3, A=0.1, U=U) - self.w2 = ct.Wall(self.r2, env, A=0.1, U=U) + self.w1 = ct.Wall(self.r1, env1, K=1e3, A=0.1, U=U) + self.w2 = ct.Wall(self.r2, env2, A=0.1, U=U) if add_mdot: - mfc1 = ct.MassFlowController(env, self.r1, mdot=0.05) - mfc2 = ct.MassFlowController(env, self.r2, mdot=0.05) + mfc1 = ct.MassFlowController(env1, self.r1, mdot=0.05) + mfc2 = ct.MassFlowController(env2, self.r2, mdot=0.05) if add_surf: self.interface1 = ct.Interface('diamond.yaml', 'diamond_100', @@ -1471,11 +1543,18 @@ def test_component_index(self): def test_component_names(self): self.create_reactors(add_surf=True) - for i in range(self.net1.n_vars): + for i in range(self.r1.n_vars): assert self.r1.component_index(self.r1.component_name(i)) == i assert self.net1.component_name(i) == '{}: {}'.format(self.r1.name, self.r1.component_name(i)) + offset = self.r1.n_vars + for i in range(self.surf1.n_vars): + assert self.surf1.component_index(self.surf1.component_name(i)) == i + name = f'{self.surf1.name}: {self.surf1.component_name(i)}' + assert self.net1.component_name(i + offset) == name + + def integrate(self, surf=False): for t in np.arange(0.5, 50, 1.0): self.net1.advance(t) @@ -1505,6 +1584,11 @@ def test_with_surface_reactions(self): self.net1.rtol = self.net2.rtol = 1e-9 self.integrate(surf=True) + def test_invalid_modification(self): + self.create_reactors(add_surf=True) + with pytest.raises(ct.CanteraError, match="Cannot add a surface"): + ct.ReactorSurface(self.interface1, self.r2, A=0.2) + def test_preconditioner_unsupported(self): self.create_reactors() self.net2.preconditioner = ct.AdaptivePreconditioner() @@ -1919,9 +2003,14 @@ def test_mixed_reactor_types(self): surf, gas = self.import_phases() r1 = ct.FlowReactor(gas) r2 = ct.IdealGasReactor(gas) - with pytest.raises(ct.CanteraError, match="Cannot mix Reactor types"): + with pytest.raises(ct.CanteraError, match="FlowReactors must be used alone"): ct.ReactorNet([r1, r2]) + r1 = ct.FlowReactor(gas) + r2 = ct.IdealGasReactor(gas) + with pytest.raises(ct.CanteraError, match="Cannot mix Reactor types"): + ct.ReactorNet([r2, r1]) + def test_unrecoverable_integrator_errors(self): surf, gas = self.import_phases() @@ -2223,17 +2312,17 @@ def make_reactors(self): self.r2 = ct.IdealGasReactor(self.gas) self.r2.volume = 0.01 - self.net = ct.ReactorNet([self.r1, self.r2]) - def test_coverages(self): self.make_reactors() + surf1 = ct.ReactorSurface(self.interface, self.r1) surf1.coverages = {'c6HH':0.3, 'c6HM':0.7} assert surf1.coverages[0] == approx(0.3) assert surf1.coverages[1] == approx(0.0) assert surf1.coverages[4] == approx(0.7) - self.net.advance(1e-5) + net = ct.ReactorNet([self.r1, self.r2]) + net.advance(1e-5) C_left = surf1.coverages self.make_reactors() @@ -2241,7 +2330,8 @@ def test_coverages(self): surf2.coverages = 'c6HH:0.3, c6HM:0.7' assert surf2.coverages[0] == approx(0.3) assert surf2.coverages[4] == approx(0.7) - self.net.advance(1e-5) + net = ct.ReactorNet([self.r1, self.r2]) + net.advance(1e-5) C_right = surf2.coverages assert sum(C_left) == approx(1.0) @@ -2263,12 +2353,14 @@ def test_coverages_regression1(self, test_data_path): surf1.coverages = C assert surf1.coverages == approx(C) + + net = ct.ReactorNet([self.r1, self.r2]) data = [] test_file = self.test_work_path / "test_coverages_regression1.csv" reference_file = test_data_path / "WallKinetics-coverages-regression1.csv" data = [] for t in np.linspace(1e-6, 1e-3): - self.net.advance(t) + net.advance(t) data.append([t, self.r1.T, self.r1.phase.P, self.r1.mass] + list(self.r1.phase.X) + list(surf1.coverages)) np.savetxt(test_file, data, delimiter=',') @@ -2288,12 +2380,14 @@ def test_coverages_regression2(self, test_data_path): surf.coverages = C assert surf.coverages == approx(C) + + net = ct.ReactorNet([self.r1, self.r2]) data = [] test_file = self.test_work_path / "test_coverages_regression2.csv" reference_file = test_data_path / "WallKinetics-coverages-regression2.csv" data = [] for t in np.linspace(1e-6, 1e-3): - self.net.advance(t) + net.advance(t) data.append([t, self.r1.T, self.r1.phase.P, self.r1.mass] + list(self.r1.phase.X) + list(surf.coverages)) np.savetxt(test_file, data, delimiter=',') @@ -2330,19 +2424,17 @@ def test_incompatible_bulk(self): gas = ct.Solution("h2o2.yaml") r = ct.Reactor(gas, clone=False) rsurf = ct.ReactorSurface(surf, r, clone=False) - net = ct.ReactorNet([r]) with pytest.raises(ct.CanteraError, match="does not have an adjacent phase named 'ohmech'"): - net.initialize() + net = ct.ReactorNet([r]) def test_mismatched_bulk(self): surf = ct.Interface("ptcombust.yaml", "Pt_surf") gas = ct.Solution("ptcombust.yaml", "gas") r = ct.Reactor(gas, clone=False) rsurf = ct.ReactorSurface(surf, [r], clone=False) - net = ct.ReactorNet([r]) with pytest.raises(ct.CanteraError, match="must be the same object"): - net.initialize() + net = ct.ReactorNet([r]) class TestReactorSensitivities: @@ -2372,17 +2464,16 @@ def test_sensitivities2(self): gas2.TPX = 900, 101325, 'H2:0.1, OH:1e-7, O2:0.1, AR:1e-5' r2 = ct.IdealGasReactor(gas2) - net = ct.ReactorNet([r1, r2]) - net.atol_sensitivity = 1e-10 - net.rtol_sensitivity = 1e-8 - surf = ct.ReactorSurface(interface, r1, A=1.5) - C = np.zeros(interface.n_species) C[0] = 0.3 C[4] = 0.7 - surf.coverages = C + + net = ct.ReactorNet([r1, r2]) + net.atol_sensitivity = 1e-10 + net.rtol_sensitivity = 1e-8 + surf.add_sensitivity_reaction(2) r2.add_sensitivity_reaction(18) @@ -2655,9 +2746,8 @@ def test_solution_reuse(): r1 = ct.Reactor(gas, clone=False) r2 = ct.ConstPressureReactor(gas, clone=False) - net = ct.ReactorNet([r1, r2]) with pytest.raises(ct.CanteraError, match="using the same Solution object"): - net.initialize() + net = ct.ReactorNet([r1, r2]) def test_interface_reuse(): @@ -2665,9 +2755,8 @@ def test_interface_reuse(): r1 = ct.Reactor(surf.adjacent["gas"], clone=False) rsurf1 = ct.ReactorSurface(surf, r1, clone=False) rsurf2 = ct.ReactorSurface(surf, r1, clone=False) - net = ct.ReactorNet([r1]) with pytest.raises(ct.CanteraError, match="using the same Solution object"): - net.initialize() + net = ct.ReactorNet([r1]) class TestCombustor: @@ -2974,6 +3063,7 @@ def setup_advance_converages_data(request): interface_phase = 'Pt_surf' request.cls.surf = ct.Interface(mechanism_file, interface_phase) request.cls.gas = request.cls.surf.adjacent["gas"] + request.cls.gas.equilibrate("HP") # get an interesting gas composition @pytest.mark.usefixtures('setup_advance_converages_data') class TestAdvanceCoverages: @@ -3006,6 +3096,19 @@ def test_different_tolerances(self): assert cov == approx(self.surf.coverages) assert any(cov != self.surf.coverages) + def test_steady_state(self): + self.surf.coverages = 'O(S):0.1, PT(S):0.5, H(S):0.4' + self.gas.TP = self.surf.TP + self.surf.advance_coverages(100.0) + cov_time = self.surf.coverages + + # Solution should be independent of initial surface coverages + self.surf.coverages = 'CO(S):0.2, PT(S):0.3, H(S):0.5' + self.surf.advance_coverages_to_steady_state() + cov_steady = self.surf.coverages + + assert cov_time == approx(cov_steady, rel=1e-8, abs=1e-12) + @pytest.fixture(scope='function') def setup_extensible_reactor_data(request): @@ -3409,16 +3512,10 @@ def test_const_pressure(self, reactor_class): ct.MassFlowController(r, downstream, mdot=160) net = ct.ReactorNet([r]) - if "Mole" not in reactor_class.__name__: - net.solve_steady() + net.solve_steady() - assert r.mass == approx(m0) - assert r.phase.T == approx(2407.35011) - else: - # Expected to raise until https://github.com/Cantera/enhancements/issues/234 - # is implemented - with pytest.raises(ct.CanteraError, match="See https://github.com"): - net.solve_steady() + assert r.mass == approx(m0) + assert r.phase.T == approx(2407.35011) @pytest.mark.parametrize("reactor_class", [ct.IdealGasReactor, ct.IdealGasMoleReactor]) @@ -3466,8 +3563,6 @@ def test_multiple_reactors(self): [ct.ConstPressureReactor, ct.IdealGasConstPressureReactor, ct.ConstPressureMoleReactor, ct.IdealGasConstPressureMoleReactor]) def test_steady_surface_disabled(self, reactor_class): - # This case demonstrated in this test should work after - # https://github.com/Cantera/enhancements/issues/234 is implemented surf = ct.Interface("methane_pox_on_pt.yaml", "Pt_surf") gas = surf.adjacent["gas"] gas.set_equivalence_ratio(0.22, "CH4:1.0", "O2:1.0, AR:3.76") @@ -3477,19 +3572,17 @@ def test_steady_surface_disabled(self, reactor_class): gas.equilibrate("HP") downstream = ct.Reservoir(gas) r = reactor_class(gas, volume=1e-2) - ct.ReactorSurface(surf, r, A=0.1) + rsurf = ct.ReactorSurface(surf, r, A=0.1) mdot = 0.2 inlet = ct.MassFlowController(upstream, r, mdot=mdot) ct.MassFlowController(r, downstream, mdot=mdot) net = ct.ReactorNet([r]) - - with pytest.raises(ct.CanteraError, match="See https://github.com"): - net.solve_steady() + net.solve_steady() # Regression values based on net.advance_to_steady_state(): - # assert r.phase.T == approx(983.7363377) - # assert r.phase.coverages[0] == approx(0.387425501) - # assert sum(r.phase.coverages) == approx(1.0) + assert r.phase.T == approx(983.7363377) + assert rsurf.phase.coverages[0] == approx(0.387425501) + assert sum(rsurf.phase.coverages) == approx(1.0) def test_jacobian(self): gas = ct.Solution("h2o2.yaml", transport_model=None) @@ -3520,7 +3613,7 @@ def test_jacobian(self): test = - mdot / mass * (1 - Y[k]) else: test = mdot / mass * Y[i] * W[k] / W[i] - assert J[i+2,k+2] == approx(test, rel=1e-4), (names[i], names[k]) + assert J[i+2,k+2] == approx(test, rel=5e-4), (names[i], names[k]) def test_logging(self, capsys): messages = [ diff --git a/test_problems/.gitignore b/test_problems/.gitignore index afbc966d1f..76771e4f53 100644 --- a/test_problems/.gitignore +++ b/test_problems/.gitignore @@ -31,5 +31,4 @@ diamondSurf/runDiamond_output.out dustyGasTransport/dustyGasTransport pureFluidTest/pureFluid surfSolverTest/surfaceSolver -surfSolverTest/surfaceSolver2 stoichSolidKinetics/stoichSolidKinetics diff --git a/test_problems/SConscript b/test_problems/SConscript index 563de6b305..633b33da4a 100644 --- a/test_problems/SConscript +++ b/test_problems/SConscript @@ -194,10 +194,6 @@ CompileAndTest('surfSolver', 'surfSolverTest', 'surfaceSolver', None, comparisons=[('results_blessed.txt', 'results.txt')], artifacts=['results.txt'], extensions=['^surfaceSolver.cpp']) -CompileAndTest('surfSolver2', 'surfSolverTest', 'surfaceSolver2', None, - comparisons=[('results2_blessed.txt', 'results2.txt')], - artifacts=['results2.txt'], - extensions=['^surfaceSolver2.cpp']) CompileAndTest('VCS-NaCl', 'VCSnonideal/NaCl_equil', 'nacl_equil', 'good_out.txt', options='-d 3', diff --git a/test_problems/diamondSurf/runDiamond.cpp b/test_problems/diamondSurf/runDiamond.cpp index 2715f6f38f..1adc3dbcb1 100644 --- a/test_problems/diamondSurf/runDiamond.cpp +++ b/test_problems/diamondSurf/runDiamond.cpp @@ -15,22 +15,23 @@ int main(int argc, char** argv) int i, k; try { - auto gas = newThermo("diamond.yaml", "gas"); + CanteraError::setStackTraceDepth(20); + auto iface = newInterface("diamond.yaml", "diamond_100"); + auto gas = iface->adjacent("gas")->thermo(); size_t nsp = gas->nSpecies(); cout.precision(4); cout << "Number of species = " << nsp << endl; - auto diamond = newThermo("diamond.yaml", "diamond"); + auto diamond = iface->adjacent("diamond")->thermo(); size_t nsp_diamond = diamond->nSpecies(); cout << "Number of species in diamond = " << nsp_diamond << endl; - auto diamond100 = newThermo("diamond.yaml", "diamond_100"); + auto diamond100 = iface->thermo(); size_t nsp_d100 = diamond100->nSpecies(); cout << "Number of species in diamond_100 = " << nsp_d100 << endl; - auto kin = newKinetics({diamond100, gas, diamond}, "diamond.yaml"); - InterfaceKinetics& ikin = dynamic_cast(*kin); - size_t nr = kin->nReactions(); + auto ikin = iface->kinetics(); + size_t nr = ikin->nReactions(); cout << "Number of reactions = " << nr << endl; double x[20]; @@ -61,7 +62,7 @@ int main(int argc, char** argv) x[0] = 1.0; diamond100->setState_TP(1200., p); - ikin.advanceCoverages(100.); + ikin->advanceCoverages(100.); // Throw some asserts in here to test that they compile AssertTrace(p == p); @@ -72,7 +73,7 @@ int main(int argc, char** argv) for (i = 0; i < 20; i++) { src[i] = 0.0; } - kin->getNetProductionRates(src); + ikin->getNetProductionRates(src); double sum = 0.0; double naH = 0.0; for (k = 0; k < 13; k++) { diff --git a/test_problems/surfSolverTest/results2_blessed.txt b/test_problems/surfSolverTest/results2_blessed.txt deleted file mode 100644 index 40306c47c4..0000000000 --- a/test_problems/surfSolverTest/results2_blessed.txt +++ /dev/null @@ -1,44 +0,0 @@ -Gas Temperature = 1.4e+03 -Gas Pressure = 1.01e+05 -Gas Phase: gas (2) - Name Conc MoleF SrcRate - (kmol/m^3) (kmol/m^2/s) - 0 H 8.69602e-05 0.00999001 -2.365619e-03 - 1 N2 0.00773858 0.889011 0.000000e+00 - 2 CH3 0 0 0.000000e+00 - 3 CH4 8.69601e-05 0.00999001 0.000000e+00 - 4 C2H2 8.69602e-05 0.00999001 -1.007791e-04 - 5 H2 0.000347841 0.03996 1.282498e-03 - 6 OH 8.69602e-07 9.99001e-05 -7.542502e-05 - 7 H2O 0.000347841 0.03996 3.880308e-05 - 8 CO 0 0 3.843405e-05 - 9 O2 8.69602e-06 0.000999001 -9.060640e-07 -Sum of gas mole fractions= 1 - -Bulk Phase: soot (11) -Bulk Temperature = 1.4e+03 -Bulk Pressure = 1.01e+05 - Name Conc MoleF SrcRate - (kmol/m^3) (kmol/m^2/s) - 0 CB-CB3 293.065 1 1.631241e-04 -Bulk Weight Growth Rate = 0.00196 kg/m^2/s -Bulk Growth Rate = 5.57e-07 m/s -Bulk Growth Rate = 2e+03 microns / hour -Density of bulk phase = 3.52e+03 kg / m^3 - = 3.52 gm / cm^3 -Sum of bulk mole fractions= 1 - -Surface Phase: soot_interface (0) -Surface Temperature = 1.4e+03 -Surface Pressure = 1.01e+05 - Name Coverage SrcRate - 0 Csoot-* 0.0184677 0.000000e+00 - 1 Csoot-H 0.981532 0.000000e+00 -Sum of coverages = 1 -Surface Phase: soot_interface (0) -Surface Temperature = 1.4e+03 -Surface Pressure = 1.01e+05 - Name Coverage SrcRate - 0 Csoot-* 0.0184677 0.000000e+00 - 1 Csoot-H 0.981532 0.000000e+00 -Sum of coverages = 1 diff --git a/test_problems/surfSolverTest/surfaceSolver.cpp b/test_problems/surfSolverTest/surfaceSolver.cpp index 1fa8aa86e3..72278fa13b 100644 --- a/test_problems/surfSolverTest/surfaceSolver.cpp +++ b/test_problems/surfSolverTest/surfaceSolver.cpp @@ -157,11 +157,10 @@ int main(int argc, char** argv) ofstream ofile("results.txt"); - iKin_ptr->setIOFlag(ioflag); /* * Solve the Equation system */ - iKin_ptr->solvePseudoSteadyStateProblem(); + iKin_ptr->solvePseudoSteadyStateProblem(ioflag); /* * Download the source terms for the rate equations @@ -192,7 +191,7 @@ int main(int argc, char** argv) gasTP->setMoleFractions(x); gasTP->setPressure(pres); - iKin_ptr->solvePseudoSteadyStateProblem(); + iKin_ptr->solvePseudoSteadyStateProblem(ioflag); iKin_ptr->getNetProductionRates(src); printGas(cout, gasTP, iKin_ptr, src); @@ -210,8 +209,10 @@ int main(int argc, char** argv) double temp = surfPhaseTP->temperature(); temp += 95; surfPhaseTP->setState_TP(temp, pres); + gasTP->setState_TP(temp, pres); + bulkPhaseTP->setState_TP(temp, pres); - iKin_ptr->solvePseudoSteadyStateProblem(); + iKin_ptr->solvePseudoSteadyStateProblem(ioflag); iKin_ptr->getNetProductionRates(src); printGas(cout, gasTP, iKin_ptr, src); @@ -227,7 +228,7 @@ int main(int argc, char** argv) /****************************************************************************/ surfPhaseTP->setState_TP(temp, pres); - iKin_ptr->solvePseudoSteadyStateProblem(); + iKin_ptr->solvePseudoSteadyStateProblem(ioflag); iKin_ptr->getNetProductionRates(src); printGas(cout, gasTP, iKin_ptr, src); diff --git a/test_problems/surfSolverTest/surfaceSolver2.cpp b/test_problems/surfSolverTest/surfaceSolver2.cpp deleted file mode 100644 index 41a5c03c27..0000000000 --- a/test_problems/surfSolverTest/surfaceSolver2.cpp +++ /dev/null @@ -1,280 +0,0 @@ -/** - * @file surfaceSolver2.cpp - * - */ - -// This file is part of Cantera. See License.txt in the top-level directory or -// at https://cantera.org/license.txt for license and copyright information. - -#define MSSIZE 200 - -#include "cantera/thermo/ThermoFactory.h" -#include "cantera/kinetics/KineticsFactory.h" -#include "cantera/kinetics/ImplicitSurfChem.h" - -#include -#include - -using namespace std; -using namespace Cantera; - -void printGas(ostream& oooo, shared_ptr gasTP, - shared_ptr iKin_ptr, double* src) -{ - double x[MSSIZE]; - double C[MSSIZE]; - oooo.precision(3); - string gasPhaseName = "gas"; - gasTP->getMoleFractions(x); - gasTP->getConcentrations(C); - double Temp = gasTP->temperature(); - double p = gasTP->pressure(); - oooo << "Gas Temperature = " << Temp << endl; - oooo << "Gas Pressure = " << p << endl; - size_t nGas = iKin_ptr->phaseIndex(gasPhaseName); - size_t kstart = iKin_ptr->kineticsSpeciesIndex(0, nGas); - oooo << "Gas Phase: " << gasPhaseName << " " - << "(" << kstart << ")" << endl; - oooo << " Name " - << " Conc MoleF SrcRate " << endl; - oooo << " " - << " (kmol/m^3) (kmol/m^2/s) " << endl; - double sum = 0.0; - size_t nspGas = gasTP->nSpecies(); - for (size_t k = 0; k < nspGas; k++) { - kstart = iKin_ptr->kineticsSpeciesIndex(k, nGas); - fmt::print(oooo, "{:4d} {:>24s} {:14g} {:14g} {:14e}\n", - k, gasTP->speciesName(k), C[k], x[k], src[kstart]); - sum += x[k]; - } - oooo << "Sum of gas mole fractions= " << sum << endl; - oooo << endl; - - -} - -void printBulk(ostream& oooo, shared_ptr bulkPhaseTP, - shared_ptr iKin_ptr, double* src) -{ - double x[MSSIZE]; - double C[MSSIZE]; - oooo.precision(3); - string bulkParticlePhaseName = bulkPhaseTP->name(); - bulkPhaseTP->getMoleFractions(x); - bulkPhaseTP->getConcentrations(C); - size_t nBulk = iKin_ptr->phaseIndex(bulkParticlePhaseName); - size_t kstart = iKin_ptr->kineticsSpeciesIndex(0, nBulk); - double dens = bulkPhaseTP->density(); - oooo << "Bulk Phase: " << bulkParticlePhaseName << " " - << "(" << kstart << ")" << endl; - double Temp = bulkPhaseTP->temperature(); - double p = bulkPhaseTP->pressure(); - oooo << "Bulk Temperature = " << Temp << endl; - oooo << "Bulk Pressure = " << p << endl; - oooo << " Name " - << " Conc MoleF SrcRate " << endl; - oooo << " " - << " (kmol/m^3) (kmol/m^2/s) " << endl; - double sum = 0.0; - double Wsum = 0.0; - const vector& molecW = bulkPhaseTP->molecularWeights(); - size_t nspBulk = bulkPhaseTP->nSpecies(); - for (size_t k = 0; k < nspBulk; k++) { - kstart = iKin_ptr->kineticsSpeciesIndex(k, nBulk); - fmt::print(oooo, "{:4d} {:>24s} {:14g} {:14g} {:14e}\n", - k, bulkPhaseTP->speciesName(k), C[k], x[k], src[kstart]); - sum += x[k]; - Wsum += src[kstart] * molecW[k]; - } - oooo << "Bulk Weight Growth Rate = " << Wsum << " kg/m^2/s" << endl; - double gr = Wsum / dens; - oooo << "Bulk Growth Rate = " << gr << " m/s" << endl; - oooo << "Bulk Growth Rate = " << gr * 1.0E6 * 3600. - << " microns / hour" << endl; - oooo << "Density of bulk phase = " << dens << " kg / m^3 "<< endl; - oooo << " = " << dens / 1.0E3 - <<" gm / cm^3 " << endl; - oooo << "Sum of bulk mole fractions= " << sum << endl; - oooo << endl; -} - -void printSurf(ostream& oooo, shared_ptr surfPhaseTP, - shared_ptr iKin_ptr, double* src) -{ - double x[MSSIZE]; - oooo.precision(3); - string surfParticlePhaseName = surfPhaseTP->name(); - surfPhaseTP->getMoleFractions(x); - size_t nSurf = iKin_ptr->phaseIndex(surfParticlePhaseName); - size_t kstart = iKin_ptr->kineticsSpeciesIndex(0, nSurf); - oooo << "Surface Phase: " << surfParticlePhaseName - << " (" << kstart << ")" << endl; - double Temp = surfPhaseTP->temperature(); - double p = surfPhaseTP->pressure(); - oooo << "Surface Temperature = " << Temp << endl; - oooo << "Surface Pressure = " << p << endl; - oooo << " Name " - << " Coverage SrcRate " << endl; - double sum = 0.0; - size_t nspSurf = surfPhaseTP->nSpecies(); - for (size_t k = 0; k < nspSurf; k++) { - kstart = iKin_ptr->kineticsSpeciesIndex(k, nSurf); - double srcK = src[kstart]; - if (fabs(srcK) < 1.0E-8) { - srcK = 0.0; - } - fmt::print(oooo, "{:4d} {:>24s} {:14g} {:14e}\n", - k, surfPhaseTP->speciesName(k), x[k], srcK); - sum += x[k]; - } - oooo << "Sum of coverages = " << sum << endl; -} - -int main(int argc, char** argv) -{ - string infile = "haca2.yaml"; - string gasPhaseName = "gas"; - string bulkParticlePhaseName = "soot"; - string surfParticlePhaseName = "soot_interface"; - int ioflag = 1; - - try { - auto gasTP = newThermo(infile, gasPhaseName); - size_t nspGas = gasTP->nSpecies(); - cout << "Number of species = " << nspGas << endl; - - auto bulkPhaseTP = newThermo(infile, bulkParticlePhaseName); - size_t nspBulk = bulkPhaseTP->nSpecies(); - cout << "Number of species in bulk phase named " << - bulkParticlePhaseName << " = " << nspBulk << endl; - - auto surfPhaseTP = newThermo(infile, surfParticlePhaseName); - size_t nsp_d100 = surfPhaseTP->nSpecies(); - cout << "Number of species in surface phase, " << surfParticlePhaseName - << " = " << nsp_d100 << endl; - - auto kin = newKinetics({surfPhaseTP, gasTP, bulkPhaseTP}, infile); - auto iKin_ptr = dynamic_pointer_cast(kin); - size_t nr = iKin_ptr->nReactions(); - cout << "Number of reactions = " << nr << endl; - - ofstream ofile("results2.txt"); - - // create a second copy of the same surface phase - // (this is a made up problem btw to check the software capability) - auto surfPhaseTP2 = newThermo(infile, surfParticlePhaseName); - size_t nsp2 = surfPhaseTP2->nSpecies(); - string pname = surfPhaseTP2->name(); - cout << "Number of species in 2nd surface phase, " << pname - << " = " << nsp2 << endl; - - // create the second InterfaceKinetics object based on the - // second surface phase. - auto kin2 = newKinetics({surfPhaseTP2, gasTP, bulkPhaseTP}, infile); - auto iKin2_ptr = dynamic_pointer_cast(kin2); - nr = iKin_ptr->nReactions(); - cout << "Number of reactions = " << nr << endl; - - double x[MSSIZE]; - - /* - * Set-up the Surface Problem - * This problem will consist of 2 identical InterfaceKinetics objects - */ - vector vecKinPtrs { iKin_ptr.get(), iKin2_ptr.get() }; - - // Create the ImplicitSurfChem problem - // Initialize it and call the pseudo steadystate capability. - ImplicitSurfChem* surfaceProb = new ImplicitSurfChem(vecKinPtrs); - surfaceProb->initialize(); - surfaceProb->setIOFlag(ioflag); - surfaceProb->solvePseudoSteadyStateProblem(); - - /* - * Download the source terms for the rate equations - */ - double src[MSSIZE]; - double src2[MSSIZE]; - iKin_ptr->getNetProductionRates(src); - iKin2_ptr->getNetProductionRates(src2); - - printGas(cout, gasTP, iKin_ptr, src); - printBulk(cout, bulkPhaseTP, iKin_ptr, src); - printSurf(cout, surfPhaseTP, iKin_ptr, src) ; - printSurf(cout, surfPhaseTP2, iKin2_ptr, src2) ; - - - printGas(ofile, gasTP, iKin_ptr, src); - printBulk(ofile, bulkPhaseTP, iKin_ptr, src); - printSurf(ofile, surfPhaseTP, iKin_ptr, src) ; - printSurf(ofile, surfPhaseTP2, iKin2_ptr, src2) ; - - /*****************************************************************************/ - /* Now Tweak the inputs and do a quick calculation */ - /****************************************************************************/ - - /* - * Set the Gas State: - * -> note that the states are set in the input file too - */ - double pres = gasTP->pressure(); - gasTP->getMoleFractions(x); - double tmp = 0.3 * std::min(x[0], x[1]); - x[0] += tmp; - x[1] -= tmp; - gasTP->setMoleFractions(x); - gasTP->setPressure(pres); - - surfaceProb->solvePseudoSteadyStateProblem(); - iKin_ptr->getNetProductionRates(src); - iKin2_ptr->getNetProductionRates(src2); - - printGas(cout, gasTP, iKin_ptr, src); - printBulk(cout, bulkPhaseTP, iKin_ptr, src); - printSurf(cout, surfPhaseTP, iKin_ptr, src) ; - printSurf(cout, surfPhaseTP2, iKin2_ptr, src2) ; - - /*****************************************************************************/ - /* Now Tweak the inputs and do a quick calculation */ - /****************************************************************************/ - - pres = surfPhaseTP->pressure(); - double temp = surfPhaseTP->temperature(); - temp += 95; - surfPhaseTP->setState_TP(temp, pres); - - surfaceProb->solvePseudoSteadyStateProblem(); - iKin_ptr->getNetProductionRates(src); - iKin2_ptr->getNetProductionRates(src2); - - printGas(cout, gasTP, iKin_ptr, src); - printBulk(cout, bulkPhaseTP, iKin_ptr, src); - printSurf(cout, surfPhaseTP, iKin_ptr, src) ; - printSurf(cout, surfPhaseTP2, iKin2_ptr, src2) ; - - /*****************************************************************************/ - /* Now Don't Tweak the inputs at all */ - /****************************************************************************/ - gasTP->setState_TP(temp, pres); - - surfaceProb->solvePseudoSteadyStateProblem(); - iKin_ptr->getNetProductionRates(src); - iKin2_ptr->getNetProductionRates(src2); - - printGas(cout, gasTP, iKin_ptr, src); - printBulk(cout, bulkPhaseTP, iKin_ptr, src); - printSurf(cout, surfPhaseTP, iKin_ptr, src) ; - printSurf(cout, surfPhaseTP2, iKin2_ptr, src2) ; - - delete surfaceProb; - surfaceProb = 0; - iKin_ptr = 0; - iKin2_ptr = 0; - appdelete(); - } catch (CanteraError& err) { - std::cout << err.what() << std::endl; - } - - return 0; -} -/***********************************************************/