Very stiff problem solver

[Ombrini]

When intra-layer reaction rate in CHR2 is too high, the initialization or the solver fails due to the small time scales.
In particular the characteristics time for reaction is 100 times less than the characteristics time of diffusion in the particle.
Any idea on how to solve the problem ?

[d-cogswell]

Does initialization succeed for a smaller applied current? How much overpotential is required?

[Ombrini]

The initialization seems working for small currents (0.1C) but the first time step is not computed by the solver ("Miscellanous exception in function: dae::solver::daeIDASolver::Solve, source file: .\ida_solver.cpp, line: 1159
Sundials IDAS solver cowardly failed to solve the system at TIME = 0.00121158; time horizon [2.30984]; IDA_LSETUP_FAIL
IDAS solver error in module 'IDAS' in function 'IDASolve': At t = 0.00121158, the linear solver setup failed unrecoverably. [IDA_LSETUP_FAIL]")
To be more clear about the code here is an example:
alpha = 0.5
sqrt_eq_K = 1e-5
kinterlayer = 0.0000001
gamma_ts = 1./((1 - c1)**2 * (1 - c2)**2)
interLayerRxn = (kinterlayer/gamma_ts)(sqrt_eq_Kact1R-act2R/sqrt_eq_K)
# interLayerRxn = 0
RxnTerm1 = -interLayerRxn
RxnTerm2 = interLayerRxn
RHS1 += stech_1RxnTerm1
RHS2 += stech_2RxnTerm2

It's a simple modification of the reaction term of the graphite, but in this case the mu_0 of the two phases is different so we have this sqr(K) that is quite high in one direction. This leads to a very high reaction rate which is not well computed if k_0 is more than 1e-7, when physically it should be at least 1e-1 or more. I can try to rethink about the physics of the system, but in general this seems a good occasion to understand the limits of the solver and possible solutions.
From what I know a backward solver should be stable also for very stiff problems, but in this case it seems differently.

[d-cogswell]

Hi @Ombrini, in your example how are the variables sqrt_eq_Kact1R, stech_1RxnTerm1, and stech_2RxnTerm2 defined? Sorry for the delay, finally getting caught up after vacation.

[Ombrini]

Hi @d-cogswell, Copying and pasting the "" was cancelled.
The equations are:
interLayerRxn = (kinterlayer/gamma_ts)(sqrt_eq_Kact1R-act2R/sqrt_eq_K)
RHS1 = stech_1*RxnTerm1 (where stech_1 = 0.5).

It seems to depend on the characteristics time of the system. In particular the ratio between the t0 given by the cathode length and 1/k0, so the reaction time. If this ratio is too big the program seems not able to reduce the stepsize efficiently when the reaction is starting. I guess it's mostly a matter of steptime tentative.
Maybe there are some daetools setting that can help for this kind of problems.
Thank you.

[d-cogswell]

So is this the model you are trying to solve below? Is k0 defined by the variable sqrt_eq_K?

            sqrt_eq_K = 1e-5
            kinterlayer = 0.0000001
            gamma_ts = 1./((1 - c1)**2 * (1 - c2)**2)
            interLayerRxn = (kinterlayer/gamma_ts)*(sqrt_eq_K*act1R-act2R/sqrt_eq_K)
            stech_1 = 0.5
            stech_2 = 0.5
            RxnTerm1 = -interLayerRxn
            RxnTerm2 = interLayerRxn
            RHS1 += stech_1*RxnTerm1
            RHS2 += stech_2*RxnTerm2

[Ombrini]

Yes it is the model,
The k0 is kinterlayer and the sort_eq_K is a constant which multiplies the activity 1 act1R and divide the activity 2 act2R.
Those lines would be posed in the line 261 of mod_electrodes of the current main.
It's basically solving the CHR2 model with a strong interlayer diffusivity (due to the fact that in this case it's an olivine structure).
To be more clear I can rewrite the line as:

interLayerRxn = ((1 - c1)**2 * (1 - c2)**2)kinterlayer(1e-5act1R-1e5act2R)

When act1R will grow a lot, and c1 will be close to 1 the interLayerRxn becomes positive, otherwise it's negative.
This basically state equilibrium between the act1 and act2.

[d-cogswell]

Using smaller error tolerances such as relTol = 1e-8 and absTol=1e-8 helps a bit. My simulation almost finishes but encounters a liner solver error.