Crashing during compilation of symbolic Jacobian function

[acoh64]

Hi everyone,

I am trying to generate a symbolic Jacobian function for solving the Cahn-Hilliard PDE in time and two spatial dimensions. For small systems (for example Nx=Ny=21), it works great. For larger systems (Nx=Ny=101), the code crashes when running the Jacobian function (I am pretty sure due to memory requirements during compilation). Is there any way to get around this? Code is below:

Δx = 0.01
Δy = 0.01
Nx = length(collect(0:Δx:1))
Ny = length(collect(0:Δy:1))
κ = 0.002
limit(a, N) = a == N+1 ? N-1 : a == 0 ? 2 : a

function CH_loop2D(du,up,p,t)
    u = reshape(up,(Nx,Ny))
    dut = copy(u)
    D = p[1]
    j_func = function (i, j)
        ip1x, im1x, ip1y, im1y = limit(i+1, Nx), limit(i-1, Nx), limit(j+1, Ny), limit(j-1, Ny)
        return log(max(1e-10,u[i,j]/(1-u[i,j]))) + 3.0*(1-2.0*u[i,j]) - κ*(u[i,ip1y]+u[ip1x,j]-4.0*u[i,j]+u[im1x,j]+u[i,im1y])/(Δx*Δy)
    end
    
    for i in 1:Nx
        for j in 1:Ny
            ip1,im1,jp1,jm1 = limit(i+1,Nx), limit(i-1,Nx), limit(j+1,Ny), limit(j-1,Ny)
            
            mu_ij = j_func(i,j)
            mu_im1j = j_func(im1,j)
            mu_ip1j = j_func(ip1,j)
            mu_ijm1 = j_func(i,jm1)
            mu_ijp1 = j_func(i,jp1)
            
            dut[i,j] = D*(((u[ip1,j]-u[im1,j])/(2.0*Δx))*((mu_ip1j-mu_im1j)/(2.0*Δx)) + ((u[i,jp1]-u[i,jm1])/(2.0*Δy))*((mu_ijp1-mu_ijm1)/(2.0*Δy)) + u[i,j]*(mu_ip1j+mu_ijp1-4.0mu_ij+mu_im1j+mu_ijm1)/(Δx*Δy))
        
        end
    end
    du .= dut[:]
    return nothing

end

p = [0.1]
tspan = (0.0,0.1)
u0 = 0.5 .+ 0.1.*rand(Nx*Ny)

@variables u_sym[1:length(u0)] p_sym[1:length(p)] t_sym
u_sym = collect(u_sym)
p_sym = collect(p_sym)
du_sym = similar(u_sym)
CH_loop2D(du_sym,u_sym,p_sym,t_sym)

J_sym = Symbolics.sparsejacobian(du_sym, u_sym)

f! = eval(build_function(du_sym, u_sym, p_sym, t_sym)[2])
J! = eval(build_function(J_sym, u_sym, p_sym, t_sym)[2])

j = Symbolics.value(substitute(J_sym, (Dict(u_sym=>rand(length(u_sym))))))

ut = rand(length(u0))
t = 0.0
p = [0.5]
@time J!(j,ut,p,t)

The kernel crashes so I don't get an error message.

My package versions are:

 Status `~/.julia/environments/v1.7/Project.toml`
  [6e4b80f9] BenchmarkTools v1.3.1
  [13f3f980] CairoMakie v0.8.3
  [41bf760c] DiffEqSensitivity v6.74.0
  [0c46a032] DifferentialEquations v7.1.0
  [7073ff75] IJulia v1.23.3
  [c3572dad] Sundials v4.9.4
  [0c5d862f] Symbolics v4.5.1
  [37e2e46d] LinearAlgebra

[xtalax]

Related issue seen here https://discourse.julialang.org/t/debugging-computer-crash-during-ode-solve/90009/13