Error trying to run LFP DFN model: "Could not find acceptable solution", "error above tolerance", "Could not find consistent states"

[ss711]

When trying to solve a DFN model with the LFP_Prada2013 cathode and Li-ion parameter sets I get the following error:
Also doesn't seem to make sense as the error 10^-15 is less than 10^-6?

CasADi - 2021-01-13 12:31:18 WARNING("roots:jac_f_z failed: NaN detected for output jac_x_x0, at nonzero index 78 (row 20, col 20).") [.../casadi/core/oracle_function.cpp:265]
---------------------------------------------------------------------------
SolverError                               Traceback (most recent call last)
~/Documents/pybamm/env/lib/python3.6/site-packages/pybamm/solvers/base_solver.py in calculate_consistent_state(self, model, time, inputs)
    486         try:
--> 487             root_sol = self.root_method._integrate(model, [time], inputs)
    488         except pybamm.SolverError as e:

~/Documents/pybamm/env/lib/python3.6/site-packages/pybamm/solvers/casadi_algebraic_solver.py in _integrate(self, model, t_eval, inputs)
    178                         """.format(
--> 179                             casadi.mmax(casadi.fabs(fun)), self.tol
    180                         )

SolverError: 
                        Could not find acceptable solution: solver terminated
                        successfully, but maximum solution error (5.01361e-16)
                        above tolerance (1e-06)
                        

During handling of the above exception, another exception occurred:

SolverError                               Traceback (most recent call last)
<ipython-input-11-39b3c3af4f00> in <module>
     24 dt = 25
     25 t_eval = np.arange(0, 36000, dt)
---> 26 solution = solver.solve(model, t_eval)

~/Documents/pybamm/env/lib/python3.6/site-packages/pybamm/solvers/base_solver.py in solve(self, model, t_eval, external_variables, inputs)
    582 
    583         # (Re-)calculate consistent initial conditions
--> 584         self._set_initial_conditions(model, ext_and_inputs, update_rhs=True)
    585 
    586         # Non-dimensionalise time

~/Documents/pybamm/env/lib/python3.6/site-packages/pybamm/solvers/base_solver.py in _set_initial_conditions(self, model, inputs, update_rhs)
    458                         y0_from_inputs[:len_rhs], y0_from_model[len_rhs:]
    459                     )
--> 460             model.y0 = self.calculate_consistent_state(model, 0, inputs)
    461 
    462     def calculate_consistent_state(self, model, time=0, inputs=None):

~/Documents/pybamm/env/lib/python3.6/site-packages/pybamm/solvers/base_solver.py in calculate_consistent_state(self, model, time, inputs)
    488         except pybamm.SolverError as e:
    489             raise pybamm.SolverError(
--> 490                 "Could not find consistent states: {}".format(e.args[0])
    491             )
    492         pybamm.logger.info("Found consistent states")

SolverError: Could not find consistent states: 
                        Could not find acceptable solution: solver terminated
                        successfully, but maximum solution error (5.01361e-16)
                        above tolerance (1e-06)

My code:

my_chemistry = pybamm.parameter_sets.Marquis2019
my_chemistry["chemistry"] = "lithium-ion"
my_chemistry['cathode'] = "LFP_Prada2013"   
my_chemistry['anode'] = 'graphite_Chen2020' 
my_chemistry['cell'] = 'A123_Lain2019'  
my_chemistry['electrolyte'] = 'lipf6_Marquis2019'
param = pybamm.ParameterValues(chemistry=my_chemistry)


model = pybamm.lithium_ion.DFN()
geometry = model.default_geometry
parameter_values = param
parameter_values.process_model(model)
parameter_values.process_geometry(geometry)
var_pts = 1000
mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)
disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
disc.process_model(model)
solver = pybamm.CasadiSolver(dt_max=25, rtol=1e-20, atol=1e-20)
dt = 25 
t_eval = np.arange(0, 36000, dt)  
solution = solver.solve(model, t_eval)

[brosaplanella]

The issue might be that the cathode and anode come from Prada 2013 and Chen 2020, respectively, but the initial concentrations are specified in the experiments entry which is left to be the default from Marquis2019.

You can try playing around with the initial concentrations to see if you can find an initial starting concentration. Alternatively, if you want to simulate an LFP & graphite battery and you know a good dataset for it, you could implement it yourself to make sure the parameters are consistent.

[valentinsulzer]

Yes, that is exactly the problem. This isn't showing up in the error message, but the root-finder is returning NaNs because in the exchange-current density c_s_max - c_s_init is negative and so its square root is a complex number.

Unless there is a particular reason for choosing this combination of parameters (is there?), I would suggest just sticking with pybamm.parameter_sets.Prada2013 as the quickest fix.

Short-term fixes we need to do are:

• Clearer error message when root-finder returns NaN
• Error message if initial conditions that are outside the bounds (in this case initial condition for the positive particle is outside its dimensionless bounds (0,1), since c_s_init > c_s_max)

Longer term fix is to provide the initial SOC instead of the initial concentrations, with automatic calculation of initial concentrations based on initial SOC and voltage limits. This is in line with some longer-term development goals but won't be done for several months at least.

[ss711]

Thank you both, I didn't realise initial concentrations were defined in experiments. I was trying to compare the LFP to the Chen and Marquis cathodes so trying to keep as much of the rest of the cell e.g. the graphite anode constant (using Chen2020). Same with electrolyte etc.