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

**URL:** <https://discourse.julialang.org/t/how-to-solve-pdaes-incorporate-mass-matrix-into-pde-solver/95301>\
**Category:** Modelling & Simulations\
**Tags:** differentialequation\
**Created:** [February 28, 2023, 3:34am UTC](https://discourse.julialang.org/t/how-to-solve-pdaes-incorporate-mass-matrix-into-pde-solver/95301 "2023-02-28T03:34:38Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![anh](https://avatars.discourse-cdn.com/v4/letter/a/67e7ee/32.png) [@anh](https://discourse.julialang.org/u/anh)\
**Post date:** [February 28, 2023, 3:34am UTC](https://discourse.julialang.org/t/how-to-solve-pdaes-incorporate-mass-matrix-into-pde-solver/95301/1 "2023-02-28T03:34:38Z")

</div>

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](https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/)), 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!

“MWE”:

```julia
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
end

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] ;
    end

end

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
    end
    return u
end

function massmx(xgrid)
    N = length(xgrid);
    M = Array{Float64}(undef,(N,4,4));
    fill!(M,0)
    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
    end
    return M
end

# 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);

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 28, 2023, 10:27am UTC](https://discourse.julialang.org/t/how-to-solve-pdaes-incorporate-mass-matrix-into-pde-solver/95301/2 "2023-02-28T10:27:49Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![anh](https://avatars.discourse-cdn.com/v4/letter/a/67e7ee/32.png) [@anh](https://discourse.julialang.org/u/anh)\
**Post date:** [February 28, 2023, 12:10pm UTC](https://discourse.julialang.org/t/how-to-solve-pdaes-incorporate-mass-matrix-into-pde-solver/95301/3 "2023-02-28T12:10:52Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 28, 2023, 12:23pm UTC](https://discourse.julialang.org/t/how-to-solve-pdaes-incorporate-mass-matrix-into-pde-solver/95301/4 "2023-02-28T12:23:06Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![anh](https://avatars.discourse-cdn.com/v4/letter/a/67e7ee/32.png) [@anh](https://discourse.julialang.org/u/anh)\
**Post date:** [February 28, 2023, 1:05pm UTC](https://discourse.julialang.org/t/how-to-solve-pdaes-incorporate-mass-matrix-into-pde-solver/95301/5 "2023-02-28T13:05:38Z")

</div>

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…

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 28, 2023, 1:23pm UTC](https://discourse.julialang.org/t/how-to-solve-pdaes-incorporate-mass-matrix-into-pde-solver/95301/6 "2023-02-28T13:23:21Z")

</div>

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?