From 0cf708ee315142abc574251952765ddbfd0e9860 Mon Sep 17 00:00:00 2001 From: sin-ha Date: Mon, 30 Mar 2020 18:30:22 +0530 Subject: [PATCH 01/24] new branch after squashing latest commits --- interfaces/cython/cantera/composite.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/interfaces/cython/cantera/composite.py b/interfaces/cython/cantera/composite.py index 2c126e08a08..c7c55b85269 100644 --- a/interfaces/cython/cantera/composite.py +++ b/interfaces/cython/cantera/composite.py @@ -509,6 +509,9 @@ def __init__(self, phase, shape=(0,), states=None, extra=None): if len(self._shape) == 1: self._indices = list(range(self._shape[0])) self._output_dummy = self._indices + elif (len(self._shape) == 0 and len(states) == 0): + self._indices = list() + self._output_dummy = list() else: self._indices = list(np.ndindex(self._shape)) self._output_dummy = self._states[..., 0] From 2aadf6bbe8ccfd6e3cf01ac928f26cd8ebbf30c9 Mon Sep 17 00:00:00 2001 From: sin-ha Date: Wed, 1 Apr 2020 19:26:48 +0530 Subject: [PATCH 02/24] Fixed slicing in SolutionArray with edge cases --- interfaces/cython/cantera/composite.py | 2 +- interfaces/cython/cantera/test/test_thermo.py | 10 ++++++++++ 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/interfaces/cython/cantera/composite.py b/interfaces/cython/cantera/composite.py index c7c55b85269..fe06226b575 100644 --- a/interfaces/cython/cantera/composite.py +++ b/interfaces/cython/cantera/composite.py @@ -509,7 +509,7 @@ def __init__(self, phase, shape=(0,), states=None, extra=None): if len(self._shape) == 1: self._indices = list(range(self._shape[0])) self._output_dummy = self._indices - elif (len(self._shape) == 0 and len(states) == 0): + elif len(self._shape) == 0 and len(states) == 0: self._indices = list() self._output_dummy = list() else: diff --git a/interfaces/cython/cantera/test/test_thermo.py b/interfaces/cython/cantera/test/test_thermo.py index 16ba5f6b507..27df9a81978 100644 --- a/interfaces/cython/cantera/test/test_thermo.py +++ b/interfaces/cython/cantera/test/test_thermo.py @@ -1803,3 +1803,13 @@ def test_slice_SolutionArray(self): soln = ct.SolutionArray(self.gas, 10) arr = soln[2:9:3] self.assertEqual(len(arr.T), 3) + + #tests for edge cases in sliced solution array + arr1 = soln[1:1] + self.assertEqual(arr1.n_elements, 3) + + arr2 = soln[slice(0)] + self.assertEqual(arr2.n_species, 9) + + arr3 = soln[None:0] + self.assertEqual(arr3.n_reactions, 28) From 911f52042a1e083ff51c11e73a4ff250525e6c3c Mon Sep 17 00:00:00 2001 From: band-a-prend Date: Sat, 29 Feb 2020 01:19:02 +0300 Subject: [PATCH 03/24] Simple fix for Sundials 5.1 compatibility The Sundials 5.1 has no any changes in CVODES API in comparison with sundials 4.0, 4.1, 5.0 therefore the compatibility configuration fix is trivial. --- SConstruct | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SConstruct b/SConstruct index 5e58d0c2da8..7416f69c6c2 100644 --- a/SConstruct +++ b/SConstruct @@ -1108,7 +1108,7 @@ if env['system_sundials'] == 'y': # Ignore the minor version, e.g. 2.4.x -> 2.4 env['sundials_version'] = '.'.join(sundials_version.split('.')[:2]) - if env['sundials_version'] not in ('2.4','2.5','2.6','2.7','3.0','3.1','3.2','4.0','4.1','5.0'): + if env['sundials_version'] not in ('2.4','2.5','2.6','2.7','3.0','3.1','3.2','4.0','4.1','5.0','5.1'): print("""ERROR: Sundials version %r is not supported.""" % env['sundials_version']) sys.exit(1) print("""INFO: Using system installation of Sundials version %s.""" % sundials_version) From bb18801c87e4c098eb6bd1a972626027df3690da Mon Sep 17 00:00:00 2001 From: band-a-prend Date: Wed, 4 Mar 2020 08:50:43 +0300 Subject: [PATCH 04/24] Update Sundials submodule sha1 (tag v5.1.0) --- ext/sundials | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/sundials b/ext/sundials index 38aa03b990e..b78e575639b 160000 --- a/ext/sundials +++ b/ext/sundials @@ -1 +1 @@ -Subproject commit 38aa03b990e616db52acd053730c31fcc3eb91f9 +Subproject commit b78e575639babe99aea9a0558d0f64732d6d729e From ceca0b0234cd285f52a3152cc07e9633c8d087b1 Mon Sep 17 00:00:00 2001 From: band-a-prend Date: Wed, 4 Mar 2020 09:21:27 +0300 Subject: [PATCH 05/24] Update Sundials submodule info to Sundials 5.1.0 --- SConstruct | 4 ++-- ext/sundials_config.h.in | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/SConstruct b/SConstruct index 7416f69c6c2..7f18c3428a6 100644 --- a/SConstruct +++ b/SConstruct @@ -1133,8 +1133,8 @@ if env['system_sundials'] == 'y': print('WARNING: External BLAS/LAPACK has been specified for Cantera ' 'but Sundials was built without this support.') else: # env['system_sundials'] == 'n' - print("""INFO: Using private installation of Sundials version 5.0.""") - env['sundials_version'] = '5.0' + print("""INFO: Using private installation of Sundials version 5.1.""") + env['sundials_version'] = '5.1' env['has_sundials_lapack'] = int(env['use_lapack']) diff --git a/ext/sundials_config.h.in b/ext/sundials_config.h.in index c53734b10fd..df25e7eef13 100644 --- a/ext/sundials_config.h.in +++ b/ext/sundials_config.h.in @@ -1,7 +1,7 @@ /* Minimal set of SUNDIALS configuration variables */ -#define SUNDIALS_VERSION "5.0.0" +#define SUNDIALS_VERSION "5.1.0" #define SUNDIALS_VERSION_MAJOR 5 -#define SUNDIALS_VERSION_MINOR 0 +#define SUNDIALS_VERSION_MINOR 1 #define SUNDIALS_VERSION_PATCH 0 #define SUNDIALS_VERSION_LABEL "" #define SUNDIALS_F77_FUNC(name,NAME) name ## _ From d33f8d42c99fd6f8778b52d5a096c7d292eb495b Mon Sep 17 00:00:00 2001 From: band-a-prend Date: Wed, 25 Mar 2020 10:24:18 +0300 Subject: [PATCH 06/24] Replace sundials version list with LooseVerion condition within version check --- SConstruct | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/SConstruct b/SConstruct index 7f18c3428a6..0afb835b204 100644 --- a/SConstruct +++ b/SConstruct @@ -1108,13 +1108,14 @@ if env['system_sundials'] == 'y': # Ignore the minor version, e.g. 2.4.x -> 2.4 env['sundials_version'] = '.'.join(sundials_version.split('.')[:2]) - if env['sundials_version'] not in ('2.4','2.5','2.6','2.7','3.0','3.1','3.2','4.0','4.1','5.0','5.1'): + sundials_ver = LooseVersion(env['sundials_version']) + if sundials_ver < LooseVersion('2.4') or sundials_ver >= LooseVersion('6.0'): print("""ERROR: Sundials version %r is not supported.""" % env['sundials_version']) sys.exit(1) print("""INFO: Using system installation of Sundials version %s.""" % sundials_version) #Determine whether or not Sundials was built with BLAS/LAPACK - if LooseVersion(env['sundials_version']) < LooseVersion('2.6'): + if sundials_ver < LooseVersion('2.6'): # In Sundials 2.4 / 2.5, SUNDIALS_BLAS_LAPACK is either 0 or 1 sundials_blas_lapack = get_expression_value(['"sundials/sundials_config.h"'], 'SUNDIALS_BLAS_LAPACK') @@ -1606,7 +1607,7 @@ linkSharedLibs = ['cantera_shared'] if env['system_sundials'] == 'y': env['sundials_libs'] = ['sundials_cvodes', 'sundials_ida', 'sundials_nvecserial'] - if env['use_lapack'] and LooseVersion(env['sundials_version']) >= LooseVersion('3.0'): + if env['use_lapack'] and sundials_ver >= LooseVersion('3.0'): if env.get('has_sundials_lapack'): env['sundials_libs'].extend(('sundials_sunlinsollapackdense', 'sundials_sunlinsollapackband')) From d945f53016d04590edd24e8e8c8ba986127ac9cd Mon Sep 17 00:00:00 2001 From: "Bryan W. Weber" Date: Fri, 27 Mar 2020 13:57:20 -0400 Subject: [PATCH 07/24] [SCons] Warn if untested SUNDIALS is used We expect that minor version number release will not introduce breaking changes, but the user should be warned that their version of SUNDIALS is untested. --- SConstruct | 3 +++ 1 file changed, 3 insertions(+) diff --git a/SConstruct b/SConstruct index 0afb835b204..5a4ac67748e 100644 --- a/SConstruct +++ b/SConstruct @@ -1112,6 +1112,9 @@ if env['system_sundials'] == 'y': if sundials_ver < LooseVersion('2.4') or sundials_ver >= LooseVersion('6.0'): print("""ERROR: Sundials version %r is not supported.""" % env['sundials_version']) sys.exit(1) + elif sundials_ver > LooseVersion('5.1'): + print("WARNING: Sundials version %r has not been tested." % env['sundials_version']) + print("""INFO: Using system installation of Sundials version %s.""" % sundials_version) #Determine whether or not Sundials was built with BLAS/LAPACK From 45e8df7a86417bfe5d4b367e4f14666a0cc189e4 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 3 Jan 2020 17:25:04 -0500 Subject: [PATCH 08/24] [Thermo] Deprecate "gas" mode of IdealSolnGasVPSS --- doc/sphinx/yaml/phases.rst | 3 +++ include/cantera/thermo/IdealSolnGasVPSS.h | 9 ++++++++- interfaces/cython/cantera/test/test_convert.py | 5 ++++- test/data/thermo-models.yaml | 2 ++ test/thermo/phaseConstructors.cpp | 4 ++++ test/thermo/thermoFromYaml.cpp | 4 ++++ test_problems/VPsilane_test/silane_equil.cpp | 3 +++ 7 files changed, 28 insertions(+), 2 deletions(-) diff --git a/doc/sphinx/yaml/phases.rst b/doc/sphinx/yaml/phases.rst index db4d7adf916..fe9f40af892 100644 --- a/doc/sphinx/yaml/phases.rst +++ b/doc/sphinx/yaml/phases.rst @@ -538,6 +538,9 @@ The ideal gas model as The ideal gas model, using variable pressure standard state methods as `described here `__. +.. deprecated:: 2.5 + Use the ``ideal-gas`` model instead. + .. _sec-yaml-ideal-molal-solution: diff --git a/include/cantera/thermo/IdealSolnGasVPSS.h b/include/cantera/thermo/IdealSolnGasVPSS.h index 4e4639a1a60..197d6f9c40e 100644 --- a/include/cantera/thermo/IdealSolnGasVPSS.h +++ b/include/cantera/thermo/IdealSolnGasVPSS.h @@ -23,6 +23,9 @@ namespace Cantera * * An ideal solution or an ideal gas approximation of a phase. Uses variable * pressure standard state methods for calculating thermodynamic properties. + * + * @deprecated "Gas" mode to be removed after Cantera 2.5. Use IdealGasPhase for + * ideal gas phases instead. */ class IdealSolnGasVPSS : public VPStandardStateTP { @@ -62,7 +65,11 @@ class IdealSolnGasVPSS : public VPStandardStateTP } //! Set this phase to represent an ideal gas - void setGasMode() { m_idealGas = true; } + void setGasMode() { + warn_deprecated("IdealSolnGasVPSS::setGasMode", + "To be removed after Cantera 2.5. Use class IdealGasPhase instead."); + m_idealGas = true; + } //! Set this phase to represent an ideal liquid or solid solution void setSolnMode() { m_idealGas = false; } diff --git a/interfaces/cython/cantera/test/test_convert.py b/interfaces/cython/cantera/test/test_convert.py index 4dac506c807..51d5e0c70bf 100644 --- a/interfaces/cython/cantera/test/test_convert.py +++ b/interfaces/cython/cantera/test/test_convert.py @@ -1187,7 +1187,9 @@ def test_vpss_and_hkft(self): Path(self.test_data_dir).joinpath("pdss_hkft.xml"), Path(self.test_work_dir).joinpath("pdss_hkft.yaml"), ) - + # @todo Remove "gas" mode test after Cantera 2.5 - "gas" mode of class + # IdealSolnGasVPSS is deprecated + ct.suppress_deprecation_warnings() for name in ["vpss_gas_pdss_hkft_phase", "vpss_soln_pdss_hkft_phase"]: ctmlPhase = ct.ThermoPhase("pdss_hkft.xml", name=name) yamlPhase = ct.ThermoPhase("pdss_hkft.yaml", name=name) @@ -1200,6 +1202,7 @@ def test_vpss_and_hkft(self): self.assertEqual(ctmlPhase.element_names, yamlPhase.element_names) self.assertEqual(ctmlPhase.species_names, yamlPhase.species_names) self.checkThermo(ctmlPhase, yamlPhase, [300, 500]) + ct.make_deprecation_warnings_fatal() def test_lattice_solid(self): ctml2yaml.convert( diff --git a/test/data/thermo-models.yaml b/test/data/thermo-models.yaml index beb17c40128..b7007583ae8 100644 --- a/test/data/thermo-models.yaml +++ b/test/data/thermo-models.yaml @@ -93,6 +93,8 @@ phases: species: [KCl(l)] thermo: Margules +# todo: Remove after Cantera 2.5 - the "gas" moode of class IdealSolnGasVPSS is +# deprecated. - name: IdealSolnGas-gas thermo: ideal-gas-VPSS species: [{gas-species: [H2O, H2, O2]}] diff --git a/test/thermo/phaseConstructors.cpp b/test/thermo/phaseConstructors.cpp index 81179f06f11..3d56d9200cf 100644 --- a/test/thermo/phaseConstructors.cpp +++ b/test/thermo/phaseConstructors.cpp @@ -344,8 +344,11 @@ TEST_F(ConstructFromScratch, RedlichKwongMFTP_missing_coeffs) EXPECT_THROW(p.setState_TP(300, 200e5), CanteraError); } +//! @todo Remove after Cantera 2.5 - "gas" mode of IdealSolnGasVPSS is +//! deprecated TEST_F(ConstructFromScratch, IdealSolnGasVPSS_gas) { + suppress_deprecation_warnings(); IdealSolnGasVPSS p; p.addSpecies(sH2O); p.addSpecies(sH2); @@ -367,6 +370,7 @@ TEST_F(ConstructFromScratch, IdealSolnGasVPSS_gas) EXPECT_NEAR(p.temperature(), 479.929, 1e-3); // based on h2o2.cti EXPECT_NEAR(p.moleFraction("H2O"), 0.01, 1e-4); EXPECT_NEAR(p.moleFraction("H2"), 0.0, 1e-4); + make_deprecation_warnings_fatal(); } TEST(PureFluidFromScratch, CarbonDioxide) diff --git a/test/thermo/thermoFromYaml.cpp b/test/thermo/thermoFromYaml.cpp index 856d96dcf7c..eb35b74037e 100644 --- a/test/thermo/thermoFromYaml.cpp +++ b/test/thermo/thermoFromYaml.cpp @@ -234,9 +234,13 @@ TEST(ThermoFromYaml, IonsFromNeutral_fromString) EXPECT_NEAR(mu[1], -2.88157316e+06, 1e-1); } +//! @todo Remove after Cantera 2.5 - "gas" mode of IdealSolnGasVPSS is +//! deprecated TEST(ThermoFromYaml, IdealSolnGas_gas) { + suppress_deprecation_warnings(); auto thermo = newThermo("thermo-models.yaml", "IdealSolnGas-gas"); + make_deprecation_warnings_fatal(); thermo->equilibrate("HP"); EXPECT_NEAR(thermo->temperature(), 479.929, 1e-3); // based on h2o2.cti EXPECT_NEAR(thermo->moleFraction("H2O"), 0.01, 1e-4); diff --git a/test_problems/VPsilane_test/silane_equil.cpp b/test_problems/VPsilane_test/silane_equil.cpp index f1cf5aa4770..98f0be48ab5 100644 --- a/test_problems/VPsilane_test/silane_equil.cpp +++ b/test_problems/VPsilane_test/silane_equil.cpp @@ -9,12 +9,15 @@ using namespace std; using namespace Cantera; +//! @todo Remove this test after Cantera 2.5 - "gas" mode of class +//! IdealSolnGasVPSS is deprecated int main(int argc, char** argv) { #if defined(_MSC_VER) && _MSC_VER < 1900 _set_output_format(_TWO_DIGIT_EXPONENT); #endif try { + suppress_deprecation_warnings(); IdealSolnGasVPSS gg("silane.xml", "silane"); ThermoPhase* g = ≫ cout.precision(4); From 931f088893e83c624514deaca68ac02630493ddd Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 3 Jan 2020 19:22:09 -0500 Subject: [PATCH 09/24] [Thermo] Deprecate class FixedChemPotSSTP This model can easily be replaced with StoichSustance using ConstCpPoly for the species thermo. The changed values in the LiSi regression tests suggest that there are some implementation errors in the FixedChemPotSSTP class which will be automatically resolved by eliminating it. --- doc/sphinx/yaml/phases.rst | 4 +++ ext/sundials | 2 +- .../cython/cantera/test/test_convert.py | 3 +++ src/thermo/FixedChemPotSSTP.cpp | 12 +++++++++ test/data/thermo-models.yaml | 1 + test/thermo/phaseConstructors.cpp | 6 ++++- test/thermo/thermoFromYaml.cpp | 3 +++ .../LatticeSolid_LiSi/Li7Si3_ls.xml | 26 ++++++++++++++++++- .../VCSnonideal/LatticeSolid_LiSi/latsol.cpp | 15 ++++++----- .../LatticeSolid_LiSi/output_blessed.txt | 13 +++++----- .../LatticeSolid_LiSi/verbose_blessed.txt | 15 +++++------ 11 files changed, 76 insertions(+), 24 deletions(-) diff --git a/doc/sphinx/yaml/phases.rst b/doc/sphinx/yaml/phases.rst index fe9f40af892..012beabdc82 100644 --- a/doc/sphinx/yaml/phases.rst +++ b/doc/sphinx/yaml/phases.rst @@ -378,6 +378,10 @@ Example:: A phase defined by a fixed value of the chemical potential, as `described here `__. +.. deprecated:: 2.5 + Use class StoichSubstance with a constant-cp species thermo model, with 'h0' + set to the desired chemical potential and 's0' set to 0. + Additional fields: ``chemical-potential`` diff --git a/ext/sundials b/ext/sundials index b78e575639b..38aa03b990e 160000 --- a/ext/sundials +++ b/ext/sundials @@ -1 +1 @@ -Subproject commit b78e575639babe99aea9a0558d0f64732d6d729e +Subproject commit 38aa03b990e616db52acd053730c31fcc3eb91f9 diff --git a/interfaces/cython/cantera/test/test_convert.py b/interfaces/cython/cantera/test/test_convert.py index 51d5e0c70bf..483befa7dfb 100644 --- a/interfaces/cython/cantera/test/test_convert.py +++ b/interfaces/cython/cantera/test/test_convert.py @@ -1054,13 +1054,16 @@ def test_fractional_stoich_coeffs(self): self.checkThermo(ctmlGas, yamlGas, [300, 500, 1300, 2000]) self.checkKinetics(ctmlGas, yamlGas, [900, 1800], [2e5, 20e5]) + # @todo Remove after Cantera 2.5 - class FixedChemPotSSTP is deprecated def test_fixed_chemical_potential_thermo(self): + ct.suppress_deprecation_warnings() ctml2yaml.convert( Path(self.test_data_dir).joinpath("LiFixed.xml"), Path(self.test_work_dir).joinpath("LiFixed.yaml"), ) ctmlGas, yamlGas = self.checkConversion("LiFixed") self.checkThermo(ctmlGas, yamlGas, [300, 500, 1300, 2000]) + ct.make_deprecation_warnings_fatal() def test_water_IAPWS95_thermo(self): ctml2yaml.convert( diff --git a/src/thermo/FixedChemPotSSTP.cpp b/src/thermo/FixedChemPotSSTP.cpp index 4d25cb035fb..b47b71b959c 100644 --- a/src/thermo/FixedChemPotSSTP.cpp +++ b/src/thermo/FixedChemPotSSTP.cpp @@ -23,22 +23,34 @@ namespace Cantera FixedChemPotSSTP::FixedChemPotSSTP() : chemPot_(0.0) { + warn_deprecated("class FixedChemPotSSTP", "To be removed after Cantera 2.5. " + "Use class StoichSubstance with a constant-cp species thermo model, " + "with 'h0' set to the desired chemical potential and 's0' set to 0."); } FixedChemPotSSTP::FixedChemPotSSTP(const std::string& infile, const std::string& id_) : chemPot_(0.0) { + warn_deprecated("class FixedChemPotSSTP", "To be removed after Cantera 2.5. " + "Use class StoichSubstance with a constant-cp species thermo model, " + "with 'h0' set to the desired chemical potential and 's0' set to 0."); initThermoFile(infile, id_); } FixedChemPotSSTP::FixedChemPotSSTP(XML_Node& xmlphase, const std::string& id_) : chemPot_(0.0) { + warn_deprecated("class FixedChemPotSSTP", "To be removed after Cantera 2.5. " + "Use class StoichSubstance with a constant-cp species thermo model, " + "with 'h0' set to the desired chemical potential and 's0' set to 0."); importPhase(xmlphase, this); } FixedChemPotSSTP::FixedChemPotSSTP(const std::string& Ename, doublereal val) : chemPot_(0.0) { + warn_deprecated("class FixedChemPotSSTP", "To be removed after Cantera 2.5. " + "Use class StoichSubstance with a constant-cp species thermo model, " + "with 'h0' set to the desired chemical potential and 's0' set to 0."); std::string pname = Ename + "Fixed"; setName(pname); setNDim(3); diff --git a/test/data/thermo-models.yaml b/test/data/thermo-models.yaml index b7007583ae8..bd5c918569d 100644 --- a/test/data/thermo-models.yaml +++ b/test/data/thermo-models.yaml @@ -15,6 +15,7 @@ phases: transport: water species: [H2O(l)] +# @todo Remove after Cantera 2.5 - class FixedChemPotSSTP is deprecated - name: Li-fixed thermo: fixed-chemical-potential species: [Li] diff --git a/test/thermo/phaseConstructors.cpp b/test/thermo/phaseConstructors.cpp index 3d56d9200cf..c4c68b59e5d 100644 --- a/test/thermo/phaseConstructors.cpp +++ b/test/thermo/phaseConstructors.cpp @@ -77,27 +77,31 @@ shared_ptr make_const_cp_species(const std::string& name, return species; } - +//! @todo Remove after Cantera 2.5 - class FixedChemPotSSTP is deprecated class FixedChemPotSstpConstructorTest : public testing::Test { }; TEST_F(FixedChemPotSstpConstructorTest, fromXML) { + suppress_deprecation_warnings(); std::unique_ptr p(newPhase("../data/LiFixed.xml")); ASSERT_EQ((int) p->nSpecies(), 1); double mu; p->getChemPotentials(&mu); ASSERT_DOUBLE_EQ(-2.3e7, mu); + make_deprecation_warnings_fatal(); } TEST_F(FixedChemPotSstpConstructorTest, SimpleConstructor) { + suppress_deprecation_warnings(); FixedChemPotSSTP p("Li", -2.3e7); ASSERT_EQ((int) p.nSpecies(), 1); double mu; p.getChemPotentials(&mu); ASSERT_DOUBLE_EQ(-2.3e7, mu); + make_deprecation_warnings_fatal(); } TEST(IonsFromNeutralConstructor, fromXML) diff --git a/test/thermo/thermoFromYaml.cpp b/test/thermo/thermoFromYaml.cpp index eb35b74037e..b3f3a3024eb 100644 --- a/test/thermo/thermoFromYaml.cpp +++ b/test/thermo/thermoFromYaml.cpp @@ -120,13 +120,16 @@ TEST(ThermoFromYaml, WaterSSTP) EXPECT_NEAR(thermo->enthalpy_mass(), -15649685.52296013, 1e-6); } +//! @todo Remove after Cantera 2.5 - class FixedChemPotSSTP is deprecated TEST(ThermoFromYaml, FixedChemPot) { + suppress_deprecation_warnings(); auto thermo = newThermo("thermo-models.yaml", "Li-fixed"); EXPECT_EQ(thermo->nSpecies(), (size_t) 1); double mu; thermo->getChemPotentials(&mu); EXPECT_DOUBLE_EQ(mu, -2.3e7); + make_deprecation_warnings_fatal(); } TEST(ThermoFromYaml, Margules) diff --git a/test_problems/VCSnonideal/LatticeSolid_LiSi/Li7Si3_ls.xml b/test_problems/VCSnonideal/LatticeSolid_LiSi/Li7Si3_ls.xml index 8b3021dd43a..7e9b8129ee6 100644 --- a/test_problems/VCSnonideal/LatticeSolid_LiSi/Li7Si3_ls.xml +++ b/test_problems/VCSnonideal/LatticeSolid_LiSi/Li7Si3_ls.xml @@ -25,7 +25,7 @@ - + @@ -55,6 +55,18 @@ + + + Li + + LiFixed + + 0.534 + + + + + @@ -148,6 +160,18 @@ + + + Li:1 + + + 298.15 + -2.3e7 + 0.0 + 0.0 + + + diff --git a/test_problems/VCSnonideal/LatticeSolid_LiSi/latsol.cpp b/test_problems/VCSnonideal/LatticeSolid_LiSi/latsol.cpp index 33c5b374592..5e3360620c3 100644 --- a/test_problems/VCSnonideal/LatticeSolid_LiSi/latsol.cpp +++ b/test_problems/VCSnonideal/LatticeSolid_LiSi/latsol.cpp @@ -2,8 +2,8 @@ #include "cantera/thermo/MargulesVPSSTP.h" #include "cantera/thermo/FixedChemPotSSTP.h" #include "cantera/thermo/LatticeSolidPhase.h" +#include "cantera/thermo/ConstCpPoly.h" #include "cantera/equil/vcs_MultiPhaseEquil.h" -#include "cantera/thermo.h" #include @@ -18,7 +18,7 @@ void testProblem(int printLvl) std::unique_ptr LiSi_solid(newPhase("Li7Si3_ls.xml", "Li7Si3_and_Interstitials(S)")); std::unique_ptr Li_liq(newPhase("Li_Liquid.xml", "Li(L)")); - FixedChemPotSSTP LiFixed("Li", -2.3E7); + std::unique_ptr LiFixed(newPhase("Li7Si3_ls.xml", "LiFixed")); MargulesVPSSTP salt("LiKCl_liquid.xml", "MoltenSalt_electrolyte"); // set states @@ -32,10 +32,10 @@ void testProblem(int printLvl) int ee = static_cast(LiSi_solid->nElements()); printf("Number of elements = %d\n", ee); - LiFixed.setState_TP(T, OneAtm); + LiFixed->setState_TP(T, OneAtm); double um[20]; - LiFixed.getChemPotentials(um); + LiFixed->getChemPotentials(um); printf(" chem pot = %g\n", um[0]); double volts = 1.635; // has some Fe in it // test suite @@ -45,13 +45,16 @@ void testProblem(int printLvl) Li_liq->getChemPotentials(um); double um_li_chempot = um[0] + dg_corr; printf("um_li_chempot = %g\n", um_li_chempot); - LiFixed.setChemicalPotential(um_li_chempot); + auto LiFixed_species = LiFixed->species(0); + auto spthermo = std::dynamic_pointer_cast(LiFixed_species->thermo); + spthermo->setParameters(298.15, um_li_chempot, 0, 0); + LiFixed->modifySpecies(0, LiFixed_species); MultiPhase mmm; mmm.addPhase(&salt, 10.); mmm.addPhase(LiSi_solid.get(), 1.); - mmm.addPhase(&LiFixed, 100.); + mmm.addPhase(LiFixed.get(), 100.); int estimateEquil = 0; diff --git a/test_problems/VCSnonideal/LatticeSolid_LiSi/output_blessed.txt b/test_problems/VCSnonideal/LatticeSolid_LiSi/output_blessed.txt index d61d5addd7a..64e19f3a984 100644 --- a/test_problems/VCSnonideal/LatticeSolid_LiSi/output_blessed.txt +++ b/test_problems/VCSnonideal/LatticeSolid_LiSi/output_blessed.txt @@ -5,7 +5,6 @@ um_li_chempot = -1.65555e+08 Trying VCS equilibrium solver Unknown Cantera EOS to VCSnonideal: 'Margules' Unknown Cantera EOS to VCSnonideal: 'LatticeSolid' -Unknown Cantera EOS to VCSnonideal: 'FixedChemPot' ================================================================================ ================ Cantera_to_vprob: START OF PROBLEM STATEMENT ==================== @@ -24,7 +23,7 @@ Unknown Cantera EOS to VCSnonideal: 'FixedChemPot' PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec TMolesInert Tmoles(kmol) MoltenSalt_electrolyte 0 0 0 UnkType: -1 2 0.000000e+00 1.000000e+01 Li7Si3_and_Interstitials(S) 1 0 0 UnkType: -1 3 0.000000e+00 1.000000e+00 - LiFixed 2 1 0 UnkType: -1 1 0.000000e+00 1.000000e+02 + LiFixed 2 1 0 Stoich Sub 1 0.000000e+00 1.000000e+02 ================================================================================ ================ Cantera_to_vprob: END OF PROBLEM STATEMENT ==================== @@ -49,7 +48,7 @@ Li7Si3_and_Interstitials(S) 1 0 0 UnkType: -1 3 0 PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec TMolesInert Tmoles(kmol) MoltenSalt_electrolyte 0 0 0 UnkType: -1 2 0.000000e+00 1.000000e+01 Li7Si3_and_Interstitials(S) 1 0 0 UnkType: -1 3 0.000000e+00 1.000000e+00 - LiFixed 2 1 0 UnkType: -1 1 0.000000e+00 1.000000e+02 + LiFixed 2 1 0 Stoich Sub 1 0.000000e+00 1.000000e+02 ================================================================================ ==================== Cantera_to_vprob: END OF PROBLEM STATEMENT ==================== @@ -78,7 +77,7 @@ Li7Si3_and_Interstitials(S) 1 0 0 UnkType: -1 3 0 PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec TMolesInert TKmoles MoltenSalt_electrolyte 0 0 0 UnkType: -1 2 0.000000e+00 1.000000e+01 Li7Si3_and_Interstitials(S) 1 0 0 UnkType: -1 3 0.000000e+00 1.000000e+00 - LiFixed 2 1 0 UnkType: -1 1 0.000000e+00 1.000000e+02 + LiFixed 2 1 0 Stoich Sub 1 0.000000e+00 1.000000e+02 Elemental Abundances: Target_kmol ElemType ElActive Li 1.105050000000E+02 0 1 @@ -145,7 +144,7 @@ VCS CALCULATION METHOD -------------------------------------------------------------------------------- Temperature = 6.3e+02 Kelvin Pressure = 1.0132e+05 Pa - total Volume = 0.35041 m**3 + total Volume = 1.6501 m**3 -------------------------------------------------------------------------------- @@ -295,14 +294,14 @@ Moles: 100.005 temperature 625.15 K pressure 1.0132e+05 Pa - density 0.001 kg/m^3 + density 534 kg/m^3 mean mol. weight 6.94 kg/kmol phase of matter unspecified 1 kg 1 kmol --------------- --------------- enthalpy -2.3855e+07 -1.6556e+08 J - internal energy -1.2399e+14 -8.6052e+14 J + internal energy -2.3855e+07 -1.6556e+08 J entropy 0 0 J/K Gibbs function -2.3855e+07 -1.6556e+08 J heat capacity c_p 0 0 J/K diff --git a/test_problems/VCSnonideal/LatticeSolid_LiSi/verbose_blessed.txt b/test_problems/VCSnonideal/LatticeSolid_LiSi/verbose_blessed.txt index 6192e0b8128..a00d6d864f7 100644 --- a/test_problems/VCSnonideal/LatticeSolid_LiSi/verbose_blessed.txt +++ b/test_problems/VCSnonideal/LatticeSolid_LiSi/verbose_blessed.txt @@ -10,7 +10,6 @@ Unknown Cantera EOS to VCSnonideal: 'LatticeSolid' vcs_Cantera_convert: Species Type 8 not known vcs_Cantera_convert: Species Type 1 not known vcs_Cantera_convert: Species Type 1 not known -Unknown Cantera EOS to VCSnonideal: 'FixedChemPot' vcs_Cantera_convert: Species Type 1 not known ================================================================================ @@ -30,7 +29,7 @@ vcs_Cantera_convert: Species Type 1 not known PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec TMolesInert Tmoles(kmol) MoltenSalt_electrolyte 0 0 0 UnkType: -1 2 0.000000e+00 1.000000e+01 Li7Si3_and_Interstitials(S) 1 0 0 UnkType: -1 3 0.000000e+00 1.000000e+00 - LiFixed 2 1 0 UnkType: -1 1 0.000000e+00 1.000000e+02 + LiFixed 2 1 0 Stoich Sub 1 0.000000e+00 1.000000e+02 ================================================================================ ================ Cantera_to_vprob: END OF PROBLEM STATEMENT ==================== @@ -55,7 +54,7 @@ Li7Si3_and_Interstitials(S) 1 0 0 UnkType: -1 3 0 PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec TMolesInert Tmoles(kmol) MoltenSalt_electrolyte 0 0 0 UnkType: -1 2 0.000000e+00 1.000000e+01 Li7Si3_and_Interstitials(S) 1 0 0 UnkType: -1 3 0.000000e+00 1.000000e+00 - LiFixed 2 1 0 UnkType: -1 1 0.000000e+00 1.000000e+02 + LiFixed 2 1 0 Stoich Sub 1 0.000000e+00 1.000000e+02 ================================================================================ ==================== Cantera_to_vprob: END OF PROBLEM STATEMENT ==================== @@ -84,7 +83,7 @@ Li7Si3_and_Interstitials(S) 1 0 0 UnkType: -1 3 0 PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec TMolesInert TKmoles MoltenSalt_electrolyte 0 0 0 UnkType: -1 2 0.000000e+00 1.000000e+01 Li7Si3_and_Interstitials(S) 1 0 0 UnkType: -1 3 0.000000e+00 1.000000e+00 - LiFixed 2 1 0 UnkType: -1 1 0.000000e+00 1.000000e+02 + LiFixed 2 1 0 Stoich Sub 1 0.000000e+00 1.000000e+02 Elemental Abundances: Target_kmol ElemType ElActive Li 1.105050000000E+02 0 1 @@ -343,7 +342,7 @@ VCS CALCULATION METHOD --- Total tentative Dimensionless Gibbs Free Energy = -4.1206857517346E+03 --- subroutine FORCE: Beginning Slope = -1.39826e-19 --- subroutine FORCE: End Slope = 9.31017e-32 - --- subroutine FORCE produced no adjustments (al = 1) + --- subroutine FORCE produced no adjustments, s2 < 0 ------------------------------------------------------------------------------------------------------- --- Summary of the Update (all species): --- Species Status Initial_KMoles Final_KMoles Initial_Mu/RT Mu/RT Init_Del_G/RT Delta_G/RT @@ -384,7 +383,7 @@ VCS CALCULATION METHOD -------------------------------------------------------------------------------- Temperature = 6.3e+02 Kelvin Pressure = 1.0132e+05 Pa - total Volume = 0.35041 m**3 + total Volume = 1.6501 m**3 -------------------------------------------------------------------------------- @@ -534,14 +533,14 @@ Moles: 100.005 temperature 625.15 K pressure 1.0132e+05 Pa - density 0.001 kg/m^3 + density 534 kg/m^3 mean mol. weight 6.94 kg/kmol phase of matter unspecified 1 kg 1 kmol --------------- --------------- enthalpy -2.3855e+07 -1.6556e+08 J - internal energy -1.2399e+14 -8.6052e+14 J + internal energy -2.3855e+07 -1.6556e+08 J entropy 0 0 J/K Gibbs function -2.3855e+07 -1.6556e+08 J heat capacity c_p 0 0 J/K From 7a080d660eab370691b81bd52dc3596e7d53e7df Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Tue, 21 Jan 2020 16:54:49 -0600 Subject: [PATCH 10/24] [numerics] Implement tabulated Func1 object - add new derived class to Func1.h that interpolates a tabulated function in C++ --- include/cantera/numerics/Func1.h | 62 ++++++++++++++++++++++++++++++-- src/numerics/Func1.cpp | 51 ++++++++++++++++++++++++++ 2 files changed, 110 insertions(+), 3 deletions(-) diff --git a/include/cantera/numerics/Func1.h b/include/cantera/numerics/Func1.h index 9c97a122edc..2af040f2d93 100644 --- a/include/cantera/numerics/Func1.h +++ b/include/cantera/numerics/Func1.h @@ -9,6 +9,7 @@ #define CT_FUNC1_H #include "cantera/base/ct_defs.h" +#include "cantera/base/ctexceptions.h" #include @@ -32,6 +33,7 @@ const int CosFuncType = 102; const int ExpFuncType = 104; const int PowFuncType = 106; const int ConstFuncType = 110; +const int TabulatedFuncType = 120; class TimesConstant1; @@ -123,6 +125,7 @@ Func1& newCompositeFunction(Func1& f1, Func1& f2); Func1& newTimesConstFunction(Func1& f1, doublereal c); Func1& newPlusConstFunction(Func1& f1, doublereal c); + //! implements the sin() function /*! * The argument to sin() is in radians @@ -165,7 +168,11 @@ class Sin1 : public Func1 virtual Func1& derivative() const; }; -/// cos + +//! implements the cos() function +/*! + * The argument to cos() is in radians + */ class Cos1 : public Func1 { public: @@ -200,7 +207,8 @@ class Cos1 : public Func1 virtual Func1& derivative() const; }; -/// exp + +//! implements the exponential function class Exp1 : public Func1 { public: @@ -233,7 +241,8 @@ class Exp1 : public Func1 virtual Func1& derivative() const; }; -/// pow + +//! implements the power function (pow) class Pow1 : public Func1 { public: @@ -265,6 +274,53 @@ class Pow1 : public Func1 virtual Func1& derivative() const; }; + +//! implements a tabulated function +class Tabulated1 : public Func1 +{ +public: + Tabulated1(const std::vector& tvec, const std::vector& fvec) : + Func1() { + if (tvec.size() != fvec.size()) { + throw CanteraError("Tabulated1::Tabulated1", + "sizes of vectors do not match ({} vs {}).", + tvec.size(), fvec.size()); + } else if (tvec.empty()) { + throw CanteraError("Tabulated1::Tabulated1", + "vectors must not be empty."); + } + m_tvec = tvec; + m_fvec = fvec; + } + + Tabulated1(const Tabulated1& b) : + Func1(b) { + } + + Tabulated1& operator=(const Tabulated1& right) { + if (&right == this) { + return *this; + } + Func1::operator=(right); + return *this; + } + + virtual std::string write(const std::string& arg) const; + virtual int ID() const { + return TabulatedFuncType; + } + virtual double eval(double t) const; + virtual Func1& duplicate() const { + return *(new Tabulated1(m_tvec, m_fvec)); + } + + virtual Func1& derivative() const; +private: + std::vector m_tvec; + std::vector m_fvec; +}; + + /** * Constant. */ diff --git a/src/numerics/Func1.cpp b/src/numerics/Func1.cpp index c4f942690e1..7feeb10fa2d 100644 --- a/src/numerics/Func1.cpp +++ b/src/numerics/Func1.cpp @@ -213,6 +213,52 @@ Func1& Pow1::derivative() const return *r; } +/******************************************************************************/ + +double Tabulated1::eval(double t) const { + size_t siz = m_tvec.size(); + if (siz) { + if (t <= m_tvec[0]) { + return m_fvec[0]; + } else if (t >= m_tvec[siz-1]) { + return m_fvec[siz-1]; + } else { + size_t ix = 0; + while (t > m_tvec[ix+1]) { + ix++; + } + double df = m_fvec[ix+1] - m_fvec[ix]; + df /= m_tvec[ix+1] - m_tvec[ix]; + df *= t - m_tvec[ix]; + return m_fvec[ix] + df; + } + } else { + return 0.; + } +} + +Func1& Tabulated1::derivative() const { + std::vector tvec; + std::vector dvec; + size_t siz = m_tvec.size(); + if (siz>1) { + for (size_t i=1; i({})", ID(), arg); @@ -233,6 +279,11 @@ string Pow1::write(const std::string& arg) const } } +string Tabulated1::write(const std::string& arg) const +{ + return fmt::format("\\mathrm{{Tabulated}}({})", arg); +} + string Const1::write(const std::string& arg) const { return fmt::format("{}", m_c); From 84d20186f143fcfc50592849451c76209534c5c1 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Tue, 21 Jan 2020 16:56:17 -0600 Subject: [PATCH 11/24] [numerics] Break out tabulated Func1 to Cython interface - the C++ defined Tabulated1 class is added as a variant of the Func1 object defined in Python - this approach is compatible with existing interfaces of various Python methods with Func1 arguments --- interfaces/cython/cantera/_cantera.pxd | 7 ++ interfaces/cython/cantera/func1.pyx | 93 +++++++++++++++++++++----- 2 files changed, 85 insertions(+), 15 deletions(-) diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index dd12c9de196..c6f59d743ed 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -38,6 +38,11 @@ cdef extern from "cantera/cython/funcWrapper.h": CxxFunc1(callback_wrapper, void*) double eval(double) except +translate_exception +cdef extern from "cantera/numerics/Func1.h": + cdef cppclass CxxTabulated1 "Cantera::Tabulated1": + CxxTabulated1(vector[double]&, vector[double]&) except +translate_exception + double eval(double) except +translate_exception + cdef extern from "cantera/base/xml.h" namespace "Cantera": cdef cppclass XML_Node: XML_Node* findByName(string) @@ -1017,6 +1022,8 @@ cdef class Func1: cdef CxxFunc1* func cdef object callable cdef object exception + cpdef void __set_callback(self, object) except * + cpdef void __set_tables(self, object, object) except * cdef class ReactorBase: cdef CxxReactorBase* rbase diff --git a/interfaces/cython/cantera/func1.pyx b/interfaces/cython/cantera/func1.pyx index 3cda2bd24be..f70a3354a8f 100644 --- a/interfaces/cython/cantera/func1.pyx +++ b/interfaces/cython/cantera/func1.pyx @@ -2,6 +2,7 @@ # at https://cantera.org/license.txt for license and copyright information. import sys +from numbers import Real as _real cdef double func_callback(double t, void* obj, void** err): """ @@ -46,36 +47,98 @@ cdef class Func1: >>> f3(6) 30.0 - For simplicity, constant functions can be defined by passing the constant + For simplicity, constant functions can be defined by passing a constant value directly:: >>> f4 = Func1(2.5) >>> f4(0.1) 2.5 + A `Func1` object representing a tabulated function is defined by sample + points and corresponding function values. Inputs are specified either by + two iterable objects containing sample point location and function values, + or a single array that concatenates those inputs in two columns. Between + sample points, values are evaluated using a linear interpolation, whereas + outside the sample interval, the value at the closest end point is returned. + Examples for tabulated `Func1` object are:: + + >>> t1 = Func1([[0, 2], [1, 1], [2, 0]]) + >>> [t1(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]] + [2.0, 2.0, 1.5, 0.5, 0.0, 0.0] + + >>> t2 = Func1([0, 1, 2], [2, 1, 0]) + >>> [t2(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]] + [2.0, 2.0, 1.5, 0.5, 0.0, 0.0] + + >>> t3 = Func1(np.array([0, 1, 2]), (2, 1, 0)) + >>> [t3(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]] + [2.0, 2.0, 1.5, 0.5, 0.0, 0.0] + Note that all methods which accept `Func1` objects will also accept the callable object and create the wrapper on their own, so it is generally unnecessary to explicitly create a `Func1` object. """ - def __cinit__(self, c): + def __cinit__(self, *args, **kwargs): self.exception = None - if hasattr(c, '__call__'): - self.callable = c + self.callable = None + + def __init__(self, *args): + if len(args) == 1: + c = args[0] + if hasattr(c, '__call__'): + # callback function + self.__set_callback(c) + else: + arr = np.array(c) + try: + if arr.ndim == 0: + # handle constants or unsized numpy arrays + k = float(c) + self.__set_callback(lambda t: k) + elif arr.size == 1: + # handle lists, tuples or numpy arrays with a single element + k = float(c[0]) + self.__set_callback(lambda t: k) + elif arr.ndim == 2: + # tabulated function (single argument) + if arr.shape[1] == 2: + time = arr[:, 0] + fval = arr[:, 1] + self.__set_tables(time, fval) + else: + raise ValueError( + "Invalid dimensions: specification of " + "tabulated function with a single array " + "requires two columns") + else: + raise TypeError + + except TypeError: + raise TypeError( + "Func1 must be constructed from a callable object, " + "a number or a numeric array with two columns") + + elif len(args) == 2: + # tabulated function (two arguments mimic C++ interface) + time, fval = args + self.__set_tables(time, fval) + else: - try: - # calling float() converts numpy arrays of size 1 to scalars - k = float(c) - except TypeError: - if hasattr(c, '__len__') and len(c) == 1: - # Handle lists or tuples with a single element - k = float(c[0]) - else: - raise TypeError('Func1 must be constructed from a number or' - ' a callable object') - self.callable = lambda t: k + raise ValueError("Invalid number of arguments") + + cpdef void __set_callback(self, c) except *: + self.callable = c self.func = new CxxFunc1(func_callback, self) + cpdef void __set_tables(self, time, fval) except *: + cdef vector[double] tvec, fvec + for t in time: + tvec.push_back(t) + for f in fval: + fvec.push_back(f) + self.func = (new CxxTabulated1(tvec, fvec)) + def __dealloc__(self): del self.func From 66ae576dc0fc0617c0dc966fdf6cccfc6305959c Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Wed, 22 Jan 2020 14:25:07 -0600 Subject: [PATCH 12/24] [numerics] Add unit tests for tabulated Func1 --- interfaces/cython/cantera/test/test_func1.py | 41 ++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/interfaces/cython/cantera/test/test_func1.py b/interfaces/cython/cantera/test/test_func1.py index b8672d45d42..e725a22073c 100644 --- a/interfaces/cython/cantera/test/test_func1.py +++ b/interfaces/cython/cantera/test/test_func1.py @@ -70,3 +70,44 @@ def test_uncopyable(self): f = ct.Func1(np.sin) with self.assertRaises(NotImplementedError): copy.copy(f) + + def test_tabulated1(self): + arr = np.array([[0, 2], [1, 1], [2, 0]]) + time = arr[:, 0] + fval = arr[:, 1] + fcn = ct.Func1(time, fval) + for t, f in zip(time, fval): + self.assertNear(f, fcn(t)) + + def test_tabulated2(self): + time = [0, 1, 2] + fval = [2, 1, 0] + fcn = ct.Func1(time, fval) + for t, f in zip(time, fval): + self.assertNear(f, fcn(t)) + + def test_tabulated3(self): + time = np.array([0, 1, 2]) + fval = np.array([2, 1, 0]) + fcn = ct.Func1(time, fval) + tt = .5*(time[1:] + time[:-1]) + ff = .5*(fval[1:] + fval[:-1]) + for t, f in zip(tt, ff): + self.assertNear(f, fcn(t)) + + def test_tabulated4(self): + time = 0, 1, 2, + fval = 2, 1, 0, + fcn = ct.Func1(time, fval) + self.assertNear(fcn(-1), fval[0]) + self.assertNear(fcn(3), fval[-1]) + + def test_failures(self): + with self.assertRaisesRegex(ValueError, 'Invalid number of arguments'): + ct.Func1(1, 2, 3) + with self.assertRaisesRegex(ValueError, 'Invalid dimensions'): + ct.Func1(np.zeros((3, 3))) + with self.assertRaisesRegex(ct.CanteraError, 'do not match'): + ct.Func1(range(2), range(3)) + with self.assertRaisesRegex(ct.CanteraError, 'must not be empty'): + ct.Func1([], []) From cb118bcbe63c5bf1e2eb8b2463ccd1d5c3e5ef48 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Thu, 23 Jan 2020 13:56:44 -0600 Subject: [PATCH 13/24] [numerics] Use Const1 objects in Python - replace Python lambda function defining constant Func1 variants by pre-existing Const1 class defined in C++ --- interfaces/cython/cantera/_cantera.pxd | 9 +++++++-- interfaces/cython/cantera/func1.pyx | 17 ++++++++++------- 2 files changed, 17 insertions(+), 9 deletions(-) diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index c6f59d743ed..56dc7e30bff 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -43,6 +43,10 @@ cdef extern from "cantera/numerics/Func1.h": CxxTabulated1(vector[double]&, vector[double]&) except +translate_exception double eval(double) except +translate_exception + cdef cppclass CxxConst1 "Cantera::Const1": + CxxConst1(double) except +translate_exception + double eval(double) except +translate_exception + cdef extern from "cantera/base/xml.h" namespace "Cantera": cdef cppclass XML_Node: XML_Node* findByName(string) @@ -1022,8 +1026,9 @@ cdef class Func1: cdef CxxFunc1* func cdef object callable cdef object exception - cpdef void __set_callback(self, object) except * - cpdef void __set_tables(self, object, object) except * + cpdef void _set_callback(self, object) except * + cpdef void _set_const(self, double) except * + cpdef void _set_tables(self, object, object) except * cdef class ReactorBase: cdef CxxReactorBase* rbase diff --git a/interfaces/cython/cantera/func1.pyx b/interfaces/cython/cantera/func1.pyx index f70a3354a8f..a8f5e1ce15d 100644 --- a/interfaces/cython/cantera/func1.pyx +++ b/interfaces/cython/cantera/func1.pyx @@ -87,24 +87,24 @@ cdef class Func1: c = args[0] if hasattr(c, '__call__'): # callback function - self.__set_callback(c) + self._set_callback(c) else: arr = np.array(c) try: if arr.ndim == 0: # handle constants or unsized numpy arrays k = float(c) - self.__set_callback(lambda t: k) + self._set_const(k) elif arr.size == 1: # handle lists, tuples or numpy arrays with a single element k = float(c[0]) - self.__set_callback(lambda t: k) + self._set_const(k) elif arr.ndim == 2: # tabulated function (single argument) if arr.shape[1] == 2: time = arr[:, 0] fval = arr[:, 1] - self.__set_tables(time, fval) + self._set_tables(time, fval) else: raise ValueError( "Invalid dimensions: specification of " @@ -121,17 +121,20 @@ cdef class Func1: elif len(args) == 2: # tabulated function (two arguments mimic C++ interface) time, fval = args - self.__set_tables(time, fval) + self._set_tables(time, fval) else: raise ValueError("Invalid number of arguments") - cpdef void __set_callback(self, c) except *: + cpdef void _set_callback(self, c) except *: self.callable = c self.func = new CxxFunc1(func_callback, self) - cpdef void __set_tables(self, time, fval) except *: + cpdef void _set_const(self, double c) except *: + self.func = (new CxxConst1(c)) + + cpdef void _set_tables(self, time, fval) except *: cdef vector[double] tvec, fvec for t in time: tvec.push_back(t) From fc32262a47060c36580b5508ff06648907470548 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sat, 25 Jan 2020 16:25:52 -0600 Subject: [PATCH 14/24] [numerics] Add previous value option to tabulated Func1 objects --- include/cantera/numerics/Func1.h | 19 ++++++++++++-- interfaces/cython/cantera/_cantera.pxd | 4 +-- interfaces/cython/cantera/func1.pyx | 29 ++++++++++++++------- src/numerics/Func1.cpp | 36 ++++++++++++++++---------- 4 files changed, 61 insertions(+), 27 deletions(-) diff --git a/include/cantera/numerics/Func1.h b/include/cantera/numerics/Func1.h index 2af040f2d93..e738aa31511 100644 --- a/include/cantera/numerics/Func1.h +++ b/include/cantera/numerics/Func1.h @@ -279,7 +279,8 @@ class Pow1 : public Func1 class Tabulated1 : public Func1 { public: - Tabulated1(const std::vector& tvec, const std::vector& fvec) : + Tabulated1(const std::vector& tvec, const std::vector& fvec, + const std::string& method = "linear") : Func1() { if (tvec.size() != fvec.size()) { throw CanteraError("Tabulated1::Tabulated1", @@ -291,6 +292,15 @@ class Tabulated1 : public Func1 } m_tvec = tvec; m_fvec = fvec; + if (method == "linear") { + m_isLinear = true; + } else if (method == "previous") { + m_isLinear = false; + } else { + throw CanteraError("Tabulated1::Tabulated1", + "interpolation method '{}' is not implemented", + method); + } } Tabulated1(const Tabulated1& b) : @@ -311,13 +321,18 @@ class Tabulated1 : public Func1 } virtual double eval(double t) const; virtual Func1& duplicate() const { - return *(new Tabulated1(m_tvec, m_fvec)); + if (m_isLinear) { + return *(new Tabulated1(m_tvec, m_fvec, "linear")); + } else { + return *(new Tabulated1(m_tvec, m_fvec, "previous")); + } } virtual Func1& derivative() const; private: std::vector m_tvec; std::vector m_fvec; + bool m_isLinear; }; diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 56dc7e30bff..48b6761c201 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -40,7 +40,7 @@ cdef extern from "cantera/cython/funcWrapper.h": cdef extern from "cantera/numerics/Func1.h": cdef cppclass CxxTabulated1 "Cantera::Tabulated1": - CxxTabulated1(vector[double]&, vector[double]&) except +translate_exception + CxxTabulated1(vector[double]&, vector[double]&, string) except +translate_exception double eval(double) except +translate_exception cdef cppclass CxxConst1 "Cantera::Const1": @@ -1028,7 +1028,7 @@ cdef class Func1: cdef object exception cpdef void _set_callback(self, object) except * cpdef void _set_const(self, double) except * - cpdef void _set_tables(self, object, object) except * + cpdef void _set_tables(self, object, object, string) except * cdef class ReactorBase: cdef CxxReactorBase* rbase diff --git a/interfaces/cython/cantera/func1.pyx b/interfaces/cython/cantera/func1.pyx index a8f5e1ce15d..9f233f2eae7 100644 --- a/interfaces/cython/cantera/func1.pyx +++ b/interfaces/cython/cantera/func1.pyx @@ -58,22 +58,29 @@ cdef class Func1: points and corresponding function values. Inputs are specified either by two iterable objects containing sample point location and function values, or a single array that concatenates those inputs in two columns. Between - sample points, values are evaluated using a linear interpolation, whereas - outside the sample interval, the value at the closest end point is returned. + sample points, values are evaluated based on the keyword argument + 'interpolation'; options are 'linear' (linear interpolation, default) or + 'previous' (nearest previous value). Outside the sample interval, the value + at the closest end point is returned. + Examples for tabulated `Func1` object are:: >>> t1 = Func1([[0, 2], [1, 1], [2, 0]]) >>> [t1(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]] [2.0, 2.0, 1.5, 0.5, 0.0, 0.0] - >>> t2 = Func1([0, 1, 2], [2, 1, 0]) + >>> t2 = Func1([[0, 2], [1, 1], [2, 0]], interpolation='previous') >>> [t2(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]] - [2.0, 2.0, 1.5, 0.5, 0.0, 0.0] + [2.0, 2.0, 2.0, 1.0, 0.0, 0.0] - >>> t3 = Func1(np.array([0, 1, 2]), (2, 1, 0)) + >>> t3 = Func1([0, 1, 2], [2, 1, 0]) >>> [t3(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]] [2.0, 2.0, 1.5, 0.5, 0.0, 0.0] + >>> t4 = Func1(np.array([0, 1, 2]), (2, 1, 0)) + >>> [t4(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]] + [2.0, 2.0, 1.5, 0.5, 0.0, 0.0] + Note that all methods which accept `Func1` objects will also accept the callable object and create the wrapper on their own, so it is generally unnecessary to explicitly create a `Func1` object. @@ -82,7 +89,7 @@ cdef class Func1: self.exception = None self.callable = None - def __init__(self, *args): + def __init__(self, *args, **kwargs): if len(args) == 1: c = args[0] if hasattr(c, '__call__'): @@ -104,7 +111,8 @@ cdef class Func1: if arr.shape[1] == 2: time = arr[:, 0] fval = arr[:, 1] - self._set_tables(time, fval) + method = kwargs.get('interpolation', 'linear') + self._set_tables(time, fval, stringify(method)) else: raise ValueError( "Invalid dimensions: specification of " @@ -121,7 +129,8 @@ cdef class Func1: elif len(args) == 2: # tabulated function (two arguments mimic C++ interface) time, fval = args - self._set_tables(time, fval) + method = kwargs.get('interpolation', 'linear') + self._set_tables(time, fval, stringify(method)) else: raise ValueError("Invalid number of arguments") @@ -134,13 +143,13 @@ cdef class Func1: cpdef void _set_const(self, double c) except *: self.func = (new CxxConst1(c)) - cpdef void _set_tables(self, time, fval) except *: + cpdef void _set_tables(self, time, fval, string method) except *: cdef vector[double] tvec, fvec for t in time: tvec.push_back(t) for f in fval: fvec.push_back(f) - self.func = (new CxxTabulated1(tvec, fvec)) + self.func = (new CxxTabulated1(tvec, fvec, method)) def __dealloc__(self): del self.func diff --git a/src/numerics/Func1.cpp b/src/numerics/Func1.cpp index 7feeb10fa2d..20530f5117a 100644 --- a/src/numerics/Func1.cpp +++ b/src/numerics/Func1.cpp @@ -227,10 +227,14 @@ double Tabulated1::eval(double t) const { while (t > m_tvec[ix+1]) { ix++; } - double df = m_fvec[ix+1] - m_fvec[ix]; - df /= m_tvec[ix+1] - m_tvec[ix]; - df *= t - m_tvec[ix]; - return m_fvec[ix] + df; + if (m_isLinear) { + double df = m_fvec[ix+1] - m_fvec[ix]; + df /= m_tvec[ix+1] - m_tvec[ix]; + df *= t - m_tvec[ix]; + return m_fvec[ix] + df; + } else { + return m_fvec[ix]; + } } } else { return 0.; @@ -241,20 +245,26 @@ Func1& Tabulated1::derivative() const { std::vector tvec; std::vector dvec; size_t siz = m_tvec.size(); - if (siz>1) { - for (size_t i=1; i1) { + for (size_t i=1; i Date: Wed, 1 Apr 2020 21:26:08 -0500 Subject: [PATCH 15/24] [numerics] Streamline and update docstrings for Tabulated1 --- include/cantera/numerics/Func1.h | 61 ++++++-------------- interfaces/cython/cantera/test/test_func1.py | 6 ++ src/numerics/Func1.cpp | 34 ++++++++++- 3 files changed, 56 insertions(+), 45 deletions(-) diff --git a/include/cantera/numerics/Func1.h b/include/cantera/numerics/Func1.h index e738aa31511..ff6036338b7 100644 --- a/include/cantera/numerics/Func1.h +++ b/include/cantera/numerics/Func1.h @@ -275,45 +275,18 @@ class Pow1 : public Func1 }; -//! implements a tabulated function +//! The Tabulated1 class implements a tabulated function class Tabulated1 : public Func1 { public: - Tabulated1(const std::vector& tvec, const std::vector& fvec, - const std::string& method = "linear") : - Func1() { - if (tvec.size() != fvec.size()) { - throw CanteraError("Tabulated1::Tabulated1", - "sizes of vectors do not match ({} vs {}).", - tvec.size(), fvec.size()); - } else if (tvec.empty()) { - throw CanteraError("Tabulated1::Tabulated1", - "vectors must not be empty."); - } - m_tvec = tvec; - m_fvec = fvec; - if (method == "linear") { - m_isLinear = true; - } else if (method == "previous") { - m_isLinear = false; - } else { - throw CanteraError("Tabulated1::Tabulated1", - "interpolation method '{}' is not implemented", - method); - } - } - - Tabulated1(const Tabulated1& b) : - Func1(b) { - } - - Tabulated1& operator=(const Tabulated1& right) { - if (&right == this) { - return *this; - } - Func1::operator=(right); - return *this; - } + //! Constructor. + /*! + * @param tvec Vector of time values + * @param fvec Vector of values + * @param method Interpolation method ('linear' or 'previous') + */ + Tabulated1(const vector_fp& tvec, const vector_fp& fvec, + const std::string& method = "linear"); virtual std::string write(const std::string& arg) const; virtual int ID() const { @@ -330,19 +303,21 @@ class Tabulated1 : public Func1 virtual Func1& derivative() const; private: - std::vector m_tvec; - std::vector m_fvec; - bool m_isLinear; + vector_fp m_tvec; //!< Vector of time values + vector_fp m_fvec; //!< Vector of function values + bool m_isLinear; //!< Boolean indicating interpolation method }; -/** - * Constant. - */ +//! The Const1 class implements a constant class Const1 : public Func1 { public: - Const1(doublereal A) : + //! Constructor. + /*! + * @param A Constant + */ + Const1(double A) : Func1() { m_c = A; } diff --git a/interfaces/cython/cantera/test/test_func1.py b/interfaces/cython/cantera/test/test_func1.py index e725a22073c..4b30f687755 100644 --- a/interfaces/cython/cantera/test/test_func1.py +++ b/interfaces/cython/cantera/test/test_func1.py @@ -102,6 +102,12 @@ def test_tabulated4(self): self.assertNear(fcn(-1), fval[0]) self.assertNear(fcn(3), fval[-1]) + def test_tabulated5(self): + time = 0, 1, 0.5, 2 + fval = 2, 1, 1, 0 + with self.assertRaisesRegex(ct.CanteraError, 'monotonically'): + ct.Func1(time, fval) + def test_failures(self): with self.assertRaisesRegex(ValueError, 'Invalid number of arguments'): ct.Func1(1, 2, 3) diff --git a/src/numerics/Func1.cpp b/src/numerics/Func1.cpp index 20530f5117a..3cf0bf9150d 100644 --- a/src/numerics/Func1.cpp +++ b/src/numerics/Func1.cpp @@ -215,6 +215,36 @@ Func1& Pow1::derivative() const /******************************************************************************/ +Tabulated1::Tabulated1(const vector_fp& tvec, const vector_fp& fvec, + const std::string& method) : + Func1() { + if (tvec.size() != fvec.size()) { + throw CanteraError("Tabulated1::Tabulated1", + "sizes of vectors do not match ({} vs {}).", + tvec.size(), fvec.size()); + } else if (tvec.empty()) { + throw CanteraError("Tabulated1::Tabulated1", + "vectors must not be empty."); + } + for (auto it = std::begin(tvec) + 1; it != std::end(tvec); it++) { + if (*(it - 1) > *it) { + throw CanteraError("Tabulated1::Tabulated1", + "time values are not increasing monotonically."); + } + } + m_tvec = tvec; + m_fvec = fvec; + if (method == "linear") { + m_isLinear = true; + } else if (method == "previous") { + m_isLinear = false; + } else { + throw CanteraError("Tabulated1::Tabulated1", + "interpolation method '{}' is not implemented", + method); + } +} + double Tabulated1::eval(double t) const { size_t siz = m_tvec.size(); if (siz) { @@ -242,8 +272,8 @@ double Tabulated1::eval(double t) const { } Func1& Tabulated1::derivative() const { - std::vector tvec; - std::vector dvec; + vector_fp tvec; + vector_fp dvec; size_t siz = m_tvec.size(); if (m_isLinear) { // piece-wise continuous derivative From 047ae463ceb09a709f39b4c89a9e1aa3d3fb5964 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Wed, 1 Apr 2020 23:07:06 -0500 Subject: [PATCH 16/24] [numerics] Improve unit tests for Tabulated1 --- interfaces/cython/cantera/test/test_func1.py | 14 +++++--- src/numerics/Func1.cpp | 35 +++++++++----------- 2 files changed, 26 insertions(+), 23 deletions(-) diff --git a/interfaces/cython/cantera/test/test_func1.py b/interfaces/cython/cantera/test/test_func1.py index 4b30f687755..101055c5ba4 100644 --- a/interfaces/cython/cantera/test/test_func1.py +++ b/interfaces/cython/cantera/test/test_func1.py @@ -103,10 +103,9 @@ def test_tabulated4(self): self.assertNear(fcn(3), fval[-1]) def test_tabulated5(self): - time = 0, 1, 0.5, 2 - fval = 2, 1, 1, 0 - with self.assertRaisesRegex(ct.CanteraError, 'monotonically'): - ct.Func1(time, fval) + fcn = ct.Func1([[0, 2], [1, 1], [2, 0]], interpolation='previous') + val = np.array([fcn(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]]) + self.assertArrayNear(val, np.array([2.0, 2.0, 2.0, 1.0, 0.0, 0.0])) def test_failures(self): with self.assertRaisesRegex(ValueError, 'Invalid number of arguments'): @@ -117,3 +116,10 @@ def test_failures(self): ct.Func1(range(2), range(3)) with self.assertRaisesRegex(ct.CanteraError, 'must not be empty'): ct.Func1([], []) + with self.assertRaisesRegex(ct.CanteraError, 'monotonically'): + ct.Func1((0, 1, 0.5, 2), (2, 1, 1, 0)) + with self.assertRaisesRegex(ct.CanteraError, 'not implemented'): + ct.Func1((0, 1, 1, 2), (2, 1, 1, 0), interpolation='quadratic') + with self.assertRaisesRegex(ct.CanteraError, 'not be empty'): + fcn = ct.Func1([], []) + diff --git a/src/numerics/Func1.cpp b/src/numerics/Func1.cpp index 3cf0bf9150d..3ffea33eec5 100644 --- a/src/numerics/Func1.cpp +++ b/src/numerics/Func1.cpp @@ -247,27 +247,24 @@ Tabulated1::Tabulated1(const vector_fp& tvec, const vector_fp& fvec, double Tabulated1::eval(double t) const { size_t siz = m_tvec.size(); - if (siz) { - if (t <= m_tvec[0]) { - return m_fvec[0]; - } else if (t >= m_tvec[siz-1]) { - return m_fvec[siz-1]; + // constructor ensures that siz > 0 + if (t <= m_tvec[0]) { + return m_fvec[0]; + } else if (t >= m_tvec[siz-1]) { + return m_fvec[siz-1]; + } else { + size_t ix = 0; + while (t > m_tvec[ix+1]) { + ix++; + } + if (m_isLinear) { + double df = m_fvec[ix+1] - m_fvec[ix]; + df /= m_tvec[ix+1] - m_tvec[ix]; + df *= t - m_tvec[ix]; + return m_fvec[ix] + df; } else { - size_t ix = 0; - while (t > m_tvec[ix+1]) { - ix++; - } - if (m_isLinear) { - double df = m_fvec[ix+1] - m_fvec[ix]; - df /= m_tvec[ix+1] - m_tvec[ix]; - df *= t - m_tvec[ix]; - return m_fvec[ix] + df; - } else { - return m_fvec[ix]; - } + return m_fvec[ix]; } - } else { - return 0.; } } From 6c9c778962dcff309cc19a9a916ffea42b927e6c Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Thu, 2 Apr 2020 21:10:01 -0500 Subject: [PATCH 17/24] [numerics] Create Tabulated1 from raw pointers --- include/cantera/numerics/Func1.h | 13 ++++++++----- interfaces/cython/cantera/_cantera.pxd | 2 +- interfaces/cython/cantera/func1.pyx | 19 +++++++++++++------ interfaces/cython/cantera/test/test_func1.py | 9 +++------ src/numerics/Func1.cpp | 20 +++++++------------- 5 files changed, 32 insertions(+), 31 deletions(-) diff --git a/include/cantera/numerics/Func1.h b/include/cantera/numerics/Func1.h index ff6036338b7..a2a8dd7b212 100644 --- a/include/cantera/numerics/Func1.h +++ b/include/cantera/numerics/Func1.h @@ -281,11 +281,12 @@ class Tabulated1 : public Func1 public: //! Constructor. /*! - * @param tvec Vector of time values - * @param fvec Vector of values + * @param n Size of tabulated value arrays + * @param tvals Pointer to time value array + * @param fvals Pointer to function value array * @param method Interpolation method ('linear' or 'previous') */ - Tabulated1(const vector_fp& tvec, const vector_fp& fvec, + Tabulated1(size_t n, const double* tvals, const double* fvals, const std::string& method = "linear"); virtual std::string write(const std::string& arg) const; @@ -295,9 +296,11 @@ class Tabulated1 : public Func1 virtual double eval(double t) const; virtual Func1& duplicate() const { if (m_isLinear) { - return *(new Tabulated1(m_tvec, m_fvec, "linear")); + return *(new Tabulated1(m_tvec.size(), &m_tvec[0], &m_fvec[0], + "linear")); } else { - return *(new Tabulated1(m_tvec, m_fvec, "previous")); + return *(new Tabulated1(m_tvec.size(), &m_tvec[0], &m_fvec[0], + "previous")); } } diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 48b6761c201..685d24a6e94 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -40,7 +40,7 @@ cdef extern from "cantera/cython/funcWrapper.h": cdef extern from "cantera/numerics/Func1.h": cdef cppclass CxxTabulated1 "Cantera::Tabulated1": - CxxTabulated1(vector[double]&, vector[double]&, string) except +translate_exception + CxxTabulated1(int, double*, double*, string) except +translate_exception double eval(double) except +translate_exception cdef cppclass CxxConst1 "Cantera::Const1": diff --git a/interfaces/cython/cantera/func1.pyx b/interfaces/cython/cantera/func1.pyx index 9f233f2eae7..cc1a92b7b3b 100644 --- a/interfaces/cython/cantera/func1.pyx +++ b/interfaces/cython/cantera/func1.pyx @@ -144,12 +144,19 @@ cdef class Func1: self.func = (new CxxConst1(c)) cpdef void _set_tables(self, time, fval, string method) except *: - cdef vector[double] tvec, fvec - for t in time: - tvec.push_back(t) - for f in fval: - fvec.push_back(f) - self.func = (new CxxTabulated1(tvec, fvec, method)) + tsiz = np.asarray(time).size + fsiz = np.asarray(fval).size + if tsiz != fsiz: + raise ValueError("Sizes of arrays do not match " + "({} vs {}).".format(tsiz, fsiz)) + elif tsiz == 0: + raise ValueError("Arrays must not be empty.") + cdef np.ndarray[np.double_t, ndim=1] tvec = np.asarray(time, + dtype=np.double) + cdef np.ndarray[np.double_t, ndim=1] fvec = np.asarray(fval, + dtype=np.double) + self.func = (new CxxTabulated1(tsiz, &tvec[0], &fvec[0], + method)) def __dealloc__(self): del self.func diff --git a/interfaces/cython/cantera/test/test_func1.py b/interfaces/cython/cantera/test/test_func1.py index 101055c5ba4..8047c6e17ff 100644 --- a/interfaces/cython/cantera/test/test_func1.py +++ b/interfaces/cython/cantera/test/test_func1.py @@ -107,19 +107,16 @@ def test_tabulated5(self): val = np.array([fcn(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]]) self.assertArrayNear(val, np.array([2.0, 2.0, 2.0, 1.0, 0.0, 0.0])) - def test_failures(self): + def test_tabulated_failures(self): with self.assertRaisesRegex(ValueError, 'Invalid number of arguments'): ct.Func1(1, 2, 3) with self.assertRaisesRegex(ValueError, 'Invalid dimensions'): ct.Func1(np.zeros((3, 3))) - with self.assertRaisesRegex(ct.CanteraError, 'do not match'): + with self.assertRaisesRegex(ValueError, 'do not match'): ct.Func1(range(2), range(3)) - with self.assertRaisesRegex(ct.CanteraError, 'must not be empty'): + with self.assertRaisesRegex(ValueError, 'must not be empty'): ct.Func1([], []) with self.assertRaisesRegex(ct.CanteraError, 'monotonically'): ct.Func1((0, 1, 0.5, 2), (2, 1, 1, 0)) with self.assertRaisesRegex(ct.CanteraError, 'not implemented'): ct.Func1((0, 1, 1, 2), (2, 1, 1, 0), interpolation='quadratic') - with self.assertRaisesRegex(ct.CanteraError, 'not be empty'): - fcn = ct.Func1([], []) - diff --git a/src/numerics/Func1.cpp b/src/numerics/Func1.cpp index 3ffea33eec5..d4f189865fa 100644 --- a/src/numerics/Func1.cpp +++ b/src/numerics/Func1.cpp @@ -215,25 +215,19 @@ Func1& Pow1::derivative() const /******************************************************************************/ -Tabulated1::Tabulated1(const vector_fp& tvec, const vector_fp& fvec, +Tabulated1::Tabulated1(size_t n, const double* tvals, const double* fvals, const std::string& method) : Func1() { - if (tvec.size() != fvec.size()) { - throw CanteraError("Tabulated1::Tabulated1", - "sizes of vectors do not match ({} vs {}).", - tvec.size(), fvec.size()); - } else if (tvec.empty()) { - throw CanteraError("Tabulated1::Tabulated1", - "vectors must not be empty."); - } - for (auto it = std::begin(tvec) + 1; it != std::end(tvec); it++) { + m_tvec.resize(n); + std::copy(tvals, tvals + n, m_tvec.begin()); + for (auto it = std::begin(m_tvec) + 1; it != std::end(m_tvec); it++) { if (*(it - 1) > *it) { throw CanteraError("Tabulated1::Tabulated1", "time values are not increasing monotonically."); } } - m_tvec = tvec; - m_fvec = fvec; + m_fvec.resize(n); + std::copy(fvals, fvals + n, m_fvec.begin()); if (method == "linear") { m_isLinear = true; } else if (method == "previous") { @@ -291,7 +285,7 @@ Func1& Tabulated1::derivative() const { dvec.push_back(0.); dvec.push_back(0.); } - return *(new Tabulated1(tvec, dvec, "previous")); + return *(new Tabulated1(tvec.size(), &tvec[0], &dvec[0], "previous")); } /******************************************************************************/ From 6ed86501c17ae390a40092ab67fd5ace8a614d86 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Thu, 2 Apr 2020 22:12:20 -0500 Subject: [PATCH 18/24] [numerics] Introduce dedicated TabulatedFunction class in Python --- interfaces/cython/cantera/_cantera.pxd | 7 +- interfaces/cython/cantera/func1.pyx | 182 ++++++++++--------- interfaces/cython/cantera/test/test_func1.py | 22 +-- 3 files changed, 107 insertions(+), 104 deletions(-) diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 685d24a6e94..aa3dd6a263d 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -43,10 +43,6 @@ cdef extern from "cantera/numerics/Func1.h": CxxTabulated1(int, double*, double*, string) except +translate_exception double eval(double) except +translate_exception - cdef cppclass CxxConst1 "Cantera::Const1": - CxxConst1(double) except +translate_exception - double eval(double) except +translate_exception - cdef extern from "cantera/base/xml.h" namespace "Cantera": cdef cppclass XML_Node: XML_Node* findByName(string) @@ -1027,7 +1023,8 @@ cdef class Func1: cdef object callable cdef object exception cpdef void _set_callback(self, object) except * - cpdef void _set_const(self, double) except * + +cdef class TabulatedFunction(Func1): cpdef void _set_tables(self, object, object, string) except * cdef class ReactorBase: diff --git a/interfaces/cython/cantera/func1.pyx b/interfaces/cython/cantera/func1.pyx index cc1a92b7b3b..c0eeaf34454 100644 --- a/interfaces/cython/cantera/func1.pyx +++ b/interfaces/cython/cantera/func1.pyx @@ -2,7 +2,6 @@ # at https://cantera.org/license.txt for license and copyright information. import sys -from numbers import Real as _real cdef double func_callback(double t, void* obj, void** err): """ @@ -54,118 +53,125 @@ cdef class Func1: >>> f4(0.1) 2.5 - A `Func1` object representing a tabulated function is defined by sample - points and corresponding function values. Inputs are specified either by - two iterable objects containing sample point location and function values, - or a single array that concatenates those inputs in two columns. Between - sample points, values are evaluated based on the keyword argument - 'interpolation'; options are 'linear' (linear interpolation, default) or - 'previous' (nearest previous value). Outside the sample interval, the value - at the closest end point is returned. + Note that all methods which accept `Func1` objects will also accept the + callable object and create the wrapper on their own, so it is often not + necessary to explicitly create a `Func1` object. + """ + def __cinit__(self, *args, **kwargs): + self.exception = None + self.callable = None + + def __init__(self, c): + if hasattr(c, '__call__'): + # callback function + self._set_callback(c) + else: + arr = np.array(c) + try: + if arr.ndim == 0: + # handle constants or unsized numpy arrays + k = float(c) + self._set_callback(lambda t: k) + elif arr.size == 1: + # handle lists, tuples or numpy arrays with a single element + k = float(c[0]) + self._set_callback(lambda t: k) + else: + raise TypeError + + except TypeError: + raise TypeError( + "'Func1' objects must be constructed from a number or " + "a callable object") from None + + cpdef void _set_callback(self, c) except *: + self.callable = c + self.func = new CxxFunc1(func_callback, self) - Examples for tabulated `Func1` object are:: + def __dealloc__(self): + del self.func + + def __call__(self, t): + return self.func.eval(t) - >>> t1 = Func1([[0, 2], [1, 1], [2, 0]]) + def __reduce__(self): + msg = "'{}' objects are not picklable".format(type(self).__name__) + raise NotImplementedError(msg) + + def __copy__(self): + msg = "'{}' objects are not copyable".format(type(self).__name__) + raise NotImplementedError(msg) + + +cdef class TabulatedFunction(Func1): + """ + A `TabulatedFunction` object representing a tabulated function is defined by + sample points and corresponding function values. Inputs are specified either + by two iterable objects containing sample point location and function + values, or a single array that concatenates those inputs in two rows or + columns. Between sample points, values are evaluated based on the optional + keyword argument 'method'; options are 'linear' (linear interpolation, + default) or 'previous' (nearest previous value). Outside the sample + interval, the value at the closest end point is returned. + + Examples for `TabulatedFunction` objects are:: + + >>> t1 = TabulatedFunction([[0, 2], [1, 1], [2, 0]]) >>> [t1(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]] [2.0, 2.0, 1.5, 0.5, 0.0, 0.0] - >>> t2 = Func1([[0, 2], [1, 1], [2, 0]], interpolation='previous') + >>> t2 = TabulatedFunction([[0, 2], [1, 1], [2, 0]], method='previous') >>> [t2(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]] [2.0, 2.0, 2.0, 1.0, 0.0, 0.0] - >>> t3 = Func1([0, 1, 2], [2, 1, 0]) + >>> t3 = TabulatedFunction([0, 1, 2], [2, 1, 0]) >>> [t3(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]] [2.0, 2.0, 1.5, 0.5, 0.0, 0.0] - >>> t4 = Func1(np.array([0, 1, 2]), (2, 1, 0)) + >>> t4 = TabulatedFunction(np.array([0, 1, 2]), (2, 1, 0)) >>> [t4(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]] [2.0, 2.0, 1.5, 0.5, 0.0, 0.0] - - Note that all methods which accept `Func1` objects will also accept the - callable object and create the wrapper on their own, so it is generally - unnecessary to explicitly create a `Func1` object. """ - def __cinit__(self, *args, **kwargs): - self.exception = None - self.callable = None - def __init__(self, *args, **kwargs): + def __init__(self, *args, method='linear'): if len(args) == 1: - c = args[0] - if hasattr(c, '__call__'): - # callback function - self._set_callback(c) + # tabulated function (single argument) + arr = np.array(args[0]) + if arr.ndim == 2: + if arr.shape[1] == 2: + time = arr[:, 0] + fval = arr[:, 1] + elif arr.shape[0] == 2: + time = arr[0, :] + fval = arr[1, :] + else: + raise ValueError("Invalid dimensions: specification of " + "tabulated function with a single array " + "requires two rows or columns") + self._set_tables(time, fval, stringify(method)) else: - arr = np.array(c) - try: - if arr.ndim == 0: - # handle constants or unsized numpy arrays - k = float(c) - self._set_const(k) - elif arr.size == 1: - # handle lists, tuples or numpy arrays with a single element - k = float(c[0]) - self._set_const(k) - elif arr.ndim == 2: - # tabulated function (single argument) - if arr.shape[1] == 2: - time = arr[:, 0] - fval = arr[:, 1] - method = kwargs.get('interpolation', 'linear') - self._set_tables(time, fval, stringify(method)) - else: - raise ValueError( - "Invalid dimensions: specification of " - "tabulated function with a single array " - "requires two columns") - else: - raise TypeError - - except TypeError: - raise TypeError( - "Func1 must be constructed from a callable object, " - "a number or a numeric array with two columns") + raise TypeError( + "'TabulatedFunction' must be constructed from " + "a numeric array with two columns") elif len(args) == 2: - # tabulated function (two arguments mimic C++ interface) + # tabulated function (two arrays mimic C++ interface) time, fval = args - method = kwargs.get('interpolation', 'linear') self._set_tables(time, fval, stringify(method)) else: - raise ValueError("Invalid number of arguments") - - - cpdef void _set_callback(self, c) except *: - self.callable = c - self.func = new CxxFunc1(func_callback, self) - - cpdef void _set_const(self, double c) except *: - self.func = (new CxxConst1(c)) + raise ValueError("Invalid number of arguments (one or two " + "arguments containing tabulated values)") cpdef void _set_tables(self, time, fval, string method) except *: - tsiz = np.asarray(time).size - fsiz = np.asarray(fval).size - if tsiz != fsiz: + tt = np.asarray(time, dtype=np.double) + ff = np.asarray(fval, dtype=np.double) + if tt.size != ff.size: raise ValueError("Sizes of arrays do not match " - "({} vs {}).".format(tsiz, fsiz)) - elif tsiz == 0: + "({} vs {})".format(tt.size, ff.size)) + elif tt.size == 0: raise ValueError("Arrays must not be empty.") - cdef np.ndarray[np.double_t, ndim=1] tvec = np.asarray(time, - dtype=np.double) - cdef np.ndarray[np.double_t, ndim=1] fvec = np.asarray(fval, - dtype=np.double) - self.func = (new CxxTabulated1(tsiz, &tvec[0], &fvec[0], + cdef np.ndarray[np.double_t, ndim=1] tvec = tt + cdef np.ndarray[np.double_t, ndim=1] fvec = ff + self.func = (new CxxTabulated1(tt.size, &tvec[0], &fvec[0], method)) - - def __dealloc__(self): - del self.func - - def __call__(self, t): - return self.func.eval(t) - - def __reduce__(self): - raise NotImplementedError('Func1 object is not picklable') - - def __copy__(self): - raise NotImplementedError('Func1 object is not copyable') diff --git a/interfaces/cython/cantera/test/test_func1.py b/interfaces/cython/cantera/test/test_func1.py index 8047c6e17ff..3c918453d88 100644 --- a/interfaces/cython/cantera/test/test_func1.py +++ b/interfaces/cython/cantera/test/test_func1.py @@ -75,21 +75,21 @@ def test_tabulated1(self): arr = np.array([[0, 2], [1, 1], [2, 0]]) time = arr[:, 0] fval = arr[:, 1] - fcn = ct.Func1(time, fval) + fcn = ct.TabulatedFunction(time, fval) for t, f in zip(time, fval): self.assertNear(f, fcn(t)) def test_tabulated2(self): time = [0, 1, 2] fval = [2, 1, 0] - fcn = ct.Func1(time, fval) + fcn = ct.TabulatedFunction(time, fval) for t, f in zip(time, fval): self.assertNear(f, fcn(t)) def test_tabulated3(self): time = np.array([0, 1, 2]) fval = np.array([2, 1, 0]) - fcn = ct.Func1(time, fval) + fcn = ct.TabulatedFunction(time, fval) tt = .5*(time[1:] + time[:-1]) ff = .5*(fval[1:] + fval[:-1]) for t, f in zip(tt, ff): @@ -98,25 +98,25 @@ def test_tabulated3(self): def test_tabulated4(self): time = 0, 1, 2, fval = 2, 1, 0, - fcn = ct.Func1(time, fval) + fcn = ct.TabulatedFunction(time, fval) self.assertNear(fcn(-1), fval[0]) self.assertNear(fcn(3), fval[-1]) def test_tabulated5(self): - fcn = ct.Func1([[0, 2], [1, 1], [2, 0]], interpolation='previous') + fcn = ct.TabulatedFunction([[0, 2], [1, 1], [2, 0]], method='previous') val = np.array([fcn(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]]) self.assertArrayNear(val, np.array([2.0, 2.0, 2.0, 1.0, 0.0, 0.0])) def test_tabulated_failures(self): with self.assertRaisesRegex(ValueError, 'Invalid number of arguments'): - ct.Func1(1, 2, 3) + ct.TabulatedFunction(1, 2, 3) with self.assertRaisesRegex(ValueError, 'Invalid dimensions'): - ct.Func1(np.zeros((3, 3))) + ct.TabulatedFunction(np.zeros((3, 3))) with self.assertRaisesRegex(ValueError, 'do not match'): - ct.Func1(range(2), range(3)) + ct.TabulatedFunction(range(2), range(3)) with self.assertRaisesRegex(ValueError, 'must not be empty'): - ct.Func1([], []) + ct.TabulatedFunction([], []) with self.assertRaisesRegex(ct.CanteraError, 'monotonically'): - ct.Func1((0, 1, 0.5, 2), (2, 1, 1, 0)) + ct.TabulatedFunction((0, 1, 0.5, 2), (2, 1, 1, 0)) with self.assertRaisesRegex(ct.CanteraError, 'not implemented'): - ct.Func1((0, 1, 1, 2), (2, 1, 1, 0), interpolation='quadratic') + ct.TabulatedFunction((0, 1, 1, 2), (2, 1, 1, 0), method='quadratic') From 95210eaf07b28141a94e9bbafa9389428c1262ec Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Fri, 3 Apr 2020 16:42:58 -0500 Subject: [PATCH 19/24] [numerics] Improve documentation for TabulatedFunction --- interfaces/cython/cantera/func1.pyx | 30 ++++++++++++-------- interfaces/cython/cantera/test/test_func1.py | 2 +- 2 files changed, 19 insertions(+), 13 deletions(-) diff --git a/interfaces/cython/cantera/func1.pyx b/interfaces/cython/cantera/func1.pyx index c0eeaf34454..9b1d8ab8a58 100644 --- a/interfaces/cython/cantera/func1.pyx +++ b/interfaces/cython/cantera/func1.pyx @@ -110,25 +110,32 @@ cdef class TabulatedFunction(Func1): by two iterable objects containing sample point location and function values, or a single array that concatenates those inputs in two rows or columns. Between sample points, values are evaluated based on the optional - keyword argument 'method'; options are 'linear' (linear interpolation, - default) or 'previous' (nearest previous value). Outside the sample - interval, the value at the closest end point is returned. + argument ``method``, which has to be supplied as a keyword; options are + ``'linear'`` (linear interpolation, default) or ``'previous'`` (nearest + previous value). Outside the sample interval, the value at the closest end + point is returned. - Examples for `TabulatedFunction` objects are:: + Examples for `TabulatedFunction` objects using a single (two-dimensional) + array as input are:: >>> t1 = TabulatedFunction([[0, 2], [1, 1], [2, 0]]) >>> [t1(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]] [2.0, 2.0, 1.5, 0.5, 0.0, 0.0] - >>> t2 = TabulatedFunction([[0, 2], [1, 1], [2, 0]], method='previous') + >>> t2 = TabulatedFunction(np.array([0, 1, 2]), (2, 1, 0)) >>> [t2(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]] - [2.0, 2.0, 2.0, 1.0, 0.0, 0.0] + [2.0, 2.0, 1.5, 0.5, 0.0, 0.0] - >>> t3 = TabulatedFunction([0, 1, 2], [2, 1, 0]) + where the optional ``method`` keyword argument changes the type of + interpolation from the ``'linear'`` default to ``'previous'``:: + + >>> t3 = TabulatedFunction([[0, 2], [1, 1], [2, 0]], method='previous') >>> [t3(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]] - [2.0, 2.0, 1.5, 0.5, 0.0, 0.0] + [2.0, 2.0, 2.0, 1.0, 0.0, 0.0] + + Alternatively, a `TabulatedFunction` can be defined using two input arrays:: - >>> t4 = TabulatedFunction(np.array([0, 1, 2]), (2, 1, 0)) + >>> t4 = TabulatedFunction([0, 1, 2], [2, 1, 0]) >>> [t4(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]] [2.0, 2.0, 1.5, 0.5, 0.0, 0.0] """ @@ -150,9 +157,8 @@ cdef class TabulatedFunction(Func1): "requires two rows or columns") self._set_tables(time, fval, stringify(method)) else: - raise TypeError( - "'TabulatedFunction' must be constructed from " - "a numeric array with two columns") + raise TypeError("'TabulatedFunction' must be constructed from " + "a numeric array with two dimensions") elif len(args) == 2: # tabulated function (two arrays mimic C++ interface) diff --git a/interfaces/cython/cantera/test/test_func1.py b/interfaces/cython/cantera/test/test_func1.py index 3c918453d88..a7d1b11fe96 100644 --- a/interfaces/cython/cantera/test/test_func1.py +++ b/interfaces/cython/cantera/test/test_func1.py @@ -119,4 +119,4 @@ def test_tabulated_failures(self): with self.assertRaisesRegex(ct.CanteraError, 'monotonically'): ct.TabulatedFunction((0, 1, 0.5, 2), (2, 1, 1, 0)) with self.assertRaisesRegex(ct.CanteraError, 'not implemented'): - ct.TabulatedFunction((0, 1, 1, 2), (2, 1, 1, 0), method='quadratic') + ct.TabulatedFunction((0, 1, 1, 2), (2, 1, 1, 0), method='not implemented') From df09a28f22a919a56b91eb62f8ad1ffdbb4f3947 Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Sat, 4 Apr 2020 20:22:47 -0500 Subject: [PATCH 20/24] [numerics] Simplify TabulatedFunction constructor --- interfaces/cython/cantera/func1.pyx | 67 +++++--------------- interfaces/cython/cantera/test/test_func1.py | 22 +++---- 2 files changed, 26 insertions(+), 63 deletions(-) diff --git a/interfaces/cython/cantera/func1.pyx b/interfaces/cython/cantera/func1.pyx index 9b1d8ab8a58..802537e3abc 100644 --- a/interfaces/cython/cantera/func1.pyx +++ b/interfaces/cython/cantera/func1.pyx @@ -106,68 +106,33 @@ cdef class Func1: cdef class TabulatedFunction(Func1): """ A `TabulatedFunction` object representing a tabulated function is defined by - sample points and corresponding function values. Inputs are specified either - by two iterable objects containing sample point location and function - values, or a single array that concatenates those inputs in two rows or - columns. Between sample points, values are evaluated based on the optional - argument ``method``, which has to be supplied as a keyword; options are - ``'linear'`` (linear interpolation, default) or ``'previous'`` (nearest - previous value). Outside the sample interval, the value at the closest end - point is returned. - - Examples for `TabulatedFunction` objects using a single (two-dimensional) - array as input are:: - - >>> t1 = TabulatedFunction([[0, 2], [1, 1], [2, 0]]) + sample points and corresponding function values. Inputs are specified by + two iterable objects containing sample point location and function values. + Between sample points, values are evaluated based on the optional argument + ``method``; options are ``'linear'`` (linear interpolation, default) or + ``'previous'`` (nearest previous value). Outside the sample interval, the + value at the closest end point is returned. + + Examples for `TabulatedFunction` objects are:: + + >>> t1 = TabulatedFunction([0, 1, 2], [2, 1, 0]) >>> [t1(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]] [2.0, 2.0, 1.5, 0.5, 0.0, 0.0] - >>> t2 = TabulatedFunction(np.array([0, 1, 2]), (2, 1, 0)) + >>> t2 = TabulatedFunction(np.array([0, 1, 2]), np.array([2, 1, 0])) >>> [t2(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]] [2.0, 2.0, 1.5, 0.5, 0.0, 0.0] - where the optional ``method`` keyword argument changes the type of - interpolation from the ``'linear'`` default to ``'previous'``:: + The optional ``method`` keyword argument changes the type of interpolation + from the ``'linear'`` default to ``'previous'``:: - >>> t3 = TabulatedFunction([[0, 2], [1, 1], [2, 0]], method='previous') + >>> t3 = TabulatedFunction([0, 1, 2], [2, 1, 0], method='previous') >>> [t3(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]] [2.0, 2.0, 2.0, 1.0, 0.0, 0.0] - - Alternatively, a `TabulatedFunction` can be defined using two input arrays:: - - >>> t4 = TabulatedFunction([0, 1, 2], [2, 1, 0]) - >>> [t4(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]] - [2.0, 2.0, 1.5, 0.5, 0.0, 0.0] """ - def __init__(self, *args, method='linear'): - if len(args) == 1: - # tabulated function (single argument) - arr = np.array(args[0]) - if arr.ndim == 2: - if arr.shape[1] == 2: - time = arr[:, 0] - fval = arr[:, 1] - elif arr.shape[0] == 2: - time = arr[0, :] - fval = arr[1, :] - else: - raise ValueError("Invalid dimensions: specification of " - "tabulated function with a single array " - "requires two rows or columns") - self._set_tables(time, fval, stringify(method)) - else: - raise TypeError("'TabulatedFunction' must be constructed from " - "a numeric array with two dimensions") - - elif len(args) == 2: - # tabulated function (two arrays mimic C++ interface) - time, fval = args - self._set_tables(time, fval, stringify(method)) - - else: - raise ValueError("Invalid number of arguments (one or two " - "arguments containing tabulated values)") + def __init__(self, time, fval, method='linear'): + self._set_tables(time, fval, stringify(method)) cpdef void _set_tables(self, time, fval, string method) except *: tt = np.asarray(time, dtype=np.double) diff --git a/interfaces/cython/cantera/test/test_func1.py b/interfaces/cython/cantera/test/test_func1.py index a7d1b11fe96..3efa0a2e324 100644 --- a/interfaces/cython/cantera/test/test_func1.py +++ b/interfaces/cython/cantera/test/test_func1.py @@ -87,6 +87,13 @@ def test_tabulated2(self): self.assertNear(f, fcn(t)) def test_tabulated3(self): + time = 0, 1, 2, + fval = 2, 1, 0, + fcn = ct.TabulatedFunction(time, fval) + self.assertNear(fcn(-1), fval[0]) + self.assertNear(fcn(3), fval[-1]) + + def test_tabulated4(self): time = np.array([0, 1, 2]) fval = np.array([2, 1, 0]) fcn = ct.TabulatedFunction(time, fval) @@ -95,23 +102,14 @@ def test_tabulated3(self): for t, f in zip(tt, ff): self.assertNear(f, fcn(t)) - def test_tabulated4(self): - time = 0, 1, 2, - fval = 2, 1, 0, - fcn = ct.TabulatedFunction(time, fval) - self.assertNear(fcn(-1), fval[0]) - self.assertNear(fcn(3), fval[-1]) - def test_tabulated5(self): - fcn = ct.TabulatedFunction([[0, 2], [1, 1], [2, 0]], method='previous') + time = [0, 1, 2] + fval = [2, 1, 0] + fcn = ct.TabulatedFunction(time, fval, method='previous') val = np.array([fcn(v) for v in [-0.5, 0, 0.5, 1.5, 2, 2.5]]) self.assertArrayNear(val, np.array([2.0, 2.0, 2.0, 1.0, 0.0, 0.0])) def test_tabulated_failures(self): - with self.assertRaisesRegex(ValueError, 'Invalid number of arguments'): - ct.TabulatedFunction(1, 2, 3) - with self.assertRaisesRegex(ValueError, 'Invalid dimensions'): - ct.TabulatedFunction(np.zeros((3, 3))) with self.assertRaisesRegex(ValueError, 'do not match'): ct.TabulatedFunction(range(2), range(3)) with self.assertRaisesRegex(ValueError, 'must not be empty'): From af38e90a951be549851be16b4e76fa47e031d8d3 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sun, 5 Apr 2020 13:17:39 -0400 Subject: [PATCH 21/24] Revert accidental change of Sundials submodule version --- ext/sundials | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/sundials b/ext/sundials index 38aa03b990e..b78e575639b 160000 --- a/ext/sundials +++ b/ext/sundials @@ -1 +1 @@ -Subproject commit 38aa03b990e616db52acd053730c31fcc3eb91f9 +Subproject commit b78e575639babe99aea9a0558d0f64732d6d729e From 60c37f9e081c572b0ab3adff6dc46d85a10cc853 Mon Sep 17 00:00:00 2001 From: "Bryan W. Weber" Date: Sun, 9 Feb 2020 14:05:19 -0500 Subject: [PATCH 22/24] [PureFluid] Use pressure criteria to decide phase Use the pressure compared with the saturation pressure to decide gas or liquid phase, rather than temperature. This change is required because for some substances, the saturation temperature for low values of the pressure was outside the limits of the equation of state, resulting in convergence errors in the saturation solver. Since phaseOfMatter can throw, we need to translate exceptions in the Cython interface as well. Add tests to prevent regression with cases identifed in #786. Fixes #786. --- interfaces/cython/cantera/_cantera.pxd | 2 +- interfaces/cython/cantera/test/test_purefluid.py | 10 +++++++++- src/thermo/PureFluidPhase.cpp | 2 +- 3 files changed, 11 insertions(+), 3 deletions(-) diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index aa3dd6a263d..7bbb1bf9fa3 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -145,7 +145,7 @@ cdef extern from "cantera/thermo/ThermoPhase.h" namespace "Cantera": # miscellaneous string type() - string phaseOfMatter() + string phaseOfMatter() except +translate_exception string report(cbool, double) except +translate_exception cbool hasPhaseTransition() cbool isPure() diff --git a/interfaces/cython/cantera/test/test_purefluid.py b/interfaces/cython/cantera/test/test_purefluid.py index 59d0983900f..c9c39b84b7c 100644 --- a/interfaces/cython/cantera/test/test_purefluid.py +++ b/interfaces/cython/cantera/test/test_purefluid.py @@ -243,6 +243,14 @@ def test_phase_of_matter(self): self.water.TQ = 300, 0.4 self.assertEqual(self.water.phase_of_matter, "liquid-gas-mix") + # These cases work after fixing GH-786 + n2 = ct.Nitrogen() + n2.TP = 100, 1000 + self.assertEqual(n2.phase_of_matter, "gas") + + co2 = ct.CarbonDioxide() + self.assertEqual(co2.phase_of_matter, "gas") + # To minimize errors when transcribing tabulated data, the input units here are: # T: K, P: MPa, rho: kg/m3, v: m3/kg, (u,h): kJ/kg, s: kJ/kg-K @@ -296,7 +304,7 @@ def a(self, T, rho): def test_has_phase_transition(self): self.assertTrue(self.fluid.has_phase_transition) - + def test_consistency_temperature(self): for state in self.states: dT = 2e-5 * state.T diff --git a/src/thermo/PureFluidPhase.cpp b/src/thermo/PureFluidPhase.cpp index 8935415829f..7d4a9ea79c9 100644 --- a/src/thermo/PureFluidPhase.cpp +++ b/src/thermo/PureFluidPhase.cpp @@ -90,7 +90,7 @@ std::string PureFluidPhase::phaseOfMatter() const return "supercritical"; } else if (m_sub->TwoPhase() == 1) { return "liquid-gas-mix"; - } else if (temperature() > satTemperature(pressure())) { + } else if (pressure() < m_sub->Ps()) { return "gas"; } else { return "liquid"; From a7d4eea7e2a8377a37fe0d8c89e0811d282052f1 Mon Sep 17 00:00:00 2001 From: sin-ha Date: Wed, 8 Apr 2020 17:48:23 +0530 Subject: [PATCH 23/24] Fixes edge cases in slicing for solutionArray --- interfaces/cython/cantera/composite.py | 5 ++--- interfaces/cython/cantera/test/test_thermo.py | 17 +++++++++-------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/interfaces/cython/cantera/composite.py b/interfaces/cython/cantera/composite.py index fe06226b575..e9bec7f5270 100644 --- a/interfaces/cython/cantera/composite.py +++ b/interfaces/cython/cantera/composite.py @@ -509,9 +509,6 @@ def __init__(self, phase, shape=(0,), states=None, extra=None): if len(self._shape) == 1: self._indices = list(range(self._shape[0])) self._output_dummy = self._indices - elif len(self._shape) == 0 and len(states) == 0: - self._indices = list() - self._output_dummy = list() else: self._indices = list(np.ndindex(self._shape)) self._output_dummy = self._states[..., 0] @@ -550,6 +547,8 @@ def __getitem__(self, index): states = self._states[index] if(isinstance(states, list)): num_rows = len(states) + if num_rows == 0: + states = None return SolutionArray(self._phase, num_rows, states) else: shape = states.shape[:-1] diff --git a/interfaces/cython/cantera/test/test_thermo.py b/interfaces/cython/cantera/test/test_thermo.py index 27df9a81978..20acaf7972c 100644 --- a/interfaces/cython/cantera/test/test_thermo.py +++ b/interfaces/cython/cantera/test/test_thermo.py @@ -1805,11 +1805,12 @@ def test_slice_SolutionArray(self): self.assertEqual(len(arr.T), 3) #tests for edge cases in sliced solution array - arr1 = soln[1:1] - self.assertEqual(arr1.n_elements, 3) - - arr2 = soln[slice(0)] - self.assertEqual(arr2.n_species, 9) - - arr3 = soln[None:0] - self.assertEqual(arr3.n_reactions, 28) + states = ct.SolutionArray(self.gas, 4) + arr1 = states[3:3] + self.assertEqual(len(arr1.T), 0) + self.assertEqual(arr1.X.shape, (0,9)) + self.assertEqual(arr1.n_reactions, 28) + + states.TP = [100,300,900,323.23], ct.one_atm + arr2 = states[slice(0)] + self.assertEqual(len(arr2.T), 0) From 1a4ed672465e5de1ac0e97734baa241404a99af7 Mon Sep 17 00:00:00 2001 From: sin-ha Date: Wed, 8 Apr 2020 18:51:17 +0530 Subject: [PATCH 24/24] Fixes slicing in solutionArray --- interfaces/cython/cantera/test/test_thermo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/interfaces/cython/cantera/test/test_thermo.py b/interfaces/cython/cantera/test/test_thermo.py index 20acaf7972c..5c5043f00ad 100644 --- a/interfaces/cython/cantera/test/test_thermo.py +++ b/interfaces/cython/cantera/test/test_thermo.py @@ -1804,7 +1804,7 @@ def test_slice_SolutionArray(self): arr = soln[2:9:3] self.assertEqual(len(arr.T), 3) - #tests for edge cases in sliced solution array + #tests for edge cases in sliced solutionArray states = ct.SolutionArray(self.gas, 4) arr1 = states[3:3] self.assertEqual(len(arr1.T), 0)