Skip to content

Commit

Permalink
Merge pull request #1020 from pybamm-team/issue-1006-process-space-var
Browse files Browse the repository at this point in the history 
Issue 1006 process space var
  • Loading branch information
@rtimms
rtimms authored May 27, 2020
2 parents 1007c87 + e86523e commit bb5f204
Show file tree
Hide file tree
Showing 5 changed files with 250 additions and 113 deletions.
2 changes: 2 additions & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@

## Features

- Allowed `ProcessedVariable` to handle cases where `len(solution.t)=1` ([#1020](https://github.com/pybamm-team/PyBaMM/pull/1020))
- Added `BackwardIndefiniteIntegral` symbol ([#1014](https://github.com/pybamm-team/PyBaMM/pull/1014))
- Added `plot` and `plot2D` to enable easy plotting of `pybamm.Array` objects ([#1008](https://github.com/pybamm-team/PyBaMM/pull/1008))
- Updated effective current collector models and added example notebook ([#1007](https://github.com/pybamm-team/PyBaMM/pull/1007))
Expand Down Expand Up @@ -41,6 +42,7 @@

## Breaking changes

- For variables discretised using finite elements the result returned by calling `ProcessedVariable` is now transposed ([#1020](https://github.com/pybamm-team/PyBaMM/pull/1020))
- Renamed "surface area density" to "surface area to volume ratio" ([#975](https://github.com/pybamm-team/PyBaMM/pull/975))
- Replaced "reaction rate" with "exchange-current density" ([#975](https://github.com/pybamm-team/PyBaMM/pull/975))
- Changed the implementation of reactions in submodels ([#948](https://github.com/pybamm-team/PyBaMM/pull/948))
Expand Down
9 changes: 9 additions & 0 deletions pybamm/solvers/base_solver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -616,6 +616,15 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
timer.format(solution.total_time),
)
)

# Raise error if solution only contains one timestep (except for algebraic
# solvers, where we may only expect one time in the solution)
if self.algebraic_solver is False and len(solution.t) == 1:
raise pybamm.SolverError(
"Solution time vector has length 1. "
"Check whether simulation terminated too early."
)

return solution

def step(
Expand Down
193 changes: 105 additions & 88 deletions pybamm/solvers/processed_variable.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,26 @@
import scipy.interpolate as interp


def make_interp2D_fun(input, interpolant):
"""
Calls and returns a 2D interpolant of the correct shape depending on the
shape of the input
"""
first_dim, second_dim, _ = input
if isinstance(first_dim, np.ndarray) and isinstance(second_dim, np.ndarray):
first_dim = first_dim[:, 0, 0]
second_dim = second_dim[:, 0]
return interpolant(second_dim, first_dim)
elif isinstance(first_dim, np.ndarray):
first_dim = first_dim[:, 0]
return interpolant(second_dim, first_dim)[:, 0]
elif isinstance(second_dim, np.ndarray):
second_dim = second_dim[:, 0]
return interpolant(second_dim, first_dim)
else:
return interpolant(second_dim, first_dim)[0]


class ProcessedVariable(object):
"""
An object that can be evaluated at arbitrary (scalars or vectors) t and x, and
Expand Down Expand Up @@ -55,48 +75,11 @@ def __init__(self, base_variable, solution, known_evals=None):
and "current collector" in self.domain
and isinstance(self.mesh[0], pybamm.ScikitSubMesh2D)
):
if len(solution.t) == 1:
# space only (steady solution)
self.initialise_2D_fixed_t_scikit_fem()
else:
self.initialise_2D_scikit_fem()
self.initialise_2D_scikit_fem()

# check variable shape
else:
if len(solution.t) == 1:
# Implementing a workaround for 0D and 1D variables. Processing
# variables that are functions of space only needs to be implemented
# properly, see #1006
if (
isinstance(self.base_eval, numbers.Number)
or len(self.base_eval.shape) == 0
or self.base_eval.shape[0] == 1
):
# Scalar value
t = self.t_sol
u = self.u_sol
inputs = {name: inp[0] for name, inp in self.inputs.items()}

entries = self.base_variable.evaluate(t, u, inputs=inputs)

def fun(t):
return entries

self._interpolation_function = fun
self.entries = entries
self.dimensions = 0
else:
# 1D function of space only
n = self.mesh[0].npts
base_shape = self.base_eval.shape[0]
if base_shape in [n, n + 1]:
self.initialise_1D(fixed_t=True)
else:
raise pybamm.SolverError(
"Solution time vector must have length > 1. Check whether "
"simulation terminated too early."
)
elif (
if (
isinstance(self.base_eval, numbers.Number)
or len(self.base_eval.shape) == 0
or self.base_eval.shape[0] == 1
Expand Down Expand Up @@ -141,9 +124,21 @@ def initialise_0D(self):
entries[idx] = self.base_variable.evaluate(t, u, inputs=inputs)

# No discretisation provided, or variable has no domain (function of t only)
self._interpolation_function = interp.interp1d(
self.t_sol, entries, kind="linear", fill_value=np.nan, bounds_error=False
)
if len(self.t_sol) == 1:
# Variable is just a scalar value, but we need to create a callable
# function to be consitent with other processed variables
def fun(t):
return entries

self._interpolation_function = fun
else:
self._interpolation_function = interp.interp1d(
self.t_sol,
entries,
kind="linear",
fill_value=np.nan,
bounds_error=False,
)

self.entries = entries
self.dimensions = 0
Expand Down Expand Up @@ -208,19 +203,22 @@ def initialise_1D(self, fixed_t=False):
self.internal_boundaries = self.mesh[0].internal_boundaries

# set up interpolation
# note that the order of 't' and 'space' is the reverse of what you'd expect
# TODO: fix processing when variable is a function of space only
if fixed_t:
# Hack for 1D space only
if len(self.t_sol) == 1:
# function of space only
interpolant = interp.interp1d(
space, entries_for_interp[:, 0], kind="linear", fill_value=np.nan
)

def interp_fun(t, z):
return interpolant(z)[:, np.newaxis]
if isinstance(z, np.ndarray):
return interpolant(z)[:, np.newaxis]
else:
return interpolant(z)

self._interpolation_function = interp_fun
else:
# function of space and time. Note that the order of 't' and 'space'
# is the reverse of what you'd expect
self._interpolation_function = interp.interp2d(
self.t_sol, space, entries_for_interp, kind="linear", fill_value=np.nan
)
Expand Down Expand Up @@ -297,40 +295,29 @@ def initialise_2D(self):
self.second_dim_pts = second_dim_pts

# set up interpolation
self._interpolation_function = interp.RegularGridInterpolator(
(first_dim_pts, second_dim_pts, self.t_sol),
entries,
method="linear",
fill_value=np.nan,
)

def initialise_2D_fixed_t_scikit_fem(self):
y_sol = self.mesh[0].edges["y"]
len_y = len(y_sol)
z_sol = self.mesh[0].edges["z"]
len_z = len(z_sol)

# Evaluate the base_variable
inputs = {name: inp[0] for name, inp in self.inputs.items()}
if len(self.t_sol) == 1:
# function of space only. Note the order of the points is the reverse
# of what you'd expect
interpolant = interp.interp2d(
second_dim_pts,
first_dim_pts,
entries[:, :, 0],
kind="linear",
fill_value=np.nan,
)

entries = np.reshape(
self.base_variable.evaluate(0, self.u_sol, inputs=inputs), [len_y, len_z]
)
def interp_fun(input):
return make_interp2D_fun(input, interpolant)

# assign attributes for reference
self.entries = entries
self.dimensions = 2
self.y_sol = y_sol
self.z_sol = z_sol
self.first_dimension = "y"
self.second_dimension = "z"
self.first_dim_pts = y_sol
self.second_dim_pts = z_sol

# set up interpolation
self._interpolation_function = interp.interp2d(
y_sol, z_sol, entries, kind="linear", fill_value=np.nan
)
self._interpolation_function = interp_fun
else:
# function of space and time.
self._interpolation_function = interp.RegularGridInterpolator(
(first_dim_pts, second_dim_pts, self.t_sol),
entries,
method="linear",
fill_value=np.nan,
)

def initialise_2D_scikit_fem(self):
y_sol = self.mesh[0].edges["y"]
Expand All @@ -349,11 +336,15 @@ def initialise_2D_scikit_fem(self):
eval_and_known_evals = self.base_variable.evaluate(
t, u, inputs=inputs, known_evals=self.known_evals[t]
)
entries[:, :, idx] = np.reshape(eval_and_known_evals[0], [len_y, len_z])
entries[:, :, idx] = np.reshape(
eval_and_known_evals[0], [len_y, len_z], order="F"
)
self.known_evals[t] = eval_and_known_evals[1]
else:
entries[:, :, idx] = np.reshape(
self.base_variable.evaluate(t, u, inputs=inputs), [len_y, len_z]
self.base_variable.evaluate(t, u, inputs=inputs),
[len_y, len_z],
order="F",
)

# assign attributes for reference
Expand All @@ -367,23 +358,49 @@ def initialise_2D_scikit_fem(self):
self.second_dim_pts = z_sol

# set up interpolation
self._interpolation_function = interp.RegularGridInterpolator(
(y_sol, z_sol, self.t_sol), entries, method="linear", fill_value=np.nan
)
if len(self.t_sol) == 1:
# function of space only. Note the order of the points is the reverse
# of what you'd expect
interpolant = interp.interp2d(
z_sol, y_sol, entries, kind="linear", fill_value=np.nan
)

def interp_fun(input):
return make_interp2D_fun(input, interpolant)

self._interpolation_function = interp_fun
else:
# function of space and time.
self._interpolation_function = interp.RegularGridInterpolator(
(y_sol, z_sol, self.t_sol), entries, method="linear", fill_value=np.nan
)

def __call__(self, t=None, x=None, r=None, y=None, z=None, warn=True):
"""
Evaluate the variable at arbitrary t (and x, r, y and/or z), using interpolation
"""
# If t is None and there is only one value of time in the soluton (i.e.
# the solution is independent of time) then we set t equal to the value
# stored in the solution. If the variable is constant (doesn't depend on
# time) evaluate arbitrarily at the first value of t. Otherwise, raise
# an error
if t is None:
if len(self.t_sol) == 1:
t = self.t_sol
elif self.base_variable.is_constant():
t = self.t_sol[0]
else:
raise ValueError(
"t cannot be None for variable {}".format(self.base_variable)
)

# Call interpolant of correct spatial dimension
if self.dimensions == 0:
out = self._interpolation_function(t)
elif self.dimensions == 1:
out = self.call_1D(t, x, r, z)
elif self.dimensions == 2:
if t is None:
out = self._interpolation_function(y, z)
else:
out = self.call_2D(t, x, r, y, z)
out = self.call_2D(t, x, r, y, z)
if warn is True and np.isnan(out).any():
pybamm.logger.warning(
"Calling variable outside interpolation range (returns 'nan')"
Expand Down
15 changes: 13 additions & 2 deletions tests/unit/test_solvers/test_base_solver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,17 @@ def test_nonmonotonic_teval(self):
):
solver.solve(model, np.array([1, 2, 3, 2]))

def test_solution_time_length_fail(self):
model = pybamm.BaseModel()
v = pybamm.Scalar(1)
model.variables = {"v": v}
solver = pybamm.DummySolver()
t_eval = np.array([0])
with self.assertRaisesRegex(
pybamm.SolverError, "Solution time vector has length 1"
):
solver.solve(model, t_eval)

def test_block_symbolic_inputs(self):
solver = pybamm.BaseSolver(rtol=1e-2, atol=1e-4)
model = pybamm.BaseModel()
Expand Down Expand Up @@ -203,13 +214,13 @@ def algebraic_eval(self, t, y, inputs):
solver.calculate_consistent_state(Model())
solver = pybamm.BaseSolver(root_method="lm")
with self.assertRaisesRegex(
pybamm.SolverError, "Could not find acceptable solution: solver terminated",
pybamm.SolverError, "Could not find acceptable solution: solver terminated"
):
solver.calculate_consistent_state(Model())
# with casadi
solver = pybamm.BaseSolver(root_method="casadi")
with self.assertRaisesRegex(
pybamm.SolverError, "Could not find acceptable solution: .../casadi",
pybamm.SolverError, "Could not find acceptable solution: .../casadi"
):
solver.calculate_consistent_state(Model())

Expand Down
Loading

0 comments on commit bb5f204

Please sign in to comment.