Skip to content

Commit

Permalink
#2418 better conditioning for algebraic equations
Browse files Browse the repository at this point in the history
  • Loading branch information
@valentinsulzer
valentinsulzer committed Nov 30, 2022
1 parent 029d0b7 commit 0cf71ad
Show file tree
Hide file tree
Showing 7 changed files with 115 additions and 81 deletions.
4 changes: 2 additions & 2 deletions examples/scripts/compare_lithium_ion.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,12 +3,12 @@
#
import pybamm

pybamm.set_logging_level("SPAM")
pybamm.set_logging_level("INFO")
# load models
models = [
# pybamm.lithium_ion.SPM(),
# pybamm.lithium_ion.SPMe(),
pybamm.lithium_ion.DFN() # {"dimensionality": 1}),
pybamm.lithium_ion.DFN({"dimensionality": 2}),
# pybamm.lithium_ion.DFN(
# {"particle": "uniform profile"}
# ),
Expand Down
2 changes: 1 addition & 1 deletion pybamm/models/submodels/electrode/ohm/full_ohm.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,7 @@ def set_algebraic(self, variables):
"interfacial current densities [A.m-3]"
]

self.algebraic[phi_s] = pybamm.div(i_s) + sum_a_j
self.algebraic[phi_s] = self.param.L_x**2 * (pybamm.div(i_s) + sum_a_j)

def set_boundary_conditions(self, variables):
Domain = self.domain.capitalize()
Expand Down
2 changes: 1 addition & 1 deletion pybamm/models/submodels/electrolyte_conductivity/full_conductivity.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ def set_algebraic(self, variables):
# Override print_name
sum_a_j.print_name = "aj"

self.algebraic = {phi_e: pybamm.div(i_e) - sum_a_j}
self.algebraic = {phi_e: self.param.L_x**2 * (pybamm.div(i_e) - sum_a_j)}

def set_initial_conditions(self, variables):
phi_e = variables["Electrolyte potential [V]"]
Expand Down
166 changes: 92 additions & 74 deletions pybamm/solvers/algebraic_solver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,82 +135,100 @@ def jac_fn(y_alg):
else:
jac_fn = None

# Methods which use least-squares are specified as either "lsq",
# which uses the default method, or with "lsq__methodname"
if self.method.startswith("lsq"):

if self.method == "lsq":
method = "trf"
else:
method = self.method[5:]
if jac_fn is None:
jac_fn = "2-point"
timer.reset()
sol = optimize.least_squares(
root_fun,
y0_alg,
method=method,
ftol=self.tol,
jac=jac_fn,
bounds=model.bounds,
**self.extra_options,
)
integration_time += timer.time()
# Methods which use minimize are specified as either "minimize",
# which uses the default method, or with "minimize__methodname"
elif self.method.startswith("minimize"):
# Adapt the root function for minimize
def root_norm(y):
return np.sum(root_fun(y) ** 2)

if jac_fn is None:
jac_norm = None
itr = 0
maxiter = 2
success = False
while not success:
# Methods which use least-squares are specified as either "lsq",
# which uses the default method, or with "lsq__methodname"
if self.method.startswith("lsq"):

if self.method == "lsq":
method = "trf"
else:
method = self.method[5:]
if jac_fn is None:
jac_fn = "2-point"
timer.reset()
sol = optimize.least_squares(
root_fun,
y0_alg,
method=method,
ftol=self.tol,
jac=jac_fn,
bounds=model.bounds,
**self.extra_options,
)
integration_time += timer.time()
# Methods which use minimize are specified as either "minimize",
# which uses the default method, or with "minimize__methodname"
elif self.method.startswith("minimize"):
# Adapt the root function for minimize
def root_norm(y):
return np.sum(root_fun(y) ** 2)

if jac_fn is None:
jac_norm = None
else:

def jac_norm(y):
return np.sum(2 * root_fun(y) * jac_fn(y), 0)

if self.method == "minimize":
method = None
else:
method = self.method[10:]
extra_options = self.extra_options
if np.any(model.bounds[0] != -np.inf) or np.any(
model.bounds[1] != np.inf
):
bounds = [
(lb, ub) for lb, ub in zip(model.bounds[0], model.bounds[1])
]
extra_options["bounds"] = bounds
timer.reset()
sol = optimize.minimize(
root_norm,
y0_alg,
method=method,
tol=self.tol,
jac=jac_norm,
**extra_options,
)
integration_time += timer.time()
else:

def jac_norm(y):
return np.sum(2 * root_fun(y) * jac_fn(y), 0)

if self.method == "minimize":
method = None
timer.reset()
sol = optimize.root(
root_fun,
y0_alg,
method=self.method,
tol=self.tol,
jac=jac_fn,
options=self.extra_options,
)
integration_time += timer.time()

if sol.success and np.all(abs(sol.fun) < self.tol):
# update initial guess for the next iteration
y0_alg = sol.x
# update solution array
y_alg[:, idx] = y0_alg
success = True
elif not sol.success:
raise pybamm.SolverError(
"Could not find acceptable solution: {}".format(sol.message)
)
else:
method = self.method[10:]
extra_options = self.extra_options
if np.any(model.bounds[0] != -np.inf) or np.any(
model.bounds[1] != np.inf
):
bounds = [
(lb, ub) for lb, ub in zip(model.bounds[0], model.bounds[1])
]
extra_options["bounds"] = bounds
timer.reset()
sol = optimize.minimize(
root_norm,
y0_alg,
method=method,
tol=self.tol,
jac=jac_norm,
**extra_options,
)
integration_time += timer.time()
else:
timer.reset()
sol = optimize.root(
root_fun,
y0_alg,
method=self.method,
tol=self.tol,
jac=jac_fn,
options=self.extra_options,
)
integration_time += timer.time()

if sol.success:
# update solution array
y_alg[:, idx] = sol.x
else:
raise pybamm.SolverError(
"Could not find acceptable solution: {}".format(sol.message)
)
y0_alg = sol.x
if itr > maxiter:
raise pybamm.SolverError(
"Could not find acceptable solution: solver terminated "
"successfully, but maximum solution error "
"({}) above tolerance ({})".format(
np.max(abs(sol.fun)), self.tol
)
)
itr += 1

# Concatenate differential part
y_diff = np.r_[[y0_diff] * len(t_eval)].T
Expand Down
14 changes: 13 additions & 1 deletion pybamm/solvers/casadi_algebraic_solver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -128,7 +128,9 @@ def _integrate(self, model, t_eval, inputs_dict=None):

# If there are no symbolic inputs, check the function is below the tol
# Skip this check if there are symbolic inputs
if success and (not any(np.isnan(fun))):
if success and (
(not any(np.isnan(fun)) and np.all(casadi.fabs(fun) < self.tol))
):
# update initial guess for the next iteration
y0_alg = y_alg_sol
y0 = casadi.vertcat(y0_diff, y0_alg)
Expand All @@ -145,6 +147,16 @@ def _integrate(self, model, t_eval, inputs_dict=None):
raise pybamm.SolverError(
"Could not find acceptable solution: solver returned NaNs"
)
else:
raise pybamm.SolverError(
"""
Could not find acceptable solution: solver terminated
successfully, but maximum solution error ({})
above tolerance ({})
""".format(
casadi.mmax(casadi.fabs(fun)), self.tol
)
)

# Concatenate differential part
y_diff = casadi.horzcat(*[y0_diff] * len(t_eval))
Expand Down
2 changes: 1 addition & 1 deletion pybamm/solvers/casadi_solver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ def __init__(
mode="safe",
rtol=1e-6,
atol=1e-6,
root_method="lm",
root_method="casadi",
root_tol=1e-6,
max_step_decrease_count=5,
dt_max=None,
Expand Down
6 changes: 5 additions & 1 deletion tests/unit/test_solvers/test_base_solver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -203,7 +203,6 @@ def __init__(self):
)
self.convert_to_format = "casadi"
self.bounds = (np.array([-np.inf]), np.array([np.inf]))
self.events = []

def rhs_eval(self, t, y, inputs):
return np.array([])
Expand All @@ -219,6 +218,11 @@ def algebraic_eval(self, t, y, inputs):
"Could not find acceptable solution: The iteration is not making",
):
solver.calculate_consistent_state(Model())
solver = pybamm.BaseSolver(root_method="lm")
with self.assertRaisesRegex(
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(
Expand Down

0 comments on commit 0cf71ad

Please sign in to comment.