Julia crashes when doing parameter estimation for a system of differential equations in VScode

I have written a battery model in Julia. This model already runs without any issues.
I want to now do parameter estimation w.r.t to testing data using DiffEqFlux.jl . However, whenever I run the diffeqflux.scimltrain command, julia crashes with this message

The terminal process “julia.exe ‘-i’, ‘–banner=no’, ‘–project=C:\Users\madmin1.julia\environments\v1.6’, ‘c:\Users\madmin1.vscode\extensions\julialang.language-julia-1.3.34\scripts\terminalserver\terminalserver.jl’, ‘\.\pipe\vsc-jl-repl-76ed54af-dfa3-4f44-b40b-ceee81f8fccd’, ‘\.\pipe\vsc-jl-cr-2d720729-fad4-47fe-9903-02a83bb0f643’, ‘USE_REVISE=true’, ‘USE_PLOTPANE=true’, ‘USE_PROGRESS=true’, ‘DEBUG_MODE=false’” terminated with exit code: 3221226356.

I have already run everything else. Its the sciml line that’s causing the Julia to crash. I have also tried restarting the Julia language server and vscode too.

I am running Windows Server 2019 and have the latest Visual studio and Julia versions respectively.
‘’’

'''
using DifferentialEquations
using Plots
using DataFrames
using DiffEqFlux, Optim
using CSV
using Flux

plotly()

global const F = 96485.0  # Faraday's constant 
global const R = 8.314462 # J.K⁻¹.mol⁻¹
T = 298 # K

function Δ!(A,N,D,R,Δr)
    """
    Creates the 2nd order FDM matrix for solid phase concentration equation with Neumann boundary conditions
    Comes from the following paper https://drive.google.com/file/d/1czf_P0tt8ViULMVcxTla_u-zaiGCsRE-/view?usp=sharing Page 41
    """
    #Neumann BCs 
    A[1,1] = -D*(R[2] + Δr)/(R[2] *Δr^2)
    A[1,2] = D*(R[2] + Δr)/(R[2] *Δr^2)
    A[N,N-1] = D*(R[N-1] - Δr)/(R[N-1] *Δr^2)
    A[N,N] = -D*(R[N-1] - Δr)/(R[N-1] *Δr^2)

    for i in 2:N-1
        A[i,i] = -2 * D/ (Δr^2)
        A[i,i-1] = D*(R[i] - Δr)/(R[i] *Δr^2)
        A[i,i+1] = D*(R[i] + Δr)/(R[i] *Δr^2)
    end  
end


function U₊_vs_c(θ)
    """
    Returns Cathode equilibrium potential. Source: https://drive.google.com/file/d/16WE-LDKcFBVcLXrl-y-Y5UyfWY1u4GvE/view?usp=sharing 
    """
    return 4.4875 - 0.8090*θ - 0.0428*tanh(18.5138*(θ-0.5542)) - 17.732*tanh(15.789*(θ-0.3117)) + 17.5842*tanh(15.9308*(θ-0.3120))
end

function U₋_vs_c(θ)
    """
    Returns Anode equilibrium potential. Source: https://drive.google.com/file/d/16WE-LDKcFBVcLXrl-y-Y5UyfWY1u4GvE/view?usp=sharing
    """
    return 0.2482 - 0.0909*tanh(29.8538*(θ-0.1234)) - 0.04478*tanh(14.9159*(θ-0.2769)) - 0.0205*tanh(30.4444*(θ-0.6103)) + 1.9793*exp(-39.3631*θ)
end


function Butler_Volmer(c⁺ₛ , c⁻ₛ , p,T=298,I=6.5)
    """
    Returns the cell voltage using Butler_Volmer equation. From paper: https://drive.google.com/file/d/1KwhI43EkDnFEZ2QGybs4o6bq7SOnPugx/view?usp=sharing
    """
    R⁺ₛ ,R⁻ₛ , ε₊ , ε₋ , A , L₊ , L₋ , D₊ , D₋ , k⁺ , k⁻ , C⁰Lᵢ , cₘₐₓ⁺ , cₘₐₓ⁻ ,Rᶠ   = p
    α⁺ = α⁻ = 0.5
    try
        i₀⁺ = k⁺ * sqrt(C⁰Lᵢ * c⁺ₛ * (cₘₐₓ⁺ - c⁺ₛ ))
        i₀⁻ = k⁻ * sqrt(C⁰Lᵢ * c⁻ₛ * (cₘₐₓ⁻ - c⁻ₛ ))
        return (R*T/(α⁺*F)) * asinh(I*R⁺ₛ/(2*ε₊*A*i₀⁺)) - (R*T/(α⁻*F)) * asinh(I*R⁻ₛ/(2*ε₋*A*i₀⁻)) + U₊_vs_c(c⁺ₛ/cₘₐₓ⁺) - U₋_vs_c(c⁻ₛ/cₘₐₓ⁻) - (Rᶠ * I)
    catch
        return 0.0# Cell_parameters["Lower cutoff voltage"]
    end
    
end


#Initial Conditions
Cₛ⁺ = ones(Float64,50) .* 17038.0 # calculated concentrations for Cathode
Cₛ⁻ = ones(Float64,50) .* 29866.0 # calculated concentrations for Anode
dCₜ = 0.0 .* ArrayPartition(Cₛ⁺,Cₛ⁻)      #Intial fluxes are zero
Cₛ = ArrayPartition(Cₛ⁺,Cₛ⁻)


function SPM!(dCₜ,Cₛ,p,t)
    Cell_R⁺ₛ ,Cell_R⁻ₛ , ε₊ , ε₋ , A , L₊ , L₋ , D₊ , D₋ , k⁺ , k⁻ , C⁰Lᵢ , cₘₐₓ⁺ , cₘₐₓ⁻ ,Rᶠ  = p
    
    I = 6.5
    N₊ = 50
    N₋ = 50

    Δr₊ = Cell_R⁺ₛ/(N₊ - 1)
    Δr₋ = Cell_R⁻ₛ/(N₋ -1)
    Rs⁺ₛ = collect(0.0:Δr₊:Cell_R⁺ₛ)
    Rs⁻ₛ = collect(0.0:Δr₋:Cell_R⁻ₛ)
    Acs⁺ = zeros(Float64,N₊,N₊) #FDM matrix for positive solid phase concentration
    Acs⁻ = zeros(Float64,N₋,N₋) #FDM matrix for negative solid phase concentration
    Bcs⁺ = zeros(Float64,N₊,1)
    Bcs⁺[N₊] = (Rs⁺ₛ[N₊ - 1] + Δr₊ )*Cell_R⁺ₛ/(Rs⁺ₛ[N₊-1] * Δr₊ * 3 * ε₊ * F * A * L₊) #Neumann BC source term for positive
    Bcs⁻ = zeros(Float64,N₋,1)
    Bcs⁻[N₋] = -(Rs⁻ₛ[N₋ - 1] + Δr₋ )*Cell_R⁻ₛ/(Rs⁻ₛ[N₋ - 1] * Δr₋ * 3 * ε₋ * F * A * L₋) #Neumann BC source term for negative
    Δ!(Acs⁺,N₊,D₊,Rs⁺ₛ,Δr₊)
    Δ!(Acs⁻,N₋,D₋,Rs⁻ₛ,Δr₋)

    dCₜ[1:N₊] = Acs⁺ * Cₛ[1:N₊] + Bcs⁺ .* I  
    dCₜ[N₊+1:N₋+N₊] = Acs⁻ * Cₛ[N₊+1:N₋+N₊] + Bcs⁻ .* I   
    
end

tspan = (0.0,2300)
p = [5.22e-6,5.86e-6,0.75,0.665,0.1027,7.56e-5,8.52e-5,4.0e-15,3.3e-14,3.42e-6,6.48e-7,1000,63104.0,33133.0,0.1]

SPM_model = ODEProblem(SPM!,Cₛ,tspan,p)
sol = solve(SPM_model,alg_hints = [:stiff])
ts_BV = collect(sol.t[1]:1:sol.t[end])
V1=[]
for t in ts_BV
    push!(V1,Butler_Volmer(sol(t)[50],sol(t)[100],p))
end


plot(ts_BV,V1,label="Initial SPM Prediction")

df = DataFrame(CSV.File("D:\\Shwetank files\\Google Drive\\Life_prediction_data_sets\\CSV_data\\SNL_18650_NCA_25C_0-100_0.5-2C_b_timeseries.csv"))
df = df[!,["Test_Time (s)","Cycle_Index","Current (A)" , "Voltage (V)"]]
dropmissing(df)
filter!(row -> row."Current (A)" < 0.0,df)

Cycle1 = filter(row -> row."Cycle_Index"==6.0,df)
Cycle1[!,"Current (A)"] = Cycle1[!,"Current (A)"].*-1
Cycle1[!,"Test_Time (s)"] = Cycle1[!,"Test_Time (s)"] .- Cycle1[!,"Test_Time (s)"][1]
plot!(Cycle1[!,"Test_Time (s)"],Cycle1[!,"Voltage (V)"],label = "SNL_1C")

function loss(p)
    tmp_prob = remake(SPM_model,p=p)
    tmp_sol = solve(tmp_prob,alg_hints = [:stiff])
    
    V1 = []
    for t in Cycle1[!,"Test_Time (s)"]
        push!(V1,Butler_Volmer(tmp_sol(t)[50],tmp_sol(t)[100],p))
    end
    l = sum(abs2,V1 - Cycle1[!,"Voltage (V)"]) |> sqrt
    return l, tmp_sol
end

pinit = p

cb= function (p,l,pred)
    println("Loss: $l")
    return false
end

res = DiffEqFlux.sciml_train(loss,pinit,ADAM(0.01),cb=cb,maxiters=10)

julia> versioninfo()
Julia Version 1.6.2
Commit 1b93d53fc4 (2021-07-14 15:36 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Xeon(R) Gold 5220 CPU @ 2.20GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, cascadelake)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS =

Try running your code in a normal REPL session outside of VSCode; you might get more info that way.

Anyways, the error code you’re getting apparently indicates heap corruption, so there’s probably not much you can do about it in your code at all.

Try running your code in a normal REPL session outside of VSCode; you might get more info that way.

Just tried and it closed before showing anything.
I will try running the code in a different PC to see if anything changes.

Edit: Nothing changed with my personal PC

Forgot put this:

(@v1.6) pkg> st
      Status `C:\Users\madmin1\.julia\environments\v1.6\Project.toml`
  [c52e3926] Atom v0.12.34
  [aae01518] BandedMatrices v0.16.10
  [a134a8b2] BlackBoxOptim v0.6.0
  [336ed68f] CSV v0.9.1
  [a93c6f00] DataFrames v1.2.2
  [aae7a2af] DiffEqFlux v1.43.0
  [9fdde737] DiffEqOperators v4.31.0
  [0c46a032] DifferentialEquations v6.18.0
  [ced4e74d] DistributionsAD v0.6.29
  [5b8099bc] DomainSets v0.5.4
  [7a1cc6ca] FFTW v1.4.3
  [587475ba] Flux v0.12.6
  [e5e0dc1b] Juno v0.8.4
  [23fbe1c1] Latexify v0.15.6
  [961ee093] ModelingToolkit v6.4.7
  [429524aa] Optim v1.4.1
  [1dea7af3] OrdinaryDiffEq v5.61.1
  [f0f68f2c] PlotlyJS v0.18.7
  [91a5bcdd] Plots v1.21.3
  [d330b81b] PyPlot v2.9.0
  [37e2e3b7] ReverseDiff v1.9.0
  [c3572dad] Sundials v4.5.3
  [0c5d862f] Symbolics v3.2.3
  [fce5fe82] Turing v0.17.4
  [fdbf4ff8] XLSX v0.7.7

I’m curious whether you solved this? I am getting the same exit code, also running fairly simple code.

Nope. I have raised it as an issue on DiffEqFlux.jl github. Link:
Julia crashes during parameter estimation for a system of differential equations

I only recently updated to Julia 1.6.2 and I am thinking it might be a VS Code thing rather than specific packages.

I did try running this in a standalone REPL and it still crashed. I don’t think its a vscode issue

Do try starting a terminal and then start Julia inside. That way you won’t lose the crash output.

1 Like

I ran this on my PC in Atom/Juno and I got a better message:

Please submit a bug report with steps to reproduce this fault, and any error messages that follow (in their entirety). Thanks.
Exception: EXCEPTION_ACCESS_VIOLATION at 0x7ff919a74408 -- RtlSetUserValueHeap at C:\Windows\SYSTEM32\ntdll.dll (unknown line)
in expression starting at C:\Users\verma\OneDrive\Desktop\Modelling\SPM.jl:140
RtlSetUserValueHeap at C:\Windows\SYSTEM32\ntdll.dll (unknown line)
RtlAllocateHeap at C:\Windows\SYSTEM32\ntdll.dll (unknown line)
RtlAllocateHeap at C:\Windows\SYSTEM32\ntdll.dll (unknown line)
malloc at C:\Windows\System32\msvcrt.dll (unknown line)

Please submit a bug report with steps to reproduce this fault, and any error messages that follow (in their entirety). Thanks.
Exception: EXCEPTION_ACCESS_VIOLATION at 0x7ff919a9e3cb --
julia> versioninfo()
Julia Version 1.6.2
Commit 1b93d53fc4 (2021-07-14 15:36 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)      
  CPU: AMD Ryzen 5 3500 6-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, znver2)
Environment:
  JULIA_EDITOR = "C:\Users\verma\AppData\Local\atom\app-1.58.0\atom.exe"  -a
  JULIA_NUM_THREADS = 6

I have finally fixed the issue!
In my ODE function, I had a collect function:

    Δr₊ = Cell_R⁺ₛ/(N₊ - 1)
    Δr₋ = Cell_R⁻ₛ/(N₋ -1)
    Rs⁺ₛ = collect(0.0:Δr₊:Cell_R⁺ₛ)
    Rs⁻ₛ = collect(0.0:Δr₋:Cell_R⁻ₛ)

I replaced the collect function with a for loop:

Δr₊ = Cell_R⁺ₛ/(N₊ - 1)
Δr₋ = Cell_R⁻ₛ/(N₋ -1)

Rs⁺ₛ = zeros(Float64,50,1)
Rs⁻ₛ = zeros(Float64,50,1)
for i in 1:50
 Rs⁺ₛ[i] = Cell_R⁺ₛ * (i-1)
 Rs⁻ₛ[i] = Cell_R⁻ₛ * (i-1)
end

And now it works!

1 Like