How to solve PDAEs (incorporate mass matrix into PDE solver)

I am trying to solve a system of 1st order linear partial differential equations with a constraint on two out of three of my dependent variables (to enforce population conservation).

I am familiar with DAEs and the associated documentation, and I have built my PDE solver by discretizing my spatial derivative using the first method outlined here (Solving Large Stiff Equations · DifferentialEquations.jl), but I have been struggling to implement the mass matrix (or the implicit DAE approach) into my PDE solver.

My first attempt at a mass-matrix approach (with the dimensions num_dependent_vars x num_dependent_vars x num_spatial_points) raised a dimensions-related error in my code (mass matrix dimensions must match the length of my initial condition “vector”, which is num_spatial_points x num_dependent_vars).

My question is: how should I incorporate an algebraic constraint on my dependent variables into my semi-discretized PDE solver? I haven’t been able to find any documentation relevant to this specific problem. I apologize in advance if I have just been missing the obvious. I don’t have much experience creating MWEs but I have attempted to simplify my code into a mostly relevant MWE which produces an error similar (although, bafflingly, not identical) to the one I stated above.

Please let me know if there is anything I can do to clarify this question. I’m also not sure if this is quite the right place to ask it but I was not sure where else to do so. Any assistance or advice would be greatly appreciated!


using DifferentialEquations, LinearAlgebra

xgrid = collect(-pi:0.2:pi)

limit(a,N) = a == N+1 ? 1 : a == 0 ? N : a

# start list of functions
function efield(t,Emax,sig,w0)
    E_t = Emax * exp(-4 * log(2) * (t/sig)^2) * cos(t * w0) ;
    return E_t

function pdesystem(dx,x,xgrid,t)

    u = @view x[:,1] ; # soc
    v = @view x[:,2] ;
    w = @view x[:,3] ;

    knum = length(xgrid)
    T2 = 500
    Eg = 1.5 / 27.21 ; # eV -> AU conversion
    Emax = 0.0014;
    sig = 2687;
    w0 = 0.006;
    Rabi = 9.5 * efield(t,Emax,sig,w0);
    dy = xgrid[3] - xgrid[1]; # for evenly spaced 1D line of x points
    for i = 1:knum
        ip1, im1 = limit(i+1,knum), limit(i-1,knum)
        dx[i,1] = ( E * (u[ip1] - u[im1])/2/dy ) - u[i] / T2  +  Eg[i] * v[i] ;
        dx[i,2] = ( E * (v[ip1] - v[im1])/2/dy ) - Eg[i] * u[i]  -  2 * Rabi * w[i] - v[i] / T2;
        dx[i,3] = ( E * (w[ip1] - w[im1])/2/dy ) + 2 * Rabi * v[i] ;
        dx[i,4] = 1 - 2 * w[i] ;


function init_sbe(xgrid)

    N = length(xgrid);
    u = Array{Float64}(undef,(N,4));
    Threads.@threads for i = 1:N
        u[i,1] = 0
        u[i,2] = 0
        u[i,3] = -1
        u[i,4] = 0
    return u

function massmx(xgrid)
    N = length(xgrid);
    M = Array{Float64}(undef,(N,4,4));
    Threads.@threads for i = 1:N
        M[i,1,1] = 1
        M[i,2,2] = 1
        M[i,3,3] = 1
        M[i,4,4] = 0
    return M

# problem set up
u0 = init_sbe(xgrid);
Nk = length(xgrid);
trange = (-12396.694214876032, 12396.694214876032) ;
M = massmx(xgrid);
f1 = ODEFunction(pdesystem,mass_matrix=M)
prob_sbe = ODEProblem(f1,u0,trange,xgrid) ;
soln1 = solve(prob_sbe,Rodas5(),reltol=1e-10);

The mass matrix should be a matrix whose properties match the vec(u) form of the equation. so for u sized Nx4, you treat this as vec(u) of length 4N and so it should be a 4N x 4N matrix. Then for your case, the top left 3N x 3N matrix will be diagonal, and the other “quadrants” would be 0 matrices

Thank you so much for your response!! This fixes my MWE! I did not realize the ODEFunction would handle the necessary dimensions reshaping automatically, this is very convenient!

When I implement this fix in my actual code, I get ERROR: Singular Exception (1) or ERROR: Singular Exception (18010), so I just want to verify: this means there is some instability in the solver, correct? And thus likely an issue with my model/solver algorithm/tolerance settings rather than a problem with the way I have set up the DAE?

Is your DAE index 1? I.e., are all of your constraints functions of the algebraic variables?

In my actual code (not the mwe), I have several dependent variables, and my constraint on three of those dependent variables is (basically): w(x,t) - u(x,t) - v(x,t) = 0 (i.e., u(x,t)+ v(x,t) = w(x,t)) which I believe is of index 1. Apologies if I am incorrect, I am rather new to all this…

One of those dependent variables is not a differential variable, i.e. it has no equation where it shows up having a derivative in time?