diff --git a/doc/sphinx/reference/onedim/governing-equations.md b/doc/sphinx/reference/onedim/governing-equations.md index 15d6adf975c..b6a04be8345 100644 --- a/doc/sphinx/reference/onedim/governing-equations.md +++ b/doc/sphinx/reference/onedim/governing-equations.md @@ -7,7 +7,7 @@ solution to reduce the three-dimensional governing equations to a single dimensi ## Axisymmetric Stagnation Flow The governing equations for a steady axisymmetric stagnation flow follow those derived -in Section 7.2 of {cite:t}`kee2017` and are implemented by class {ct}`StFlow`. +in Section 7.2 of {cite:t}`kee2017` and are implemented by class {ct}`Flow1D`. *Continuity*: @@ -99,7 +99,7 @@ $$ where $D_{ki}$ is the multicomponent diffusion coefficient and $D_k^T$ is the Soret diffusion coefficient. Inclusion of the Soret calculation must be explicitly enabled when setting up the simulation, on top of specifying a multicomponent transport model, -for example by using the {ct}`StFlow::enableSoret` method (C++) or setting the +for example by using the {ct}`Flow1D::enableSoret` method (C++) or setting the {py:attr}`~cantera.FlameBase.soret_enabled` property (Python). ## Boundary Conditions @@ -295,4 +295,4 @@ $$ $$ \dot{m}_\t{in,right} = - \rho(z_{\t{in,right}}) U_o(z_{\t{in,right}}) -$$ \ No newline at end of file +$$ diff --git a/include/cantera/clib/ctonedim.h b/include/cantera/clib/ctonedim.h index 438aa7df64d..b5f3023a5e1 100644 --- a/include/cantera/clib/ctonedim.h +++ b/include/cantera/clib/ctonedim.h @@ -59,6 +59,16 @@ extern "C" { CANTERA_CAPI int inlet_setSpreadRate(int i, double v); + CANTERA_CAPI int flow1D_new(int iph, int ikin, int itr, int itype); + CANTERA_CAPI int flow1D_setTransport(int i, int itr); + CANTERA_CAPI int flow1D_enableSoret(int i, int iSoret); + CANTERA_CAPI int flow1D_setPressure(int i, double p); + CANTERA_CAPI double flow1D_pressure(int i); + CANTERA_CAPI int flow1D_setFixedTempProfile(int i, size_t n, const double* pos, + size_t m, const double* temp); + CANTERA_CAPI int flow1D_solveEnergyEqn(int i, int flag); + + //! @todo: Remove all functions with `stflow` prefix after Cantera 3.1 CANTERA_CAPI int stflow_new(int iph, int ikin, int itr, int itype); CANTERA_CAPI int stflow_setTransport(int i, int itr); CANTERA_CAPI int stflow_enableSoret(int i, int iSoret); diff --git a/include/cantera/oneD/Boundary1D.h b/include/cantera/oneD/Boundary1D.h index da39a495baf..3c285bb3741 100644 --- a/include/cantera/oneD/Boundary1D.h +++ b/include/cantera/oneD/Boundary1D.h @@ -13,7 +13,7 @@ #include "Domain1D.h" #include "cantera/thermo/SurfPhase.h" #include "cantera/kinetics/InterfaceKinetics.h" -#include "StFlow.h" +#include "Flow1D.h" namespace Cantera { @@ -109,8 +109,8 @@ class Boundary1D : public Domain1D protected: void _init(size_t n); - StFlow* m_flow_left = nullptr; - StFlow* m_flow_right = nullptr; + Flow1D* m_flow_left = nullptr; + Flow1D* m_flow_right = nullptr; size_t m_ilr = 0; size_t m_left_nv = 0; size_t m_right_nv = 0; @@ -176,7 +176,7 @@ class Inlet1D : public Boundary1D size_t m_nsp = 0; vector m_yin; string m_xstr; - StFlow* m_flow = nullptr; + Flow1D* m_flow = nullptr; }; /** @@ -293,7 +293,7 @@ class OutletRes1D : public Boundary1D size_t m_nsp = 0; vector m_yres; string m_xstr; - StFlow* m_flow = nullptr; + Flow1D* m_flow = nullptr; }; /** diff --git a/include/cantera/oneD/Flow1D.h b/include/cantera/oneD/Flow1D.h new file mode 100644 index 00000000000..aa48218b8e8 --- /dev/null +++ b/include/cantera/oneD/Flow1D.h @@ -0,0 +1,816 @@ +//! @file Flow1D.h + +// 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_FLOW1D_H +#define CT_FLOW1D_H + +#include "Domain1D.h" +#include "cantera/base/Array.h" +#include "cantera/base/Solution.h" +#include "cantera/thermo/ThermoPhase.h" +#include "cantera/kinetics/Kinetics.h" + +namespace Cantera +{ + +//------------------------------------------ +// constants +//------------------------------------------ + +//! Offsets of solution components in the 1D solution array. +enum offset +{ + c_offset_U //! axial velocity [m/s] + , c_offset_V //! strain rate + , c_offset_T //! temperature [kelvin] + , c_offset_L //! (1/r)dP/dr + , c_offset_E //! electric field + , c_offset_Uo //! oxidizer axial velocity [m/s] + , c_offset_Y //! mass fractions +}; + +class Transport; + +//! @defgroup flowGroup Flow Domains +//! One-dimensional flow domains. +//! @ingroup onedGroup + +/** + * This class represents 1D flow domains that satisfy the one-dimensional + * similarity solution for chemically-reacting, axisymmetric flows. + * @ingroup flowGroup + */ +class Flow1D : public Domain1D +{ +public: + //-------------------------------- + // construction and destruction + //-------------------------------- + + //! Create a new flow domain. + //! @param ph Object representing the gas phase. This object will be used + //! to evaluate all thermodynamic, kinetic, and transport properties. + //! @param nsp Number of species. + //! @param points Initial number of grid points + Flow1D(ThermoPhase* ph = 0, size_t nsp = 1, size_t points = 1); + + //! Delegating constructor + Flow1D(shared_ptr th, size_t nsp = 1, size_t points = 1); + + //! Create a new flow domain. + //! @param sol Solution object used to evaluate all thermodynamic, kinetic, and + //! transport properties + //! @param id name of flow domain + //! @param points initial number of grid points + Flow1D(shared_ptr sol, const string& id="", size_t points=1); + + ~Flow1D(); + + string domainType() const override; + + //! @name Problem Specification + //! @{ + + void setupGrid(size_t n, const double* z) override; + + void resetBadValues(double* xg) override; + + ThermoPhase& phase() { + return *m_thermo; + } + + Kinetics& kinetics() { + return *m_kin; + } + + void setKinetics(shared_ptr kin) override; + + void setTransport(shared_ptr trans) override; + + //! Set the transport model + //! @since New in %Cantera 3.0. + void setTransportModel(const string& trans); + + //! Retrieve transport model + //! @since New in %Cantera 3.0. + string transportModel() const; + + //! Enable thermal diffusion, also known as Soret diffusion. + //! Requires that multicomponent transport properties be + //! enabled to carry out calculations. + void enableSoret(bool withSoret) { + m_do_soret = withSoret; + } + bool withSoret() const { + return m_do_soret; + } + + //! Compute species diffusive fluxes with respect to + //! their mass fraction gradients (fluxGradientBasis = ThermoBasis::mass) + //! or mole fraction gradients (fluxGradientBasis = ThermoBasis::molar, default) + //! when using the mixture-averaged transport model. + //! @param fluxGradientBasis set flux computation to mass or mole basis + //! @since New in %Cantera 3.1. + void setFluxGradientBasis(ThermoBasis fluxGradientBasis); + + //! Compute species diffusive fluxes with respect to + //! their mass fraction gradients (fluxGradientBasis = ThermoBasis::mass) + //! or mole fraction gradients (fluxGradientBasis = ThermoBasis::molar, default) + //! when using the mixture-averaged transport model. + //! @return the basis used for flux computation (mass or mole fraction gradients) + //! @since New in %Cantera 3.1. + ThermoBasis fluxGradientBasis() const { + return m_fluxGradientBasis; + } + + //! Set the pressure. Since the flow equations are for the limit of small + //! Mach number, the pressure is very nearly constant throughout the flow. + void setPressure(double p) { + m_press = p; + } + + //! The current pressure [Pa]. + double pressure() const { + return m_press; + } + + //! Write the initial solution estimate into array x. + void _getInitialSoln(double* x) override; + + void _finalize(const double* x) override; + + //! Sometimes it is desired to carry out the simulation using a specified + //! temperature profile, rather than computing it by solving the energy + //! equation. This method specifies this profile. + void setFixedTempProfile(vector& zfixed, vector& tfixed) { + m_zfix = zfixed; + m_tfix = tfixed; + } + + /** + * Set the temperature fixed point at grid point j, and disable the energy + * equation so that the solution will be held to this value. + */ + void setTemperature(size_t j, double t) { + m_fixedtemp[j] = t; + m_do_energy[j] = false; + } + + //! The fixed temperature value at point j. + double T_fixed(size_t j) const { + return m_fixedtemp[j]; + } + + //! @} + + string componentName(size_t n) const override; + + size_t componentIndex(const string& name) const override; + + //! Returns true if the specified component is an active part of the solver state + virtual bool componentActive(size_t n) const; + + void show(const double* x) override; + + shared_ptr asArray(const double* soln) const override; + void fromArray(SolutionArray& arr, double* soln) override; + + //! Set flow configuration for freely-propagating flames, using an internal point + //! with a fixed temperature as the condition to determine the inlet mass flux. + void setFreeFlow() { + m_dovisc = false; + m_isFree = true; + m_usesLambda = false; + } + + //! Set flow configuration for axisymmetric counterflow flames, using specified + //! inlet mass fluxes. + void setAxisymmetricFlow() { + m_dovisc = true; + m_isFree = false; + m_usesLambda = true; + } + + //! Set flow configuration for burner-stabilized flames, using specified inlet mass + //! fluxes. + void setUnstrainedFlow() { + m_dovisc = false; + m_isFree = false; + m_usesLambda = false; + } + + void solveEnergyEqn(size_t j=npos); + + //! Get the solving stage (used by IonFlow specialization) + //! @since New in %Cantera 3.0 + virtual size_t getSolvingStage() const; + + //! Solving stage mode for handling ionized species (used by IonFlow specialization) + //! - @c stage=1: the fluxes of charged species are set to zero + //! - @c stage=2: the electric field equation is solved, and the drift flux for + //! ionized species is evaluated + virtual void setSolvingStage(const size_t stage); + + //! Set to solve electric field in a point (used by IonFlow specialization) + virtual void solveElectricField(size_t j=npos); + + //! Set to fix voltage in a point (used by IonFlow specialization) + virtual void fixElectricField(size_t j=npos); + + //! Retrieve flag indicating whether electric field is solved or not (used by + //! IonFlow specialization) + virtual bool doElectricField(size_t j) const; + + //! Turn radiation on / off. + void enableRadiation(bool doRadiation) { + m_do_radiation = doRadiation; + } + + //! Returns `true` if the radiation term in the energy equation is enabled + bool radiationEnabled() const { + return m_do_radiation; + } + + //! Return radiative heat loss at grid point j + double radiativeHeatLoss(size_t j) const { + return m_qdotRadiation[j]; + } + + //! Set the emissivities for the boundary values + /*! + * Reads the emissivities for the left and right boundary values in the + * radiative term and writes them into the variables, which are used for the + * calculation. + */ + void setBoundaryEmissivities(double e_left, double e_right); + + //! Return emissivity at left boundary + double leftEmissivity() const { + return m_epsilon_left; + } + + //! Return emissivity at right boundary + double rightEmissivity() const { + return m_epsilon_right; + } + + void fixTemperature(size_t j=npos); + + /** + * @name Two-Point control method + * + * In this method two control points are designated in the 1D domain, and the value + * of the temperature at these points is fixed. The values of the control points are + * imposed and thus serve as a boundary condition that affects the solution of the + * governing equations in the 1D domain. The imposition of fixed points in the + * domain means that the original set of governing equations' boundary conditions + * would over-determine the problem. Thus, the boundary conditions are changed to + * reflect the fact that the control points are serving as internal boundary + * conditions. + * + * The imposition of the two internal boundary conditions requires that two other + * boundary conditions be changed. The first is the boundary condition for the + * continuity equation at the left boundary, which is changed to be a value that is + * derived from the solution at the left boundary. The second is the continuity + * boundary condition at the right boundary, which is also determined from the flow + * solution by using the oxidizer axial velocity equation variable to compute the + * mass flux at the right boundary. + * + * This method is based on the work of Nishioka et al. @cite nishioka1996 . + */ + //! @{ + + //! Returns the temperature at the left control point + double leftControlPointTemperature() const; + + //! Returns the z-coordinate of the left control point + double leftControlPointCoordinate() const; + + //! Sets the temperature of the left control point + void setLeftControlPointTemperature(double temperature); + + //! Sets the coordinate of the left control point + void setLeftControlPointCoordinate(double z_left); + + //! Returns the temperature at the right control point + double rightControlPointTemperature() const; + + //! Returns the z-coordinate of the right control point + double rightControlPointCoordinate() const; + + //! Sets the temperature of the right control point + void setRightControlPointTemperature(double temperature); + + //! Sets the coordinate of the right control point + void setRightControlPointCoordinate(double z_right); + + //! Sets the status of the two-point control + void enableTwoPointControl(bool twoPointControl); + + //! Returns the status of the two-point control + bool twoPointControlEnabled() const { + return m_twoPointControl; + } + //! @} + + bool doEnergy(size_t j) { + return m_do_energy[j]; + } + + //! Change the grid size. Called after grid refinement. + void resize(size_t components, size_t points) override; + + //! Set the gas object state to be consistent with the solution at point j. + void setGas(const double* x, size_t j); + + //! Set the gas state to be consistent with the solution at the midpoint + //! between j and j + 1. + void setGasAtMidpoint(const double* x, size_t j); + + double density(size_t j) const { + return m_rho[j]; + } + + /** + * Retrieve flag indicating whether flow is freely propagating. + * The flow is unstrained and the axial mass flow rate is not specified. + * For free flame propagation, the axial velocity is determined by the solver. + * @since New in %Cantera 3.0 + */ + bool isFree() const { + return m_isFree; + } + + /** + * Retrieve flag indicating whether flow uses radial momentum. + * If `true`, radial momentum equation for @f$ V @f$ as well as + * @f$ d\Lambda/dz = 0 @f$ are solved; if `false`, @f$ \Lambda(z) = 0 @f$ and + * @f$ V(z) = 0 @f$ by definition. + * @since New in %Cantera 3.0 + */ + bool isStrained() const { + return m_usesLambda; + } + + void setViscosityFlag(bool dovisc) { + m_dovisc = dovisc; + } + + /** + * Evaluate the residual functions for axisymmetric stagnation flow. + * If jGlobal == npos, the residual function is evaluated at all grid points. + * Otherwise, the residual function is only evaluated at grid points j-1, j, + * and j+1. This option is used to efficiently evaluate the Jacobian numerically. + * + * These residuals at all the boundary grid points are evaluated using a default + * boundary condition that may be modified by a boundary object that is attached + * to the domain. The boundary object connected will modify these equations by + * subtracting the boundary object's values for V, T, mdot, etc. As a result, + * these residual equations will force the solution variables to the values of + * the connected boundary object. + * + * @param jGlobal Global grid point at which to update the residual + * @param[in] xGlobal Global state vector + * @param[out] rsdGlobal Global residual vector + * @param[out] diagGlobal Global boolean mask indicating whether each solution + * component has a time derivative (1) or not (0). + * @param[in] rdt Reciprocal of the timestep (`rdt=0` implies steady-state.) + */ + void eval(size_t jGlobal, double* xGlobal, double* rsdGlobal, + integer* diagGlobal, double rdt) override; + + //! Index of the species on the left boundary with the largest mass fraction + size_t leftExcessSpecies() const { + return m_kExcessLeft; + } + + //! Index of the species on the right boundary with the largest mass fraction + size_t rightExcessSpecies() const { + return m_kExcessRight; + } + +protected: + AnyMap getMeta() const override; + void setMeta(const AnyMap& state) override; + + //! @name Updates of cached properties + //! These methods are called by eval() to update cached properties and data that are + //! used for the evaluation of the governing equations. + //! @{ + + /** + * Update the thermodynamic properties from point j0 to point j1 + * (inclusive), based on solution x. + * + * The gas state is set to be consistent with the solution at the + * points from j0 to j1. + * + * Properties that are computed and cached are: + * * #m_rho (density) + * * #m_wtm (mean molecular weight) + * * #m_cp (specific heat capacity) + * * #m_hk (species specific enthalpies) + * * #m_wdot (species production rates) + */ + void updateThermo(const double* x, size_t j0, size_t j1) { + for (size_t j = j0; j <= j1; j++) { + setGas(x,j); + m_rho[j] = m_thermo->density(); + m_wtm[j] = m_thermo->meanMolecularWeight(); + m_cp[j] = m_thermo->cp_mass(); + m_thermo->getPartialMolarEnthalpies(&m_hk(0, j)); + m_kin->getNetProductionRates(&m_wdot(0, j)); + } + } + + //! Update the transport properties at grid points in the range from `j0` + //! to `j1`, based on solution `x`. + virtual void updateTransport(double* x, size_t j0, size_t j1); + + //! Update the diffusive mass fluxes. + virtual void updateDiffFluxes(const double* x, size_t j0, size_t j1); + + //! Update the properties (thermo, transport, and diffusion flux). + //! This function is called in eval after the points which need + //! to be updated are defined. + virtual void updateProperties(size_t jg, double* x, size_t jmin, size_t jmax); + + /** + * Computes the radiative heat loss vector over points jmin to jmax and stores + * the data in the qdotRadiation variable. + * + * The simple radiation model used was established by Liu and Rogg + * @cite liu1991. This model considers the radiation of CO2 and H2O. + * + * This model uses the optically thin limit and the gray-gas approximation to + * simply calculate a volume specified heat flux out of the Planck absorption + * coefficients, the boundary emissivities and the temperature. Polynomial lines + * calculate the species Planck coefficients for H2O and CO2. The data for the + * lines are taken from the RADCAL program @cite RADCAL. + * The coefficients for the polynomials are taken from + * [TNF Workshop](https://tnfworkshop.org/radiation/) material. + */ + void computeRadiation(double* x, size_t jmin, size_t jmax); + + //! @} + + //! @name Governing Equations + //! Methods called by eval() to calculate residuals for individual governing + //! equations. + + /** + * Evaluate the continuity equation residual. + * + * This function calculates the residual of the continuity equation + * @f[ + * \frac{d(\rho u)}{dz} + 2\rho V = 0 + * @f] + * + * Axisymmetric flame: + * The continuity equation propagates information from right-to-left. + * The @f$ \rho u @f$ at point 0 is dependent on @f$ \rho u @f$ at point 1, + * but not on @f$ \dot{m} @f$ from the inlet. + * + * Freely-propagating flame: + * The continuity equation propagates information away from a fixed temperature + * point that is set in the domain. + * + * Unstrained flame: + * A specified mass flux; the main example being burner-stabilized flames. + * + * The default boundary condition for the continuity equation is + * (@f$ u = 0 @f$) at the left and right boundary. + * + * @param[in] x Local domain state vector, includes variables like temperature, + * density, etc. + * @param[out] rsd Local domain residual vector that stores the continuity + * equation residuals. + * @param[out] diag Local domain diagonal matrix that controls whether an entry + * has a time-derivative (used by the solver). + * @param[in] rdt Reciprocal of the timestep. + * @param[in] jmin The index for the starting point in the local domain grid. + * @param[in] jmax The index for the ending point in the local domain grid. + */ + virtual void evalContinuity(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax); + + /** + * Evaluate the momentum equation residual. + * + * The function calculates the radial momentum equation defined as + * @f[ + * \rho u \frac{dV}{dz} + \rho V^2 = + * \frac{d}{dz}\left( \mu \frac{dV}{dz} \right) - \Lambda + * @f] + * + * The radial momentum equation is used for axisymmetric flows, and incorporates + * terms for time and spatial variations of radial velocity (@f$ V @f$). The + * default boundary condition is zero radial velocity (@f$ V @f$) at the left + * and right boundary. + * + * For argument explanation, see evalContinuity(). + */ + virtual void evalMomentum(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax); + + /** + * Evaluate the lambda equation residual. + * + * The function calculates the lambda equation as + * @f[ + * \frac{d\Lambda}{dz} = 0 + * @f] + * + * The lambda equation serves as an eigenvalue that allows the momentum + * equation and continuity equations to be simultaneously satisfied in + * axisymmetric flows. The lambda equation propagates information from + * left-to-right. The default boundary condition is @f$ \Lambda = 0 @f$ + * at the left and zero flux at the right boundary. + * + * For argument explanation, see evalContinuity(). + */ + virtual void evalLambda(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax); + + /** + * Evaluate the energy equation residual. + * + * The function calculates the energy equation: + * @f[ + * \rho c_p u \frac{dT}{dz} = + * \frac{d}{dz}\left( \lambda \frac{dT}{dz} \right) + * - \sum_k h_kW_k\dot{\omega}_k + * - \sum_k j_k \frac{dh_k}{dz} + * @f] + * + * The energy equation includes contributions from + * chemical reactions and diffusion. Default is zero temperature (@f$ T @f$) + * at the left and right boundaries. These boundary values are updated by the + * specific boundary object connected to the domain. + * + * For argument explanation, see evalContinuity(). + */ + virtual void evalEnergy(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax); + + /** + * Evaluate the species equations' residuals. + * + * The function calculates the species equations as + * @f[ + * \rho u \frac{dY_k}{dz} + \frac{dj_k}{dz} = W_k\dot{\omega}_k + * @f] + * + * The species equations include terms for temporal and spatial variations + * of species mass fractions (@f$ Y_k @f$). The default boundary condition is zero + * flux for species at the left and right boundary. + * + * For argument explanation, see evalContinuity(). + */ + virtual void evalSpecies(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax); + + /** + * Evaluate the electric field equation residual to be zero everywhere. + * + * The electric field equation is implemented in the IonFlow class. The default + * boundary condition is zero electric field (@f$ E @f$) at the boundary, + * and @f$ E @f$ is zero within the domain. + * + * For argument explanation, see evalContinuity(). + */ + virtual void evalElectricField(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax); + + //! @} End of Governing Equations + + //! Alternate version of evalContinuity with legacy signature. + //! Implemented by StFlow; included here to prevent compiler warnings about shadowed + //! virtual functions. + //! @deprecated To be removed after Cantera 3.1. + virtual void evalContinuity(size_t j, double* x, double* r, int* diag, double rdt); + + /** + * Evaluate the oxidizer axial velocity equation residual. + * + * The function calculates the oxidizer axial velocity equation as + * @f[ + * \frac{dU_{o}}{dz} = 0 + * @f] + * + * This equation serves as a dummy equation that is used only in the context of + * two-point flame control, and serves as the way for two interior control points to + * be specified while maintaining block tridiagonal structure. The default boundary + * condition is @f$ U_o = 0 @f$ at the right and zero flux at the left boundary. + * + * For argument explanation, see evalContinuity(). + */ + virtual void evalUo(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax); + + //! @name Solution components + //! @{ + + double T(const double* x, size_t j) const { + return x[index(c_offset_T, j)]; + } + double& T(double* x, size_t j) { + return x[index(c_offset_T, j)]; + } + double T_prev(size_t j) const { + return prevSoln(c_offset_T, j); + } + + double rho_u(const double* x, size_t j) const { + return m_rho[j]*x[index(c_offset_U, j)]; + } + + double u(const double* x, size_t j) const { + return x[index(c_offset_U, j)]; + } + + double V(const double* x, size_t j) const { + return x[index(c_offset_V, j)]; + } + double V_prev(size_t j) const { + return prevSoln(c_offset_V, j); + } + + double lambda(const double* x, size_t j) const { + return x[index(c_offset_L, j)]; + } + + //! Solution component for oxidizer velocity, @see evalUo + double Uo(const double* x, size_t j) const { + return x[index(c_offset_Uo, j)]; + } + + double Y(const double* x, size_t k, size_t j) const { + return x[index(c_offset_Y + k, j)]; + } + + double& Y(double* x, size_t k, size_t j) { + return x[index(c_offset_Y + k, j)]; + } + + double Y_prev(size_t k, size_t j) const { + return prevSoln(c_offset_Y + k, j); + } + + double X(const double* x, size_t k, size_t j) const { + return m_wtm[j]*Y(x,k,j)/m_wt[k]; + } + + double flux(size_t k, size_t j) const { + return m_flux(k, j); + } + //! @} + + //! @name Convective spatial derivatives + //! + //! These use upwind differencing, assuming u(z) is negative + //! @{ + double dVdz(const double* x, size_t j) const { + size_t jloc = (u(x,j) > 0.0 ? j : j + 1); + return (V(x,jloc) - V(x,jloc-1))/m_dz[jloc-1]; + } + + double dYdz(const double* x, size_t k, size_t j) const { + size_t jloc = (u(x,j) > 0.0 ? j : j + 1); + return (Y(x,k,jloc) - Y(x,k,jloc-1))/m_dz[jloc-1]; + } + + double dTdz(const double* x, size_t j) const { + size_t jloc = (u(x,j) > 0.0 ? j : j + 1); + return (T(x,jloc) - T(x,jloc-1))/m_dz[jloc-1]; + } + //! @} + + double shear(const double* x, size_t j) const { + double c1 = m_visc[j-1]*(V(x,j) - V(x,j-1)); + double c2 = m_visc[j]*(V(x,j+1) - V(x,j)); + return 2.0*(c2/(z(j+1) - z(j)) - c1/(z(j) - z(j-1)))/(z(j+1) - z(j-1)); + } + + double divHeatFlux(const double* x, size_t j) const { + double c1 = m_tcon[j-1]*(T(x,j) - T(x,j-1)); + double c2 = m_tcon[j]*(T(x,j+1) - T(x,j)); + return -2.0*(c2/(z(j+1) - z(j)) - c1/(z(j) - z(j-1)))/(z(j+1) - z(j-1)); + } + + size_t mindex(size_t k, size_t j, size_t m) { + return m*m_nsp*m_nsp + m_nsp*j + k; + } + + //! Get the gradient of species specific molar enthalpies + virtual void grad_hk(const double* x, size_t j); + + //--------------------------------------------------------- + // member data + //--------------------------------------------------------- + + double m_press = -1.0; // pressure + + // grid parameters + vector m_dz; + + // mixture thermo properties + vector m_rho; //!< Vector of size #m_nsp to cache densities + vector m_wtm; //!< Vector of size #m_nsp to cache mean molecular weights + + // species thermo properties + vector m_wt; + vector m_cp; //!< Vector of size #m_nsp to cache specific heat capacities + + // transport properties + vector m_visc; + vector m_tcon; + //! Array of size #m_nsp by #m_points for saving density times diffusion + //! coefficient times species molar mass divided by mean molecular weight + vector m_diff; + vector m_multidiff; + Array2D m_dthermal; + Array2D m_flux; + + //! Array of size #m_nsp by #m_points for saving molar enthalpies + Array2D m_hk; + + //! Array of size #m_nsp by #m_points-1 for saving enthalpy fluxes + Array2D m_dhk_dz; + + //! Array of size #m_nsp by #m_points for saving species production rates + Array2D m_wdot; + + size_t m_nsp; //!< Number of species in the mechanism + + ThermoPhase* m_thermo = nullptr; + Kinetics* m_kin = nullptr; + Transport* m_trans = nullptr; + + // boundary emissivities for the radiation calculations + double m_epsilon_left = 0.0; + double m_epsilon_right = 0.0; + + //! Indices within the ThermoPhase of the radiating species. First index is + //! for CO2, second is for H2O. + vector m_kRadiating; + + // flags + vector m_do_energy; + bool m_do_soret = false; + ThermoBasis m_fluxGradientBasis = ThermoBasis::molar; + vector m_do_species; + bool m_do_multicomponent = false; + + //! flag for the radiative heat loss + bool m_do_radiation = false; + + //! radiative heat loss vector + vector m_qdotRadiation; + + // fixed T and Y values + vector m_fixedtemp; + vector m_zfix; + vector m_tfix; + + //! Index of species with a large mass fraction at each boundary, for which + //! the mass fraction may be calculated as 1 minus the sum of the other mass + //! fractions + size_t m_kExcessLeft = 0; + size_t m_kExcessRight = 0; + + bool m_dovisc; + bool m_isFree; + bool m_usesLambda; + + //! Flag for activating two-point flame control + bool m_twoPointControl = false; + + //! Location of the left control point when two-point control is enabled + double m_zLeft = Undef; + + //! Temperature of the left control point when two-point control is enabled + double m_tLeft = Undef; + + //! Location of the right control point when two-point control is enabled + double m_zRight = Undef; + + //! Temperature of the right control point when two-point control is enabled + double m_tRight = Undef; + +public: + //! Location of the point where temperature is fixed + double m_zfixed = Undef; + + //! Temperature at the point used to fix the flame location + double m_tfixed = -1.0; + +private: + vector m_ybar; +}; + +} + +#endif diff --git a/include/cantera/oneD/IonFlow.h b/include/cantera/oneD/IonFlow.h index 484bdb51080..5392d9616dc 100644 --- a/include/cantera/oneD/IonFlow.h +++ b/include/cantera/oneD/IonFlow.h @@ -6,7 +6,7 @@ #ifndef CT_IONFLOW_H #define CT_IONFLOW_H -#include "cantera/oneD/StFlow.h" +#include "cantera/oneD/Flow1D.h" namespace Cantera { @@ -25,7 +25,7 @@ namespace Cantera * * @ingroup flowGroup */ -class IonFlow : public StFlow +class IonFlow : public Flow1D { public: IonFlow(ThermoPhase* ph = 0, size_t nsp = 1, size_t points = 1); diff --git a/include/cantera/oneD/StFlow.h b/include/cantera/oneD/StFlow.h index 8f96b13b0a6..40e826a7fcb 100644 --- a/include/cantera/oneD/StFlow.h +++ b/include/cantera/oneD/StFlow.h @@ -6,49 +6,21 @@ #ifndef CT_STFLOW_H #define CT_STFLOW_H -#include "Domain1D.h" -#include "cantera/base/Array.h" -#include "cantera/base/Solution.h" -#include "cantera/thermo/ThermoPhase.h" -#include "cantera/kinetics/Kinetics.h" +#include "Flow1D.h" namespace Cantera { -//------------------------------------------ -// constants -//------------------------------------------ - -//! Offsets of solution components in the 1D solution array. -enum offset -{ - c_offset_U //! axial velocity [m/s] - , c_offset_V //! strain rate - , c_offset_T //! temperature [kelvin] - , c_offset_L //! (1/r)dP/dr - , c_offset_E //! electric field - , c_offset_Uo //! oxidizer axial velocity [m/s] - , c_offset_Y //! mass fractions -}; - -class Transport; - -//! @defgroup flowGroup Flow Domains -//! One-dimensional flow domains. -//! @ingroup onedGroup - /** * This class represents 1D flow domains that satisfy the one-dimensional * similarity solution for chemically-reacting, axisymmetric flows. + * + * @deprecated To be removed after Cantera 3.1; replaced by Flow1D. * @ingroup flowGroup */ -class StFlow : public Domain1D +class StFlow : public Flow1D { public: - //-------------------------------- - // construction and destruction - //-------------------------------- - //! Create a new flow domain. //! @param ph Object representing the gas phase. This object will be used //! to evaluate all thermodynamic, kinetic, and transport properties. @@ -66,735 +38,28 @@ class StFlow : public Domain1D //! @param points initial number of grid points StFlow(shared_ptr sol, const string& id="", size_t points=1); - ~StFlow(); - - string domainType() const override; - - //! @name Problem Specification - //! @{ - - void setupGrid(size_t n, const double* z) override; - - void resetBadValues(double* xg) override; - - ThermoPhase& phase() { - return *m_thermo; - } - - Kinetics& kinetics() { - return *m_kin; - } - - void setKinetics(shared_ptr kin) override; - - void setTransport(shared_ptr trans) override; - - //! Set the transport model - //! @since New in %Cantera 3.0. - void setTransportModel(const string& trans); - - //! Retrieve transport model - //! @since New in %Cantera 3.0. - string transportModel() const; - - //! Enable thermal diffusion, also known as Soret diffusion. - //! Requires that multicomponent transport properties be - //! enabled to carry out calculations. - void enableSoret(bool withSoret) { - m_do_soret = withSoret; - } - bool withSoret() const { - return m_do_soret; - } - - //! Compute species diffusive fluxes with respect to - //! their mass fraction gradients (fluxGradientBasis = ThermoBasis::mass) - //! or mole fraction gradients (fluxGradientBasis = ThermoBasis::molar, default) - //! when using the mixture-averaged transport model. - //! @param fluxGradientBasis set flux computation to mass or mole basis - //! @since New in %Cantera 3.1. - void setFluxGradientBasis(ThermoBasis fluxGradientBasis); - - //! Compute species diffusive fluxes with respect to - //! their mass fraction gradients (fluxGradientBasis = ThermoBasis::mass) - //! or mole fraction gradients (fluxGradientBasis = ThermoBasis::molar, default) - //! when using the mixture-averaged transport model. - //! @return the basis used for flux computation (mass or mole fraction gradients) - //! @since New in %Cantera 3.1. - ThermoBasis fluxGradientBasis() const { - return m_fluxGradientBasis; - } - - //! Set the pressure. Since the flow equations are for the limit of small - //! Mach number, the pressure is very nearly constant throughout the flow. - void setPressure(double p) { - m_press = p; - } - - //! The current pressure [Pa]. - double pressure() const { - return m_press; - } - - //! Write the initial solution estimate into array x. - void _getInitialSoln(double* x) override; - - void _finalize(const double* x) override; - - //! Sometimes it is desired to carry out the simulation using a specified - //! temperature profile, rather than computing it by solving the energy - //! equation. This method specifies this profile. - void setFixedTempProfile(vector& zfixed, vector& tfixed) { - m_zfix = zfixed; - m_tfix = tfixed; - } - - /** - * Set the temperature fixed point at grid point j, and disable the energy - * equation so that the solution will be held to this value. - */ - void setTemperature(size_t j, double t) { - m_fixedtemp[j] = t; - m_do_energy[j] = false; - } - - //! The fixed temperature value at point j. - double T_fixed(size_t j) const { - return m_fixedtemp[j]; - } - - //! @} - - string componentName(size_t n) const override; - - size_t componentIndex(const string& name) const override; - - //! Returns true if the specified component is an active part of the solver state - virtual bool componentActive(size_t n) const; - - //! Print the solution. - void show(const double* x) override; - - shared_ptr asArray(const double* soln) const override; - void fromArray(SolutionArray& arr, double* soln) override; - - //! Set flow configuration for freely-propagating flames, using an internal point - //! with a fixed temperature as the condition to determine the inlet mass flux. - void setFreeFlow() { - m_dovisc = false; - m_isFree = true; - m_usesLambda = false; - } - - //! Set flow configuration for axisymmetric counterflow flames, using specified - //! inlet mass fluxes. - void setAxisymmetricFlow() { - m_dovisc = true; - m_isFree = false; - m_usesLambda = true; - } - - //! Set flow configuration for burner-stabilized flames, using specified inlet mass - //! fluxes. - void setUnstrainedFlow() { - m_dovisc = false; - m_isFree = false; - m_usesLambda = false; - } - - void solveEnergyEqn(size_t j=npos); - - //! Get the solving stage (used by IonFlow specialization) - //! @since New in %Cantera 3.0 - virtual size_t getSolvingStage() const; - - //! Solving stage mode for handling ionized species (used by IonFlow specialization) - //! - @c stage=1: the fluxes of charged species are set to zero - //! - @c stage=2: the electric field equation is solved, and the drift flux for - //! ionized species is evaluated - virtual void setSolvingStage(const size_t stage); - - //! Set to solve electric field in a point (used by IonFlow specialization) - virtual void solveElectricField(size_t j=npos); - - //! Set to fix voltage in a point (used by IonFlow specialization) - virtual void fixElectricField(size_t j=npos); - - //! Retrieve flag indicating whether electric field is solved or not (used by - //! IonFlow specialization) - virtual bool doElectricField(size_t j) const; - - //! Turn radiation on / off. - void enableRadiation(bool doRadiation) { - m_do_radiation = doRadiation; - } - - //! Returns `true` if the radiation term in the energy equation is enabled - bool radiationEnabled() const { - return m_do_radiation; - } - - //! Return radiative heat loss at grid point j - double radiativeHeatLoss(size_t j) const { - return m_qdotRadiation[j]; - } - - //! Set the emissivities for the boundary values - /*! - * Reads the emissivities for the left and right boundary values in the - * radiative term and writes them into the variables, which are used for the - * calculation. - */ - void setBoundaryEmissivities(double e_left, double e_right); - - //! Return emissivity at left boundary - double leftEmissivity() const { - return m_epsilon_left; - } - - //! Return emissivity at right boundary - double rightEmissivity() const { - return m_epsilon_right; - } + void eval(size_t j, double* x, double* r, integer* mask, double rdt) override; - void fixTemperature(size_t j=npos); + //! Evaluate all residual components at the right boundary. + virtual void evalRightBoundary(double* x, double* res, int* diag, double rdt); - /** - * @name Two-Point control method - * - * In this method two control points are designated in the 1D domain, and the value - * of the temperature at these points is fixed. The values of the control points are - * imposed and thus serve as a boundary condition that affects the solution of the - * governing equations in the 1D domain. The imposition of fixed points in the - * domain means that the original set of governing equations' boundary conditions - * would over-determine the problem. Thus, the boundary conditions are changed to - * reflect the fact that the control points are serving as internal boundary - * conditions. - * - * The imposition of the two internal boundary conditions requires that two other - * boundary conditions be changed. The first is the boundary condition for the - * continuity equation at the left boundary, which is changed to be a value that is - * derived from the solution at the left boundary. The second is the continuity - * boundary condition at the right boundary, which is also determined from the flow - * solution by using the oxidizer axial velocity equation variable to compute the - * mass flux at the right boundary. - * - * This method is based on the work of Nishioka et al. @cite nishioka1996 . - */ - //! @{ - - //! Returns the temperature at the left control point - double leftControlPointTemperature() const; - - //! Returns the z-coordinate of the left control point - double leftControlPointCoordinate() const; - - //! Sets the temperature of the left control point - void setLeftControlPointTemperature(double temperature); - - //! Sets the coordinate of the left control point - void setLeftControlPointCoordinate(double z_left); - - //! Returns the temperature at the right control point - double rightControlPointTemperature() const; - - //! Returns the z-coordinate of the right control point - double rightControlPointCoordinate() const; - - //! Sets the temperature of the right control point - void setRightControlPointTemperature(double temperature); - - //! Sets the coordinate of the right control point - void setRightControlPointCoordinate(double z_right); - - //! Sets the status of the two-point control - void enableTwoPointControl(bool twoPointControl); - - //! Returns the status of the two-point control - bool twoPointControlEnabled() const { - return m_twoPointControl; - } - //! @} - - bool doEnergy(size_t j) { - return m_do_energy[j]; - } - - //! Change the grid size. Called after grid refinement. - void resize(size_t components, size_t points) override; - - //! Set the gas object state to be consistent with the solution at point j. - void setGas(const double* x, size_t j); - - //! Set the gas state to be consistent with the solution at the midpoint - //! between j and j + 1. - void setGasAtMidpoint(const double* x, size_t j); - - double density(size_t j) const { - return m_rho[j]; - } - - /** - * Retrieve flag indicating whether flow is freely propagating. - * The flow is unstrained and the axial mass flow rate is not specified. - * For free flame propagation, the axial velocity is determined by the solver. - * @since New in %Cantera 3.0 - */ - bool isFree() const { - return m_isFree; - } - - /** - * Retrieve flag indicating whether flow uses radial momentum. - * If `true`, radial momentum equation for @f$ V @f$ as well as - * @f$ d\Lambda/dz = 0 @f$ are solved; if `false`, @f$ \Lambda(z) = 0 @f$ and - * @f$ V(z) = 0 @f$ by definition. - * @since New in %Cantera 3.0 - */ - bool isStrained() const { - return m_usesLambda; - } - - void setViscosityFlag(bool dovisc) { - m_dovisc = dovisc; - } - - /** - * Evaluate the residual functions for axisymmetric stagnation flow. - * If jGlobal == npos, the residual function is evaluated at all grid points. - * Otherwise, the residual function is only evaluated at grid points j-1, j, - * and j+1. This option is used to efficiently evaluate the Jacobian numerically. - * - * These residuals at all the boundary grid points are evaluated using a default - * boundary condition that may be modified by a boundary object that is attached - * to the domain. The boundary object connected will modify these equations by - * subtracting the boundary object's values for V, T, mdot, etc. As a result, - * these residual equations will force the solution variables to the values of - * the connected boundary object. - * - * @param jGlobal Global grid point at which to update the residual - * @param[in] xGlobal Global state vector - * @param[out] rsdGlobal Global residual vector - * @param[out] diagGlobal Global boolean mask indicating whether each solution - * component has a time derivative (1) or not (0). - * @param[in] rdt Reciprocal of the timestep (`rdt=0` implies steady-state.) - */ - void eval(size_t jGlobal, double* xGlobal, double* rsdGlobal, - integer* diagGlobal, double rdt) override; - - //! Index of the species on the left boundary with the largest mass fraction - size_t leftExcessSpecies() const { - return m_kExcessLeft; - } - - //! Index of the species on the right boundary with the largest mass fraction - size_t rightExcessSpecies() const { - return m_kExcessRight; - } + void evalContinuity(size_t j, double* x, double* r, int* diag, double rdt) override; protected: - AnyMap getMeta() const override; - void setMeta(const AnyMap& state) override; - double wdot(size_t k, size_t j) const { return m_wdot(k,j); } - //! Update the properties (thermo, transport, and diffusion flux). - //! This function is called in eval after the points which need - //! to be updated are defined. - virtual void updateProperties(size_t jg, double* x, size_t jmin, size_t jmax); - - /** - * Computes the radiative heat loss vector over points jmin to jmax and stores - * the data in the qdotRadiation variable. - * - * The simple radiation model used was established by Liu and Rogg - * @cite liu1991. This model considers the radiation of CO2 and H2O. - * - * This model uses the optically thin limit and the gray-gas approximation to - * simply calculate a volume specified heat flux out of the Planck absorption - * coefficients, the boundary emissivities and the temperature. Polynomial lines - * calculate the species Planck coefficients for H2O and CO2. The data for the - * lines are taken from the RADCAL program @cite RADCAL. - * The coefficients for the polynomials are taken from - * [TNF Workshop](https://tnfworkshop.org/radiation/) material. - */ - void computeRadiation(double* x, size_t jmin, size_t jmax); - - /** - * Evaluate the continuity equation residual. - * - * This function calculates the residual of the continuity equation - * @f[ - * \frac{d(\rho u)}{dz} + 2\rho V = 0 - * @f] - * - * Axisymmetric flame: - * The continuity equation propagates information from right-to-left. - * The @f$ \rho u @f$ at point 0 is dependent on @f$ \rho u @f$ at point 1, - * but not on @f$ \dot{m} @f$ from the inlet. - * - * Freely-propagating flame: - * The continuity equation propagates information away from a fixed temperature - * point that is set in the domain. - * - * Unstrained flame: - * A specified mass flux; the main example being burner-stabilized flames. - * - * The default boundary condition for the continuity equation is - * (@f$ u = 0 @f$) at the left and right boundary. - * - * @param[in] x Local domain state vector, includes variables like temperature, - * density, etc. - * @param[out] rsd Local domain residual vector that stores the continuity - * equation residuals. - * @param[out] diag Local domain diagonal matrix that controls whether an entry - * has a time-derivative (used by the solver). - * @param[in] rdt Reciprocal of the timestep. - * @param[in] jmin The index for the starting point in the local domain grid. - * @param[in] jmax The index for the ending point in the local domain grid. - */ - virtual void evalContinuity(double* x, double* rsd, int* diag, - double rdt, size_t jmin, size_t jmax); - - /** - * Evaluate the momentum equation residual. - * - * The function calculates the radial momentum equation defined as - * @f[ - * \rho u \frac{dV}{dz} + \rho V^2 = - * \frac{d}{dz}\left( \mu \frac{dV}{dz} \right) - \Lambda - * @f] - * - * The radial momentum equation is used for axisymmetric flows, and incorporates - * terms for time and spatial variations of radial velocity (@f$ V @f$). The - * default boundary condition is zero radial velocity (@f$ V @f$) at the left - * and right boundary. - * - * For argument explanation, see evalContinuity(). - */ - virtual void evalMomentum(double* x, double* rsd, int* diag, - double rdt, size_t jmin, size_t jmax); - - /** - * Evaluate the lambda equation residual. - * - * The function calculates the lambda equation as - * @f[ - * \frac{d\Lambda}{dz} = 0 - * @f] - * - * The lambda equation serves as an eigenvalue that allows the momentum - * equation and continuity equations to be simultaneously satisfied in - * axisymmetric flows. The lambda equation propagates information from - * left-to-right. The default boundary condition is @f$ \Lambda = 0 @f$ - * at the left and zero flux at the right boundary. - * - * For argument explanation, see evalContinuity(). - */ - virtual void evalLambda(double* x, double* rsd, int* diag, - double rdt, size_t jmin, size_t jmax); - - /** - * Evaluate the energy equation residual. - * - * The function calculates the energy equation: - * @f[ - * \rho c_p u \frac{dT}{dz} = - * \frac{d}{dz}\left( \lambda \frac{dT}{dz} \right) - * - \sum_k h_kW_k\dot{\omega}_k - * - \sum_k j_k \frac{dh_k}{dz} - * @f] - * - * The energy equation includes contributions from - * chemical reactions and diffusion. Default is zero temperature (@f$ T @f$) - * at the left and right boundaries. These boundary values are updated by the - * specific boundary object connected to the domain. - * - * For argument explanation, see evalContinuity(). - */ - virtual void evalEnergy(double* x, double* rsd, int* diag, - double rdt, size_t jmin, size_t jmax); - - /** - * Evaluate the species equations' residuals. - * - * The function calculates the species equations as - * @f[ - * \rho u \frac{dY_k}{dz} + \frac{dj_k}{dz} = W_k\dot{\omega}_k - * @f] - * - * The species equations include terms for temporal and spatial variations - * of species mass fractions (@f$ Y_k @f$). The default boundary condition is zero - * flux for species at the left and right boundary. - * - * For argument explanation, see evalContinuity(). - */ - virtual void evalSpecies(double* x, double* rsd, int* diag, - double rdt, size_t jmin, size_t jmax); - - /** - * Evaluate the electric field equation residual to be zero everywhere. - * - * The electric field equation is implemented in the IonFlow class. The default - * boundary condition is zero electric field (@f$ E @f$) at the boundary, - * and @f$ E @f$ is zero within the domain. - * - * For argument explanation, see evalContinuity(). - */ - virtual void evalElectricField(double* x, double* rsd, int* diag, - double rdt, size_t jmin, size_t jmax); - - /** - * Evaluate the oxidizer axial velocity equation residual. - * - * The function calculates the oxidizer axial velocity equation as - * @f[ - * \frac{dU_{o}}{dz} = 0 - * @f] - * - * This equation serves as a dummy equation that is used only in the context of - * two-point flame control, and serves as the way for two interior control points to - * be specified while maintaining block tridiagonal structure. The default boundary - * condition is @f$ U_o = 0 @f$ at the right and zero flux at the left boundary. - * - * For argument explanation, see evalContinuity(). - */ - virtual void evalUo(double* x, double* rsd, int* diag, - double rdt, size_t jmin, size_t jmax); - - /** - * Update the thermodynamic properties from point j0 to point j1 - * (inclusive), based on solution x. - * - * The gas state is set to be consistent with the solution at the - * points from j0 to j1. - * - * Properties that are computed and cached are: - * * #m_rho (density) - * * #m_wtm (mean molecular weight) - * * #m_cp (specific heat capacity) - * * #m_hk (species specific enthalpies) - * * #m_wdot (species production rates) - */ - void updateThermo(const double* x, size_t j0, size_t j1) { - for (size_t j = j0; j <= j1; j++) { - setGas(x,j); - m_rho[j] = m_thermo->density(); - m_wtm[j] = m_thermo->meanMolecularWeight(); - m_cp[j] = m_thermo->cp_mass(); - m_thermo->getPartialMolarEnthalpies(&m_hk(0, j)); - m_kin->getNetProductionRates(&m_wdot(0, j)); - } - } - - //! @name Solution components - //! @{ - - double T(const double* x, size_t j) const { - return x[index(c_offset_T, j)]; - } - double& T(double* x, size_t j) { - return x[index(c_offset_T, j)]; - } - double T_prev(size_t j) const { - return prevSoln(c_offset_T, j); - } - - double rho_u(const double* x, size_t j) const { - return m_rho[j]*x[index(c_offset_U, j)]; - } - - double u(const double* x, size_t j) const { - return x[index(c_offset_U, j)]; - } - - double V(const double* x, size_t j) const { - return x[index(c_offset_V, j)]; - } - double V_prev(size_t j) const { - return prevSoln(c_offset_V, j); - } - - double lambda(const double* x, size_t j) const { - return x[index(c_offset_L, j)]; - } - - //! Solution component for oxidizer velocity, @see evalUo - double Uo(const double* x, size_t j) const { - return x[index(c_offset_Uo, j)]; - } - - double Y(const double* x, size_t k, size_t j) const { - return x[index(c_offset_Y + k, j)]; - } - - double& Y(double* x, size_t k, size_t j) { - return x[index(c_offset_Y + k, j)]; - } - - double Y_prev(size_t k, size_t j) const { - return prevSoln(c_offset_Y + k, j); - } - - double X(const double* x, size_t k, size_t j) const { - return m_wtm[j]*Y(x,k,j)/m_wt[k]; - } - - double flux(size_t k, size_t j) const { - return m_flux(k, j); + //! Write the net production rates at point `j` into array `m_wdot` + void getWdot(double* x, size_t j) { + setGas(x,j); + m_kin->getNetProductionRates(&m_wdot(0,j)); } - //! @} - //! @name convective spatial derivatives. - //! - //! These use upwind differencing, assuming u(z) is negative - //! @{ - double dVdz(const double* x, size_t j) const { - size_t jloc = (u(x,j) > 0.0 ? j : j + 1); - return (V(x,jloc) - V(x,jloc-1))/m_dz[jloc-1]; - } - - double dYdz(const double* x, size_t k, size_t j) const { - size_t jloc = (u(x,j) > 0.0 ? j : j + 1); - return (Y(x,k,jloc) - Y(x,k,jloc-1))/m_dz[jloc-1]; - } - - double dTdz(const double* x, size_t j) const { - size_t jloc = (u(x,j) > 0.0 ? j : j + 1); - return (T(x,jloc) - T(x,jloc-1))/m_dz[jloc-1]; - } - //! @} - - double shear(const double* x, size_t j) const { - double c1 = m_visc[j-1]*(V(x,j) - V(x,j-1)); - double c2 = m_visc[j]*(V(x,j+1) - V(x,j)); - return 2.0*(c2/(z(j+1) - z(j)) - c1/(z(j) - z(j-1)))/(z(j+1) - z(j-1)); - } - - double divHeatFlux(const double* x, size_t j) const { - double c1 = m_tcon[j-1]*(T(x,j) - T(x,j-1)); - double c2 = m_tcon[j]*(T(x,j+1) - T(x,j)); - return -2.0*(c2/(z(j+1) - z(j)) - c1/(z(j) - z(j-1)))/(z(j+1) - z(j-1)); - } - - size_t mindex(size_t k, size_t j, size_t m) { - return m*m_nsp*m_nsp + m_nsp*j + k; - } - - //! Update the diffusive mass fluxes. - virtual void updateDiffFluxes(const double* x, size_t j0, size_t j1); - - //! Get the gradient of species specific molar enthalpies - virtual void grad_hk(const double* x, size_t j); - - //--------------------------------------------------------- - // member data - //--------------------------------------------------------- - - double m_press = -1.0; // pressure - - // grid parameters - vector m_dz; - - // mixture thermo properties - vector m_rho; - vector m_wtm; - - // species thermo properties - vector m_wt; - vector m_cp; - - // transport properties - vector m_visc; - vector m_tcon; - //! Array of size #m_nsp by #m_points for saving density times diffusion - //! coefficient times species molar mass divided by mean molecular weight - vector m_diff; - vector m_multidiff; - Array2D m_dthermal; - Array2D m_flux; - - //! Array of size #m_nsp by #m_points for saving molar enthalpies - Array2D m_hk; - - //! Array of size #m_nsp by #m_points-1 for saving enthalpy fluxes - Array2D m_dhk_dz; - - // production rates - Array2D m_wdot; - - size_t m_nsp; //!< Number of species in the mechanism - - ThermoPhase* m_thermo = nullptr; - Kinetics* m_kin = nullptr; - Transport* m_trans = nullptr; - - // boundary emissivities for the radiation calculations - double m_epsilon_left = 0.0; - double m_epsilon_right = 0.0; - - //! Indices within the ThermoPhase of the radiating species. First index is - //! for CO2, second is for H2O. - vector m_kRadiating; - - // flags - vector m_do_energy; - bool m_do_soret = false; - ThermoBasis m_fluxGradientBasis = ThermoBasis::molar; - vector m_do_species; - bool m_do_multicomponent = false; - - //! flag for the radiative heat loss - bool m_do_radiation = false; - - //! radiative heat loss vector - vector m_qdotRadiation; - - // fixed T and Y values - vector m_fixedtemp; - vector m_zfix; - vector m_tfix; - - //! Index of species with a large mass fraction at each boundary, for which - //! the mass fraction may be calculated as 1 minus the sum of the other mass - //! fractions - size_t m_kExcessLeft = 0; - size_t m_kExcessRight = 0; - - bool m_dovisc; - bool m_isFree; - bool m_usesLambda; - - //! Update the transport properties at grid points in the range from `j0` - //! to `j1`, based on solution `x`. - virtual void updateTransport(double* x, size_t j0, size_t j1); - - //! Flag for activating two-point flame control - bool m_twoPointControl = false; - - //! Location of the left control point when two-point control is enabled - double m_zLeft = Undef; - - //! Temperature of the left control point when two-point control is enabled - double m_tLeft = Undef; - - //! Location of the right control point when two-point control is enabled - double m_zRight = Undef; - - //! Temperature of the right control point when two-point control is enabled - double m_tRight = Undef; - -public: - //! Location of the point where temperature is fixed - double m_zfixed = Undef; - - //! Temperature at the point used to fix the flame location - double m_tfixed = -1.0; - -private: - vector m_ybar; + //! Evaluate the residual function. This function is called in eval + //! after updateProperties is called. + virtual void evalResidual(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax); }; } diff --git a/include/cantera/onedim.h b/include/cantera/onedim.h index a74c0c3eda4..eef1579b0c8 100644 --- a/include/cantera/onedim.h +++ b/include/cantera/onedim.h @@ -17,7 +17,7 @@ #include "oneD/Sim1D.h" #include "oneD/Domain1D.h" #include "oneD/Boundary1D.h" -#include "oneD/StFlow.h" +#include "oneD/Flow1D.h" #include "oneD/refine.h" #endif diff --git a/interfaces/cython/cantera/_onedim.pxd b/interfaces/cython/cantera/_onedim.pxd index cf2800bf30f..3863e088096 100644 --- a/interfaces/cython/cantera/_onedim.pxd +++ b/interfaces/cython/cantera/_onedim.pxd @@ -64,8 +64,8 @@ cdef extern from "cantera/oneD/Boundary1D.h": cbool coverageEnabled() -cdef extern from "cantera/oneD/StFlow.h": - cdef cppclass CxxStFlow "Cantera::StFlow" (CxxDomain1D): +cdef extern from "cantera/oneD/Flow1D.h": + cdef cppclass CxxFlow1D "Cantera::Flow1D" (CxxDomain1D): void setTransportModel(const string&) except +translate_exception string type() string transportModel() @@ -173,7 +173,7 @@ cdef class ReactingSurface1D(Boundary1D): cdef public Kinetics surface cdef class FlowBase(Domain1D): - cdef CxxStFlow* flow + cdef CxxFlow1D* flow cdef class Sim1D: cdef CxxSim1D* sim diff --git a/interfaces/cython/cantera/_onedim.pyx b/interfaces/cython/cantera/_onedim.pyx index 4da961939b7..e7335ab0a1a 100644 --- a/interfaces/cython/cantera/_onedim.pyx +++ b/interfaces/cython/cantera/_onedim.pyx @@ -449,7 +449,7 @@ cdef class FlowBase(Domain1D): self._domain = CxxNewDomain1D( stringify(self._domain_type), phase._base, stringify(name)) self.domain = self._domain.get() - self.flow = self.domain + self.flow = self.domain def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) diff --git a/interfaces/matlab_experimental/OneDim/Domain1D.m b/interfaces/matlab_experimental/OneDim/Domain1D.m index 0a992756ca5..658e8bb5a8b 100644 --- a/interfaces/matlab_experimental/OneDim/Domain1D.m +++ b/interfaces/matlab_experimental/OneDim/Domain1D.m @@ -108,12 +108,12 @@ function delete(d) function set.energyEnabled(d, flag) d.energyEnabled = flag; - ctFunc('stflow_solveEnergyEqn', d.domainID, int8(flag)); + ctFunc('flow1D_solveEnergyEqn', d.domainID, int8(flag)); end function set.soretEnabled(d, flag) d.soretEnabled = flag; - ctFunc('stflow_enableSoret', d.domainID, int8(flag)); + ctFunc('flow1D_enableSoret', d.domainID, int8(flag)); end %% Domain Get Methods @@ -359,7 +359,7 @@ function setTransientTolerances(d, component, rtol, atol) end function set.transport(d, itr) - ctFunc('stflow_setTransport', d.domainID, itr); + ctFunc('flow1D_setTransport', d.domainID, itr); end function setupGrid(d, grid) diff --git a/interfaces/matlab_experimental/OneDim/Flow1D.m b/interfaces/matlab_experimental/OneDim/Flow1D.m index 253699ed815..1da4c7fc1fd 100644 --- a/interfaces/matlab_experimental/OneDim/Flow1D.m +++ b/interfaces/matlab_experimental/OneDim/Flow1D.m @@ -31,11 +31,11 @@ %% Flow1D Class Methods function pressure = get.P(d) - pressure = ctFunc('stflow_pressure', d.domainID); + pressure = ctFunc('flow1D_pressure', d.domainID); end function set.P(d, p) - ctFunc('stflow_setPressure', d.domainID, p); + ctFunc('flow1D_setPressure', d.domainID, p); end function setFixedTempProfile(d, profile) @@ -57,11 +57,11 @@ function setFixedTempProfile(d, profile) if sz(1) == 2 l = length(profile(1, :)); - ctFunc('stflow_setFixedTempProfile', d.domainID, ... + ctFunc('flow1D_setFixedTempProfile', d.domainID, ... l, profile(1, :), l, profile(2, :)); elseif sz(2) == 2 l = length(profile(:, 1)); - ctFunc('stflow_setFixedTempProfile', d.domainID, ... + ctFunc('flow1D_setFixedTempProfile', d.domainID, ... l, profile(:, 1), l, profile(:, 2)); else error('Wrong temperature profile array shape.'); @@ -71,4 +71,4 @@ function setFixedTempProfile(d, profile) end -end \ No newline at end of file +end diff --git a/samples/cxx/flamespeed/flamespeed.cpp b/samples/cxx/flamespeed/flamespeed.cpp index 09e7749ee05..2859038d613 100644 --- a/samples/cxx/flamespeed/flamespeed.cpp +++ b/samples/cxx/flamespeed/flamespeed.cpp @@ -8,7 +8,7 @@ * * flamespeed [equivalence_ratio] [refine_grid] [loglevel] * - * Requires: cantera >= 3.0 + * Requires: cantera >= 3.1 * * .. tags:: C++, combustion, 1D flow, premixed flame, flame speed, saving output */ @@ -56,7 +56,7 @@ int flamespeed(double phi, bool refine_grid, int loglevel) //-------- step 1: create the flow ------------- - auto flow = newDomain("gas-flow", sol, "flow"); + auto flow = newDomain("gas-flow", sol, "flow"); flow->setFreeFlow(); // create an initial grid diff --git a/src/clib/ctonedim.cpp b/src/clib/ctonedim.cpp index 5bc04656108..4896bb9781f 100644 --- a/src/clib/ctonedim.cpp +++ b/src/clib/ctonedim.cpp @@ -398,13 +398,13 @@ extern "C" { } } - //------------------ stagnation flow domains -------------------- + //------------------ flow domains -------------------- - int stflow_new(int iph, int ikin, int itr, int itype) + int flow1D_new(int iph, int ikin, int itr, int itype) { try { auto ph = ThermoCabinet::at(iph); - auto x = make_shared(ph, ph->nSpecies(), 2); + auto x = make_shared(ph, ph->nSpecies(), 2); if (itype == 1) { x->setAxisymmetricFlow(); } else if (itype == 2) { @@ -420,47 +420,47 @@ extern "C" { } } - int stflow_setTransport(int i, int itr) + int flow1D_setTransport(int i, int itr) { try { - DomainCabinet::get(i).setTransport(TransportCabinet::at(itr)); + DomainCabinet::get(i).setTransport(TransportCabinet::at(itr)); return 0; } catch (...) { return handleAllExceptions(-1, ERR); } } - int stflow_enableSoret(int i, int iSoret) + int flow1D_enableSoret(int i, int iSoret) { try { bool withSoret = (iSoret > 0); - DomainCabinet::get(i).enableSoret(withSoret); + DomainCabinet::get(i).enableSoret(withSoret); return 0; } catch (...) { return handleAllExceptions(-1, ERR); } } - int stflow_setPressure(int i, double p) + int flow1D_setPressure(int i, double p) { try { - DomainCabinet::get(i).setPressure(p); + DomainCabinet::get(i).setPressure(p); return 0; } catch (...) { return handleAllExceptions(-1, ERR); } } - double stflow_pressure(int i) + double flow1D_pressure(int i) { try { - return DomainCabinet::get(i).pressure(); + return DomainCabinet::get(i).pressure(); } catch (...) { return handleAllExceptions(DERR, DERR); } } - int stflow_setFixedTempProfile(int i, size_t n, const double* pos, + int flow1D_setFixedTempProfile(int i, size_t n, const double* pos, size_t m, const double* temp) { try { @@ -469,20 +469,20 @@ extern "C" { vpos[j] = pos[j]; vtemp[j] = temp[j]; } - DomainCabinet::get(i).setFixedTempProfile(vpos, vtemp); + DomainCabinet::get(i).setFixedTempProfile(vpos, vtemp); return 0; } catch (...) { return handleAllExceptions(-1, ERR); } } - int stflow_solveEnergyEqn(int i, int flag) + int flow1D_solveEnergyEqn(int i, int flag) { try { if (flag > 0) { - DomainCabinet::get(i).solveEnergyEqn(npos); + DomainCabinet::get(i).solveEnergyEqn(npos); } else { - DomainCabinet::get(i).fixTemperature(npos); + DomainCabinet::get(i).fixTemperature(npos); } return 0; } catch (...) { @@ -490,6 +490,77 @@ extern "C" { } } + int stflow_new(int iph, int ikin, int itr, int itype) + { + try { + throw NotImplementedError("stflow_new", + "Function replaced by 'flow1D_new'."); + } catch (...) { + return handleAllExceptions(-1, ERR); + } + } + + int stflow_setTransport(int i, int itr) + { + try { + throw NotImplementedError("stflow_setTransport", + "Function replaced by 'flow1D_setTransport'."); + } catch (...) { + return handleAllExceptions(-1, ERR); + } + } + + int stflow_enableSoret(int i, int iSoret) + { + try { + throw NotImplementedError("stflow_enableSoret", + "Function replaced by 'flow1D_enableSoret'."); + } catch (...) { + return handleAllExceptions(-1, ERR); + } + } + + int stflow_setPressure(int i, double p) + { + try { + throw NotImplementedError("stflow_setPressure", + "Function replaced by 'flow1D_setPressure'."); + } catch (...) { + return handleAllExceptions(-1, ERR); + } + } + + double stflow_pressure(int i) + { + try { + throw NotImplementedError("stflow_pressure", + "Function replaced by 'flow1D_pressure'."); + } catch (...) { + return handleAllExceptions(DERR, DERR); + } + } + + int stflow_setFixedTempProfile(int i, size_t n, const double* pos, + size_t m, const double* temp) + { + try { + throw NotImplementedError("stflow_setFixedTempProfile", + "Function replaced by 'flow1D_setFixedTempProfile'."); + } catch (...) { + return handleAllExceptions(-1, ERR); + } + } + + int stflow_solveEnergyEqn(int i, int flag) + { + try { + throw NotImplementedError("stflow_solveEnergyEqn", + "Function replaced by 'flow1D_solveEnergyEqn'."); + } catch (...) { + return handleAllExceptions(-1, ERR); + } + } + //------------------- Sim1D -------------------------------------- int sim1D_new(size_t nd, const int* domains) diff --git a/src/oneD/Boundary1D.cpp b/src/oneD/Boundary1D.cpp index 5b8cc7d8185..f698e88d14e 100644 --- a/src/oneD/Boundary1D.cpp +++ b/src/oneD/Boundary1D.cpp @@ -6,7 +6,7 @@ #include "cantera/base/SolutionArray.h" #include "cantera/oneD/Boundary1D.h" #include "cantera/oneD/OneDim.h" -#include "cantera/oneD/StFlow.h" +#include "cantera/oneD/Flow1D.h" using namespace std; @@ -42,7 +42,7 @@ void Boundary1D::_init(size_t n) } m_left_loc = container().start(m_index-1); m_left_points = r.nPoints(); - m_flow_left = dynamic_cast(&r); + m_flow_left = dynamic_cast(&r); if (m_flow_left != nullptr) { m_phase_left = &m_flow_left->phase(); } @@ -64,7 +64,7 @@ void Boundary1D::_init(size_t n) m_right_nsp = 0; } m_right_loc = container().start(m_index+1); - m_flow_right = dynamic_cast(&r); + m_flow_right = dynamic_cast(&r); if (m_flow_right != nullptr) { m_phase_right = &m_flow_right->phase(); } @@ -212,7 +212,7 @@ void Inlet1D::eval(size_t jg, double* xg, double* rg, // When using two-point control, the mass flow rate at the left inlet is // not specified. Instead, the mass flow rate is dictated by the // velocity at the left inlet, which comes from the U variable. The - // default boundary condition specified in the StFlow.cpp file already + // default boundary condition specified in the Flow1D.cpp file already // handles this case. We only need to update the stored value of m_mdot // so that other equations that use the quantity are consistent. m_mdot = m_flow->density(0)*xb[c_offset_U]; diff --git a/src/oneD/DomainFactory.cpp b/src/oneD/DomainFactory.cpp index e11957e2259..de164e3a32b 100644 --- a/src/oneD/DomainFactory.cpp +++ b/src/oneD/DomainFactory.cpp @@ -5,8 +5,9 @@ #include "cantera/oneD/DomainFactory.h" #include "cantera/oneD/Boundary1D.h" -#include "cantera/oneD/StFlow.h" +#include "cantera/oneD/Flow1D.h" #include "cantera/oneD/IonFlow.h" +#include "cantera/oneD/StFlow.h" #include "cantera/transport/Transport.h" namespace Cantera @@ -39,37 +40,40 @@ DomainFactory::DomainFactory() return new ReactingSurf1D(solution, id); }); reg("gas-flow", [](shared_ptr solution, const string& id) { + return new Flow1D(solution, id); + }); + reg("legacy-flow", [](shared_ptr solution, const string& id) { return new StFlow(solution, id); }); reg("ion-flow", [](shared_ptr solution, const string& id) { return new IonFlow(solution, id); }); reg("free-flow", [](shared_ptr solution, const string& id) { - StFlow* ret; + Flow1D* ret; if (solution->transport()->transportModel() == "ionized-gas") { ret = new IonFlow(solution, id); } else { - ret = new StFlow(solution, id); + ret = new Flow1D(solution, id); } ret->setFreeFlow(); return ret; }); reg("axisymmetric-flow", [](shared_ptr solution, const string& id) { - StFlow* ret; + Flow1D* ret; if (solution->transport()->transportModel() == "ionized-gas") { ret = new IonFlow(solution, id); } else { - ret = new StFlow(solution, id); + ret = new Flow1D(solution, id); } ret->setAxisymmetricFlow(); return ret; }); reg("unstrained-flow", [](shared_ptr solution, const string& id) { - StFlow* ret; + Flow1D* ret; if (solution->transport()->transportModel() == "ionized-gas") { ret = new IonFlow(solution, id); } else { - ret = new StFlow(solution, id); + ret = new Flow1D(solution, id); } ret->setUnstrainedFlow(); return ret; diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp new file mode 100644 index 00000000000..cb82b164d13 --- /dev/null +++ b/src/oneD/Flow1D.cpp @@ -0,0 +1,1288 @@ +//! @file Flow1D.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. + +#include "cantera/base/SolutionArray.h" +#include "cantera/oneD/Flow1D.h" +#include "cantera/oneD/refine.h" +#include "cantera/transport/Transport.h" +#include "cantera/transport/TransportFactory.h" +#include "cantera/numerics/funcs.h" +#include "cantera/base/global.h" + +using namespace std; + +namespace Cantera +{ + +Flow1D::Flow1D(ThermoPhase* ph, size_t nsp, size_t points) : + Domain1D(nsp+c_offset_Y, points), + m_nsp(nsp) +{ + m_points = points; + + if (ph == 0) { + return; // used to create a dummy object + } + m_thermo = ph; + + size_t nsp2 = m_thermo->nSpecies(); + if (nsp2 != m_nsp) { + m_nsp = nsp2; + Domain1D::resize(m_nsp+c_offset_Y, points); + } + + // make a local copy of the species molecular weight vector + m_wt = m_thermo->molecularWeights(); + + // set pressure based on associated thermo object + setPressure(m_thermo->pressure()); + + // the species mass fractions are the last components in the solution + // vector, so the total number of components is the number of species + // plus the offset of the first mass fraction. + m_nv = c_offset_Y + m_nsp; + + // enable all species equations by default + m_do_species.resize(m_nsp, true); + + // but turn off the energy equation at all points + m_do_energy.resize(m_points,false); + + m_diff.resize(m_nsp*m_points); + m_multidiff.resize(m_nsp*m_nsp*m_points); + m_flux.resize(m_nsp,m_points); + m_wdot.resize(m_nsp,m_points, 0.0); + m_hk.resize(m_nsp, m_points, 0.0); + m_dhk_dz.resize(m_nsp, m_points - 1, 0.0); + m_ybar.resize(m_nsp); + m_qdotRadiation.resize(m_points, 0.0); + + //-------------- default solution bounds -------------------- + setBounds(c_offset_U, -1e20, 1e20); // no bounds on u + setBounds(c_offset_V, -1e20, 1e20); // no bounds on V + setBounds(c_offset_T, 200.0, 2*m_thermo->maxTemp()); // temperature bounds + setBounds(c_offset_L, -1e20, 1e20); // lambda should be negative + setBounds(c_offset_E, -1e20, 1e20); // no bounds on electric field + setBounds(c_offset_Uo, -1e20, 1e20); // no bounds on Uo + // mass fraction bounds + for (size_t k = 0; k < m_nsp; k++) { + setBounds(c_offset_Y+k, -1.0e-7, 1.0e5); + } + + //-------------------- grid refinement ------------------------- + m_refiner->setActive(c_offset_U, false); + m_refiner->setActive(c_offset_V, false); + m_refiner->setActive(c_offset_T, false); + m_refiner->setActive(c_offset_L, false); + m_refiner->setActive(c_offset_Uo, false); + + vector gr; + for (size_t ng = 0; ng < m_points; ng++) { + gr.push_back(1.0*ng/m_points); + } + setupGrid(m_points, gr.data()); + + // Find indices for radiating species + m_kRadiating.resize(2, npos); + m_kRadiating[0] = m_thermo->speciesIndex("CO2"); + m_kRadiating[1] = m_thermo->speciesIndex("H2O"); +} + +Flow1D::Flow1D(shared_ptr th, size_t nsp, size_t points) + : Flow1D(th.get(), nsp, points) +{ + auto sol = Solution::create(); + sol->setThermo(th); + setSolution(sol); +} + +Flow1D::Flow1D(shared_ptr sol, const string& id, size_t points) + : Flow1D(sol->thermo().get(), sol->thermo()->nSpecies(), points) +{ + setSolution(sol); + m_id = id; + m_kin = m_solution->kinetics().get(); + m_trans = m_solution->transport().get(); + if (m_trans->transportModel() == "none") { + throw CanteraError("Flow1D::Flow1D", + "An appropriate transport model\nshould be set when instantiating the " + "Solution ('gas') object."); + } + m_solution->registerChangedCallback(this, [this]() { + setKinetics(m_solution->kinetics()); + setTransport(m_solution->transport()); + }); +} + +Flow1D::~Flow1D() +{ + if (m_solution) { + m_solution->removeChangedCallback(this); + } +} + +string Flow1D::domainType() const { + if (m_isFree) { + return "free-flow"; + } + if (m_usesLambda) { + return "axisymmetric-flow"; + } + return "unstrained-flow"; +} + +void Flow1D::setKinetics(shared_ptr kin) +{ + m_kin = kin.get(); + m_solution->setKinetics(kin); +} + +void Flow1D::setTransport(shared_ptr trans) +{ + if (!trans) { + throw CanteraError("Flow1D::setTransport", "Unable to set empty transport."); + } + m_trans = trans.get(); + if (m_trans->transportModel() == "none") { + throw CanteraError("Flow1D::setTransport", "Invalid Transport model 'none'."); + } + m_do_multicomponent = (m_trans->transportModel() == "multicomponent" || + m_trans->transportModel() == "multicomponent-CK"); + + m_diff.resize(m_nsp * m_points); + if (m_do_multicomponent) { + m_multidiff.resize(m_nsp * m_nsp*m_points); + m_dthermal.resize(m_nsp, m_points, 0.0); + } + m_solution->setTransport(trans); +} + +void Flow1D::resize(size_t ncomponents, size_t points) +{ + Domain1D::resize(ncomponents, points); + m_rho.resize(m_points, 0.0); + m_wtm.resize(m_points, 0.0); + m_cp.resize(m_points, 0.0); + m_visc.resize(m_points, 0.0); + m_tcon.resize(m_points, 0.0); + + m_diff.resize(m_nsp*m_points); + if (m_do_multicomponent) { + m_multidiff.resize(m_nsp*m_nsp*m_points); + m_dthermal.resize(m_nsp, m_points, 0.0); + } + m_flux.resize(m_nsp,m_points); + m_wdot.resize(m_nsp,m_points, 0.0); + m_hk.resize(m_nsp, m_points, 0.0); + m_dhk_dz.resize(m_nsp, m_points - 1, 0.0); + m_do_energy.resize(m_points,false); + m_qdotRadiation.resize(m_points, 0.0); + m_fixedtemp.resize(m_points); + + m_dz.resize(m_points-1); + m_z.resize(m_points); +} + +void Flow1D::setupGrid(size_t n, const double* z) +{ + resize(m_nv, n); + + m_z[0] = z[0]; + for (size_t j = 1; j < m_points; j++) { + if (z[j] <= z[j-1]) { + throw CanteraError("Flow1D::setupGrid", + "grid points must be monotonically increasing"); + } + m_z[j] = z[j]; + m_dz[j-1] = m_z[j] - m_z[j-1]; + } +} + +void Flow1D::resetBadValues(double* xg) +{ + double* x = xg + loc(); + for (size_t j = 0; j < m_points; j++) { + double* Y = x + m_nv*j + c_offset_Y; + m_thermo->setMassFractions(Y); + m_thermo->getMassFractions(Y); + } +} + +void Flow1D::setTransportModel(const string& trans) +{ + m_solution->setTransportModel(trans); +} + +string Flow1D::transportModel() const { + return m_trans->transportModel(); +} + +void Flow1D::setFluxGradientBasis(ThermoBasis fluxGradientBasis) { + m_fluxGradientBasis = fluxGradientBasis; + if (transportModel() != "mixture-averaged-CK" && + transportModel() != "mixture-averaged") { + warn_user("Flow1D::setFluxGradientBasis", + "Setting fluxGradientBasis only affects " + "the mixture-averaged diffusion model."); + } +} + +void Flow1D::_getInitialSoln(double* x) +{ + for (size_t j = 0; j < m_points; j++) { + T(x,j) = m_thermo->temperature(); + m_thermo->getMassFractions(&Y(x, 0, j)); + m_rho[j] = m_thermo->density(); + } +} + +void Flow1D::setGas(const double* x, size_t j) +{ + m_thermo->setTemperature(T(x,j)); + const double* yy = x + m_nv*j + c_offset_Y; + m_thermo->setMassFractions_NoNorm(yy); + m_thermo->setPressure(m_press); +} + +void Flow1D::setGasAtMidpoint(const double* x, size_t j) +{ + m_thermo->setTemperature(0.5*(T(x,j)+T(x,j+1))); + const double* yyj = x + m_nv*j + c_offset_Y; + const double* yyjp = x + m_nv*(j+1) + c_offset_Y; + for (size_t k = 0; k < m_nsp; k++) { + m_ybar[k] = 0.5*(yyj[k] + yyjp[k]); + } + m_thermo->setMassFractions_NoNorm(m_ybar.data()); + m_thermo->setPressure(m_press); +} + +void Flow1D::_finalize(const double* x) +{ + if (!m_do_multicomponent && m_do_soret) { + throw CanteraError("Flow1D::_finalize", + "Thermal diffusion (the Soret effect) is enabled, and requires " + "using a multicomponent transport model."); + } + + size_t nz = m_zfix.size(); + bool e = m_do_energy[0]; + for (size_t j = 0; j < m_points; j++) { + if (e || nz == 0) { + m_fixedtemp[j] = T(x, j); + } else { + double zz = (z(j) - z(0))/(z(m_points - 1) - z(0)); + double tt = linearInterp(zz, m_zfix, m_tfix); + m_fixedtemp[j] = tt; + } + } + if (e) { + solveEnergyEqn(); + } + + if (m_isFree) { + // If the domain contains the temperature fixed point, make sure that it + // is correctly set. This may be necessary when the grid has been modified + // externally. + if (m_tfixed != Undef) { + for (size_t j = 0; j < m_points; j++) { + if (z(j) == m_zfixed) { + return; // fixed point is already set correctly + } + } + + for (size_t j = 0; j < m_points - 1; j++) { + // Find where the temperature profile crosses the current + // fixed temperature. + if ((T(x, j) - m_tfixed) * (T(x, j+1) - m_tfixed) <= 0.0) { + m_tfixed = T(x, j+1); + m_zfixed = z(j+1); + return; + } + } + } + } +} + +void Flow1D::eval(size_t jGlobal, double* xGlobal, double* rsdGlobal, + integer* diagGlobal, double rdt) +{ + // If evaluating a Jacobian, and the global point is outside the domain of + // influence for this domain, then skip evaluating the residual + if (jGlobal != npos && (jGlobal + 1 < firstPoint() || jGlobal > lastPoint() + 1)) { + return; + } + + // start of local part of global arrays + double* x = xGlobal + loc(); + double* rsd = rsdGlobal + loc(); + integer* diag = diagGlobal + loc(); + + size_t jmin, jmax; + if (jGlobal == npos) { // evaluate all points + jmin = 0; + jmax = m_points - 1; + } else { // evaluate points for Jacobian + size_t jpt = (jGlobal == 0) ? 0 : jGlobal - firstPoint(); + jmin = std::max(jpt, 1) - 1; + jmax = std::min(jpt+1,m_points-1); + } + + updateProperties(jGlobal, x, jmin, jmax); + + if (m_do_radiation) { // Calculation of qdotRadiation + computeRadiation(x, jmin, jmax); + } + + evalContinuity(x, rsd, diag, rdt, jmin, jmax); + evalMomentum(x, rsd, diag, rdt, jmin, jmax); + evalEnergy(x, rsd, diag, rdt, jmin, jmax); + evalLambda(x, rsd, diag, rdt, jmin, jmax); + evalElectricField(x, rsd, diag, rdt, jmin, jmax); + evalUo(x, rsd, diag, rdt, jmin, jmax); + evalSpecies(x, rsd, diag, rdt, jmin, jmax); +} + +void Flow1D::updateProperties(size_t jg, double* x, size_t jmin, size_t jmax) +{ + // properties are computed for grid points from j0 to j1 + size_t j0 = std::max(jmin, 1) - 1; + size_t j1 = std::min(jmax+1,m_points-1); + + updateThermo(x, j0, j1); + if (jg == npos || m_force_full_update) { + // update transport properties only if a Jacobian is not being + // evaluated, or if specifically requested + updateTransport(x, j0, j1); + } + if (jg == npos) { + double* Yleft = x + index(c_offset_Y, jmin); + m_kExcessLeft = distance(Yleft, max_element(Yleft, Yleft + m_nsp)); + double* Yright = x + index(c_offset_Y, jmax); + m_kExcessRight = distance(Yright, max_element(Yright, Yright + m_nsp)); + } + + // update the species diffusive mass fluxes whether or not a + // Jacobian is being evaluated + updateDiffFluxes(x, j0, j1); +} + +void Flow1D::updateTransport(double* x, size_t j0, size_t j1) +{ + if (m_do_multicomponent) { + for (size_t j = j0; j < j1; j++) { + setGasAtMidpoint(x,j); + double wtm = m_thermo->meanMolecularWeight(); + double rho = m_thermo->density(); + m_visc[j] = (m_dovisc ? m_trans->viscosity() : 0.0); + m_trans->getMultiDiffCoeffs(m_nsp, &m_multidiff[mindex(0,0,j)]); + + // Use m_diff as storage for the factor outside the summation + for (size_t k = 0; k < m_nsp; k++) { + m_diff[k+j*m_nsp] = m_wt[k] * rho / (wtm*wtm); + } + + m_tcon[j] = m_trans->thermalConductivity(); + if (m_do_soret) { + m_trans->getThermalDiffCoeffs(m_dthermal.ptrColumn(0) + j*m_nsp); + } + } + } else { // mixture averaged transport + for (size_t j = j0; j < j1; j++) { + setGasAtMidpoint(x,j); + m_visc[j] = (m_dovisc ? m_trans->viscosity() : 0.0); + + if (m_fluxGradientBasis == ThermoBasis::molar) { + m_trans->getMixDiffCoeffs(&m_diff[j*m_nsp]); + } else { + m_trans->getMixDiffCoeffsMass(&m_diff[j*m_nsp]); + } + + double rho = m_thermo->density(); + + if (m_fluxGradientBasis == ThermoBasis::molar) { + double wtm = m_thermo->meanMolecularWeight(); + for (size_t k=0; k < m_nsp; k++) { + m_diff[k+j*m_nsp] *= m_wt[k] * rho / wtm; + } + } else { + for (size_t k=0; k < m_nsp; k++) { + m_diff[k+j*m_nsp] *= rho; + } + } + m_tcon[j] = m_trans->thermalConductivity(); + } + } +} + +void Flow1D::updateDiffFluxes(const double* x, size_t j0, size_t j1) +{ + if (m_do_multicomponent) { + for (size_t j = j0; j < j1; j++) { + double dz = z(j+1) - z(j); + for (size_t k = 0; k < m_nsp; k++) { + double sum = 0.0; + for (size_t m = 0; m < m_nsp; m++) { + sum += m_wt[m] * m_multidiff[mindex(k,m,j)] * (X(x,m,j+1)-X(x,m,j)); + } + m_flux(k,j) = sum * m_diff[k+j*m_nsp] / dz; + } + } + } else { + for (size_t j = j0; j < j1; j++) { + double sum = 0.0; + double dz = z(j+1) - z(j); + if (m_fluxGradientBasis == ThermoBasis::molar) { + for (size_t k = 0; k < m_nsp; k++) { + m_flux(k,j) = m_diff[k+m_nsp*j] * (X(x,k,j) - X(x,k,j+1))/dz; + sum -= m_flux(k,j); + } + } else { + for (size_t k = 0; k < m_nsp; k++) { + m_flux(k,j) = m_diff[k+m_nsp*j] * (Y(x,k,j) - Y(x,k,j+1))/dz; + sum -= m_flux(k,j); + } + } + // correction flux to insure that \sum_k Y_k V_k = 0. + for (size_t k = 0; k < m_nsp; k++) { + m_flux(k,j) += sum*Y(x,k,j); + } + } + } + + if (m_do_soret) { + for (size_t m = j0; m < j1; m++) { + double gradlogT = 2.0 * (T(x,m+1) - T(x,m)) / + ((T(x,m+1) + T(x,m)) * (z(m+1) - z(m))); + for (size_t k = 0; k < m_nsp; k++) { + m_flux(k,m) -= m_dthermal(k,m)*gradlogT; + } + } + } +} + +void Flow1D::computeRadiation(double* x, size_t jmin, size_t jmax) +{ + // Variable definitions for the Planck absorption coefficient and the + // radiation calculation: + double k_P_ref = 1.0*OneAtm; + + // Polynomial coefficients: + const double c_H2O[6] = {-0.23093, -1.12390, 9.41530, -2.99880, + 0.51382, -1.86840e-5}; + const double c_CO2[6] = {18.741, -121.310, 273.500, -194.050, + 56.310, -5.8169}; + + // Calculation of the two boundary values + double boundary_Rad_left = m_epsilon_left * StefanBoltz * pow(T(x, 0), 4); + double boundary_Rad_right = m_epsilon_right * StefanBoltz * pow(T(x, m_points - 1), 4); + + for (size_t j = jmin; j < jmax; j++) { + // calculation of the mean Planck absorption coefficient + double k_P = 0; + // Absorption coefficient for H2O + if (m_kRadiating[1] != npos) { + double k_P_H2O = 0; + for (size_t n = 0; n <= 5; n++) { + k_P_H2O += c_H2O[n] * pow(1000 / T(x, j), (double) n); + } + k_P_H2O /= k_P_ref; + k_P += m_press * X(x, m_kRadiating[1], j) * k_P_H2O; + } + // Absorption coefficient for CO2 + if (m_kRadiating[0] != npos) { + double k_P_CO2 = 0; + for (size_t n = 0; n <= 5; n++) { + k_P_CO2 += c_CO2[n] * pow(1000 / T(x, j), (double) n); + } + k_P_CO2 /= k_P_ref; + k_P += m_press * X(x, m_kRadiating[0], j) * k_P_CO2; + } + + // Calculation of the radiative heat loss term + double radiative_heat_loss = 2 * k_P *(2 * StefanBoltz * pow(T(x, j), 4) + - boundary_Rad_left - boundary_Rad_right); + + // set the radiative heat loss vector + m_qdotRadiation[j] = radiative_heat_loss; + } +} + +void Flow1D::evalContinuity(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax) +{ + // The left boundary has the same form for all cases. + if (jmin == 0) { // left boundary + rsd[index(c_offset_U,jmin)] = -(rho_u(x,jmin + 1) - rho_u(x,jmin))/m_dz[jmin] + -(density(jmin + 1)*V(x,jmin + 1) + + density(jmin)*V(x,jmin)); + + diag[index(c_offset_U,jmin)] = 0; // Algebraic constraint + } + + if (jmax == m_points - 1) { // right boundary + if (m_usesLambda) { // axisymmetric flow + rsd[index(c_offset_U, jmax)] = rho_u(x, jmax); + } else { // right boundary (same for unstrained/free-flow) + rsd[index(c_offset_U, jmax)] = rho_u(x, jmax) - rho_u(x, jmax - 1); + } + diag[index(c_offset_U, jmax)] = 0; // Algebraic constraint + } + + // j0 and j1 are constrained to only interior points + size_t j0 = std::max(jmin, 1); + size_t j1 = std::min(jmax, m_points - 2); + if (m_usesLambda) { // "axisymmetric-flow" + for (size_t j = j0; j <= j1; j++) { // interior points + // For "axisymmetric-flow", the continuity equation propagates the + // mass flow rate information to the left (j+1 -> j) from the value + // specified at the right boundary. The lambda information propagates + // in the opposite direction. + rsd[index(c_offset_U,j)] = -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j] + -(density(j+1)*V(x,j+1) + density(j)*V(x,j)); + + diag[index(c_offset_U, j)] = 0; // Algebraic constraint + } + } else if (m_isFree) { // "free-flow" + for (size_t j = j0; j <= j1; j++) { + // terms involving V are zero as V=0 by definition + if (grid(j) > m_zfixed) { + rsd[index(c_offset_U,j)] = -(rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1]; + } else if (grid(j) == m_zfixed) { + if (m_do_energy[j]) { + rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed); + } else { + rsd[index(c_offset_U,j)] = (rho_u(x,j) - m_rho[0]*0.3); // why 0.3? + } + } else { // grid(j < m_zfixed + rsd[index(c_offset_U,j)] = -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j]; + } + diag[index(c_offset_U, j)] = 0; // Algebraic constraint + } + } else { // "unstrained-flow" (fixed mass flow rate) + for (size_t j = j0; j <= j1; j++) { + rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j - 1); + diag[index(c_offset_U, j)] = 0; // Algebraic constraint + } + } +} + +void Flow1D::evalMomentum(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax) +{ + if (!m_usesLambda) { //disable this equation + for (size_t j = jmin; j <= jmax; j++) { + rsd[index(c_offset_V, j)] = V(x, j); + diag[index(c_offset_V, j)] = 0; + } + return; + } + + if (jmin == 0) { // left boundary + rsd[index(c_offset_V,jmin)] = V(x,jmin); + } + + if (jmax == m_points - 1) { // right boundary + rsd[index(c_offset_V,jmax)] = V(x,jmax); + diag[index(c_offset_V,jmax)] = 0; + } + + // j0 and j1 are constrained to only interior points + size_t j0 = std::max(jmin, 1); + size_t j1 = std::min(jmax, m_points - 2); + for (size_t j = j0; j <= j1; j++) { // interior points + rsd[index(c_offset_V,j)] = (shear(x, j) - lambda(x, j) + - rho_u(x, j) * dVdz(x, j) + - m_rho[j] * V(x, j) * V(x, j)) / m_rho[j] + - rdt * (V(x, j) - V_prev(j)); + diag[index(c_offset_V, j)] = 1; + } +} + +void Flow1D::evalLambda(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax) +{ + if (!m_usesLambda) { // disable this equation + for (size_t j = jmin; j <= jmax; j++) { + rsd[index(c_offset_L, j)] = lambda(x, j); + diag[index(c_offset_L, j)] = 0; + } + return; + } + + if (jmin == 0) { // left boundary + if (m_twoPointControl) { + rsd[index(c_offset_L, jmin)] = lambda(x,jmin+1) - lambda(x,jmin); + } else { + rsd[index(c_offset_L, jmin)] = -rho_u(x, jmin); + } + } + + if (jmax == m_points - 1) { // right boundary + rsd[index(c_offset_L, jmax)] = lambda(x, jmax) - lambda(x, jmax-1); + diag[index(c_offset_L, jmax)] = 0; + } + + // j0 and j1 are constrained to only interior points + size_t j0 = std::max(jmin, 1); + size_t j1 = std::min(jmax, m_points - 2); + for (size_t j = j0; j <= j1; j++) { // interior points + if (m_twoPointControl) { + if (grid(j) == m_zLeft) { + rsd[index(c_offset_L, j)] = T(x,j) - m_tLeft; + } else if (grid(j) > m_zLeft) { + rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1); + } else if (grid(j) < m_zLeft) { + rsd[index(c_offset_L, j)] = lambda(x,j+1) - lambda(x,j); + } + } else { + rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1); + } + diag[index(c_offset_L, j)] = 0; + } +} + +void Flow1D::evalEnergy(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax) +{ + if (jmin == 0) { // left boundary + rsd[index(c_offset_T,jmin)] = T(x,jmin); + } + + if (jmax == m_points - 1) { // right boundary + rsd[index(c_offset_T, jmax)] = T(x,jmax); + } + + // j0 and j1 are constrained to only interior points + size_t j0 = std::max(jmin, 1); + size_t j1 = std::min(jmax, m_points - 2); + for (size_t j = j0; j <= j1; j++) { + if (m_do_energy[j]) { + grad_hk(x, j); + double sum = 0.0; + for (size_t k = 0; k < m_nsp; k++) { + double flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j)); + sum += m_wdot(k,j)*m_hk(k,j); + sum += flxk * m_dhk_dz(k,j) / m_wt[k]; + } + + rsd[index(c_offset_T, j)] = - m_cp[j]*rho_u(x,j)*dTdz(x,j) + - divHeatFlux(x,j) - sum; + rsd[index(c_offset_T, j)] /= (m_rho[j]*m_cp[j]); + rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j)); + rsd[index(c_offset_T, j)] -= (m_qdotRadiation[j] / (m_rho[j] * m_cp[j])); + diag[index(c_offset_T, j)] = 1; + } else { + // residual equations if the energy equation is disabled + rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j); + diag[index(c_offset_T, j)] = 0; + } + } +} + +void Flow1D::evalUo(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax) +{ + if (!m_twoPointControl) { // disable this equation + for (size_t j = jmin; j <= jmax; j++) { + rsd[index(c_offset_Uo, j)] = Uo(x, j); + diag[index(c_offset_Uo, j)] = 0; + } + return; + } + + if (jmin == 0) { // left boundary + rsd[index(c_offset_Uo, jmin)] = Uo(x,jmin+1) - Uo(x,jmin); + diag[index(c_offset_Uo, jmin)] = 0; + } + + if (jmax == m_points - 1) { // right boundary + if(m_twoPointControl) { + rsd[index(c_offset_Uo, jmax)] = Uo(x,jmax) - Uo(x,jmax-1); + } + diag[index(c_offset_Uo, jmax)] = 0; + } + + // j0 and j1 are constrained to only interior points + size_t j0 = std::max(jmin, 1); + size_t j1 = std::min(jmax, m_points - 2); + for (size_t j = j0; j <= j1; j++) { // interior points + if (m_twoPointControl) { + if (grid(j) == m_zRight) { + rsd[index(c_offset_Uo, j)] = T(x,j) - m_tRight; + } else if (grid(j) > m_zRight) { + rsd[index(c_offset_Uo, j)] = Uo(x,j) - Uo(x,j-1); + } else if (grid(j) < m_zRight) { + rsd[index(c_offset_Uo, j)] = Uo(x,j+1) - Uo(x,j); + } + } else { + rsd[index(c_offset_Uo, j)] = Uo(x,j+1) - Uo(x,j); + } + diag[index(c_offset_Uo, j)] = 0; + } +} + +void Flow1D::evalSpecies(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax) +{ + if (jmin == 0) { // left boundary + double sum = 0.0; + for (size_t k = 0; k < m_nsp; k++) { + sum += Y(x,k,jmin); + rsd[index(c_offset_Y + k, jmin)] = -(m_flux(k,jmin) + + rho_u(x,jmin) * Y(x,k,jmin)); + } + rsd[index(c_offset_Y + leftExcessSpecies(), jmin)] = 1.0 - sum; + } + + if (jmax == m_points - 1) { // right boundary + double sum = 0.0; + for (size_t k = 0; k < m_nsp; k++) { + sum += Y(x,k,jmax); + rsd[index(k+c_offset_Y,jmax)] = m_flux(k,jmax-1) + + rho_u(x,jmax)*Y(x,k,jmax); + } + rsd[index(c_offset_Y + rightExcessSpecies(), jmax)] = 1.0 - sum; + diag[index(c_offset_Y + rightExcessSpecies(), jmax)] = 0; + } + + // j0 and j1 are constrained to only interior points + size_t j0 = std::max(jmin, 1); + size_t j1 = std::min(jmax, m_points - 2); + for (size_t j = j0; j <= j1; j++) { + for (size_t k = 0; k < m_nsp; k++) { + double convec = rho_u(x,j)*dYdz(x,k,j); + double diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1)) / (z(j+1) - z(j-1)); + rsd[index(c_offset_Y + k, j)] = (m_wt[k]*m_wdot(k,j) + - convec - diffus)/m_rho[j] + - rdt*(Y(x,k,j) - Y_prev(k,j)); + diag[index(c_offset_Y + k, j)] = 1; + } + } +} + +void Flow1D::evalElectricField(double* x, double* rsd, int* diag, + double rdt, size_t jmin, size_t jmax) +{ + for (size_t j = jmin; j <= jmax; j++) { + // The same value is used for left/right/interior points + rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)]; + } +} + +void Flow1D::evalContinuity(size_t j, double* x, double* rsd, int* diag, double rdt) +{ + throw CanteraError("Flow1D::evalContinuity", + "Overloaded by StFlow; to be removed after Cantera 3.1"); +} + +void Flow1D::show(const double* x) +{ + writelog(" Pressure: {:10.4g} Pa\n", m_press); + + Domain1D::show(x); + + if (m_do_radiation) { + writeline('-', 79, false, true); + writelog("\n z radiative heat loss"); + writeline('-', 79, false, true); + for (size_t j = 0; j < m_points; j++) { + writelog("\n {:10.4g} {:10.4g}", m_z[j], m_qdotRadiation[j]); + } + writelog("\n"); + } +} + +string Flow1D::componentName(size_t n) const +{ + switch (n) { + case c_offset_U: + return "velocity"; + case c_offset_V: + return "spread_rate"; + case c_offset_T: + return "T"; + case c_offset_L: + return "lambda"; + case c_offset_E: + return "eField"; + case c_offset_Uo: + return "Uo"; + default: + if (n >= c_offset_Y && n < (c_offset_Y + m_nsp)) { + return m_thermo->speciesName(n - c_offset_Y); + } else { + return ""; + } + } +} + +size_t Flow1D::componentIndex(const string& name) const +{ + if (name=="velocity") { + return c_offset_U; + } else if (name=="spread_rate") { + return c_offset_V; + } else if (name=="T") { + return c_offset_T; + } else if (name=="lambda") { + return c_offset_L; + } else if (name == "eField") { + return c_offset_E; + } else if (name == "Uo") { + return c_offset_Uo; + } else { + for (size_t n=c_offset_Y; ntransportModel(); + + state["phase"]["name"] = m_thermo->name(); + AnyValue source = m_thermo->input().getMetadata("filename"); + state["phase"]["source"] = source.empty() ? "" : source.asString(); + + state["radiation-enabled"] = m_do_radiation; + if (m_do_radiation) { + state["emissivity-left"] = m_epsilon_left; + state["emissivity-right"] = m_epsilon_right; + } + + set energy_flags(m_do_energy.begin(), m_do_energy.end()); + if (energy_flags.size() == 1) { + state["energy-enabled"] = m_do_energy[0]; + } else { + state["energy-enabled"] = m_do_energy; + } + + state["Soret-enabled"] = m_do_soret; + + state["flux-gradient-basis"] = static_cast(m_fluxGradientBasis); + + set species_flags(m_do_species.begin(), m_do_species.end()); + if (species_flags.size() == 1) { + state["species-enabled"] = m_do_species[0]; + } else { + for (size_t k = 0; k < m_nsp; k++) { + state["species-enabled"][m_thermo->speciesName(k)] = m_do_species[k]; + } + } + + state["refine-criteria"]["ratio"] = m_refiner->maxRatio(); + state["refine-criteria"]["slope"] = m_refiner->maxDelta(); + state["refine-criteria"]["curve"] = m_refiner->maxSlope(); + state["refine-criteria"]["prune"] = m_refiner->prune(); + state["refine-criteria"]["grid-min"] = m_refiner->gridMin(); + state["refine-criteria"]["max-points"] = + static_cast(m_refiner->maxPoints()); + + if (m_zfixed != Undef) { + state["fixed-point"]["location"] = m_zfixed; + state["fixed-point"]["temperature"] = m_tfixed; + } + + // Two-point control meta data + if (m_twoPointControl) { + state["continuation-method"]["type"] = "two-point"; + state["continuation-method"]["left-location"] = m_zLeft; + state["continuation-method"]["right-location"] = m_zRight; + state["continuation-method"]["left-temperature"] = m_tLeft; + state["continuation-method"]["right-temperature"] = m_tRight; + } + + return state; +} + +shared_ptr Flow1D::asArray(const double* soln) const +{ + auto arr = SolutionArray::create( + m_solution, static_cast(nPoints()), getMeta()); + arr->addExtra("grid", false); // leading entry + AnyValue value; + value = m_z; + arr->setComponent("grid", value); + vector data(nPoints()); + for (size_t i = 0; i < nComponents(); i++) { + if (componentActive(i)) { + auto name = componentName(i); + for (size_t j = 0; j < nPoints(); j++) { + data[j] = soln[index(i, j)]; + } + if (!arr->hasComponent(name)) { + arr->addExtra(name, componentIndex(name) > c_offset_Y); + } + value = data; + arr->setComponent(name, value); + } + } + value = m_rho; + arr->setComponent("D", value); // use density rather than pressure + + if (m_do_radiation) { + arr->addExtra("radiative-heat-loss", true); // add at end + value = m_qdotRadiation; + arr->setComponent("radiative-heat-loss", value); + } + + return arr; +} + +void Flow1D::fromArray(SolutionArray& arr, double* soln) +{ + Domain1D::setMeta(arr.meta()); + arr.setLoc(0); + auto phase = arr.thermo(); + m_press = phase->pressure(); + + const auto grid = arr.getComponent("grid").as>(); + setupGrid(nPoints(), &grid[0]); + + for (size_t i = 0; i < nComponents(); i++) { + if (!componentActive(i)) { + continue; + } + string name = componentName(i); + if (arr.hasComponent(name)) { + const vector data = arr.getComponent(name).as>(); + for (size_t j = 0; j < nPoints(); j++) { + soln[index(i,j)] = data[j]; + } + } else { + warn_user("Flow1D::fromArray", "Saved state does not contain values for " + "component '{}' in domain '{}'.", name, id()); + } + } + + updateProperties(npos, soln + loc(), 0, m_points - 1); + setMeta(arr.meta()); +} + +void Flow1D::setMeta(const AnyMap& state) +{ + if (state.hasKey("energy-enabled")) { + const AnyValue& ee = state["energy-enabled"]; + if (ee.isScalar()) { + m_do_energy.assign(nPoints(), ee.asBool()); + } else { + m_do_energy = ee.asVector(nPoints()); + } + } + + setTransportModel(state.getString("transport-model", "mixture-averaged")); + + if (state.hasKey("Soret-enabled")) { + m_do_soret = state["Soret-enabled"].asBool(); + } + + if (state.hasKey("flux-gradient-basis")) { + m_fluxGradientBasis = static_cast( + state["flux-gradient-basis"].asInt()); + } + + if (state.hasKey("species-enabled")) { + const AnyValue& se = state["species-enabled"]; + if (se.isScalar()) { + m_do_species.assign(m_thermo->nSpecies(), se.asBool()); + } else { + m_do_species = se.asVector(m_thermo->nSpecies()); + } + } + + if (state.hasKey("radiation-enabled")) { + m_do_radiation = state["radiation-enabled"].asBool(); + if (m_do_radiation) { + m_epsilon_left = state["emissivity-left"].asDouble(); + m_epsilon_right = state["emissivity-right"].asDouble(); + } + } + + if (state.hasKey("refine-criteria")) { + const AnyMap& criteria = state["refine-criteria"].as(); + double ratio = criteria.getDouble("ratio", m_refiner->maxRatio()); + double slope = criteria.getDouble("slope", m_refiner->maxDelta()); + double curve = criteria.getDouble("curve", m_refiner->maxSlope()); + double prune = criteria.getDouble("prune", m_refiner->prune()); + m_refiner->setCriteria(ratio, slope, curve, prune); + + if (criteria.hasKey("grid-min")) { + m_refiner->setGridMin(criteria["grid-min"].asDouble()); + } + if (criteria.hasKey("max-points")) { + m_refiner->setMaxPoints(criteria["max-points"].asInt()); + } + } + + if (state.hasKey("fixed-point")) { + m_zfixed = state["fixed-point"]["location"].asDouble(); + m_tfixed = state["fixed-point"]["temperature"].asDouble(); + } + + // Two-point control meta data + if (state.hasKey("continuation-method")) { + const AnyMap& cm = state["continuation-method"].as(); + if (cm["type"] == "two-point") { + m_twoPointControl = true; + m_zLeft = cm["left-location"].asDouble(); + m_zRight = cm["right-location"].asDouble(); + m_tLeft = cm["left-temperature"].asDouble(); + m_tRight = cm["right-temperature"].asDouble(); + } else { + warn_user("Flow1D::setMeta", "Unknown continuation method '{}'.", + cm["type"].asString()); + } + } +} + +void Flow1D::solveEnergyEqn(size_t j) +{ + bool changed = false; + if (j == npos) { + for (size_t i = 0; i < m_points; i++) { + if (!m_do_energy[i]) { + changed = true; + } + m_do_energy[i] = true; + } + } else { + if (!m_do_energy[j]) { + changed = true; + } + m_do_energy[j] = true; + } + m_refiner->setActive(c_offset_U, true); + m_refiner->setActive(c_offset_V, true); + m_refiner->setActive(c_offset_T, true); + if (changed) { + needJacUpdate(); + } +} + +size_t Flow1D::getSolvingStage() const +{ + throw NotImplementedError("Flow1D::getSolvingStage", + "Not used by '{}' objects.", type()); +} + +void Flow1D::setSolvingStage(const size_t stage) +{ + throw NotImplementedError("Flow1D::setSolvingStage", + "Not used by '{}' objects.", type()); +} + +void Flow1D::solveElectricField(size_t j) +{ + throw NotImplementedError("Flow1D::solveElectricField", + "Not used by '{}' objects.", type()); +} + +void Flow1D::fixElectricField(size_t j) +{ + throw NotImplementedError("Flow1D::fixElectricField", + "Not used by '{}' objects.", type()); +} + +bool Flow1D::doElectricField(size_t j) const +{ + throw NotImplementedError("Flow1D::doElectricField", + "Not used by '{}' objects.", type()); +} + +void Flow1D::setBoundaryEmissivities(double e_left, double e_right) +{ + if (e_left < 0 || e_left > 1) { + throw CanteraError("Flow1D::setBoundaryEmissivities", + "The left boundary emissivity must be between 0.0 and 1.0!"); + } else if (e_right < 0 || e_right > 1) { + throw CanteraError("Flow1D::setBoundaryEmissivities", + "The right boundary emissivity must be between 0.0 and 1.0!"); + } else { + m_epsilon_left = e_left; + m_epsilon_right = e_right; + } +} + +void Flow1D::fixTemperature(size_t j) +{ + bool changed = false; + if (j == npos) { + for (size_t i = 0; i < m_points; i++) { + if (m_do_energy[i]) { + changed = true; + } + m_do_energy[i] = false; + } + } else { + if (m_do_energy[j]) { + changed = true; + } + m_do_energy[j] = false; + } + m_refiner->setActive(c_offset_U, false); + m_refiner->setActive(c_offset_V, false); + m_refiner->setActive(c_offset_T, false); + if (changed) { + needJacUpdate(); + } +} + +void Flow1D::grad_hk(const double* x, size_t j) +{ + for(size_t k = 0; k < m_nsp; k++) { + if (u(x, j) > 0.0) { + m_dhk_dz(k,j) = (m_hk(k,j) - m_hk(k,j-1)) / m_dz[j - 1]; + } + else { + m_dhk_dz(k,j) = (m_hk(k,j+1) - m_hk(k,j)) / m_dz[j]; + } + } +} + +// Two-point control functions +double Flow1D::leftControlPointTemperature() const +{ + if (m_twoPointControl) { + if (m_zLeft != Undef) { + return m_tLeft; + } else { + throw CanteraError("Flow1D::leftControlPointTemperature", + "Invalid operation: left control point location is not set"); + } + } else { + throw CanteraError("Flow1D::leftControlPointTemperature", + "Invalid operation: two-point control is not enabled."); + } +} + +double Flow1D::leftControlPointCoordinate() const +{ + if (m_twoPointControl) { + if (m_zLeft != Undef) { + return m_zLeft; + } else { + throw CanteraError("Flow1D::leftControlPointCoordinate", + "Invalid operation: left control point location is not set"); + } + } else { + throw CanteraError("Flow1D::leftControlPointCoordinate", + "Invalid operation: two-point control is not enabled."); + } +} + +void Flow1D::setLeftControlPointTemperature(double temperature) +{ + if (m_twoPointControl) { + if (m_zLeft != Undef) { + m_tLeft = temperature; + } else { + throw CanteraError("Flow1D::setLeftControlPointTemperature", + "Invalid operation: left control point location is not set"); + } + } else { + throw CanteraError("Flow1D::setLeftControlPointTemperature", + "Invalid operation: two-point control is not enabled."); + } +} + +void Flow1D::setLeftControlPointCoordinate(double z_left) +{ + if (m_twoPointControl) { + m_zLeft = z_left; + } else { + throw CanteraError("Flow1D::setLeftControlPointCoordinate", + "Invalid operation: two-point control is not enabled."); + } +} + +double Flow1D::rightControlPointTemperature() const +{ + if (m_twoPointControl) { + if (m_zRight != Undef) { + return m_tRight; + } else { + throw CanteraError("Flow1D::rightControlPointTemperature", + "Invalid operation: right control point location is not set"); + } + } else { + throw CanteraError("Flow1D::rightControlPointTemperature", + "Invalid operation: two-point control is not enabled."); + } +} + +double Flow1D::rightControlPointCoordinate() const +{ + if (m_twoPointControl) { + if (m_zRight != Undef) { + return m_zRight; + } else { + throw CanteraError("Flow1D::rightControlPointCoordinate", + "Invalid operation: right control point location is not set"); + } + } else { + throw CanteraError("Flow1D::rightControlPointCoordinate", + "Invalid operation: two-point control is not enabled."); + } +} + +void Flow1D::setRightControlPointTemperature(double temperature) +{ + if (m_twoPointControl) { + if (m_zRight != Undef) { + m_tRight = temperature; + } else { + throw CanteraError("Flow1D::setRightControlPointTemperature", + "Invalid operation: right control point location is not set"); + } + } else { + throw CanteraError("Flow1D::setRightControlPointTemperature", + "Invalid operation: two-point control is not enabled."); + } +} + +void Flow1D::setRightControlPointCoordinate(double z_right) +{ + if (m_twoPointControl) { + m_zRight = z_right; + } else { + throw CanteraError("Flow1D::setRightControlPointCoordinate", + "Invalid operation: two-point control is not enabled."); + } +} + +void Flow1D::enableTwoPointControl(bool twoPointControl) +{ + if (isStrained()) { + m_twoPointControl = twoPointControl; + } else { + throw CanteraError("Flow1D::enableTwoPointControl", + "Invalid operation: two-point control can only be used" + "with axisymmetric flames."); + } +} + +} // namespace diff --git a/src/oneD/IonFlow.cpp b/src/oneD/IonFlow.cpp index c1363266d14..9630bd6f19c 100644 --- a/src/oneD/IonFlow.cpp +++ b/src/oneD/IonFlow.cpp @@ -4,7 +4,7 @@ // at https://cantera.org/license.txt for license and copyright information. #include "cantera/oneD/IonFlow.h" -#include "cantera/oneD/StFlow.h" +#include "cantera/oneD/Flow1D.h" #include "cantera/oneD/refine.h" #include "cantera/transport/Transport.h" #include "cantera/numerics/funcs.h" @@ -16,7 +16,7 @@ namespace Cantera { IonFlow::IonFlow(ThermoPhase* ph, size_t nsp, size_t points) : - StFlow(ph, nsp, points) + Flow1D(ph, nsp, points) { // make a local copy of species charge for (size_t k = 0; k < m_nsp; k++) { @@ -82,7 +82,7 @@ string IonFlow::domainType() const { } void IonFlow::resize(size_t components, size_t points){ - StFlow::resize(components, points); + Flow1D::resize(components, points); m_mobility.resize(m_nsp*m_points); m_do_species.resize(m_nsp,true); m_do_electric_field.resize(m_points,false); @@ -93,13 +93,13 @@ bool IonFlow::componentActive(size_t n) const if (n == c_offset_E) { return true; } else { - return StFlow::componentActive(n); + return Flow1D::componentActive(n); } } void IonFlow::updateTransport(double* x, size_t j0, size_t j1) { - StFlow::updateTransport(x,j0,j1); + Flow1D::updateTransport(x,j0,j1); for (size_t j = j0; j < j1; j++) { setGasAtMidpoint(x,j); m_trans->getMobilities(&m_mobility[j*m_nsp]); @@ -202,7 +202,7 @@ void IonFlow::setSolvingStage(const size_t stage) void IonFlow::evalElectricField(double* x, double* rsd, int* diag, double rdt, size_t jmin, size_t jmax) { - StFlow::evalElectricField(x, rsd, diag, rdt, jmin, jmax); + Flow1D::evalElectricField(x, rsd, diag, rdt, jmin, jmax); if (m_stage != 2) { return; } @@ -227,7 +227,7 @@ void IonFlow::evalElectricField(double* x, double* rsd, int* diag, void IonFlow::evalSpecies(double* x, double* rsd, int* diag, double rdt, size_t jmin, size_t jmax) { - StFlow::evalSpecies(x, rsd, diag, rdt, jmin, jmax); + Flow1D::evalSpecies(x, rsd, diag, rdt, jmin, jmax); if (m_stage != 2) { return; } @@ -311,7 +311,7 @@ void IonFlow::setElectronTransport(vector& tfix, vector& diff_e, void IonFlow::_finalize(const double* x) { - StFlow::_finalize(x); + Flow1D::_finalize(x); bool p = m_do_electric_field[0]; if (p) { diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp index 555001dde59..3d03688afee 100644 --- a/src/oneD/Sim1D.cpp +++ b/src/oneD/Sim1D.cpp @@ -7,7 +7,7 @@ #include "cantera/oneD/Sim1D.h" #include "cantera/oneD/MultiJac.h" -#include "cantera/oneD/StFlow.h" +#include "cantera/oneD/Flow1D.h" #include "cantera/oneD/MultiNewton.h" #include "cantera/oneD/refine.h" #include "cantera/numerics/funcs.h" @@ -622,7 +622,7 @@ int Sim1D::setFixedTemperature(double t) size_t mfixed = npos; // loop over current grid to determine where new point is needed - StFlow* d_free = dynamic_cast(&domain(n)); + Flow1D* d_free = dynamic_cast(&domain(n)); size_t npnow = d.nPoints(); size_t nstart = znew.size(); if (d_free && d_free->isFree()) { @@ -702,7 +702,7 @@ double Sim1D::fixedTemperature() { double t_fixed = std::numeric_limits::quiet_NaN(); for (size_t n = 0; n < nDomains(); n++) { - StFlow* d = dynamic_cast(&domain(n)); + Flow1D* d = dynamic_cast(&domain(n)); if (d && d->isFree() && d->m_tfixed > 0) { t_fixed = d->m_tfixed; break; @@ -715,7 +715,7 @@ double Sim1D::fixedTemperatureLocation() { double z_fixed = std::numeric_limits::quiet_NaN(); for (size_t n = 0; n < nDomains(); n++) { - StFlow* d = dynamic_cast(&domain(n)); + Flow1D* d = dynamic_cast(&domain(n)); if (d && d->isFree() && d->m_tfixed > 0) { z_fixed = d->m_zfixed; break; @@ -735,7 +735,7 @@ void Sim1D::setLeftControlPoint(double temperature) continue; } - StFlow& d_axis = dynamic_cast(domain(n)); + Flow1D& d_axis = dynamic_cast(domain(n)); size_t np = d_axis.nPoints(); // Skip if two-point control is not enabled @@ -786,7 +786,7 @@ void Sim1D::setRightControlPoint(double temperature) continue; } - StFlow& d_axis = dynamic_cast(domain(n)); + Flow1D& d_axis = dynamic_cast(domain(n)); size_t np = d_axis.nPoints(); // Skip if two-point control is not enabled diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp index 0ab9cd50e44..f6631441aa5 100644 --- a/src/oneD/StFlow.cpp +++ b/src/oneD/StFlow.cpp @@ -3,12 +3,7 @@ // 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/base/SolutionArray.h" #include "cantera/oneD/StFlow.h" -#include "cantera/oneD/refine.h" -#include "cantera/transport/Transport.h" -#include "cantera/transport/TransportFactory.h" -#include "cantera/numerics/funcs.h" #include "cantera/base/global.h" using namespace std; @@ -17,1265 +12,303 @@ namespace Cantera { StFlow::StFlow(ThermoPhase* ph, size_t nsp, size_t points) : - Domain1D(nsp+c_offset_Y, points), - m_nsp(nsp) + Flow1D(ph, nsp, points) { - m_points = points; - - if (ph == 0) { - return; // used to create a dummy object - } - m_thermo = ph; - - size_t nsp2 = m_thermo->nSpecies(); - if (nsp2 != m_nsp) { - m_nsp = nsp2; - Domain1D::resize(m_nsp+c_offset_Y, points); - } - - // make a local copy of the species molecular weight vector - m_wt = m_thermo->molecularWeights(); - - // set pressure based on associated thermo object - setPressure(m_thermo->pressure()); - - // the species mass fractions are the last components in the solution - // vector, so the total number of components is the number of species - // plus the offset of the first mass fraction. - m_nv = c_offset_Y + m_nsp; - - // enable all species equations by default - m_do_species.resize(m_nsp, true); - - // but turn off the energy equation at all points - m_do_energy.resize(m_points,false); - - m_diff.resize(m_nsp*m_points); - m_multidiff.resize(m_nsp*m_nsp*m_points); - m_flux.resize(m_nsp,m_points); - m_wdot.resize(m_nsp,m_points, 0.0); - m_hk.resize(m_nsp, m_points, 0.0); - m_dhk_dz.resize(m_nsp, m_points - 1, 0.0); - m_ybar.resize(m_nsp); - m_qdotRadiation.resize(m_points, 0.0); - - //-------------- default solution bounds -------------------- - setBounds(c_offset_U, -1e20, 1e20); // no bounds on u - setBounds(c_offset_V, -1e20, 1e20); // no bounds on V - setBounds(c_offset_T, 200.0, 2*m_thermo->maxTemp()); // temperature bounds - setBounds(c_offset_L, -1e20, 1e20); // lambda should be negative - setBounds(c_offset_E, -1e20, 1e20); // no bounds on electric field - setBounds(c_offset_Uo, -1e20, 1e20); // no bounds on Uo - // mass fraction bounds - for (size_t k = 0; k < m_nsp; k++) { - setBounds(c_offset_Y+k, -1.0e-7, 1.0e5); - } - - //-------------------- grid refinement ------------------------- - m_refiner->setActive(c_offset_U, false); - m_refiner->setActive(c_offset_V, false); - m_refiner->setActive(c_offset_T, false); - m_refiner->setActive(c_offset_L, false); - m_refiner->setActive(c_offset_Uo, false); - - vector gr; - for (size_t ng = 0; ng < m_points; ng++) { - gr.push_back(1.0*ng/m_points); - } - setupGrid(m_points, gr.data()); - - // Find indices for radiating species - m_kRadiating.resize(2, npos); - m_kRadiating[0] = m_thermo->speciesIndex("CO2"); - m_kRadiating[1] = m_thermo->speciesIndex("H2O"); + warn_deprecated("StFlow::StFlow", + "To be removed after Cantera 3.1. Class replaced by Flow1D."); } StFlow::StFlow(shared_ptr th, size_t nsp, size_t points) - : StFlow(th.get(), nsp, points) + : Flow1D(th, nsp, points) { - auto sol = Solution::create(); - sol->setThermo(th); - setSolution(sol); + warn_deprecated("StFlow::StFlow", + "To be removed after Cantera 3.1. Class replaced by Flow1D."); } StFlow::StFlow(shared_ptr sol, const string& id, size_t points) - : StFlow(sol->thermo().get(), sol->thermo()->nSpecies(), points) + : Flow1D(sol, id, points) { - setSolution(sol); - m_id = id; - m_kin = m_solution->kinetics().get(); - m_trans = m_solution->transport().get(); - if (m_trans->transportModel() == "none") { - throw CanteraError("StFlow::StFlow", - "An appropriate transport model\nshould be set when instantiating the " - "Solution ('gas') object."); - } - m_solution->registerChangedCallback(this, [this]() { - setKinetics(m_solution->kinetics()); - setTransport(m_solution->transport()); - }); + warn_deprecated("StFlow::StFlow", + "To be removed after Cantera 3.1. Class replaced by Flow1D."); } -StFlow::~StFlow() +void StFlow::eval(size_t jg, double* xg, double* rg, integer* diagg, double rdt) { - if (m_solution) { - m_solution->removeChangedCallback(this); - } -} - -string StFlow::domainType() const { - if (m_isFree) { - return "free-flow"; - } - if (m_usesLambda) { - return "axisymmetric-flow"; - } - return "unstrained-flow"; -} - -void StFlow::setKinetics(shared_ptr kin) -{ - m_kin = kin.get(); - m_solution->setKinetics(kin); -} - -void StFlow::setTransport(shared_ptr trans) -{ - if (!trans) { - throw CanteraError("StFlow::setTransport", "Unable to set empty transport."); - } - m_trans = trans.get(); - if (m_trans->transportModel() == "none") { - throw CanteraError("StFlow::setTransport", "Invalid Transport model 'none'."); - } - m_do_multicomponent = (m_trans->transportModel() == "multicomponent" || - m_trans->transportModel() == "multicomponent-CK"); - - m_diff.resize(m_nsp * m_points); - if (m_do_multicomponent) { - m_multidiff.resize(m_nsp * m_nsp*m_points); - m_dthermal.resize(m_nsp, m_points, 0.0); - } - m_solution->setTransport(trans); -} - -void StFlow::resize(size_t ncomponents, size_t points) -{ - Domain1D::resize(ncomponents, points); - m_rho.resize(m_points, 0.0); - m_wtm.resize(m_points, 0.0); - m_cp.resize(m_points, 0.0); - m_visc.resize(m_points, 0.0); - m_tcon.resize(m_points, 0.0); - - m_diff.resize(m_nsp*m_points); - if (m_do_multicomponent) { - m_multidiff.resize(m_nsp*m_nsp*m_points); - m_dthermal.resize(m_nsp, m_points, 0.0); - } - m_flux.resize(m_nsp,m_points); - m_wdot.resize(m_nsp,m_points, 0.0); - m_hk.resize(m_nsp, m_points, 0.0); - m_dhk_dz.resize(m_nsp, m_points - 1, 0.0); - m_do_energy.resize(m_points,false); - m_qdotRadiation.resize(m_points, 0.0); - m_fixedtemp.resize(m_points); - - m_dz.resize(m_points-1); - m_z.resize(m_points); -} - -void StFlow::setupGrid(size_t n, const double* z) -{ - resize(m_nv, n); - - m_z[0] = z[0]; - for (size_t j = 1; j < m_points; j++) { - if (z[j] <= z[j-1]) { - throw CanteraError("StFlow::setupGrid", - "grid points must be monotonically increasing"); - } - m_z[j] = z[j]; - m_dz[j-1] = m_z[j] - m_z[j-1]; - } -} - -void StFlow::resetBadValues(double* xg) -{ - double* x = xg + loc(); - for (size_t j = 0; j < m_points; j++) { - double* Y = x + m_nv*j + c_offset_Y; - m_thermo->setMassFractions(Y); - m_thermo->getMassFractions(Y); - } -} - -void StFlow::setTransportModel(const string& trans) -{ - m_solution->setTransportModel(trans); -} - -string StFlow::transportModel() const { - return m_trans->transportModel(); -} - -void StFlow::setFluxGradientBasis(ThermoBasis fluxGradientBasis) { - m_fluxGradientBasis = fluxGradientBasis; - if (transportModel() != "mixture-averaged-CK" && - transportModel() != "mixture-averaged") { - warn_user("StFlow::setFluxGradientBasis", - "Setting fluxGradientBasis only affects " - "the mixture-averaged diffusion model."); - } -} - -void StFlow::_getInitialSoln(double* x) -{ - for (size_t j = 0; j < m_points; j++) { - T(x,j) = m_thermo->temperature(); - m_thermo->getMassFractions(&Y(x, 0, j)); - m_rho[j] = m_thermo->density(); - } -} - -void StFlow::setGas(const double* x, size_t j) -{ - m_thermo->setTemperature(T(x,j)); - const double* yy = x + m_nv*j + c_offset_Y; - m_thermo->setMassFractions_NoNorm(yy); - m_thermo->setPressure(m_press); -} - -void StFlow::setGasAtMidpoint(const double* x, size_t j) -{ - m_thermo->setTemperature(0.5*(T(x,j)+T(x,j+1))); - const double* yyj = x + m_nv*j + c_offset_Y; - const double* yyjp = x + m_nv*(j+1) + c_offset_Y; - for (size_t k = 0; k < m_nsp; k++) { - m_ybar[k] = 0.5*(yyj[k] + yyjp[k]); - } - m_thermo->setMassFractions_NoNorm(m_ybar.data()); - m_thermo->setPressure(m_press); -} - -void StFlow::_finalize(const double* x) -{ - if (!m_do_multicomponent && m_do_soret) { - throw CanteraError("StFlow::_finalize", - "Thermal diffusion (the Soret effect) is enabled, and requires " - "using a multicomponent transport model."); - } - - size_t nz = m_zfix.size(); - bool e = m_do_energy[0]; - for (size_t j = 0; j < m_points; j++) { - if (e || nz == 0) { - m_fixedtemp[j] = T(x, j); - } else { - double zz = (z(j) - z(0))/(z(m_points - 1) - z(0)); - double tt = linearInterp(zz, m_zfix, m_tfix); - m_fixedtemp[j] = tt; - } - } - if (e) { - solveEnergyEqn(); - } - - if (m_isFree) { - // If the domain contains the temperature fixed point, make sure that it - // is correctly set. This may be necessary when the grid has been modified - // externally. - if (m_tfixed != Undef) { - for (size_t j = 0; j < m_points; j++) { - if (z(j) == m_zfixed) { - return; // fixed point is already set correctly - } - } - - for (size_t j = 0; j < m_points - 1; j++) { - // Find where the temperature profile crosses the current - // fixed temperature. - if ((T(x, j) - m_tfixed) * (T(x, j+1) - m_tfixed) <= 0.0) { - m_tfixed = T(x, j+1); - m_zfixed = z(j+1); - return; - } - } - } - } -} - -void StFlow::eval(size_t jGlobal, double* xGlobal, double* rsdGlobal, - integer* diagGlobal, double rdt) -{ - // If evaluating a Jacobian, and the global point is outside the domain of + // if evaluating a Jacobian, and the global point is outside the domain of // influence for this domain, then skip evaluating the residual - if (jGlobal != npos && (jGlobal + 1 < firstPoint() || jGlobal > lastPoint() + 1)) { + if (jg != npos && (jg + 1 < firstPoint() || jg > lastPoint() + 1)) { return; } // start of local part of global arrays - double* x = xGlobal + loc(); - double* rsd = rsdGlobal + loc(); - integer* diag = diagGlobal + loc(); + double* x = xg + loc(); + double* rsd = rg + loc(); + integer* diag = diagg + loc(); size_t jmin, jmax; - if (jGlobal == npos) { // evaluate all points + if (jg == npos) { // evaluate all points jmin = 0; jmax = m_points - 1; } else { // evaluate points for Jacobian - size_t jpt = (jGlobal == 0) ? 0 : jGlobal - firstPoint(); + size_t jpt = (jg == 0) ? 0 : jg - firstPoint(); jmin = std::max(jpt, 1) - 1; jmax = std::min(jpt+1,m_points-1); } - updateProperties(jGlobal, x, jmin, jmax); - - if (m_do_radiation) { // Calculation of qdotRadiation - computeRadiation(x, jmin, jmax); - } - - evalContinuity(x, rsd, diag, rdt, jmin, jmax); - evalMomentum(x, rsd, diag, rdt, jmin, jmax); - evalEnergy(x, rsd, diag, rdt, jmin, jmax); - evalLambda(x, rsd, diag, rdt, jmin, jmax); - evalElectricField(x, rsd, diag, rdt, jmin, jmax); + updateProperties(jg, x, jmin, jmax); + evalResidual(x, rsd, diag, rdt, jmin, jmax); evalUo(x, rsd, diag, rdt, jmin, jmax); - evalSpecies(x, rsd, diag, rdt, jmin, jmax); -} - -void StFlow::updateProperties(size_t jg, double* x, size_t jmin, size_t jmax) -{ - // properties are computed for grid points from j0 to j1 - size_t j0 = std::max(jmin, 1) - 1; - size_t j1 = std::min(jmax+1,m_points-1); - - updateThermo(x, j0, j1); - if (jg == npos || m_force_full_update) { - // update transport properties only if a Jacobian is not being - // evaluated, or if specifically requested - updateTransport(x, j0, j1); - } - if (jg == npos) { - double* Yleft = x + index(c_offset_Y, jmin); - m_kExcessLeft = distance(Yleft, max_element(Yleft, Yleft + m_nsp)); - double* Yright = x + index(c_offset_Y, jmax); - m_kExcessRight = distance(Yright, max_element(Yright, Yright + m_nsp)); - } - - // update the species diffusive mass fluxes whether or not a - // Jacobian is being evaluated - updateDiffFluxes(x, j0, j1); } -void StFlow::computeRadiation(double* x, size_t jmin, size_t jmax) -{ - // Variable definitions for the Planck absorption coefficient and the - // radiation calculation: - double k_P_ref = 1.0*OneAtm; - - // Polynomial coefficients: - const double c_H2O[6] = {-0.23093, -1.12390, 9.41530, -2.99880, - 0.51382, -1.86840e-5}; - const double c_CO2[6] = {18.741, -121.310, 273.500, -194.050, - 56.310, -5.8169}; - - // Calculation of the two boundary values - double boundary_Rad_left = m_epsilon_left * StefanBoltz * pow(T(x, 0), 4); - double boundary_Rad_right = m_epsilon_right * StefanBoltz * pow(T(x, m_points - 1), 4); - - for (size_t j = jmin; j < jmax; j++) { - // calculation of the mean Planck absorption coefficient - double k_P = 0; - // Absorption coefficient for H2O - if (m_kRadiating[1] != npos) { - double k_P_H2O = 0; - for (size_t n = 0; n <= 5; n++) { - k_P_H2O += c_H2O[n] * pow(1000 / T(x, j), (double) n); - } - k_P_H2O /= k_P_ref; - k_P += m_press * X(x, m_kRadiating[1], j) * k_P_H2O; - } - // Absorption coefficient for CO2 - if (m_kRadiating[0] != npos) { - double k_P_CO2 = 0; - for (size_t n = 0; n <= 5; n++) { - k_P_CO2 += c_CO2[n] * pow(1000 / T(x, j), (double) n); - } - k_P_CO2 /= k_P_ref; - k_P += m_press * X(x, m_kRadiating[0], j) * k_P_CO2; - } - - // Calculation of the radiative heat loss term - double radiative_heat_loss = 2 * k_P *(2 * StefanBoltz * pow(T(x, j), 4) - - boundary_Rad_left - boundary_Rad_right); - - // set the radiative heat loss vector - m_qdotRadiation[j] = radiative_heat_loss; - } -} - -void StFlow::evalContinuity(double* x, double* rsd, int* diag, - double rdt, size_t jmin, size_t jmax) -{ - // The left boundary has the same form for all cases. - if (jmin == 0) { // left boundary - rsd[index(c_offset_U,jmin)] = -(rho_u(x,jmin + 1) - rho_u(x,jmin))/m_dz[jmin] - -(density(jmin + 1)*V(x,jmin + 1) - + density(jmin)*V(x,jmin)); - - diag[index(c_offset_U,jmin)] = 0; // Algebraic constraint - } - - if (jmax == m_points - 1) { // right boundary - if (m_usesLambda) { // axisymmetric flow - rsd[index(c_offset_U, jmax)] = rho_u(x, jmax); - } else { // right boundary (same for unstrained/free-flow) - rsd[index(c_offset_U, jmax)] = rho_u(x, jmax) - rho_u(x, jmax - 1); - } - diag[index(c_offset_U, jmax)] = 0; // Algebraic constraint - } - - // j0 and j1 are constrained to only interior points - size_t j0 = std::max(jmin, 1); - size_t j1 = std::min(jmax, m_points - 2); - if (m_usesLambda) { // "axisymmetric-flow" - for (size_t j = j0; j <= j1; j++) { // interior points - // For "axisymmetric-flow", the continuity equation propagates the - // mass flow rate information to the left (j+1 -> j) from the value - // specified at the right boundary. The lambda information propagates - // in the opposite direction. - rsd[index(c_offset_U,j)] = -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j] - -(density(j+1)*V(x,j+1) + density(j)*V(x,j)); - - diag[index(c_offset_U, j)] = 0; // Algebraic constraint - } - } else if (m_isFree) { // "free-flow" - for (size_t j = j0; j <= j1; j++) { - // terms involving V are zero as V=0 by definition - if (grid(j) > m_zfixed) { - rsd[index(c_offset_U,j)] = -(rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1]; - } else if (grid(j) == m_zfixed) { - if (m_do_energy[j]) { - rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed); - } else { - rsd[index(c_offset_U,j)] = (rho_u(x,j) - m_rho[0]*0.3); // why 0.3? - } - } else { // grid(j < m_zfixed - rsd[index(c_offset_U,j)] = -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j]; - } - diag[index(c_offset_U, j)] = 0; // Algebraic constraint - } - } else { // "unstrained-flow" (fixed mass flow rate) - for (size_t j = j0; j <= j1; j++) { - rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j - 1); - diag[index(c_offset_U, j)] = 0; // Algebraic constraint - } - } -} - -void StFlow::evalMomentum(double* x, double* rsd, int* diag, +void StFlow::evalResidual(double* x, double* rsd, int* diag, double rdt, size_t jmin, size_t jmax) { - if (!m_usesLambda) { //disable this equation - for (size_t j = jmin; j <= jmax; j++) { - rsd[index(c_offset_V, j)] = V(x, j); - diag[index(c_offset_V, j)] = 0; - } - return; - } - - if (jmin == 0) { // left boundary - rsd[index(c_offset_V,jmin)] = V(x,jmin); - } - - if (jmax == m_points - 1) { // right boundary - rsd[index(c_offset_V,jmax)] = V(x,jmax); - diag[index(c_offset_V,jmax)] = 0; - } - - // j0 and j1 are constrained to only interior points - size_t j0 = std::max(jmin, 1); - size_t j1 = std::min(jmax, m_points - 2); - for (size_t j = j0; j <= j1; j++) { // interior points - rsd[index(c_offset_V,j)] = (shear(x, j) - lambda(x, j) - - rho_u(x, j) * dVdz(x, j) - - m_rho[j] * V(x, j) * V(x, j)) / m_rho[j] - - rdt * (V(x, j) - V_prev(j)); - diag[index(c_offset_V, j)] = 1; - } -} + //---------------------------------------------------- + // evaluate the residual equations at all required + // grid points + //---------------------------------------------------- -void StFlow::evalLambda(double* x, double* rsd, int* diag, - double rdt, size_t jmin, size_t jmax) -{ - if (!m_usesLambda) { // disable this equation - for (size_t j = jmin; j <= jmax; j++) { - rsd[index(c_offset_L, j)] = lambda(x, j); - diag[index(c_offset_L, j)] = 0; - } - return; - } - - if (jmin == 0) { // left boundary - if (m_twoPointControl) { - rsd[index(c_offset_L, jmin)] = lambda(x,jmin+1) - lambda(x,jmin); - } else { - rsd[index(c_offset_L, jmin)] = -rho_u(x, jmin); - } - } - - if (jmax == m_points - 1) { // right boundary - rsd[index(c_offset_L, jmax)] = lambda(x, jmax) - lambda(x, jmax-1); - diag[index(c_offset_L, jmax)] = 0; - } - - // j0 and j1 are constrained to only interior points - size_t j0 = std::max(jmin, 1); - size_t j1 = std::min(jmax, m_points - 2); - for (size_t j = j0; j <= j1; j++) { // interior points - if (m_twoPointControl) { - if (grid(j) == m_zLeft) { - rsd[index(c_offset_L, j)] = T(x,j) - m_tLeft; - } else if (grid(j) > m_zLeft) { - rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1); - } else if (grid(j) < m_zLeft) { - rsd[index(c_offset_L, j)] = lambda(x,j+1) - lambda(x,j); - } - } else { - rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1); - } - diag[index(c_offset_L, j)] = 0; - } -} - -void StFlow::evalEnergy(double* x, double* rsd, int* diag, - double rdt, size_t jmin, size_t jmax) -{ - if (jmin == 0) { // left boundary - rsd[index(c_offset_T,jmin)] = T(x,jmin); - } - - if (jmax == m_points - 1) { // right boundary - rsd[index(c_offset_T, jmax)] = T(x,jmax); - } - - // j0 and j1 are constrained to only interior points - size_t j0 = std::max(jmin, 1); - size_t j1 = std::min(jmax, m_points - 2); - for (size_t j = j0; j <= j1; j++) { - if (m_do_energy[j]) { - grad_hk(x, j); - double sum = 0.0; - for (size_t k = 0; k < m_nsp; k++) { - double flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j)); - sum += wdot(k,j)*m_hk(k,j); - sum += flxk * m_dhk_dz(k,j) / m_wt[k]; + // calculation of qdotRadiation (see docstring of enableRadiation) + if (m_do_radiation) { + // variable definitions for the Planck absorption coefficient and the + // radiation calculation: + double k_P_ref = 1.0*OneAtm; + + // polynomial coefficients: + const double c_H2O[6] = {-0.23093, -1.12390, 9.41530, -2.99880, + 0.51382, -1.86840e-5}; + const double c_CO2[6] = {18.741, -121.310, 273.500, -194.050, + 56.310, -5.8169}; + + // calculation of the two boundary values + double boundary_Rad_left = m_epsilon_left * StefanBoltz * pow(T(x, 0), 4); + double boundary_Rad_right = m_epsilon_right * StefanBoltz * pow(T(x, m_points - 1), 4); + + // loop over all grid points + for (size_t j = jmin; j < jmax; j++) { + // helping variable for the calculation + double radiative_heat_loss = 0; + + // calculation of the mean Planck absorption coefficient + double k_P = 0; + // absorption coefficient for H2O + if (m_kRadiating[1] != npos) { + double k_P_H2O = 0; + for (size_t n = 0; n <= 5; n++) { + k_P_H2O += c_H2O[n] * pow(1000 / T(x, j), (double) n); + } + k_P_H2O /= k_P_ref; + k_P += m_press * X(x, m_kRadiating[1], j) * k_P_H2O; } - - rsd[index(c_offset_T, j)] = - m_cp[j]*rho_u(x,j)*dTdz(x,j) - - divHeatFlux(x,j) - sum; - rsd[index(c_offset_T, j)] /= (m_rho[j]*m_cp[j]); - rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j)); - rsd[index(c_offset_T, j)] -= (m_qdotRadiation[j] / (m_rho[j] * m_cp[j])); - diag[index(c_offset_T, j)] = 1; - } else { - // residual equations if the energy equation is disabled - rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j); - diag[index(c_offset_T, j)] = 0; - } - } -} - -void StFlow::evalUo(double* x, double* rsd, int* diag, - double rdt, size_t jmin, size_t jmax) -{ - if (!m_twoPointControl) { // disable this equation - for (size_t j = jmin; j <= jmax; j++) { - rsd[index(c_offset_Uo, j)] = Uo(x, j); - diag[index(c_offset_Uo, j)] = 0; - } - return; - } - - if (jmin == 0) { // left boundary - rsd[index(c_offset_Uo, jmin)] = Uo(x,jmin+1) - Uo(x,jmin); - diag[index(c_offset_Uo, jmin)] = 0; - } - - if (jmax == m_points - 1) { // right boundary - if(m_twoPointControl) { - rsd[index(c_offset_Uo, jmax)] = Uo(x,jmax) - Uo(x,jmax-1); - } - diag[index(c_offset_Uo, jmax)] = 0; - } - - // j0 and j1 are constrained to only interior points - size_t j0 = std::max(jmin, 1); - size_t j1 = std::min(jmax, m_points - 2); - for (size_t j = j0; j <= j1; j++) { // interior points - if (m_twoPointControl) { - if (grid(j) == m_zRight) { - rsd[index(c_offset_Uo, j)] = T(x,j) - m_tRight; - } else if (grid(j) > m_zRight) { - rsd[index(c_offset_Uo, j)] = Uo(x,j) - Uo(x,j-1); - } else if (grid(j) < m_zRight) { - rsd[index(c_offset_Uo, j)] = Uo(x,j+1) - Uo(x,j); + // absorption coefficient for CO2 + if (m_kRadiating[0] != npos) { + double k_P_CO2 = 0; + for (size_t n = 0; n <= 5; n++) { + k_P_CO2 += c_CO2[n] * pow(1000 / T(x, j), (double) n); + } + k_P_CO2 /= k_P_ref; + k_P += m_press * X(x, m_kRadiating[0], j) * k_P_CO2; } - } else { - rsd[index(c_offset_Uo, j)] = Uo(x,j+1) - Uo(x,j); - } - diag[index(c_offset_Uo, j)] = 0; - } -} -void StFlow::evalSpecies(double* x, double* rsd, int* diag, - double rdt, size_t jmin, size_t jmax) -{ - if (jmin == 0) { // left boundary - double sum = 0.0; - for (size_t k = 0; k < m_nsp; k++) { - sum += Y(x,k,jmin); - rsd[index(c_offset_Y + k, jmin)] = -(m_flux(k,jmin) + - rho_u(x,jmin) * Y(x,k,jmin)); - } - rsd[index(c_offset_Y + leftExcessSpecies(), jmin)] = 1.0 - sum; - } + // calculation of the radiative heat loss term + radiative_heat_loss = 2 * k_P *(2 * StefanBoltz * pow(T(x, j), 4) + - boundary_Rad_left - boundary_Rad_right); - if (jmax == m_points - 1) { // right boundary - double sum = 0.0; - for (size_t k = 0; k < m_nsp; k++) { - sum += Y(x,k,jmax); - rsd[index(k+c_offset_Y,jmax)] = m_flux(k,jmax-1) + - rho_u(x,jmax)*Y(x,k,jmax); + // set the radiative heat loss vector + m_qdotRadiation[j] = radiative_heat_loss; } - rsd[index(c_offset_Y + rightExcessSpecies(), jmax)] = 1.0 - sum; - diag[index(c_offset_Y + rightExcessSpecies(), jmax)] = 0; } - // j0 and j1 are constrained to only interior points - size_t j0 = std::max(jmin, 1); - size_t j1 = std::min(jmax, m_points - 2); - for (size_t j = j0; j <= j1; j++) { - for (size_t k = 0; k < m_nsp; k++) { - double convec = rho_u(x,j)*dYdz(x,k,j); - double diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1)) / (z(j+1) - z(j-1)); - rsd[index(c_offset_Y + k, j)] = (m_wt[k]*(wdot(k,j)) - - convec - diffus)/m_rho[j] - - rdt*(Y(x,k,j) - Y_prev(k,j)); - diag[index(c_offset_Y + k, j)] = 1; - } - } -} - -void StFlow::evalElectricField(double* x, double* rsd, int* diag, - double rdt, size_t jmin, size_t jmax) -{ for (size_t j = jmin; j <= jmax; j++) { - // The same value is used for left/right/interior points - rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)]; - } -} - -void StFlow::updateTransport(double* x, size_t j0, size_t j1) -{ - if (m_do_multicomponent) { - for (size_t j = j0; j < j1; j++) { - setGasAtMidpoint(x,j); - double wtm = m_thermo->meanMolecularWeight(); - double rho = m_thermo->density(); - m_visc[j] = (m_dovisc ? m_trans->viscosity() : 0.0); - m_trans->getMultiDiffCoeffs(m_nsp, &m_multidiff[mindex(0,0,j)]); - - // Use m_diff as storage for the factor outside the summation - for (size_t k = 0; k < m_nsp; k++) { - m_diff[k+j*m_nsp] = m_wt[k] * rho / (wtm*wtm); - } - - m_tcon[j] = m_trans->thermalConductivity(); - if (m_do_soret) { - m_trans->getThermalDiffCoeffs(m_dthermal.ptrColumn(0) + j*m_nsp); - } - } - } else { // mixture averaged transport - for (size_t j = j0; j < j1; j++) { - setGasAtMidpoint(x,j); - m_visc[j] = (m_dovisc ? m_trans->viscosity() : 0.0); - - if (m_fluxGradientBasis == ThermoBasis::molar) { - m_trans->getMixDiffCoeffs(&m_diff[j*m_nsp]); + //---------------------------------------------- + // left boundary + //---------------------------------------------- + + if (j == 0) { + // these may be modified by a boundary object + + // Continuity. This propagates information right-to-left, since + // rho_u at point 0 is dependent on rho_u at point 1, but not on + // mdot from the inlet. + rsd[index(c_offset_U,0)] = + -(rho_u(x,1) - rho_u(x,0))/m_dz[0] + -(density(1)*V(x,1) + density(0)*V(x,0)); + + // the inlet (or other) object connected to this one will modify + // these equations by subtracting its values for V, T, and mdot. As + // a result, these residual equations will force the solution + // variables to the values for the boundary object + rsd[index(c_offset_V,0)] = V(x,0); + rsd[index(c_offset_T,0)] = T(x,0); + if (m_usesLambda) { + rsd[index(c_offset_L, 0)] = -rho_u(x, 0); } else { - m_trans->getMixDiffCoeffsMass(&m_diff[j*m_nsp]); - } - - double rho = m_thermo->density(); - - if (m_fluxGradientBasis == ThermoBasis::molar) { - double wtm = m_thermo->meanMolecularWeight(); - for (size_t k=0; k < m_nsp; k++) { - m_diff[k+j*m_nsp] *= m_wt[k] * rho / wtm; - } - } else { - for (size_t k=0; k < m_nsp; k++) { - m_diff[k+j*m_nsp] *= rho; - } + rsd[index(c_offset_L, 0)] = lambda(x, 0); + diag[index(c_offset_L, 0)] = 0; } - m_tcon[j] = m_trans->thermalConductivity(); - } - } -} - -void StFlow::show(const double* x) -{ - writelog(" Pressure: {:10.4g} Pa\n", m_press); - - Domain1D::show(x); - if (m_do_radiation) { - writeline('-', 79, false, true); - writelog("\n z radiative heat loss"); - writeline('-', 79, false, true); - for (size_t j = 0; j < m_points; j++) { - writelog("\n {:10.4g} {:10.4g}", m_z[j], m_qdotRadiation[j]); - } - writelog("\n"); - } -} - -void StFlow::updateDiffFluxes(const double* x, size_t j0, size_t j1) -{ - if (m_do_multicomponent) { - for (size_t j = j0; j < j1; j++) { - double dz = z(j+1) - z(j); + // The default boundary condition for species is zero flux. However, + // the boundary object may modify this. + double sum = 0.0; for (size_t k = 0; k < m_nsp; k++) { - double sum = 0.0; - for (size_t m = 0; m < m_nsp; m++) { - sum += m_wt[m] * m_multidiff[mindex(k,m,j)] * (X(x,m,j+1)-X(x,m,j)); - } - m_flux(k,j) = sum * m_diff[k+j*m_nsp] / dz; + sum += Y(x,k,0); + rsd[index(c_offset_Y + k, 0)] = + -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0)); } - } - } else { - for (size_t j = j0; j < j1; j++) { - double sum = 0.0; - double dz = z(j+1) - z(j); - if (m_fluxGradientBasis == ThermoBasis::molar) { - for (size_t k = 0; k < m_nsp; k++) { - m_flux(k,j) = m_diff[k+m_nsp*j] * (X(x,k,j) - X(x,k,j+1))/dz; - sum -= m_flux(k,j); - } + rsd[index(c_offset_Y + leftExcessSpecies(), 0)] = 1.0 - sum; + + // set residual of poisson's equ to zero + rsd[index(c_offset_E, 0)] = x[index(c_offset_E, j)]; + } else if (j == m_points - 1) { + evalRightBoundary(x, rsd, diag, rdt); + } else { // interior points + evalContinuity(j, x, rsd, diag, rdt); + // set residual of poisson's equ to zero + rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)]; + + //------------------------------------------------ + // Radial momentum equation + // + // \rho dV/dt + \rho u dV/dz + \rho V^2 + // = d(\mu dV/dz)/dz - lambda + //------------------------------------------------- + if (m_usesLambda) { + rsd[index(c_offset_V,j)] = + (shear(x, j) - lambda(x, j) - rho_u(x, j) * dVdz(x, j) + - m_rho[j] * V(x, j) * V(x, j)) / m_rho[j] + - rdt * (V(x, j) - V_prev(j)); + diag[index(c_offset_V, j)] = 1; } else { - for (size_t k = 0; k < m_nsp; k++) { - m_flux(k,j) = m_diff[k+m_nsp*j] * (Y(x,k,j) - Y(x,k,j+1))/dz; - sum -= m_flux(k,j); - } + rsd[index(c_offset_V, j)] = V(x, j); + diag[index(c_offset_V, j)] = 0; } - // correction flux to insure that \sum_k Y_k V_k = 0. - for (size_t k = 0; k < m_nsp; k++) { - m_flux(k,j) += sum*Y(x,k,j); - } - } - } - if (m_do_soret) { - for (size_t m = j0; m < j1; m++) { - double gradlogT = 2.0 * (T(x,m+1) - T(x,m)) / - ((T(x,m+1) + T(x,m)) * (z(m+1) - z(m))); + //------------------------------------------------- + // Species equations + // + // \rho dY_k/dt + \rho u dY_k/dz + dJ_k/dz + // = M_k\omega_k + //------------------------------------------------- + getWdot(x,j); for (size_t k = 0; k < m_nsp; k++) { - m_flux(k,m) -= m_dthermal(k,m)*gradlogT; + double convec = rho_u(x,j)*dYdz(x,k,j); + double diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1)) + / (z(j+1) - z(j-1)); + rsd[index(c_offset_Y + k, j)] + = (m_wt[k]*(wdot(k,j)) + - convec - diffus)/m_rho[j] + - rdt*(Y(x,k,j) - Y_prev(k,j)); + diag[index(c_offset_Y + k, j)] = 1; } - } - } -} -string StFlow::componentName(size_t n) const -{ - switch (n) { - case c_offset_U: - return "velocity"; - case c_offset_V: - return "spread_rate"; - case c_offset_T: - return "T"; - case c_offset_L: - return "lambda"; - case c_offset_E: - return "eField"; - case c_offset_Uo: - return "Uo"; - default: - if (n >= c_offset_Y && n < (c_offset_Y + m_nsp)) { - return m_thermo->speciesName(n - c_offset_Y); - } else { - return ""; - } - } -} - -size_t StFlow::componentIndex(const string& name) const -{ - if (name=="velocity") { - return c_offset_U; - } else if (name=="spread_rate") { - return c_offset_V; - } else if (name=="T") { - return c_offset_T; - } else if (name=="lambda") { - return c_offset_L; - } else if (name == "eField") { - return c_offset_E; - } else if (name == "Uo") { - return c_offset_Uo; - } else { - for (size_t n=c_offset_Y; ntransportModel(); - - state["phase"]["name"] = m_thermo->name(); - AnyValue source = m_thermo->input().getMetadata("filename"); - state["phase"]["source"] = source.empty() ? "" : source.asString(); - - state["radiation-enabled"] = m_do_radiation; - if (m_do_radiation) { - state["emissivity-left"] = m_epsilon_left; - state["emissivity-right"] = m_epsilon_right; - } - - set energy_flags(m_do_energy.begin(), m_do_energy.end()); - if (energy_flags.size() == 1) { - state["energy-enabled"] = m_do_energy[0]; - } else { - state["energy-enabled"] = m_do_energy; - } - - state["Soret-enabled"] = m_do_soret; - - state["flux-gradient-basis"] = static_cast(m_fluxGradientBasis); - - set species_flags(m_do_species.begin(), m_do_species.end()); - if (species_flags.size() == 1) { - state["species-enabled"] = m_do_species[0]; - } else { - for (size_t k = 0; k < m_nsp; k++) { - state["species-enabled"][m_thermo->speciesName(k)] = m_do_species[k]; - } - } - - state["refine-criteria"]["ratio"] = m_refiner->maxRatio(); - state["refine-criteria"]["slope"] = m_refiner->maxDelta(); - state["refine-criteria"]["curve"] = m_refiner->maxSlope(); - state["refine-criteria"]["prune"] = m_refiner->prune(); - state["refine-criteria"]["grid-min"] = m_refiner->gridMin(); - state["refine-criteria"]["max-points"] = - static_cast(m_refiner->maxPoints()); - - if (m_zfixed != Undef) { - state["fixed-point"]["location"] = m_zfixed; - state["fixed-point"]["temperature"] = m_tfixed; - } - - // Two-point control meta data - if (m_twoPointControl) { - state["continuation-method"]["type"] = "two-point"; - state["continuation-method"]["left-location"] = m_zLeft; - state["continuation-method"]["right-location"] = m_zRight; - state["continuation-method"]["left-temperature"] = m_tLeft; - state["continuation-method"]["right-temperature"] = m_tRight; - } + //----------------------------------------------- + // energy equation + // + // \rho c_p dT/dt + \rho c_p u dT/dz + // = d(k dT/dz)/dz + // - sum_k(\omega_k h_k_ref) + // - sum_k(J_k c_p_k / M_k) dT/dz + //----------------------------------------------- + if (m_do_energy[j]) { + + setGas(x,j); + double dtdzj = dTdz(x,j); + double sum = 0.0; - return state; -} + grad_hk(x, j); + for (size_t k = 0; k < m_nsp; k++) { + double flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j)); + sum += wdot(k,j)*m_hk(k,j); + sum += flxk * m_dhk_dz(k,j) / m_wt[k]; + } -shared_ptr StFlow::asArray(const double* soln) const -{ - auto arr = SolutionArray::create( - m_solution, static_cast(nPoints()), getMeta()); - arr->addExtra("grid", false); // leading entry - AnyValue value; - value = m_z; - arr->setComponent("grid", value); - vector data(nPoints()); - for (size_t i = 0; i < nComponents(); i++) { - if (componentActive(i)) { - auto name = componentName(i); - for (size_t j = 0; j < nPoints(); j++) { - data[j] = soln[index(i, j)]; - } - if (!arr->hasComponent(name)) { - arr->addExtra(name, componentIndex(name) > c_offset_Y); + rsd[index(c_offset_T, j)] = - m_cp[j]*rho_u(x,j)*dtdzj + - divHeatFlux(x,j) - sum; + rsd[index(c_offset_T, j)] /= (m_rho[j]*m_cp[j]); + rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j)); + rsd[index(c_offset_T, j)] -= (m_qdotRadiation[j] / (m_rho[j] * m_cp[j])); + diag[index(c_offset_T, j)] = 1; + } else { + // residual equations if the energy equation is disabled + rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j); + diag[index(c_offset_T, j)] = 0; } - value = data; - arr->setComponent(name, value); - } - } - value = m_rho; - arr->setComponent("D", value); // use density rather than pressure - - if (m_do_radiation) { - arr->addExtra("radiative-heat-loss", true); // add at end - value = m_qdotRadiation; - arr->setComponent("radiative-heat-loss", value); - } - return arr; -} - -void StFlow::fromArray(SolutionArray& arr, double* soln) -{ - Domain1D::setMeta(arr.meta()); - arr.setLoc(0); - auto phase = arr.thermo(); - m_press = phase->pressure(); - - const auto grid = arr.getComponent("grid").as>(); - setupGrid(nPoints(), &grid[0]); - - for (size_t i = 0; i < nComponents(); i++) { - if (!componentActive(i)) { - continue; - } - string name = componentName(i); - if (arr.hasComponent(name)) { - const vector data = arr.getComponent(name).as>(); - for (size_t j = 0; j < nPoints(); j++) { - soln[index(i,j)] = data[j]; + if (m_usesLambda) { + rsd[index(c_offset_L, j)] = lambda(x, j) - lambda(x, j - 1); + } else { + rsd[index(c_offset_L, j)] = lambda(x, j); } - } else { - warn_user("StFlow::fromArray", "Saved state does not contain values for " - "component '{}' in domain '{}'.", name, id()); + diag[index(c_offset_L, j)] = 0; } } - - updateProperties(npos, soln + loc(), 0, m_points - 1); - setMeta(arr.meta()); } -void StFlow::setMeta(const AnyMap& state) +void StFlow::evalRightBoundary(double* x, double* rsd, int* diag, double rdt) { - if (state.hasKey("energy-enabled")) { - const AnyValue& ee = state["energy-enabled"]; - if (ee.isScalar()) { - m_do_energy.assign(nPoints(), ee.asBool()); - } else { - m_do_energy = ee.asVector(nPoints()); - } - } - - setTransportModel(state.getString("transport-model", "mixture-averaged")); + size_t j = m_points - 1; - if (state.hasKey("Soret-enabled")) { - m_do_soret = state["Soret-enabled"].asBool(); - } - - if (state.hasKey("flux-gradient-basis")) { - m_fluxGradientBasis = static_cast( - state["flux-gradient-basis"].asInt()); - } - - if (state.hasKey("species-enabled")) { - const AnyValue& se = state["species-enabled"]; - if (se.isScalar()) { - m_do_species.assign(m_thermo->nSpecies(), se.asBool()); - } else { - m_do_species = se.asVector(m_thermo->nSpecies()); - } - } + // the boundary object connected to the right of this one may modify or + // replace these equations. The default boundary conditions are zero u, V, + // and T, and zero diffusive flux for all species. - if (state.hasKey("radiation-enabled")) { - m_do_radiation = state["radiation-enabled"].asBool(); - if (m_do_radiation) { - m_epsilon_left = state["emissivity-left"].asDouble(); - m_epsilon_right = state["emissivity-right"].asDouble(); - } - } - - if (state.hasKey("refine-criteria")) { - const AnyMap& criteria = state["refine-criteria"].as(); - double ratio = criteria.getDouble("ratio", m_refiner->maxRatio()); - double slope = criteria.getDouble("slope", m_refiner->maxDelta()); - double curve = criteria.getDouble("curve", m_refiner->maxSlope()); - double prune = criteria.getDouble("prune", m_refiner->prune()); - m_refiner->setCriteria(ratio, slope, curve, prune); - - if (criteria.hasKey("grid-min")) { - m_refiner->setGridMin(criteria["grid-min"].asDouble()); - } - if (criteria.hasKey("max-points")) { - m_refiner->setMaxPoints(criteria["max-points"].asInt()); - } - } - - if (state.hasKey("fixed-point")) { - m_zfixed = state["fixed-point"]["location"].asDouble(); - m_tfixed = state["fixed-point"]["temperature"].asDouble(); - } - - // Two-point control meta data - if (state.hasKey("continuation-method")) { - const AnyMap& cm = state["continuation-method"].as(); - if (cm["type"] == "two-point") { - m_twoPointControl = true; - m_zLeft = cm["left-location"].asDouble(); - m_zRight = cm["right-location"].asDouble(); - m_tLeft = cm["left-temperature"].asDouble(); - m_tRight = cm["right-temperature"].asDouble(); - } else { - warn_user("StFlow::setMeta", "Unknown continuation method '{}'.", - cm["type"].asString()); - } + rsd[index(c_offset_V,j)] = V(x,j); + diag[index(c_offset_V,j)] = 0; + double sum = 0.0; + // set residual of poisson's equ to zero + rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)]; + for (size_t k = 0; k < m_nsp; k++) { + sum += Y(x,k,j); + rsd[index(k+c_offset_Y,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j); } -} - -void StFlow::solveEnergyEqn(size_t j) -{ - bool changed = false; - if (j == npos) { - for (size_t i = 0; i < m_points; i++) { - if (!m_do_energy[i]) { - changed = true; - } - m_do_energy[i] = true; - } + rsd[index(c_offset_Y + rightExcessSpecies(), j)] = 1.0 - sum; + diag[index(c_offset_Y + rightExcessSpecies(), j)] = 0; + if (m_usesLambda) { + rsd[index(c_offset_U, j)] = rho_u(x, j); } else { - if (!m_do_energy[j]) { - changed = true; - } - m_do_energy[j] = true; + rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j-1); } - m_refiner->setActive(c_offset_U, true); - m_refiner->setActive(c_offset_V, true); - m_refiner->setActive(c_offset_T, true); - if (changed) { - needJacUpdate(); - } -} - -size_t StFlow::getSolvingStage() const -{ - throw NotImplementedError("StFlow::getSolvingStage", - "Not used by '{}' objects.", type()); -} - -void StFlow::setSolvingStage(const size_t stage) -{ - throw NotImplementedError("StFlow::setSolvingStage", - "Not used by '{}' objects.", type()); -} -void StFlow::solveElectricField(size_t j) -{ - throw NotImplementedError("StFlow::solveElectricField", - "Not used by '{}' objects.", type()); + rsd[index(c_offset_L, j)] = lambda(x, j) - lambda(x, j-1); + diag[index(c_offset_L, j)] = 0; + rsd[index(c_offset_T, j)] = T(x, j); } -void StFlow::fixElectricField(size_t j) +void StFlow::evalContinuity(size_t j, double* x, double* rsd, int* diag, double rdt) { - throw NotImplementedError("StFlow::fixElectricField", - "Not used by '{}' objects.", type()); -} - -bool StFlow::doElectricField(size_t j) const -{ - throw NotImplementedError("StFlow::doElectricField", - "Not used by '{}' objects.", type()); -} - -void StFlow::setBoundaryEmissivities(double e_left, double e_right) -{ - if (e_left < 0 || e_left > 1) { - throw CanteraError("StFlow::setBoundaryEmissivities", - "The left boundary emissivity must be between 0.0 and 1.0!"); - } else if (e_right < 0 || e_right > 1) { - throw CanteraError("StFlow::setBoundaryEmissivities", - "The right boundary emissivity must be between 0.0 and 1.0!"); - } else { - m_epsilon_left = e_left; - m_epsilon_right = e_right; - } -} - -void StFlow::fixTemperature(size_t j) -{ - bool changed = false; - if (j == npos) { - for (size_t i = 0; i < m_points; i++) { - if (m_do_energy[i]) { - changed = true; + //algebraic constraint + diag[index(c_offset_U, j)] = 0; + //---------------------------------------------- + // Continuity equation + // + // d(\rho u)/dz + 2\rho V = 0 + //---------------------------------------------- + if (m_usesLambda) { + // Note that this propagates the mass flow rate information to the left + // (j+1 -> j) from the value specified at the right boundary. The + // lambda information propagates in the opposite direction. + rsd[index(c_offset_U,j)] = + -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j] + -(density(j+1)*V(x,j+1) + density(j)*V(x,j)); + } else if (m_isFree) { + // terms involving V are zero as V=0 by definition + if (grid(j) > m_zfixed) { + rsd[index(c_offset_U,j)] = + - (rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1]; + } else if (grid(j) == m_zfixed) { + if (m_do_energy[j]) { + rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed); + } else { + rsd[index(c_offset_U,j)] = (rho_u(x,j) + - m_rho[0]*0.3); // why 0.3? } - m_do_energy[i] = false; - } - } else { - if (m_do_energy[j]) { - changed = true; - } - m_do_energy[j] = false; - } - m_refiner->setActive(c_offset_U, false); - m_refiner->setActive(c_offset_V, false); - m_refiner->setActive(c_offset_T, false); - if (changed) { - needJacUpdate(); - } -} - -void StFlow::grad_hk(const double* x, size_t j) -{ - for(size_t k = 0; k < m_nsp; k++) { - if (u(x, j) > 0.0) { - m_dhk_dz(k,j) = (m_hk(k,j) - m_hk(k,j-1)) / m_dz[j - 1]; - } - else { - m_dhk_dz(k,j) = (m_hk(k,j+1) - m_hk(k,j)) / m_dz[j]; + } else if (grid(j) < m_zfixed) { + rsd[index(c_offset_U,j)] = + - (rho_u(x,j+1) - rho_u(x,j))/m_dz[j]; } - } -} - -// Two-point control functions -double StFlow::leftControlPointTemperature() const -{ - if (m_twoPointControl) { - if (m_zLeft != Undef) { - return m_tLeft; - } else { - throw CanteraError("StFlow::leftControlPointTemperature", - "Invalid operation: left control point location is not set"); - } - } else { - throw CanteraError("StFlow::leftControlPointTemperature", - "Invalid operation: two-point control is not enabled."); - } -} - -double StFlow::leftControlPointCoordinate() const -{ - if (m_twoPointControl) { - if (m_zLeft != Undef) { - return m_zLeft; - } else { - throw CanteraError("StFlow::leftControlPointCoordinate", - "Invalid operation: left control point location is not set"); - } - } else { - throw CanteraError("StFlow::leftControlPointCoordinate", - "Invalid operation: two-point control is not enabled."); - } -} - -void StFlow::setLeftControlPointTemperature(double temperature) -{ - if (m_twoPointControl) { - if (m_zLeft != Undef) { - m_tLeft = temperature; - } else { - throw CanteraError("StFlow::setLeftControlPointTemperature", - "Invalid operation: left control point location is not set"); - } - } else { - throw CanteraError("StFlow::setLeftControlPointTemperature", - "Invalid operation: two-point control is not enabled."); - } -} - -void StFlow::setLeftControlPointCoordinate(double z_left) -{ - if (m_twoPointControl) { - m_zLeft = z_left; - } else { - throw CanteraError("StFlow::setLeftControlPointCoordinate", - "Invalid operation: two-point control is not enabled."); - } -} - -double StFlow::rightControlPointTemperature() const -{ - if (m_twoPointControl) { - if (m_zRight != Undef) { - return m_tRight; - } else { - throw CanteraError("StFlow::rightControlPointTemperature", - "Invalid operation: right control point location is not set"); - } - } else { - throw CanteraError("StFlow::rightControlPointTemperature", - "Invalid operation: two-point control is not enabled."); - } -} - -double StFlow::rightControlPointCoordinate() const -{ - if (m_twoPointControl) { - if (m_zRight != Undef) { - return m_zRight; - } else { - throw CanteraError("StFlow::rightControlPointCoordinate", - "Invalid operation: right control point location is not set"); - } - } else { - throw CanteraError("StFlow::rightControlPointCoordinate", - "Invalid operation: two-point control is not enabled."); - } -} - -void StFlow::setRightControlPointTemperature(double temperature) -{ - if (m_twoPointControl) { - if (m_zRight != Undef) { - m_tRight = temperature; - } else { - throw CanteraError("StFlow::setRightControlPointTemperature", - "Invalid operation: right control point location is not set"); - } - } else { - throw CanteraError("StFlow::setRightControlPointTemperature", - "Invalid operation: two-point control is not enabled."); - } -} - -void StFlow::setRightControlPointCoordinate(double z_right) -{ - if (m_twoPointControl) { - m_zRight = z_right; - } else { - throw CanteraError("StFlow::setRightControlPointCoordinate", - "Invalid operation: two-point control is not enabled."); - } -} - -void StFlow::enableTwoPointControl(bool twoPointControl) -{ - if (isStrained()) { - m_twoPointControl = twoPointControl; } else { - throw CanteraError("StFlow::enableTwoPointControl", - "Invalid operation: two-point control can only be used" - "with axisymmetric flames."); + // unstrained with fixed mass flow rate + rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j - 1); } } diff --git a/src/oneD/refine.cpp b/src/oneD/refine.cpp index 719d26c47b2..5364e7e25e3 100644 --- a/src/oneD/refine.cpp +++ b/src/oneD/refine.cpp @@ -4,7 +4,7 @@ // at https://cantera.org/license.txt for license and copyright information. #include "cantera/oneD/refine.h" -#include "cantera/oneD/StFlow.h" +#include "cantera/oneD/Flow1D.h" #include "cantera/base/global.h" using namespace std; @@ -144,7 +144,7 @@ int Refiner::analyze(size_t n, const double* z, const double* x) } } - StFlow* fflame = dynamic_cast(m_domain); + Flow1D* fflame = dynamic_cast(m_domain); // Refine based on properties of the grid itself for (size_t j = 1; j < n-1; j++) { diff --git a/test/clib/test_ctonedim.cpp b/test/clib/test_ctonedim.cpp index 5a7a178b7be..cf4752b20d3 100644 --- a/test/clib/test_ctonedim.cpp +++ b/test/clib/test_ctonedim.cpp @@ -26,7 +26,7 @@ TEST(ctonedim, freeflow) int flow = domain_new("free-flow", sol, "flow"); ASSERT_GE(flow, 0); domain_setID(flow, "flow"); - ASSERT_NEAR(stflow_pressure(flow), P, 1e-5); + ASSERT_NEAR(flow1D_pressure(flow), P, 1e-5); int buflen = domain_type3(flow, 0, 0); char* buf = new char[buflen]; @@ -56,10 +56,10 @@ TEST(ctonedim, freeflow_from_parts) thermo_setPressure(ph, P); int itype = 2; // free flow - int flow = stflow_new(ph, kin, tr, itype); + int flow = flow1D_new(ph, kin, tr, itype); ASSERT_GE(flow, 0); domain_setID(flow, "flow"); - ASSERT_NEAR(stflow_pressure(flow), P, 1e-5); + ASSERT_NEAR(flow1D_pressure(flow), P, 1e-5); int buflen = domain_type3(flow, 0, 0); char* buf = new char[buflen]; @@ -124,7 +124,7 @@ TEST(ctonedim, catcomb_stack) // flow int itype = 1; // free flow - int flow = stflow_new(gas, gas_kin, trans, itype); + int flow = flow1D_new(gas, gas_kin, trans, itype); domain_setID(flow, "flow"); // reacting surface @@ -173,7 +173,7 @@ TEST(ctonedim, freeflame_from_parts) // flow int itype = 2; // free flow - int flow = stflow_new(ph, kin, tr, itype); + int flow = flow1D_new(ph, kin, tr, itype); domain_setID(flow, "flow"); // grid @@ -232,7 +232,7 @@ TEST(ctonedim, freeflame_from_parts) sim1D_setFixedTemperature(flame, 0.85 * T + .15 * Tad); // solve and save - stflow_solveEnergyEqn(flow, 1); + flow1D_solveEnergyEqn(flow, 1); bool refine_grid = false; int loglevel = 0; sim1D_solve(flame, loglevel, refine_grid); @@ -254,3 +254,25 @@ TEST(ctonedim, freeflame_from_parts) Tprev = T; } } + +TEST(ctonedim, stflow_tests) +{ + //! @todo: To be removed after Cantera 3.1 + ct_resetStorage(); + auto gas = newThermo("h2o2.yaml", "ohmech"); + + int sol = soln_newSolution("h2o2.yaml", "ohmech", "default"); + int ph = soln_thermo(sol); + int kin = soln_kinetics(sol); + int tr = soln_transport(sol); + + // spot check some errors + int itype = 2; // free flow + int ret = stflow_new(ph, kin, tr, itype); + ASSERT_EQ(ret, -1); // -1 is an error code + + int flow = flow1D_new(ph, kin, tr, itype); + ASSERT_EQ(stflow_setTransport(flow, tr), -1); + ASSERT_EQ(stflow_pressure(flow), DERR); // DERR is an error code + ASSERT_EQ(stflow_setPressure(flow, OneAtm), -1); +} diff --git a/test/clib/utils.h b/test/clib/utils.h new file mode 100644 index 00000000000..e69de29bb2d diff --git a/test/oneD/test_oneD.cpp b/test/oneD/test_oneD.cpp index 47de95d47b7..0e5e1d72ebd 100644 --- a/test/oneD/test_oneD.cpp +++ b/test/oneD/test_oneD.cpp @@ -5,6 +5,7 @@ #include "cantera/core.h" #include "cantera/onedim.h" #include "cantera/oneD/DomainFactory.h" +#include "cantera/oneD/StFlow.h" #include "cantera/oneD/IonFlow.h" using namespace Cantera; @@ -35,7 +36,7 @@ TEST(onedim, freeflame) double Tad = gas->temperature(); // flow - auto flow = newDomain("free-flow", sol, "flow"); + auto flow = newDomain("free-flow", sol, "flow"); // grid int nz = 21; @@ -102,15 +103,107 @@ TEST(onedim, freeflame) } } +TEST(onedim, legacy) +{ + //! @todo: Remove after Cantera 3.1 + auto sol = newSolution("h2o2.yaml", "ohmech", "mixture-averaged"); + auto gas = sol->thermo(); + size_t nsp = gas->nSpecies(); + + // reactants + double uin = .3; + double T = 300; + double P = 101325; + string X = "H2:0.65, O2:0.5, AR:2"; + gas->setState_TPX(T, P, X); + double rho_in = gas->density(); + vector yin(nsp); + gas->getMassFractions(&yin[0]); + + // product estimate + gas->equilibrate("HP"); + vector yout(nsp); + gas->getMassFractions(&yout[0]); + double rho_out = gas->density(); + double Tad = gas->temperature(); + + // flow + suppress_deprecation_warnings(); + auto flow = newDomain("legacy-flow", sol); + make_deprecation_warnings_fatal(); + flow->setID("flow"); + flow->setFreeFlow(); + + // grid + int nz = 21; + double lz = 0.02; + vector z(nz); + double dz = lz; + dz /= (double)(nz - 1); + for (int iz = 0; iz < nz; iz++) { + z[iz] = iz * dz; + } + flow->setupGrid(nz, &z[0]); + + // inlet + auto inlet = newDomain("inlet", sol); + inlet->setMoleFractions(X); + inlet->setMdot(uin * rho_in); + inlet->setTemperature(T); + + // outlet + auto outlet = newDomain("outlet", sol); + double uout = inlet->mdot() / rho_out; + + // set up simulation + vector> domains { inlet, flow, outlet }; + Sim1D flame(domains); + int dom = static_cast(flame.domainIndex("flow")); + ASSERT_EQ(dom, 1); + + // set up initial guess + vector locs{0.0, 0.3, 0.7, 1.0}; + vector value{uin, uin, uout, uout}; + flame.setInitialGuess("velocity", locs, value); + value = {T, T, Tad, Tad}; + flame.setInitialGuess("T", locs, value); + for (size_t i = 0; i < nsp; i++) { + value = {yin[i], yin[i], yout[i], yout[i]}; + flame.setInitialGuess(gas->speciesName(i), locs, value); + } + + // simulation settings + double ratio = 15.0; + double slope = 0.3; + double curve = 0.5; + flame.setRefineCriteria(dom, ratio, slope, curve); + flame.setFixedTemperature(0.85 * T + .15 * Tad); + + // solve + flow->solveEnergyEqn(); + bool refine_grid = false; + int loglevel = 0; + flame.solve(loglevel, refine_grid); + + ASSERT_EQ(flow->nPoints(), static_cast(nz + 1)); + size_t comp = flow->componentIndex("T"); + double Tprev = flame.value(dom, comp, 0); + for (size_t n = 0; n < flow->nPoints(); n++) { + T = flame.value(dom, comp, n); + ASSERT_GE(T, Tprev); + Tprev = T; + } +} + TEST(onedim, flame_types) { auto sol = newSolution("h2o2.yaml", "ohmech", "mixture-averaged"); - auto free = newDomain("free-flow", sol, "flow"); + auto free = newDomain("free-flow", sol, "flow"); ASSERT_EQ(free->type(), "free-flow"); - auto symm = newDomain("axisymmetric-flow", sol, "flow"); + auto symm = newDomain("axisymmetric-flow", sol, "flow"); ASSERT_EQ(symm->type(), "axisymmetric-flow"); - auto burner = newDomain("unstrained-flow", sol, "flow"); + auto burner = newDomain("unstrained-flow", sol, "flow"); ASSERT_EQ(burner->type(), "unstrained-flow"); }