Skip to content

Commit

Permalink
#191 try again to remove ida import
Browse files Browse the repository at this point in the history
  • Loading branch information
martinjrobins committed May 8, 2019
1 parent ea7eb8c commit a4f553c
Showing 1 changed file with 15 additions and 34 deletions.
49 changes: 15 additions & 34 deletions pybamm/solvers/scikits_dae_solver.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,16 +7,14 @@

import numpy as np
import importlib
from scipy.sparse import issparse
import scipy.sparse as sparse

scikits_odes_spec = importlib.util.find_spec("scikits")
if scikits_odes_spec is not None:
scikits_odes_spec = importlib.util.find_spec("scikits.odes")
if scikits_odes_spec is not None:
scikits_odes = importlib.util.module_from_spec(scikits_odes_spec)
scikits_odes_spec.loader.exec_module(scikits_odes)
else:
jac_class = object


class ScikitsDaeSolver(pybamm.DaeSolver):
Expand Down Expand Up @@ -81,10 +79,20 @@ def rootfn(t, y, ydot, return_root):
extra_options = {"old_api": False, "rtol": self.tol, "atol": self.tol}

if jacobian:
# Put the user-supplied Jacobian into the SUNDIALS Class
jacfn = JacobianFunctionIDA()
jacfn.set_jacobian(mass_matrix, jacobian)
extra_options.update({"jacfn": jacfn})
jac_y0_t0 = jacobian(t_eval[0], y0)
if sparse.issparse(jac_y0_t0):
def jacfn(t, y, ydot, residuals, cj, J):
jac_eval = jacobian(t, y) - cj * mass_matrix
J[:][:] = jac_eval.toarray()

else:
def jacfn(t, y, ydot, residuals, cj, J):
jac_eval = jacobian(t, y) - cj * mass_matrix
J[:][:] = jac_eval

extra_options.update({
"jacfn": jacfn
})

if events:
extra_options.update({"rootfn": rootfn, "nr_rootfns": len(events)})
Expand All @@ -98,30 +106,3 @@ def rootfn(t, y, ydot, return_root):

# return solution, we need to tranpose y to match scipy's interface
return sol.values.t, np.transpose(sol.values.y)


class JacobianFunctionIDA(jac_class):
def set_jacobian(self, mass_matrix, jacobian):
"""
Sets the user supplied mass matrix and Jacobian function for the DAE model.
Parameters
----------
mass_matrix : array_like
The (sparse) mass matrix for the chosen spatial method.
jacobian : method
A function that takes in t and y and returns the Jacobian.
"""
self.mass_matrix = mass_matrix
self.jacobian = jacobian

def evaluate(self, t, y, ydot, residuals, cj, return_jacobian):
# scikits_odes requires the full (dense) jacobian
jac_eval = self.jacobian(t, y) - cj * self.mass_matrix
if issparse(jac_eval):
return_jacobian[:][:] = jac_eval.toarray()
else:
return_jacobian[:][:] = jac_eval

return 0

0 comments on commit a4f553c

Please sign in to comment.