NDSOLVE-Error ("For the method \!$$\"IDA\"$$, only machine real code is available. "

I am trying to solve mass transfer equations in an adsorption column to obtain the breakthrough curves. Previously, I asked a question on this and received great help:NDSolve for system of PDEs -Error: fewer dependent variables!

Now, I have changed the equations a little bit to account for the properties of another reactor and material. This is the updated code:

(*From Vitaliy Kaurov for nice display of operators*)
pdConv[f_] :=
f /. Derivative[inds__][g_][vars__] :>
Apply[Defer[D[g[vars], ##]] &,
Transpose[{{vars}, {inds}}] /. {{var_, 0} :>
Sequence[], {var_, 1} :> {var}}]]
(*PDEModels/tutorial/HeatTransfer/HeatTransferVerificationTests#\
463435833*)
ClearAll[HeatTransferModel]
HeatTransferModel[T_, X_List, k_, \[Rho]_, Cp_, Velocity_, Source_] :=
Module[{V, Q, a = k},
V = If[Velocity === "NoFlow",
Q = If[Source === "NoSource", 0, Source];
If[FreeQ[a, _?VectorQ], a = a*IdentityMatrix[Length[X]]];
If[VectorQ[a], a = DiagonalMatrix[a]];
a = PiecewiseExpand[Piecewise[{{-a, True}}]];
Inactive[Div][a.Inactive[Grad][T, X], X] + V - Q]
TimeHeatTransferModel[T_, TimeVar_, X_List, k_, \[Rho]_, Cp_,
Velocity_, Source_] := \[Rho]*Cp*D[T, {TimeVar, 1}] +
HeatTransferModel[T, X, k, \[Rho], Cp, Velocity, Source]
(*Create Parametric PDE Operators for Fluid and Solid*)
(*Include Small Diffusive Term to Fluid Otherwise MMA Might Complain \
parmfop =
TimeHeatTransferModel[c[t, x], t, {x}, {d}, 1,
1, {u}, -(\[Rho] (1 - \[Epsilon])/\[Epsilon]) Q];
parmsop =
TimeHeatTransferModel[q[t, x], t, {x}, {0}, 1, 1, "NoFlow", Q];
parmfop // pdConv
parmsop // pdConv

eps = 0.43;
qflow = 5* (10^-7); (*m^3/s*)
dint = 0.022;(*m*)
Across = \[Pi] * (dint^2)/4;(*m^2*)
usup = qflow/Across;
ueps = usup/eps;
Tgas = 301;(*K*)
Rgas = 8.3145;(*J/mol/K*)
Pgas = 1.02;(*bar*)
yf = 0.2;
c0 = yf*Pgas*101325/(Rgas*Tgas);
KL = 0.226*60;(*1/min*)
rho = 1228.5;
qm = 5.09;(*mol/kg*)
nm = 0.429;
K0 = 4.31* (10^-4);(*1/bar*)
dH = -29.38;(*kJ/mol*)
Keq = (K0*Exp[-dH*1000/(Rgas*Tgas)])/101325;(*1/Pascal*)
tend = 10;(*min*)
xend = 0.171;(*m*)
Qsource = -KL * (q[t,
x] - (qm*
Keq*(c[t, x]*Rgas*
Tgas)/((1 + (Keq*(c[t, x]*Rgas*Tgas))^nm)^(1/nm))));
pdef = (parmfop == 0) /. {\[Epsilon] -> eps, u -> ueps, \[Rho] -> rho,
d -> 1, Q -> Qsource};
pdes = (parmsop == 0) /. {Q -> Qsource};
dc = DirichletCondition[c[t, x] == c0 (1 - Exp[-40 t]), x == 0];
icf = c[0, x] == 0;
ics = q[0, x] == 0;
{cifn, qifn} =
NDSolveValue[{pdef, pdes, dc, icf, ics}, {c, q}, {t, 0, tend}, {x,
0, xend},
Method -> {"MethodOfLines", "TemporalVariable" -> t,
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.5}}}];


However, after running the code I receive this error: NDSolveValue::mconly: For the method IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions.

I checked the equations and I think the main issue is most probably in the Qsource equation, as the code works fine when I replace it with a simple equation. On the other hand, that equation is obtained by fitting isotherm models to experimental data and should be used as it is.

• It's probably (Keq*(c[t, x]*Rgas*Tgas))^nm, because the derivative of c is negative driving c negative: D[c[t, x], t] /. c -> cifn /. {t -> 0, x -> 0.1} Commented Sep 17, 2020 at 18:01
• My first guess is that there's some problem with the equations. I figure you'd know better than me whether they are what they are supposed to be. You should be able to inspect the FEM solution for t == 0. For instance, you can plot the initial values of the time derivative with Plot[ Evaluate[D[c[t, x], t] /. c -> cifn /. {t -> 0}], {x, 0, xend}]. You might try a smaller MaxCellMeasure, if the negative values are from overshoot by the interpolation. Commented Sep 18, 2020 at 16:11