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!
“MWE”:
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);