Sensitivity analysis for DAEs

[acoh64]

Hi everyone,

I am trying to solve forward sensitivity equations for the 1D Cahn-Hilliard equation (with the eventual goal of moving to 2D). I am getting a StackOverflowError with the following code:

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

function CH_loop(du,u,p,t)

    D = p[1]
    for i in 1:Nx
        ip1,im1 = limit(i+1,Nx), limit(i-1,Nx)
        du[i] = D*(((u[ip1]-u[im1])/(2.0*Δx))*((u[ip1+Nx]-u[im1+Nx])/(2.0*Δx)) + u[i]*(u[ip1+Nx]-2.0*u[i+Nx]+u[im1+Nx])/Δx^2)
        du[i+Nx] =  u[i+Nx] - (log(u[i]/(1-u[i])) + 3.0*(1-2.0*u[i]) - κ*(u[ip1]-2.0*u[i]+u[im1])/Δx^2)
    end

end

function chem_pot_init(c,p)

    D = p[1]
    res = zeros(Nx)
    for i in 1:Nx
        ip1,im1 = limit(i+1,Nx), limit(i-1,Nx)
        res[i] =  log(c[i]/(1-c[i])) + 3.0*(1-2.0*c[i]) - κ*(c[ip1]-2.0*c[i]+c[im1])/Δx^2
    end

    return res

end

p = [0.1]
u0 = 0.5 .+ 0.1.*rand(Nx)
mu0 = chem_pot_init(u0,p)
init = [u0;mu0]
tspan = (0.0,1.0)
M = Diagonal([ones(Nx);zeros(Nx)])

du0 = copy(init)
jac_sparsity = Symbolics.jacobian_sparsity((du,u)->CH_loop(du,u,p,0.0),du0,init)

func1 = ODEFunction(CH_loop,mass_matrix=M,jac_prototype=float.(jac_sparsity))
prob1 = ODEProblem(func1, init, tspan, p)

sol1 = solve(prob1,saveat=collect(tspan[1]:0.1:tspan[2]))
sense_prob = ODEForwardSensitivityProblem(func1,init,tspan,p)

solve(sense_prob,save_everystep=false)

The error is:

ERROR: StackOverflowError:

These are my package versions:

  [41bf760c] DiffEqSensitivity v6.74.0
  [0c46a032] DifferentialEquations v7.1.0
  [961ee093] ModelingToolkit v8.11.3
  [0c5d862f] Symbolics v4.5.1

[ChrisRackauckas]

@YingboMa can you take a quick look at this?

Is there a reason you're using ODEForwardSensitivityProblem rather than duals?

[acoh64]

Thanks! When I have been solving these equations for large systems, implicit solvers have been much faster, so I thought continuous forward or adjoint sensitivity analysis would be faster. Also, I am generally only optimizing ~10-20 parameters so I believe FSA should be most efficient. I have tried using duals as well and get different errors.

[ChrisRackauckas]

res = zeros(Nx) that's why you get errors with dual numbers. Make that res = similar(c) or something so it matches eltypes and you'll be fine.

[YingboMa]

julia> solve(sense_prob,save_everystep=false)
retcode: Success
Interpolation: 1st order linear
t: 2-element Vector{Float64}:
 0.0
 1.0
u: 2-element Vector{Vector{Float64}}:
 [0.5435269377237694, 0.5150956772555779, 0.501326117349611, 0.5890160312636273, 0.572724604965598, 0.595889016206594, 0.5572862102374384, 0.5665907426421062, 0.5349023180669545, 0.5002292722210575  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 [0.1225288248099461, 0.13027393794145153, 0.1546890611139187, 0.1984491465587749, 0.2635202229625326, 0.3487703544745312, 0.4487633159649244, 0.5544124766090672, 0.6552038023623508, 0.7420449312556077  …  0.07085716146817547, 0.10014869004452133, 0.12499537251100165, 0.14417535699731207, 0.1570171914981712, 0.16406148156383946, 0.16699013182361092, 0.16772230320764003, 0.1676638527937975, 0.16756845765983977]

(@v1.7) pkg> st DiffEqSensitivity
      Status `~/.julia/environments/v1.7/Project.toml`
  [41bf760c] DiffEqSensitivity v6.75.0

it works for me.

[acoh64]

I was using Julia v1.6.2 and the problem was fixed after updating to v1.7.2