Solve piecewise PDE system

I am trying to solve a system of PDEs describing the diffusion and reaction of 4 different substances.

The reaction is as follows:

The concentration of E does not matter and is not calculated.

My eqn. system is set up as

as the system is symmetric,

The initial conditions are defined piecewise, so that there are three different regions with differing concentrations initially.

The boundary condition is no flux over the edges of the system:

My code is

clear vars; close all;

DH = 9.31e-5; % cm^2 s^-1
DOH = 5.3e-5; % cm^2 s^-1
DDMP = 2.25e-5; % cm^2 s^-1
DAcetone = 1.28e-5; % cm^2 s^-1

T = 298; % K
global R
R = 3.144; % J K^-1 mol^-1

% concentration guesses
scale_factor = 3;

cDMP_A = 0.4*scale_factor; % mol L^-1
cDMP_B = 0*scale_factor;
cDMP_C = 0*scale_factor;

cOH_A = 0.35*scale_factor; % mol L^-1
cH_A = 1e-14/cOH_A;
cOH_B = 1e-7; % mol L^-1
cH_B = 1e-7;
cH_C = 0.251*scale_factor;
cOH_C = 1e-14/cH_C;

% define mesh
x = 0:0.001:1.379;
t = 0:0.001:60;

sol = pdepe(0, @(x, t, u, dudx)pdefunc(x, t, u, dudx, [DDMP; DAcetone; DH; DOH], T), @(x)pdeic(x, 0.1895, 1.1895, 1.379, cDMP_A, cH_A, cOH_A, cH_C, cOH_C), @pdebc, x, t);

% functions
function [c, f, s] = pdefunc(x, t, u, dudx, D, T)
    
    global R

    c = ones(4, 1);
    f = D.*dudx; 
    s = [-7.32e10.*exp(-46200/(R.*T)).*u(1).*u(3); ...
         7.32e10.*exp(-46200/(R.*T)).*u(1).*u(3); ...
         -1e20.*u(3).*u(4);
         -1e20.*u(3).*u(4)];

end

function u0 = pdeic(x, xa, xb, xc, cDMP_A, cH_A, cOH_A, cH_C, cOH_C)

    if x >= 0 && x <= xa
        u0 = [cDMP_A; ...
              0; ...
              cH_A; ...
              cOH_A];
    elseif x > xa && x <= xb
        u0 = [0; ...
              0; ...
              1e-7; ...
              1e-7];
    elseif x > xb && x <= xc
        u0 = [0; ...
              0; ...
              cH_C; 
              cOH_C];
    end
end

function [pl, ql, pr, qr] = pdebc(xl, ul, xr, ur, t)

    pl = [0; 0; 0; 0];
    pr = [0; 0; 0; 0];
    ql = [1; 1; 1; 1];
    qr = [1; 1; 1; 1];

end

The output I get is

Warning: Failure at t=7.476760e-13.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed
(1.615587e-27) at time t. 

Expected output: Solved PDE.

How can I change my code to let it calculate my concentrations successfully?

[J. Alex Lee]

maybe it's not the issue, but wouldn't be out of the question: can you scale the problem differnetly so that you don't have these enormous numbers as prefactors?

[Aaron Löwenstein]

I think I actually found out that the issue is the large difference in prefactors of diffusion and reaction constants. However, I don't know how to scale them in a way that this is not a problem anymore.

I thought about removing the reaction term for the production of E, as this process is considered to be instantanenous and instead "alter" the solution of C and D in real-time. I would do it in a way that I update the solution for each position after each time step by following the following equation, which I derived for my system. This equation takes into account the instantaneous nature of the reaction of C and D.

A very similar equation holds for the concentration of D, where denotes the updated concentration and K is a constant.

How could I, for every time step, recompute the solution based on the solution the pde solver gave me, to then feed it back to the solver?

[J. Alex Lee]

i think you are now also starting to think also about disparate time scales (stiffness). to focus first on the scale, you'd typically attempt to formulate the dimensionless problem; this also has the advantage of reducing all your physical parameters into fewer dimensionless parameters.

For example you have the Arrhenius term, where

-46200/(R.*T)

should be dimensionless (whatever constant 46200 must have units of J/mol), so you could define the dimensionless temperature as

Tdimless = -46200/(R.*T)

you can try to apply same kind of analysis to your actual equations and variables. For example you can scale your variables by their boundary conditions, and then the constants you have as boundary conditions.

Hopefully at the end of this, your dimensionless variables will be of order unity, or at least of similar orders.

[Aaron Löwenstein]

Thank you, this helped a bit. By making all variables dimensionless, I was able to get all constants closer together in terms of magnitude so that the solver finished without error. I was able to do this by introducing

However, when inspecting the solution, approx. 25% of all concentrations are negative, which is unphysical. Could this be due to the solver? I'm quite certain my equations are correct.

[J. Alex Lee]

Just to make sure you are doing dimensional analysis correctly, when you define , are you replacing every instance of c in all your equations with ?

[Aaron Löwenstein]

Yes, I replaced everything correctly. However I noticed that the solution seems to oscillate, creating high concentrations on either end in alternating fashion.

[J. Alex Lee]

spurious oscillations might indicate discretization problems...with pdepe, do you need to worry about things like CFL condition, or is pdepe supposed to deal with it for you?

[Aaron Löwenstein]

I never heard about the CFL condition before and just read about it. How I understand pdepe is that it uses different timesteps internally in order to try and yield stable solutions and then interpolates the solution at the time steps that one provides. This means that pdepe is supposed to deal with this condition, if I am not mistaken.

[J. Alex Lee]

agree...i also wondered at first if your very sharp steps (in space) in the initial conditions will be problematic...what if you try smoothing out the IC into more sigmoidal shapes...do you get the same oscillations?