Solving PDEs with mass conservation

I want to solve PDEs with mass conservation but have encountered difficulties on implementing mass conservation.

Suppose we have two interacting and diffusing species, y and z, each of which exists in either active or inactive form, denoted by subscripts A and I respectively. Mass conservation ensures that the total amounts Y_Tot := y_A + y_I and Z_Tot := z_A + z_I are preserved.

Therefore, we describe the dynamics of y_A and z_A only, and we express y_I and z_I as y_I = Y_Tot - y_A and z_A = Z_Tot - z_I, where y_A (z_A) represents the overall amount of y_A (z_A) present within the unidimensional domain.

My question is the following: how to integrate the concentration vector over the whole spatial domain (to compute y_A and z_A) in the function pdefun that defines the coefficients of the PDE? Do you have any suggestion on how can the whole concentration vector be accessible at each iteration step within the pdefun function?

Do you have any reference (github, textbook, exercise book, paper) where a similar problem has been solved in Matlab?

Please, find below the code. Thank you in advance for any help!

Main:

%% main_pde_mass_conservation
%  PDEs modelling diffusion with mass conservation

clear all
close all
clc

%% Parameters appearing in the source term of the PDEs
parameters.alpha = 1;
parameters.gamma = 0.7;
parameters.k = 0.5;

% Total concentrations
parameters.Y_tot = 50;
parameters.Z_tot = 20;

% Diffusion coefficients
parameters.D_y = 0.1;
parameters.D_z = 0.1;


%% Simulation

% Time span of integration
T_end = 100;
tspan = linspace(0,T_end,100)';

% Spatial mesh
L_domain = 1;
lmesh = linspace(0,L_domain,100);
parameters.lmesh = lmesh;

% External step input
input_max = 0.1;
ext_input_signal = input_max .* [ones(1,length(lmesh)/2) zeros(1,length(lmesh)/2)];
parameters.ext_input_signal = ext_input_signal;

% Solve equation
m = 0;
sol = pdepe(m,@(l,t,conc,dconc_dl)pdeFun(l,t,conc,dconc_dl,parameters),@pdeIC,@pdeBC,lmesh,tspan);

y_sol = sol(:,:,1);
z_sol = sol(:,:,2);


%% Plot

% Plot simulation results in time and space
figure('units','normalized','outerposition',[0 0 1 1])
Fig1 = figure(1);
surf(lmesh,tspan,y_sol,'FaceColor','interp','EdgeColor','interp','FaceLighting','flat')
xlabel('space')
ylabel('time')
zlabel('y(l,t)')
view([150 25])
% set direction for x axis
ax = gca;
ax.XDir = 'reverse';


figure('units','normalized','outerposition',[0 0 1 1])
Fig2 = figure(2);
surf(lmesh,tspan,z_sol,'FaceColor','interp','EdgeColor','interp')
xlabel('space')
ylabel('time')
zlabel('z(l,t)')
view([150 25])
% set direction for x axis
ax = gca;
ax.XDir = 'reverse';


% Plot species at t=T_end (end of simulation)
figure('units','normalized','outerposition',[0 0 1 1])
Fig3 = figure(3);
p(1) = plot(lmesh,y_sol(end,:),'color','b','Linewidth',1.5);
hold on
p(2) = plot(lmesh,z_sol(end,:),'color','m','Linewidth',1.5);
xlabel('space')
labels = ['y_A'; 'z_A'];
lgd = legend([p(1) p(2)],labels);


%% Function pdefun that defines the coefficients of the PDE
function [coupling_pde,flux_term,source_term] = pdeFun(l,t,conc,dconc_dl,parameters)

    % Parametres
    alpha = parameters.alpha;
    gamma = parameters.gamma;
    k = parameters.k;

    Y_tot = parameters.Y_tot;
    Z_tot = parameters.Z_tot;

    D_y = parameters.D_y;
    D_z = parameters.D_z;

    % External input signal
    ext_input_signal = interp1(parameters.lmesh,parameters.ext_input_signal,l);


    % State variables
    y_A = conc(1);
    z_A = conc(2);

    %% HOW TO COMPUTE y_I and z_I?
    % The implementation below is NOT correct as y_A and z_A are the
    % concentrations in the specific point of the mesh, while I need to
    % integrate the concentration vector over the whole spatial domain
    y_I = Y_tot - y_A;
    z_I= Z_tot - z_A;



    % PDE
    coupling_pde = [1; 1];

    flux_term = [D_y; D_z] .* dconc_dl;

    source_term = [ext_input_signal*y_A + k*y_I*z_I - gamma*z_A*y_A; ... 
        -ext_input_signal*z_A + k*z_I + alpha*y_A];

end


%% Function pdeIC that defines the initial conditions
function conc0 = pdeIC(space)

    if space <= 0.5
        conc0 = [0.5; 0.3];
    else
        conc0 = [0; 0];
    end

end


%% Function pdeBC that defines the boundary conditions
function [p_left,q_left,p_right,q_right] = pdeBC(l_left,conc_left,l_right,conc_right,t)

    p_left = [0; 0];
    q_left = [1; 1];

    p_right = [0; 0];
    q_right = [1; 1];

end

Answer 1

[Bill Greene]

I don't know how to integrate the dependent variables using pdepe. However, I have written a PDE solver with an input syntax similar to pdepe that does support this. It has an option, "vectorized", that results in a single call to the pdefun for all points in the spatial domain. This PDE solver, pde1dM, can be downloaded here.

Here is my version of your MATLAB script that uses pde1dM.

function matlabAnswers12_19_2020
%% main_pde_mass_conservation
%  PDEs modelling diffusion with mass conservation
% clear all
% close all
% clc
%% Parameters appearing in the source term of the PDEs
parameters.alpha = 1;
parameters.gamma = 0.7;
parameters.k = 0.5;
% Total concentrations
parameters.Y_tot = 50;
parameters.Z_tot = 20;
% Diffusion coefficients
parameters.D_y = 0.1;
parameters.D_z = 0.1;
%% Simulation
% Time span of integration
T_end = 100;
tspan = linspace(0,T_end,30)';
% Spatial mesh
L_domain = 1;
lmesh = linspace(0,L_domain,100);
parameters.lmesh = lmesh;
% External step input
input_max = 0.1;
ext_input_signal = input_max .* [ones(1,length(lmesh)/2) zeros(1,length(lmesh)/2)];
parameters.ext_input_signal = ext_input_signal;
% Solve equation
m = 0;
pdef = @(l,t,conc,dconc_dl)pdeFun(l,t,conc,dconc_dl,parameters);
if 0
  sol = pdepe(m,pdef,@pdeIC,@pdeBC,lmesh,tspan);
else
  opts.vectorized='on';
  sol = pde1dM(m,pdef,@pdeIC,@pdeBC,lmesh,tspan,opts);
end
y_sol = sol(:,:,1);
z_sol = sol(:,:,2);
%% Plot
% Plot simulation results in time and space
%figure('units','normalized','outerposition',[0 0 1 1])
figure;
surf(lmesh,tspan,y_sol,'FaceColor','interp','EdgeColor','interp','FaceLighting','flat')
xlabel('space')
ylabel('time')
zlabel('y(l,t)')
view([150 25])
% set direction for x axis
ax = gca;
ax.XDir = 'reverse';
%figure('units','normalized','outerposition',[0 0 1 1])
figure
surf(lmesh,tspan,z_sol,'FaceColor','interp','EdgeColor','interp')
xlabel('space')
ylabel('time')
zlabel('z(l,t)')
view([150 25])
% set direction for x axis
ax = gca;
ax.XDir = 'reverse';
% Plot species at t=T_end (end of simulation)
%figure('units','normalized','outerposition',[0 0 1 1])
figure
p(1) = plot(lmesh,y_sol(end,:),'color','b','Linewidth',1.5);
hold on
p(2) = plot(lmesh,z_sol(end,:),'color','m','Linewidth',1.5);
xlabel('space')
labels = ['y_A'; 'z_A'];
legend([p(1) p(2)],labels);
end

%% Function pdefun that defines the coefficients of the PDE
function [coupling_pde,flux_term,source_term] = pdeFun(l,t,conc,dconc_dl,parameters)
% Parametres
alpha = parameters.alpha;
gamma = parameters.gamma;
k = parameters.k;
Y_tot = parameters.Y_tot;
Z_tot = parameters.Z_tot;
D_y = parameters.D_y;
D_z = parameters.D_z;
% External input signal
ext_input_signal = interp1(parameters.lmesh,parameters.ext_input_signal,l);
% State variables
if 0
  y_A = conc(1,:);
  z_A = conc(2,:);
  %% HOW TO COMPUTE y_I and z_I?
  % The implementation below is NOT correct as y_A and z_A are the
  % concentrations in the specific point of the mesh, while I need to
  % integrate the concentration vector over the whole spatial domain
else
    % integrate the dependent variables over the spatial domain
  y_A=trapz(l, conc(1,:));
  z_A=trapz(l, conc(2,:));
end
y_I = Y_tot - y_A;
z_I= Z_tot - z_A;
nl=length(l); % number of spatial points to evaluate
% PDE
coupling_pde = ones(2,nl);
flux_term = [D_y; D_z] .* dconc_dl;
source_term = [ext_input_signal.*y_A + k*y_I.*z_I - gamma*z_A.*y_A; ...
  -ext_input_signal.*z_A + k*z_I + alpha*y_A];
end
%% Function pdeIC that defines the initial conditions
function conc0 = pdeIC(space)
if space <= 0.5
  conc0 = [0.5; 0.3];
else
  conc0 = [0; 0];
end
end
%% Function pdeBC that defines the boundary conditions
function [p_left,q_left,p_right,q_right] = pdeBC(l_left,conc_left,l_right,conc_right,t)
    p_left = [0; 0];
    q_left = [1; 1];
    p_right = [0; 0];
    q_right = [1; 1];
end

[Irene Zorzan]

Thank you very much for your help, pde1dM seems to pefectly fit my needs!