From 6c5850adaf2b3bc6e284bf219f4cea0fb21713de Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Fri, 30 Oct 2020 15:49:01 -0400 Subject: [PATCH 1/5] #1220 only instantiate parameter classes once --- CHANGELOG.md | 1 + pybamm/__init__.py | 9 +++-- pybamm/expression_tree/symbol.py | 4 +-- pybamm/expression_tree/unary_operators.py | 8 ++--- pybamm/geometry/battery_geometry.py | 2 +- pybamm/parameters/electrical_parameters.py | 5 ++- pybamm/parameters/geometric_parameters.py | 3 ++ pybamm/parameters/lead_acid_parameters.py | 6 ++-- pybamm/parameters/lithium_ion_parameters.py | 6 ++-- pybamm/parameters/thermal_parameters.py | 5 ++- .../test_models/standard_output_tests.py | 36 +++++++++---------- .../test_unary_operators.py | 2 +- .../test_geometry/test_battery_geometry.py | 4 +-- .../test_parameters/test_current_functions.py | 6 ++-- .../test_electrical_parameters.py | 2 +- .../test_geometric_parameters.py | 2 +- .../test_lithium_ion_parameters.py | 2 +- 17 files changed, 58 insertions(+), 45 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f9c492050e..be712ed80d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ ## Optimizations +- Only instantiate the geometric, electrical and thermal parameter classes once ## Bug fixes diff --git a/pybamm/__init__.py b/pybamm/__init__.py index 8fb1d06f41..c452e5b487 100644 --- a/pybamm/__init__.py +++ b/pybamm/__init__.py @@ -164,9 +164,12 @@ def version(formatted=False): # from .parameters.parameter_values import ParameterValues from .parameters import constants -from .parameters.geometric_parameters import GeometricParameters -from .parameters.electrical_parameters import ElectricalParameters -from .parameters.thermal_parameters import ThermalParameters +from .parameters.geometric_parameters import geometric_parameters, GeometricParameters +from .parameters.electrical_parameters import ( + electrical_parameters, + ElectricalParameters, +) +from .parameters.thermal_parameters import thermal_parameters, ThermalParameters from .parameters.lithium_ion_parameters import LithiumIonParameters from .parameters.lead_acid_parameters import LeadAcidParameters from .parameters import parameter_sets diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py index e0db67aa98..4e9335adfe 100644 --- a/pybamm/expression_tree/symbol.py +++ b/pybamm/expression_tree/symbol.py @@ -297,8 +297,8 @@ def visualise(self, filename): pybamm.logger.error("Please install graphviz>=2.42.2 to use dot exporter") def relabel_tree(self, symbol, counter): - """ Finds all children of a symbol and assigns them a new id so that they can be - visualised properly using the graphviz output + """Finds all children of a symbol and assigns them a new id so that they can be + visualised properly using the graphviz output """ name = symbol.name if name == "div": diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py index ac244248d5..47ae30917a 100644 --- a/pybamm/expression_tree/unary_operators.py +++ b/pybamm/expression_tree/unary_operators.py @@ -1147,14 +1147,14 @@ def x_average(symbol): if a.id == b.id == c.id: return a else: - geo = pybamm.GeometricParameters() + geo = pybamm.geometric_parameters l_n = geo.l_n l_s = geo.l_s l_p = geo.l_p return (l_n * a + l_s * b + l_p * c) / (l_n + l_s + l_p) # Otherwise, use Integral to calculate average value else: - geo = pybamm.GeometricParameters() + geo = pybamm.geometric_parameters if symbol.domain == ["negative electrode"]: x = pybamm.standard_spatial_vars.x_n l = geo.l_n @@ -1214,7 +1214,7 @@ def z_average(symbol): return symbol.orphans[0] # Otherwise, use Integral to calculate average value else: - geo = pybamm.GeometricParameters() + geo = pybamm.geometric_parameters z = pybamm.standard_spatial_vars.z l_z = geo.l_z return Integral(symbol, z) / l_z @@ -1251,7 +1251,7 @@ def yz_average(symbol): return symbol.orphans[0] # Otherwise, use Integral to calculate average value else: - geo = pybamm.GeometricParameters() + geo = pybamm.geometric_parameters y = pybamm.standard_spatial_vars.y z = pybamm.standard_spatial_vars.z l_y = geo.l_y diff --git a/pybamm/geometry/battery_geometry.py b/pybamm/geometry/battery_geometry.py index d3e0b68901..7b05b7ea72 100644 --- a/pybamm/geometry/battery_geometry.py +++ b/pybamm/geometry/battery_geometry.py @@ -22,7 +22,7 @@ def battery_geometry(include_particles=True, current_collector_dimension=0): """ var = pybamm.standard_spatial_vars - geo = pybamm.GeometricParameters() + geo = pybamm.geometric_parameters l_n = geo.l_n l_s = geo.l_s diff --git a/pybamm/parameters/electrical_parameters.py b/pybamm/parameters/electrical_parameters.py index c233a4c622..6e82c9e55c 100644 --- a/pybamm/parameters/electrical_parameters.py +++ b/pybamm/parameters/electrical_parameters.py @@ -17,7 +17,7 @@ class ElectricalParameters: def __init__(self): # Get geometric parameters - self.geo = pybamm.GeometricParameters() + self.geo = pybamm.geometric_parameters # Set parameters self._set_dimensional_parameters() @@ -64,3 +64,6 @@ def _set_dimensionless_parameters(self): / self.I_typ * pybamm.Function(np.sign, self.I_typ) ) + + +electrical_parameters = ElectricalParameters() diff --git a/pybamm/parameters/geometric_parameters.py b/pybamm/parameters/geometric_parameters.py index 98d6f254c1..7efae1228c 100644 --- a/pybamm/parameters/geometric_parameters.py +++ b/pybamm/parameters/geometric_parameters.py @@ -101,3 +101,6 @@ def _set_dimensionless_parameters(self): self.l_tab_p = self.L_tab_p / self.L_z self.centre_y_tab_p = self.Centre_y_tab_p / self.L_z self.centre_z_tab_p = self.Centre_z_tab_p / self.L_z + + +geometric_parameters = GeometricParameters() diff --git a/pybamm/parameters/lead_acid_parameters.py b/pybamm/parameters/lead_acid_parameters.py index 3ccd0736ad..162793aef5 100644 --- a/pybamm/parameters/lead_acid_parameters.py +++ b/pybamm/parameters/lead_acid_parameters.py @@ -22,9 +22,9 @@ class LeadAcidParameters: def __init__(self): # Get geometric, electrical and thermal parameters - self.geo = pybamm.GeometricParameters() - self.elec = pybamm.ElectricalParameters() - self.therm = pybamm.ThermalParameters() + self.geo = pybamm.geometric_parameters + self.elec = pybamm.electrical_parameters + self.therm = pybamm.thermal_parameters # Set parameters and scales self._set_dimensional_parameters() diff --git a/pybamm/parameters/lithium_ion_parameters.py b/pybamm/parameters/lithium_ion_parameters.py index 9ec595e44d..29bc22c841 100644 --- a/pybamm/parameters/lithium_ion_parameters.py +++ b/pybamm/parameters/lithium_ion_parameters.py @@ -36,9 +36,9 @@ def __init__(self, options=None): self.options = options # Get geometric, electrical and thermal parameters - self.geo = pybamm.GeometricParameters() - self.elec = pybamm.ElectricalParameters() - self.therm = pybamm.ThermalParameters() + self.geo = pybamm.geometric_parameters + self.elec = pybamm.electrical_parameters + self.therm = pybamm.thermal_parameters # Set parameters and scales self._set_dimensional_parameters() diff --git a/pybamm/parameters/thermal_parameters.py b/pybamm/parameters/thermal_parameters.py index 81c9880274..69251ccd96 100644 --- a/pybamm/parameters/thermal_parameters.py +++ b/pybamm/parameters/thermal_parameters.py @@ -18,7 +18,7 @@ class ThermalParameters: def __init__(self): # Get geometric parameters - self.geo = pybamm.GeometricParameters() + self.geo = pybamm.geometric_parameters # Set parameters self._set_dimensional_parameters() @@ -178,3 +178,6 @@ def _set_dimensionless_parameters(self): def T_amb(self, t): "Dimensionless ambient temperature" return (self.T_amb_dim(t) - self.T_ref) / self.Delta_T + + +thermal_parameters = ThermalParameters() diff --git a/tests/integration/test_models/standard_output_tests.py b/tests/integration/test_models/standard_output_tests.py index d186fc4596..8530aae667 100644 --- a/tests/integration/test_models/standard_output_tests.py +++ b/tests/integration/test_models/standard_output_tests.py @@ -67,7 +67,7 @@ def __init__(self, model, param, disc, solution, operating_condition): # Use dimensional time and space self.t = solution.t * model.timescale_eval - geo = pybamm.GeometricParameters() + geo = pybamm.geometric_parameters L_x = param.evaluate(geo.L_x) self.x_n = disc.mesh["negative electrode"].nodes * L_x @@ -127,10 +127,10 @@ def __init__(self, model, param, disc, solution, operating_condition): def test_each_reaction_overpotential(self): """Testing that: - - discharge: eta_r_n > 0, eta_r_p < 0 - - charge: eta_r_n < 0, eta_r_p > 0 - - off: eta_r_n == 0, eta_r_p == 0 - """ + - discharge: eta_r_n > 0, eta_r_p < 0 + - charge: eta_r_n < 0, eta_r_p > 0 + - off: eta_r_n == 0, eta_r_p == 0 + """ tol = 0.01 t, x_n, x_p = self.t, self.x_n, self.x_p if self.operating_condition == "discharge": @@ -145,9 +145,9 @@ def test_each_reaction_overpotential(self): def test_overpotentials(self): """Testing that all are: - - discharge: . < 0 - - charge: . > 0 - - off: . == 0 + - discharge: . < 0 + - charge: . > 0 + - off: . == 0 """ tol = 0.001 if self.operating_condition == "discharge": @@ -168,10 +168,10 @@ def test_overpotentials(self): ) def test_ocps(self): - """ Testing that: - - discharge: ocp_n increases, ocp_p decreases - - charge: ocp_n decreases, ocp_p increases - - off: ocp_n, ocp_p constant + """Testing that: + - discharge: ocp_n increases, ocp_p decreases + - charge: ocp_n decreases, ocp_p increases + - off: ocp_n, ocp_p constant """ neg_end_vs_start = self.ocp_n_av(self.t[-1]) - self.ocp_n_av(self.t[1]) pos_end_vs_start = self.ocp_p_av(self.t[-1]) - self.ocp_p_av(self.t[1]) @@ -187,9 +187,9 @@ def test_ocps(self): def test_ocv(self): """Testing that: - - discharge: ocv decreases - - charge: ocv increases - - off: ocv constant + - discharge: ocv decreases + - charge: ocv increases + - off: ocv constant """ end_vs_start = self.ocv_av(self.t[-1]) - self.ocv_av(self.t[1]) @@ -203,9 +203,9 @@ def test_ocv(self): def test_voltage(self): """Testing that: - - discharge: voltage decreases - - charge: voltage increases - - off: voltage constant + - discharge: voltage decreases + - charge: voltage increases + - off: voltage constant """ end_vs_start = self.voltage(self.t[-1]) - self.voltage(self.t[1]) diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py index b5d6f6bdb0..9b4f6a218f 100644 --- a/tests/unit/test_expression_tree/test_unary_operators.py +++ b/tests/unit/test_expression_tree/test_unary_operators.py @@ -460,7 +460,7 @@ def test_x_average(self): pybamm.x_average(symbol_on_edges) # Particle domains - geo = pybamm.GeometricParameters() + geo = pybamm.geometric_parameters l_n = geo.l_n l_p = geo.l_p diff --git a/tests/unit/test_geometry/test_battery_geometry.py b/tests/unit/test_geometry/test_battery_geometry.py index 668a4b2234..e2b11a064d 100644 --- a/tests/unit/test_geometry/test_battery_geometry.py +++ b/tests/unit/test_geometry/test_battery_geometry.py @@ -17,7 +17,7 @@ def test_geometry_keys(self): def test_geometry(self): var = pybamm.standard_spatial_vars - geo = pybamm.GeometricParameters() + geo = pybamm.geometric_parameters for cc_dimension in [0, 1, 2]: geometry = pybamm.battery_geometry(current_collector_dimension=cc_dimension) self.assertIsInstance(geometry, pybamm.Geometry) @@ -42,7 +42,7 @@ class TestReadParameters(unittest.TestCase): # This is the most complicated geometry and should test the parameters are # all returned for the deepest dict def test_read_parameters(self): - geo = pybamm.GeometricParameters() + geo = pybamm.geometric_parameters L_n = geo.L_n L_s = geo.L_s L_p = geo.L_p diff --git a/tests/unit/test_parameters/test_current_functions.py b/tests/unit/test_parameters/test_current_functions.py index b640842d69..2222d3d3b8 100644 --- a/tests/unit/test_parameters/test_current_functions.py +++ b/tests/unit/test_parameters/test_current_functions.py @@ -10,7 +10,7 @@ class TestCurrentFunctions(unittest.TestCase): def test_constant_current(self): # test simplify - param = pybamm.ElectricalParameters() + param = pybamm.electrical_parameters current = param.current_with_time parameter_values = pybamm.ParameterValues( { @@ -24,7 +24,7 @@ def test_constant_current(self): def test_get_current_data(self): # test process parameters - param = pybamm.ElectricalParameters() + param = pybamm.electrical_parameters dimensional_current = param.dimensional_current_with_time parameter_values = pybamm.ParameterValues( { @@ -47,7 +47,7 @@ def my_fun(t, A, omega): return A * pybamm.sin(2 * np.pi * omega * t) # choose amplitude and frequency - param = pybamm.ElectricalParameters() + param = pybamm.electrical_parameters A = param.I_typ omega = pybamm.Parameter("omega") diff --git a/tests/unit/test_parameters/test_electrical_parameters.py b/tests/unit/test_parameters/test_electrical_parameters.py index 0a51f7bb9e..dd513d1d42 100644 --- a/tests/unit/test_parameters/test_electrical_parameters.py +++ b/tests/unit/test_parameters/test_electrical_parameters.py @@ -9,7 +9,7 @@ class TestElectricalParameters(unittest.TestCase): def test_current_functions(self): # create current functions - param = pybamm.ElectricalParameters() + param = pybamm.electrical_parameters dimensional_current = param.dimensional_current_with_time dimensional_current_density = param.dimensional_current_density_with_time dimensionless_current = param.current_with_time diff --git a/tests/unit/test_parameters/test_geometric_parameters.py b/tests/unit/test_parameters/test_geometric_parameters.py index 63df878fac..69c769830c 100644 --- a/tests/unit/test_parameters/test_geometric_parameters.py +++ b/tests/unit/test_parameters/test_geometric_parameters.py @@ -8,7 +8,7 @@ class TestGeometricParameters(unittest.TestCase): def test_macroscale_parameters(self): - geo = pybamm.GeometricParameters() + geo = pybamm.geometric_parameters L_n = geo.L_n L_s = geo.L_s L_p = geo.L_p diff --git a/tests/unit/test_parameters/test_lithium_ion_parameters.py b/tests/unit/test_parameters/test_lithium_ion_parameters.py index 8d48a84b10..965e604056 100644 --- a/tests/unit/test_parameters/test_lithium_ion_parameters.py +++ b/tests/unit/test_parameters/test_lithium_ion_parameters.py @@ -194,7 +194,7 @@ def test_thermal_parameters(self): # values.evaluate(param.tau_th_yz), 1.4762 * 10 ** (3), 2 # ) - # thermal = pybamm.ThermalParameters() + # thermal = pybamm.thermal_parameters # np.testing.assert_almost_equal( # values.evaluate(thermal.rho_eff_dim), 1.8116 * 10 ** (6), 2 # ) From 5de58b789d5aca9c31e41142f9f1fbb57aae1502 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Fri, 30 Oct 2020 16:25:15 -0400 Subject: [PATCH 2/5] #1229 rewrite is_constant to be faster --- pybamm/expression_tree/array.py | 4 ++++ pybamm/expression_tree/binary_operators.py | 4 ++++ pybamm/expression_tree/concatenations.py | 4 ++++ pybamm/expression_tree/functions.py | 4 ++++ pybamm/expression_tree/scalar.py | 8 +++++--- pybamm/expression_tree/symbol.py | 15 +++------------ pybamm/expression_tree/unary_operators.py | 4 ++++ tests/unit/test_expression_tree/test_symbol.py | 5 +---- 8 files changed, 29 insertions(+), 19 deletions(-) diff --git a/pybamm/expression_tree/array.py b/pybamm/expression_tree/array.py index 95622571b4..cd068d0b0d 100644 --- a/pybamm/expression_tree/array.py +++ b/pybamm/expression_tree/array.py @@ -105,6 +105,10 @@ def _base_evaluate(self, t=None, y=None, y_dot=None, inputs=None): """ See :meth:`pybamm.Symbol._base_evaluate()`. """ return self._entries + def is_constant(self): + """ See :meth:`pybamm.Symbol.is_constant()`. """ + return True + def linspace(start, stop, num=50, **kwargs): """ diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py index ef49658ad5..8011246806 100644 --- a/pybamm/expression_tree/binary_operators.py +++ b/pybamm/expression_tree/binary_operators.py @@ -218,6 +218,10 @@ def evaluates_on_edges(self, dimension): dimension ) + def is_constant(self): + """ See :meth:`pybamm.Symbol.is_constant()`. """ + return self.left.is_constant() and self.right.is_constant() + class Power(BinaryOperator): """A node in the expression tree representing a `**` power operator diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py index 456029314a..f9f9e627bb 100644 --- a/pybamm/expression_tree/concatenations.py +++ b/pybamm/expression_tree/concatenations.py @@ -103,6 +103,10 @@ def _evaluate_for_shape(self): [child.evaluate_for_shape() for child in self.children] ) + def is_constant(self): + """ See :meth:`pybamm.Symbol.is_constant()`. """ + return all(child.is_constant() for child in self.children) + class NumpyConcatenation(Concatenation): """A node in the expression tree representing a concatenation of equations, when we diff --git a/pybamm/expression_tree/functions.py b/pybamm/expression_tree/functions.py index 1d901d16ef..aad86d272d 100644 --- a/pybamm/expression_tree/functions.py +++ b/pybamm/expression_tree/functions.py @@ -174,6 +174,10 @@ def evaluates_on_edges(self, dimension): """ See :meth:`pybamm.Symbol.evaluates_on_edges()`. """ return any(child.evaluates_on_edges(dimension) for child in self.children) + def is_constant(self): + """ See :meth:`pybamm.Symbol.is_constant()`. """ + return all(child.is_constant() for child in self.children) + def _evaluate_for_shape(self): """ Default behaviour: has same shape as all child diff --git a/pybamm/expression_tree/scalar.py b/pybamm/expression_tree/scalar.py index fffc461960..4f571a2a3b 100644 --- a/pybamm/expression_tree/scalar.py +++ b/pybamm/expression_tree/scalar.py @@ -24,9 +24,7 @@ class Scalar(pybamm.Symbol): """ def __init__(self, value, name=None, domain=[]): - """ - - """ + """""" # set default name if not provided self.value = value if name is None: @@ -62,3 +60,7 @@ def _jac(self, variable): def new_copy(self): """ See :meth:`pybamm.Symbol.new_copy()`. """ return Scalar(self.value, self.name, self.domain) + + def is_constant(self): + """ See :meth:`pybamm.Symbol.is_constant()`. """ + return True diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py index 4e9335adfe..e1242a85ff 100644 --- a/pybamm/expression_tree/symbol.py +++ b/pybamm/expression_tree/symbol.py @@ -602,24 +602,15 @@ def _evaluate_for_shape(self): def is_constant(self): """returns true if evaluating the expression is not dependent on `t` or `y` - or `u` + or `inputs` See Also -------- evaluate : evaluate the expression """ - # if any of the nodes are instances of any of these types, then the whole - # expression depends on either t or y or u - search_types = ( - pybamm.Variable, - pybamm.StateVector, - pybamm.Time, - pybamm.InputParameter, - ) - - # do the search, return true if no relevent nodes are found - return not self.has_symbol_of_classes(search_types) + # Default behaviour is False + return False def evaluate_ignoring_errors(self, t=0): """ diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py index 47ae30917a..51f6f3e4f4 100644 --- a/pybamm/expression_tree/unary_operators.py +++ b/pybamm/expression_tree/unary_operators.py @@ -87,6 +87,10 @@ def evaluates_on_edges(self, dimension): """ See :meth:`pybamm.Symbol.evaluates_on_edges()`. """ return self.child.evaluates_on_edges(dimension) + def is_constant(self): + """ See :meth:`pybamm.Symbol.is_constant()`. """ + return self.child.is_constant() + class Negate(UnaryOperator): """A node in the expression tree representing a `-` negation operator diff --git a/tests/unit/test_expression_tree/test_symbol.py b/tests/unit/test_expression_tree/test_symbol.py index 2c23e09fc7..6153711e66 100644 --- a/tests/unit/test_expression_tree/test_symbol.py +++ b/tests/unit/test_expression_tree/test_symbol.py @@ -185,14 +185,11 @@ def test_symbol_is_constant(self): self.assertFalse(a.is_constant()) a = pybamm.Parameter("a") - self.assertTrue(a.is_constant()) + self.assertFalse(a.is_constant()) a = pybamm.Scalar(1) * pybamm.Variable("a") self.assertFalse(a.is_constant()) - a = pybamm.Scalar(1) * pybamm.Parameter("a") - self.assertTrue(a.is_constant()) - a = pybamm.Scalar(1) * pybamm.StateVector(slice(10)) self.assertFalse(a.is_constant()) From 9426665e274bfa4836b94f3a91889e8180a66890 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Fri, 30 Oct 2020 19:17:59 -0400 Subject: [PATCH 3/5] #1220 some more optimizations --- pybamm/expression_tree/input_parameter.py | 4 ++ pybamm/expression_tree/symbol.py | 66 ++++++++++++------- pybamm/meshes/meshes.py | 13 +++- pybamm/solvers/base_solver.py | 10 +-- pybamm/solvers/casadi_algebraic_solver.py | 9 +++ pybamm/solvers/casadi_solver.py | 15 +++-- pybamm/util.py | 10 ++- .../test_meshes/test_scikit_fem_submesh.py | 5 +- tests/unit/test_timer.py | 25 ++++--- 9 files changed, 102 insertions(+), 55 deletions(-) diff --git a/pybamm/expression_tree/input_parameter.py b/pybamm/expression_tree/input_parameter.py index e63242cf81..bbc78b7a8f 100644 --- a/pybamm/expression_tree/input_parameter.py +++ b/pybamm/expression_tree/input_parameter.py @@ -36,6 +36,10 @@ def set_expected_size(self, size): "Specify the size that the input parameter should be" self._expected_size = size + # We also need to update the saved size and shape + self._saved_size = size + self._saved_shape = (size, 1) + def _evaluate_for_shape(self): """ Returns the scalar 'NaN' to represent the shape of a parameter. diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py index e1242a85ff..68bda1a61f 100644 --- a/pybamm/expression_tree/symbol.py +++ b/pybamm/expression_tree/symbol.py @@ -174,8 +174,8 @@ def domain(self, domain): @property def auxiliary_domains(self): - "Returns domains that are not the primary domain" - return {k: v for k, v in self._domains.items() if k != "primary"} + "Returns auxiliary domains" + return self._auxiliary_domains @auxiliary_domains.setter def auxiliary_domains(self, auxiliary_domains): @@ -193,6 +193,7 @@ def auxiliary_domains(self, auxiliary_domains): if len(set(values)) != len(values): raise pybamm.DomainError("All auxiliary domains must be different") + self._auxiliary_domains = auxiliary_domains self._domains.update(auxiliary_domains) @property @@ -203,10 +204,16 @@ def secondary_domain(self): def copy_domains(self, symbol): "Copy the domains from a given symbol, bypassing checks" self._domains = symbol.domains + self._domain = self._domains["primary"] + self._auxiliary_domains = { + k: v for k, v in self._domains.items() if k != "primary" + } def clear_domains(self): "Clear domains, bypassing checks" self._domains = {"primary": []} + self._domain = [] + self._auxiliary_domains = {} def get_children_auxiliary_domains(self, children): "Combine auxiliary domains from children, at all levels" @@ -727,37 +734,46 @@ def size(self): """ Size of an object, found by evaluating it with appropriate t and y """ - return np.prod(self.shape) + try: + return self._saved_size + except AttributeError: + self._saved_size = np.prod(self.shape) + return self._saved_size @property def shape(self): """ Shape of an object, found by evaluating it with appropriate t and y. """ - # Default behaviour is to try to evaluate the object directly - # Try with some large y, to avoid having to unpack (slow) try: - y = np.nan * np.ones((1000, 1)) - evaluated_self = self.evaluate(0, y, y, inputs="shape test") - # If that fails, fall back to calculating how big y should really be - except ValueError: - unpacker = pybamm.SymbolUnpacker(pybamm.StateVector) - state_vectors_in_node = unpacker.unpack_symbol(self).values() - if state_vectors_in_node == []: - y = None + return self._saved_shape + except AttributeError: + # Default behaviour is to try to evaluate the object directly + # Try with some large y, to avoid having to unpack (slow) + try: + y = np.nan * np.ones((1000, 1)) + evaluated_self = self.evaluate(0, y, y, inputs="shape test") + # If that fails, fall back to calculating how big y should really be + except ValueError: + unpacker = pybamm.SymbolUnpacker(pybamm.StateVector) + state_vectors_in_node = unpacker.unpack_symbol(self).values() + if state_vectors_in_node == []: + y = None + else: + min_y_size = max( + len(x._evaluation_array) for x in state_vectors_in_node + ) + # Pick a y that won't cause RuntimeWarnings + y = np.nan * np.ones((min_y_size, 1)) + evaluated_self = self.evaluate(0, y, y, inputs="shape test") + + # Return shape of evaluated object + if isinstance(evaluated_self, numbers.Number): + self._saved_shape = () else: - min_y_size = max( - len(x._evaluation_array) for x in state_vectors_in_node - ) - # Pick a y that won't cause RuntimeWarnings - y = np.nan * np.ones((min_y_size, 1)) - evaluated_self = self.evaluate(0, y, y, inputs="shape test") - - # Return shape of evaluated object - if isinstance(evaluated_self, numbers.Number): - return () - else: - return evaluated_self.shape + self._saved_shape = evaluated_self.shape + + return self._saved_shape @property def size_for_testing(self): diff --git a/pybamm/meshes/meshes.py b/pybamm/meshes/meshes.py index 6f273afcc6..4b70f8ca6b 100644 --- a/pybamm/meshes/meshes.py +++ b/pybamm/meshes/meshes.py @@ -94,7 +94,18 @@ def __init__(self, geometry, submesh_types, var_pts): else: for lim, sym in spatial_limits.items(): if isinstance(sym, pybamm.Symbol): - sym_eval = sym.evaluate() + try: + sym_eval = sym.evaluate() + except NotImplementedError as e: + if sym.has_symbol_of_classes(pybamm.Parameter): + raise pybamm.DiscretisationError( + "Parameter values have not yet been set for " + "geometry. Make sure that something like " + "`param.process_geometry(geometry)` has been " + "run." + ) + else: + raise e elif isinstance(sym, numbers.Number): sym_eval = sym geometry[domain][spatial_variable][lim] = sym_eval diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py index fb1d1c9c09..981b53d2f2 100644 --- a/pybamm/solvers/base_solver.py +++ b/pybamm/solvers/base_solver.py @@ -563,16 +563,13 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None): # Set up (if not done already) if model not in self.models_set_up: self.set_up(model, ext_and_inputs, t_eval) - set_up_time = timer.time() self.models_set_up.update( {model: {"initial conditions": model.concatenated_initial_conditions}} ) else: ics_set_up = self.models_set_up[model]["initial conditions"] # Check that initial conditions have not been updated - if ics_set_up.id == model.concatenated_initial_conditions.id: - set_up_time = 0 - else: + if ics_set_up.id != model.concatenated_initial_conditions.id: # If the new initial conditions are different, set up again # Doing the whole setup again might be slow, but no need to prematurely # optimize this @@ -580,7 +577,7 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None): self.models_set_up[model][ "initial conditions" ] = model.concatenated_initial_conditions - set_up_time = timer.time() + set_up_time = timer.time() # (Re-)calculate consistent initial conditions self._set_initial_conditions(model, ext_and_inputs, update_rhs=True) @@ -806,12 +803,11 @@ def step( ) self.set_up(model, ext_and_inputs) t = 0.0 - set_up_time = timer.time() else: # initialize with old solution t = old_solution.t[-1] model.y0 = old_solution.y[:, -1] - set_up_time = 0 + set_up_time = timer.time() # (Re-)calculate consistent initial conditions self._set_initial_conditions(model, ext_and_inputs, update_rhs=False) diff --git a/pybamm/solvers/casadi_algebraic_solver.py b/pybamm/solvers/casadi_algebraic_solver.py index 25fe6aff1d..0b484edfbe 100644 --- a/pybamm/solvers/casadi_algebraic_solver.py +++ b/pybamm/solvers/casadi_algebraic_solver.py @@ -66,6 +66,15 @@ def _integrate(self, model, t_eval, inputs=None): inputs = casadi.vertcat(*[v for v in inputs.values()]) y0 = model.y0 + + # If y0 already satisfies the tolerance for all t then keep it + if has_symbolic_inputs is False and all( + np.all(abs(model.casadi_algebraic(t, y0, inputs).full()) < self.tol) + for t in t_eval + ): + pybamm.logger.debug("Keeping same solution at all times") + return pybamm.Solution(t_eval, y0, termination="success") + # The casadi algebraic solver can read rhs equations, but leaves them unchanged # i.e. the part of the solution vector that corresponds to the differential # equations will be equal to the initial condition provided. This allows this diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index d13af9785d..6f7272b889 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -310,14 +310,19 @@ def create_integrator(self, model, inputs, t_eval=None): if model in self.integrators: # If we're not using the grid, we don't need to change the integrator if use_grid is False: - return self.integrators[model] + return self.integrators[model][0] # Otherwise, create new integrator with an updated grid + # We don't need to update the grid if reusing the same t_eval else: method, problem, options = self.integrator_specs[model] - options["grid"] = t_eval - integrator = casadi.integrator("F", method, problem, options) - self.integrators[model] = (integrator, use_grid) - return integrator + t_eval_old = options["grid"] + if np.array_equal(t_eval_old, t_eval): + return self.integrators[model][0] + else: + options["grid"] = t_eval + integrator = casadi.integrator("F", method, problem, options) + self.integrators[model] = (integrator, use_grid) + return integrator else: y0 = model.y0 rhs = model.casadi_rhs diff --git a/pybamm/util.py b/pybamm/util.py index f58e04280f..46625e0d63 100644 --- a/pybamm/util.py +++ b/pybamm/util.py @@ -154,10 +154,14 @@ def format(self, time=None): """ if time is None: time = self.time() - if time < 1e-2: - return str(time) + " seconds" + if time < 1e-6: + return "{:.3f} ns".format(time * 1e9) + if time < 1e-3: + return "{:.3f} us".format(time * 1e6) + if time < 1: + return "{:.3f} ms".format(time * 1e3) elif time < 60: - return str(round(time, 2)) + " seconds" + return "{:.3f} s".format(time) output = [] time = int(round(time)) units = [(604800, "week"), (86400, "day"), (3600, "hour"), (60, "minute")] diff --git a/tests/unit/test_meshes/test_scikit_fem_submesh.py b/tests/unit/test_meshes/test_scikit_fem_submesh.py index ade3e4c760..207016bf1c 100644 --- a/tests/unit/test_meshes/test_scikit_fem_submesh.py +++ b/tests/unit/test_meshes/test_scikit_fem_submesh.py @@ -77,7 +77,10 @@ def test_init_failure(self): var = pybamm.standard_spatial_vars var_pts = {var.x_n: 10, var.x_s: 10, var.x_p: 10, var.y: 10, var.z: 10} # there are parameters in the variables that need to be processed - with self.assertRaises(NotImplementedError): + with self.assertRaisesRegex( + pybamm.DiscretisationError, + "Parameter values have not yet been set for geometry", + ): pybamm.Mesh(geometry, submesh_types, var_pts) lims = {var.x_n: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}} diff --git a/tests/unit/test_timer.py b/tests/unit/test_timer.py index 74d8eabcd6..531479881f 100644 --- a/tests/unit/test_timer.py +++ b/tests/unit/test_timer.py @@ -29,20 +29,19 @@ def test_timing(self): self.assertLess(b, a) def test_timer_format(self): - import sys - t = pybamm.Timer() - self.assertEqual(t.format(1e-3), "0.001 seconds") - self.assertEqual(t.format(0.000123456789), "0.000123456789 seconds") - self.assertEqual(t.format(0.123456789), "0.12 seconds") - if sys.hexversion < 0x3000000: - self.assertEqual(t.format(2), "2.0 seconds") - else: - self.assertEqual(t.format(2), "2 seconds") - self.assertEqual(t.format(2.5), "2.5 seconds") - self.assertEqual(t.format(12.5), "12.5 seconds") - self.assertEqual(t.format(59.41), "59.41 seconds") - self.assertEqual(t.format(59.4126347547), "59.41 seconds") + self.assertEqual(t.format(1e-9), "1.000 ns") + self.assertEqual(t.format(0.000000123456789), "123.457 ns") + self.assertEqual(t.format(1e-6), "1.000 us") + self.assertEqual(t.format(0.000123456789), "123.457 us") + self.assertEqual(t.format(0.999e-3), "999.000 us") + self.assertEqual(t.format(1e-3), "1.000 ms") + self.assertEqual(t.format(0.123456789), "123.457 ms") + self.assertEqual(t.format(2), "2.000 s") + self.assertEqual(t.format(2.5), "2.500 s") + self.assertEqual(t.format(12.5), "12.500 s") + self.assertEqual(t.format(59.41), "59.410 s") + self.assertEqual(t.format(59.4126347547), "59.413 s") self.assertEqual(t.format(60.2), "1 minute, 0 seconds") self.assertEqual(t.format(61), "1 minute, 1 second") self.assertEqual(t.format(121), "2 minutes, 1 second") From 327790583aad45a22d6a613725b6b2245a8c3879 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Fri, 30 Oct 2020 19:57:51 -0400 Subject: [PATCH 4/5] #1220 update changelog and fix formatting --- CHANGELOG.md | 5 ++++- pybamm/expression_tree/scalar.py | 1 - pybamm/meshes/meshes.py | 4 ++-- 3 files changed, 6 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index be712ed80d..efcc93217e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,7 +9,10 @@ ## Optimizations -- Only instantiate the geometric, electrical and thermal parameter classes once +- Avoid unnecessary repeated computations in the solvers ([#1222](https://github.com/pybamm-team/PyBaMM/pull/1222)) +- Rewrite `Symbol.is_constant` to be more efficient ([#1222](https://github.com/pybamm-team/PyBaMM/pull/1222)) +- Cache shape and size calculations ([#1222](https://github.com/pybamm-team/PyBaMM/pull/1222)) +- Only instantiate the geometric, electrical and thermal parameter classes once ([#1222](https://github.com/pybamm-team/PyBaMM/pull/1222)) ## Bug fixes diff --git a/pybamm/expression_tree/scalar.py b/pybamm/expression_tree/scalar.py index 4f571a2a3b..5bde1f147e 100644 --- a/pybamm/expression_tree/scalar.py +++ b/pybamm/expression_tree/scalar.py @@ -24,7 +24,6 @@ class Scalar(pybamm.Symbol): """ def __init__(self, value, name=None, domain=[]): - """""" # set default name if not provided self.value = value if name is None: diff --git a/pybamm/meshes/meshes.py b/pybamm/meshes/meshes.py index 4b70f8ca6b..3eb88a1d2b 100644 --- a/pybamm/meshes/meshes.py +++ b/pybamm/meshes/meshes.py @@ -96,7 +96,7 @@ def __init__(self, geometry, submesh_types, var_pts): if isinstance(sym, pybamm.Symbol): try: sym_eval = sym.evaluate() - except NotImplementedError as e: + except NotImplementedError as error: if sym.has_symbol_of_classes(pybamm.Parameter): raise pybamm.DiscretisationError( "Parameter values have not yet been set for " @@ -105,7 +105,7 @@ def __init__(self, geometry, submesh_types, var_pts): "run." ) else: - raise e + raise error elif isinstance(sym, numbers.Number): sym_eval = sym geometry[domain][spatial_variable][lim] = sym_eval From ebc0eabb0d3b459ce7c0fd3750f9d72fe39e63d5 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Sun, 1 Nov 2020 23:12:20 -0500 Subject: [PATCH 5/5] #1220 coverage --- pybamm/expression_tree/parameter.py | 9 ++++++- pybamm/expression_tree/symbol.py | 15 +++++------ pybamm/solvers/casadi_algebraic_solver.py | 1 + .../unit/test_expression_tree/test_symbol.py | 5 +++- tests/unit/test_meshes/test_meshes.py | 26 +++++++++++++++++++ tests/unit/test_models/test_base_model.py | 6 +++++ .../test_casadi_algebraic_solver.py | 20 ++++++++++++++ 7 files changed, 71 insertions(+), 11 deletions(-) diff --git a/pybamm/expression_tree/parameter.py b/pybamm/expression_tree/parameter.py index eca440ebd2..83ce18a0dc 100644 --- a/pybamm/expression_tree/parameter.py +++ b/pybamm/expression_tree/parameter.py @@ -35,6 +35,10 @@ def _evaluate_for_shape(self): """ return np.nan + def is_constant(self): + """ See :meth:`pybamm.Symbol.is_constant()`. """ + return True + class FunctionParameter(pybamm.Symbol): """A node in the expression tree representing a function parameter @@ -60,7 +64,10 @@ class FunctionParameter(pybamm.Symbol): """ def __init__( - self, name, inputs, diff_variable=None, + self, + name, + inputs, + diff_variable=None, ): # assign diff variable self.diff_variable = diff_variable diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py index 68bda1a61f..51b0db4af6 100644 --- a/pybamm/expression_tree/symbol.py +++ b/pybamm/expression_tree/symbol.py @@ -299,7 +299,7 @@ def visualise(self, filename): DotExporter( new_node, nodeattrfunc=lambda node: 'label="{}"'.format(node.label) ).to_picture(filename) - except FileNotFoundError: + except FileNotFoundError: # pragma: no cover # raise error but only through logger so that test passes pybamm.logger.error("Please install graphviz>=2.42.2 to use dot exporter") @@ -757,14 +757,11 @@ def shape(self): except ValueError: unpacker = pybamm.SymbolUnpacker(pybamm.StateVector) state_vectors_in_node = unpacker.unpack_symbol(self).values() - if state_vectors_in_node == []: - y = None - else: - min_y_size = max( - len(x._evaluation_array) for x in state_vectors_in_node - ) - # Pick a y that won't cause RuntimeWarnings - y = np.nan * np.ones((min_y_size, 1)) + min_y_size = max( + max(len(x._evaluation_array) for x in state_vectors_in_node), 1 + ) + # Pick a y that won't cause RuntimeWarnings + y = np.nan * np.ones((min_y_size, 1)) evaluated_self = self.evaluate(0, y, y, inputs="shape test") # Return shape of evaluated object diff --git a/pybamm/solvers/casadi_algebraic_solver.py b/pybamm/solvers/casadi_algebraic_solver.py index 0b484edfbe..5c33b39795 100644 --- a/pybamm/solvers/casadi_algebraic_solver.py +++ b/pybamm/solvers/casadi_algebraic_solver.py @@ -155,6 +155,7 @@ def _integrate(self, model, t_eval, inputs=None): ): # update initial guess for the next iteration y0_alg = y_alg_sol + y0 = casadi.vertcat(y0_diff, y0_alg) # update solution array if y_alg is None: y_alg = y_alg_sol diff --git a/tests/unit/test_expression_tree/test_symbol.py b/tests/unit/test_expression_tree/test_symbol.py index 6153711e66..8b2cc02019 100644 --- a/tests/unit/test_expression_tree/test_symbol.py +++ b/tests/unit/test_expression_tree/test_symbol.py @@ -185,7 +185,10 @@ def test_symbol_is_constant(self): self.assertFalse(a.is_constant()) a = pybamm.Parameter("a") - self.assertFalse(a.is_constant()) + self.assertTrue(a.is_constant()) + + a = pybamm.Scalar(1) * pybamm.Parameter("a") + self.assertTrue(a.is_constant()) a = pybamm.Scalar(1) * pybamm.Variable("a") self.assertFalse(a.is_constant()) diff --git a/tests/unit/test_meshes/test_meshes.py b/tests/unit/test_meshes/test_meshes.py index f8a3c50acd..ff71555293 100644 --- a/tests/unit/test_meshes/test_meshes.py +++ b/tests/unit/test_meshes/test_meshes.py @@ -107,6 +107,32 @@ def test_init_failure(self): with self.assertRaisesRegex(KeyError, "Points not given"): pybamm.Mesh(geometry, submesh_types, var_pts) + # Not processing geometry parameters + geometry = pybamm.battery_geometry() + + var = pybamm.standard_spatial_vars + var_pts = {var.x_n: 10, var.x_s: 10, var.x_p: 12, var.r_n: 5, var.r_p: 6} + + submesh_types = { + "negative electrode": pybamm.MeshGenerator(pybamm.Uniform1DSubMesh), + "separator": pybamm.MeshGenerator(pybamm.Uniform1DSubMesh), + "positive electrode": pybamm.MeshGenerator(pybamm.Uniform1DSubMesh), + "negative particle": pybamm.MeshGenerator(pybamm.Uniform1DSubMesh), + "positive particle": pybamm.MeshGenerator(pybamm.Uniform1DSubMesh), + "current collector": pybamm.MeshGenerator(pybamm.SubMesh0D), + } + + with self.assertRaisesRegex(pybamm.DiscretisationError, "Parameter values"): + pybamm.Mesh(geometry, submesh_types, var_pts) + + # Geometry has an unrecognized variable type + var = pybamm.standard_spatial_vars + geometry["negative electrode"] = { + var.x_n: {"min": 0, "max": pybamm.Variable("var")} + } + with self.assertRaisesRegex(NotImplementedError, "for symbol var"): + pybamm.Mesh(geometry, submesh_types, var_pts) + def test_mesh_sizes(self): param = pybamm.ParameterValues( values={ diff --git a/tests/unit/test_models/test_base_model.py b/tests/unit/test_models/test_base_model.py index d301ca0449..a054067539 100644 --- a/tests/unit/test_models/test_base_model.py +++ b/tests/unit/test_models/test_base_model.py @@ -577,6 +577,12 @@ def test_export_casadi(self): variable_names = ["Volume-averaged cell temperature"] out = sim.built_model.export_casadi_objects(variable_names) + # Test fails if not discretised + with self.assertRaisesRegex( + pybamm.DiscretisationError, "Cannot automatically discretise model" + ): + model.export_casadi_objects(["Electrolyte concentration"]) + @unittest.skipIf(platform.system() == "Windows", "Skipped for Windows") def test_generate_casadi(self): model = pybamm.BaseModel() diff --git a/tests/unit/test_solvers/test_casadi_algebraic_solver.py b/tests/unit/test_solvers/test_casadi_algebraic_solver.py index 1d8422318f..862ea298ca 100644 --- a/tests/unit/test_solvers/test_casadi_algebraic_solver.py +++ b/tests/unit/test_solvers/test_casadi_algebraic_solver.py @@ -105,6 +105,26 @@ def test_model_solver_with_time(self): sol[1, :], ) + def test_model_solver_with_time_not_changing(self): + # Create model + model = pybamm.BaseModel() + var1 = pybamm.Variable("var1") + var2 = pybamm.Variable("var2") + model.algebraic = {var1: var1 - 3, var2: 2 * var1 - var2} + model.initial_conditions = {var1: pybamm.Scalar(1), var2: pybamm.Scalar(4)} + model.variables = {"var1": var1, "var2": var2} + + disc = pybamm.Discretisation() + disc.process_model(model) + + # Solve + t_eval = np.linspace(0, 1) + solver = pybamm.CasadiAlgebraicSolver() + solution = solver.solve(model, t_eval) + + sol = np.vstack((3 + 0 * t_eval, 6 + 0 * t_eval)) + np.testing.assert_array_almost_equal(solution.y, sol) + def test_model_solver_with_bounds(self): # Note: we need a better test case to test this functionality properly # Create model