Convergence problem when using Driftfusion code. "Unable to meet integration tolerances" or "Spatial discretization has failed".

I wrote a code where called Driftfusion code. My code works perfectly for a certain input. But when I change the input, making the problem more complex, the code did not work, giving me the following message:

Warning: Failure at t=4.101531e+02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (9.094947e-13) at time t.

> In ode15s (line 717)

In pdepe (line 289)

In df (line 161)

In jumptoV (line 88)

In AA_20230404_Xray_jvst_vsx_Evst (line 136)

Warning: Time integration has failed. Solution is available at requested time points up to t=3.872039e+02.

> In pdepe (line 303)

In df (line 161)

In jumptoV (line 88)

In AA_20230404_Xray_jvst_vsx_Evst (line 136)

I make the code work by increasing the x-mesh. But this "solution" failed for a certain difficulty of the problem to be solved. Increasing the x-mesh for that difficulty gave the following error:

Error using pdepe (line 293)

Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative.

Error in df (line 161)

u = pdepe(par.m,@dfpde,@dfic,@dfbc,x,t,options);

Error in equilibrate (line 69)

sol = df(sol, par);

Error in AA_20230404_Xray_jvst_vsx_Evst (line 79)

soleq_cell1 = equilibrate(par_cell1);

The x-mesh is calculated with the following code, but I have compared the x-mesh that works and the one that doesn't, and I haven't found anything suspicious:

% Error function for each layer
d = par.dcell;

Unable to resolve the name 'par.dcell'.

p = par.layer_points;
dcum0 = par.dcum0;
% For backwards compatibility
if length(par.xmesh_coeff) < length(p)
    xmesh_coeff = 0.7*ones(1, length(p));
else
    xmesh_coeff = par.xmesh_coeff;
end
xcell = cell(1,length(p));
for i = 1:length(p)
    if any(strcmp(par.layer_type{1,i}, {'layer', 'active'}))
        parr = -0.5 : (1/p(i)) : 0.5;
        x_layer = erf(2*pi*xmesh_coeff(i)*parr);
        x_layer = x_layer-x_layer(1); % Subtract base to get zero
        x_layer = x_layer./max(x_layer); % Normalise the funciton
        x_layer = dcum0(i) + x_layer * d(i);
        xcell{i} = x_layer(1:end-1);
    elseif any(strcmp(par.layer_type{1,i}, {'junction', 'interface'}))
        x_layer = linspace(dcum0(i), dcum0(i+1)-(d(i)/p(i)), p(i));
        xcell{i} = x_layer;
    end
end
x = [xcell{:}];
x = [x, dcum0(end)];

[Torsten]

Nothing can be said without executable code that reproduces one of your error messages.

[Agustin]

Can I provide you with the complete code?

[Torsten]

Is this a threat ? :-)

Seriously: The forum usually helps with MATLAB problems, not with numerical problems that arise from a certain physical model. Usually, the above error about the failed spatial discretization results from the boundary condition function you provide for pdepe. So if your code is long and hard to understand, it would be helpful if you could narrow down the cause for failure before posting it.

[Agustin]

Thank you very much for your reply! I understand, I am writing on this forum because it is clearly not a problem with the physical model. The code works perfectly until I increase the discretisation. Or alternatively, "ode15s" does not converge when the solution of the differential equation has a jump too abrupt. Both cases are MATLAB problems.

[Torsten]

No, these are numerical issues that are related to your "jump", not MATLAB problems.

But don't let us split hairs. If you are allowed to include executable code here that reproduces the error quite fast, try it.