version 0.2.0 update

Replaced the LU linear solve in initialization with a KLU linear solver

Added version control which reevaluates the model files for breaking PETLION versions

Minor changes with the SEI model

Project.toml

@@ -1,14 +1,15 @@
 name = "PETLION"
 uuid = "5e0a28e4-193c-47fa-bbb8-9c901cc1ac2c"
 authors = ["Marc D. Berliner", "Richard D. Braatz"]
-version = "0.1.6"
+version = "0.2.0"
 
 [deps]
 BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 GeneralizedGenerated = "6b9d7cbe-bcb9-11e9-073f-15a7a543e2eb"
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+KLU = "ef3ab10e-7fda-4108-b977-705223b18434"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
@@ -17,7 +18,6 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
-SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
@@ -26,6 +26,7 @@ BSON = "0.3"
 Dierckx = "0.5"
 GeneralizedGenerated = "0.3"
 IfElse = "0.1.0"
+KLU = "0.2"
 RecipesBase = "1"
 RecursiveArrayTools = "2"
 SciMLBase = "1"
@@ -34,7 +35,7 @@ SpecialFunctions = "1"
 StatsBase = "0.3 - 0.33"
 Sundials = "4"
 Symbolics = "0.1, 1, 2, 3"
-julia = "1, 1.5, 1.6"
+julia = "1, 1.5, 1.6, 1.7"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/PETLION.jl

@@ -5,7 +5,8 @@ using SciMLBase: DAEFunction, DAEProblem, step!, init
 using Dierckx: Spline1D
 using GeneralizedGenerated: mk_function, RuntimeFn
 using LinearAlgebra: diagind, Tridiagonal, norm
-using SparseArrays: sparse, findnz, SparseMatrixCSC, spzeros
+using KLU: klu, klu!, KLUFactorization
+using SparseArrays: sparse, findnz, SparseMatrixCSC, spzeros, spdiagm
 using SparseDiffTools: matrix_colors, ForwardColorJacCache, forwarddiff_color_jacobian!
 using RecursiveArrayTools: VectorOfArray
 using Symbolics: @variables, Num, gradient, jacobian_sparsity, expand_derivatives, Differential, get_variables, sparsejacobian, substitute, simplify, build_function, IfElse
@@ -15,7 +16,6 @@ using SHA: sha1
 
 import Sundials
 import LinearAlgebra
-import SuiteSparse
 
 # Must be loaded last
 using BSON: @load, @save

src/generate_functions.jl

@@ -1,4 +1,4 @@
-const PETLION_VERSION = (0,1,6)
+const PETLION_VERSION = (0,2,0)
 const options = Dict{Symbol,Any}(
     :SAVE_SYMBOLIC_FUNCTIONS => true,
     :FILE_DIRECTORY => pwd(),
@@ -24,11 +24,49 @@ function load_functions(p::AbstractParam)
     return funcs
 end
 
+function get_saved_model_version(p::AbstractParam)
+    """
+    Gets the version of PETLION used to generate the saved model
+    """
+    str = strings_directory_func(p)
+    str *= "/info.txt"
+
+    out = readline(str)
+    out = replace(out, "PETLION version: v" => "")
+    # convert this to a tuple
+    out = (Meta.parse.(split(out,"."))...,)
+end
+
+function remove_model_files(p::AbstractParam)
+    """
+    Removes symbolic files for the model
+    """
+    str = strings_directory_func(p) * "/"
+    for x in readdir(str)
+        rm(str*x)
+    end
+end
+
 function load_functions_symbolic(p::AbstractParam)
     dir = strings_directory_func(p) * "/"
     files_exist = isdir(dir)
 
-    if files_exist && options[:SAVE_SYMBOLIC_FUNCTIONS]
+    if files_exist
+        file_version = get_saved_model_version(p)
+
+        # have there been any breaking changes since creating the functions?
+        no_breaking_changes = (file_version[1] == PETLION_VERSION[1]) && (file_version[2] == PETLION_VERSION[2])
+
+        if !no_breaking_changes
+            @warn "Breaking updates encountered: re-evaluating model..."
+            remove_model_files(p)
+        end
+
+    else
+        no_breaking_changes = false
+    end
+
+    if files_exist && no_breaking_changes && options[:SAVE_SYMBOLIC_FUNCTIONS]
         ## residuals
         initial_guess! = include(dir * "initial_guess.jl")
         f_alg!         = include(dir * "f_alg.jl")

src/model_evaluation.jl

@@ -367,14 +367,14 @@ end
     YP   .= 0.0
     J     = J_alg.sp
     γ     = 0.0
-    L     = J_alg.L
+    L     = J_alg.L # KLU factorization
 
     # starting loop for Newton's method
     @inbounds for iter in 1:itermax
         # updating res, Y, and J
         R_alg(res,t,Y,YP,p,run)
         J_alg(t,Y,YP,γ,p,run)
-        LinearAlgebra.lu!(L, J)
+        klu!(L, J)
 
         Y_old .= Y_new
         Y_new .-= L\res

src/params.jl

@@ -404,7 +404,7 @@ function system_NMC_LiC6(θ, funcs, cathode, anode;
     # Effective electrolyte conductivity function
     K_eff = K_eff,
     # Thermodynamic factor, ∂ln(f)/∂ln(c_e)
-    thermodynamic_factor = thermodynamic_factor,
+    thermodynamic_factor = thermodynamic_factor_linear,
     # By default, this
     ## Custon functions
     # Reaction rate equation will use the reaction defined by the cathode

src/physics_equations/numerical_tools.jl

@@ -199,21 +199,21 @@ function interpolate_electrolyte_concetration_fluxes(c_e, p::AbstractParam)
     Δx = Δx_values(p.N)
 
     # Fluxes within the positive electrode
-    @inbounds @views c_e_flux_p = (c_e[2:p.N.p] .- c_e[1:p.N.p-1])/(Δx.p*p.θ[:l_p])
+    @inbounds @views ∂ₓc_e_p = (c_e[2:p.N.p] .- c_e[1:p.N.p-1])/(Δx.p*p.θ[:l_p])
 
     # Fluxes at the separator-positive interface
-    @inbounds @views c_e_flux_ps = (c_e[p.N.p+1] .- c_e[p.N.p]) / ((Δx.p*p.θ[:l_p]/2+Δx.s*p.θ[:l_s]/2))
+    @inbounds @views ∂ₓc_e_ps = (c_e[p.N.p+1] .- c_e[p.N.p]) / ((Δx.p*p.θ[:l_p]/2+Δx.s*p.θ[:l_s]/2))
 
     # Fluxes within the separator
-    @inbounds @views c_e_flux_s = (c_e[p.N.p+2:p.N.p+p.N.s] .- c_e[p.N.p+1:p.N.p+p.N.s-1])/(Δx.s*p.θ[:l_s])
+    @inbounds @views ∂ₓc_e_s = (c_e[p.N.p+2:p.N.p+p.N.s] .- c_e[p.N.p+1:p.N.p+p.N.s-1])/(Δx.s*p.θ[:l_s])
 
     # Fluxes at the separator-negative interface
-    @inbounds @views c_e_flux_sn = (c_e[p.N.p+p.N.s+1] .- c_e[p.N.p+p.N.s]) / ((Δx.n*p.θ[:l_n]/2+Δx.s*p.θ[:l_s]/2))
+    @inbounds @views ∂ₓc_e_sn = (c_e[p.N.p+p.N.s+1] .- c_e[p.N.p+p.N.s]) / ((Δx.n*p.θ[:l_n]/2+Δx.s*p.θ[:l_s]/2))
 
     # Fluxes within the negative electrode
-    @inbounds @views c_e_flux_n = (c_e[p.N.p+p.N.s+2:end] .- c_e[p.N.p+p.N.s+1:end-1])/(Δx.n*p.θ[:l_n])
+    @inbounds @views ∂ₓc_e_n = (c_e[p.N.p+p.N.s+2:end] .- c_e[p.N.p+p.N.s+1:end-1])/(Δx.n*p.θ[:l_n])
 
-    return [c_e_flux_p; c_e_flux_ps], [c_e_flux_s; c_e_flux_sn], c_e_flux_n
+    return [∂ₓc_e_p; ∂ₓc_e_ps], [∂ₓc_e_s; ∂ₓc_e_sn], ∂ₓc_e_n
 end
 
 Δx_values(N) = (p=1/N.p, s=1/N.s, n=1/N.n, a=1/N.a, z=1/N.z)

src/physics_equations/residuals.jl

@@ -249,24 +249,24 @@ function residuals_c_s_avg!(res, states, ∂states, p::T) where {jac,temp,T<:Abs
             c_s = @inbounds @views c_s_avg[ind]
 
             # First order derivatives matrix multiplication
-            ∂ₓc_s_avg = deriv[1](c_s)
+            ∂ᵣc_s = deriv[1](c_s)
 
             # Boundary condition at r = 1
-            ∂ₓc_s_avg[N_r] = -j[i]/D_s_eff[i]*Rp
+            ∂ᵣc_s[N_r] = -j[i]/D_s_eff[i]*Rp
 
             # Boundary condition at r = 0
-            ∂ₓc_s_avg[1] = 0
+            ∂ᵣc_s[1] = 0
 
             # Second order derivatives matrix multiplication
-            ∂ₓₓc_s_avg = deriv[2](c_s)
+            ∂ᵣᵣc_s = deriv[2](c_s)
 
             # Neumann BC at r = 1
-            ∂ₓₓc_s_avg[end] += 50deriv[2].Δx*∂ₓc_s_avg[N_r]*(deriv[2].coeff)
+            ∂ᵣᵣc_s[end] += 50deriv[2].Δx*∂ᵣc_s[N_r]*(deriv[2].coeff)
 
             # Make the RHS vector for this particle
             @inbounds rhsCs[ind] .= (D_s_eff[i]./Rp^2) .* [
-                3∂ₓₓc_s_avg[1]
-                ∂ₓₓc_s_avg[2:end]+2.0./range(1/(N_r-1), 1, length=N_r-1).*∂ₓc_s_avg[2:end]
+                3∂ᵣᵣc_s[1]
+                ∂ᵣᵣc_s[2:end]+2.0./range(1/(N_r-1), 1, length=N_r-1).*∂ᵣc_s[2:end]
                 ]
 
         end
@@ -314,13 +314,13 @@ function residuals_c_s_avg!(res, states, ∂states, p::T) where {jac,temp,T<:Abs
             ind = ind = (1:N_r) .+ N_r*(i-1)
             c_s = c_s_avg[ind]
 
-            ∂ₓc_s_avg = diffusion_matrix*reverse(c_s)
-            ∂ₓc_s_avg[1] = -j[i]*Rp*0.5/D_s_eff[i] # modified BC value due to cheb scheme
-            ∂ₓc_s_avg[end] = 0 # no flux BC
+            ∂ᵣc_s = diffusion_matrix*reverse(c_s)
+            ∂ᵣc_s[1] = -j[i]*Rp*0.5/D_s_eff[i] # modified BC value due to cheb scheme
+            ∂ᵣc_s[end] = 0 # no flux BC
 
-            rhs_numerator = reverse(diffusion_matrix*(4*D_s_eff[i]*((radial_position .+ 1).^2).*∂ₓc_s_avg/(Rp^2)))
+            rhs_numerator = reverse(diffusion_matrix*(4*D_s_eff[i]*((radial_position .+ 1).^2).*∂ᵣc_s/(Rp^2)))
 
-            rhs_limit_vector = (4*D_s_eff[i]/Rp^2)*3*(diffusion_matrix*∂ₓc_s_avg) # limit at r_tilde tends to -1 (at center)
+            rhs_limit_vector = (4*D_s_eff[i]/Rp^2)*3*(diffusion_matrix*∂ᵣc_s) # limit at r_tilde tends to -1 (at center)
 
             @inbounds rhsCs[ind] .= [
                 rhs_limit_vector[end] # L'hopital's rule at the center of the particle
@@ -422,7 +422,11 @@ function residuals_T!(res, states, ∂states, p)
     T_BC_sx =  p.θ[:h_cell]*(p.θ[:T_amb]-T[1])/(Δx.a*p.θ[:l_a])
     T_BC_dx = -p.θ[:h_cell]*(T[end]-p.θ[:T_amb])/(Δx.z*p.θ[:l_z])
 
-    block_tridiag(N) = sparse(Matrix(Tridiagonal{eltype(I_density)}(ones(N-1),-[1;2ones(N-1)],ones(N-1))))
+    block_tridiag(N) = block_tridiag(N) = spdiagm(
+        -1 => ones(eltype(I_density),N-1),
+        0 => -[1;2ones(eltype(I_density),N-2);1],
+        +1 => ones(eltype(I_density),N-1),
+        )
 
     # Positive current collector
     A_a = p.θ[:λ_a].*block_tridiag(p.N.a)
@@ -649,7 +653,7 @@ function residuals_j_s!(res, states, p::AbstractParam)
     η_s = Φ_s.n .- Φ_e.n .- p.θ[:Uref_s] .- F.*j_s.*R_film
     α = 0.5
 
-    j_s_calc = (p.θ[:i_0_jside].*(I_density/I1C)^p.θ[:w]./F).*(exp.((1-α).*F./(R.*T.n).*η_s)-exp.(-α.*F./(R.*T.n).*η_s))
+    j_s_calc = -abs.((p.θ[:i_0_jside].*(I_density/I1C)^p.θ[:w]./F).*(-exp.(-α.*F./(R.*T.n).*η_s)))
 
     # Only activate the side reaction during charge
     j_s_calc .= [IfElse.ifelse(I_density > 0, x, 0) for x in j_s_calc]
@@ -730,7 +734,7 @@ function residuals_Φ_e!(res, states, p::AbstractParam)
     ## Electrolyte fluxes
     # Evaluate the interpolation of the electrolyte concentration fluxes at the
     # edges of the control volumes.
-    c_e_flux_p, c_e_flux_s, c_e_flux_n = PETLION.interpolate_electrolyte_concetration_fluxes(c_e, p)
+    ∂ₓc_e_p, ∂ₓc_e_s, ∂ₓc_e_n = PETLION.interpolate_electrolyte_concetration_fluxes(c_e, p)
 
     ## RHS arrays
     ν_p,ν_s,ν_n = p.numerics.thermodynamic_factor(c_e.p, c_e.s, c_e.n, T.p, T.s, T.n, p)
@@ -739,9 +743,9 @@ function residuals_Φ_e!(res, states, p::AbstractParam)
     K = 2R.*(1-p.θ[:t₊]).*ν[1:end-1]/F
 
     prod_tot = [
-        K̂_eff_p.*T̄_p.*c_e_flux_p./c̄_e_p # p
-        K̂_eff_s.*T̄_s.*c_e_flux_s./c̄_e_s # s
-        K̂_eff_n[1:end-1].*T̄_n.*c_e_flux_n./c̄_e_n # n
+        K̂_eff_p.*T̄_p.*∂ₓc_e_p./c̄_e_p # p
+        K̂_eff_s.*T̄_s.*∂ₓc_e_s./c̄_e_s # s
+        K̂_eff_n[1:end-1].*T̄_n.*∂ₓc_e_n./c̄_e_n # n
     ]
 
     prod_tot[2:end] .-= prod_tot[1:end-1]
@@ -791,7 +795,11 @@ function residuals_Φ_s!(res, states, p::AbstractParam)
     f_p[1]   += -I_density*(Δx.p*p.θ[:l_p]/σ_eff_p)
     f_n[end] += +I_density*(Δx.n*p.θ[:l_n]/σ_eff_n)
 
-    block_tridiag(N) = sparse(Matrix(Tridiagonal{eltype(j.p)}(ones(N-1),-[1;2ones(N-2);1],ones(N-1))))
+    block_tridiag(N) = spdiagm(
+        -1 => ones(eltype(j.p),N-1),
+        0 => -[1;2ones(eltype(j.p),N-2);1],
+        +1 => ones(eltype(j.p),N-1),
+        )
 
     A_p = block_tridiag(p.N.p)
     A_n = block_tridiag(p.N.n)

src/physics_equations/scalar_residual.jl

@@ -13,6 +13,8 @@
 
 @inline calc_T(Y::Vector{<:Number}, p::AbstractParamTemp{true}) = @inbounds @views Y[p.ind.T]
 @inline calc_T(::Vector{<:Number}, p::AbstractParamTemp{false}) = repeat([p.θ[:T₀]], p.N.a+p.N.p+p.N.s+p.N.n+p.N.z)
+@inline calc_K_eff(Y::Vector{<:Number}, p::AbstractParamTemp{true})  = @inbounds @views p.numerics.K_eff(Y[p.ind.c_e.p], Y[p.ind.c_e.s], Y[p.ind.c_e.n], Y[p.ind.T.p], Y[p.ind.T.s], Y[p.ind.T.n], p)
+@inline calc_K_eff(Y::Vector{<:Number}, p::AbstractParamTemp{false}) = @inbounds @views p.numerics.K_eff(Y[p.ind.c_e.p], Y[p.ind.c_e.s], Y[p.ind.c_e.n], repeat([p.θ[:T₀]], p.N.p), repeat([p.θ[:T₀]], p.N.s), repeat([p.θ[:T₀]], p.N.n), p)
 
 @inline method_I(Y, p)   = calc_I(Y,p)
 @inline method_V(Y, p)   = calc_V(Y,p)
@@ -316,9 +318,9 @@ function _get_jacobian_combined(J_sp_base,J_base_func,J_sp_scalar,J_scalar_func,
     J_sp = [J_sp_base; J_sp_scalar']
 
     if lu_decomposition
-        L = LinearAlgebra.lu(J_sp)
+        L = klu(J_sp)
     else
-        L = LinearAlgebra.lu(sparse([1],[1],[1.0]))
+        L = klu(sparse([1],[1],[1.0]))
     end
 
     ind_base   = findall(J_sp.rowval .< size(J_sp,1))

src/structures.jl

@@ -115,7 +115,7 @@ struct jacobian_combined{
     J_scalar::SubArray{Float64, 1, Vector{Float64}, Tuple{Vector{Int64}}, false}
     θ_tot::Vector{Float64}
     θ_keys::Vector{Symbol}
-    L::SuiteSparse.UMFPACK.UmfpackLU{Float64, Int64}
+    L::KLUFactorization{Float64, Int64}
 end
 Base.getindex(J::jacobian_combined,i...) = getindex(J.sp,i...)
 Base.axes(J::jacobian_combined,i...) = axes(J.sp,i...)
@@ -576,6 +576,7 @@ function Base.show(io::IO, run::T) where {T<:AbstractRun}
 
     print(io, str)
 end
+(funcs::model_funcs)(model::model_output) = (@assert !isempty(model); (@inbounds funcs(model.results[end].run)))
 Base.show(io::IO, funcs::model_funcs) = println(io,"PETLION model functions")
 function Base.show(io::IO, model::model_output)
     results = model.results
@@ -650,6 +651,7 @@ end
     print(io, str[1:end-1])
 end
 
+#=
 function _MTK_MatVecProd(A, x; simple::Bool = true)
     """
     Change matrix-vector multiplication in Symbolics to ignore 0s in the matrix.
@@ -682,6 +684,7 @@ end
 # overloads of * to use _MTK_MatVecProd when appropriate
 Base.:*(A::AbstractMatrix, x::AbstractVector{Num}; simple::Bool = true) = _MTK_MatVecProd(A, x; simple=simple)
 Base.:*(A::Union{Array{T,2}, Array{T,2}, Array{T,2}, Array{T,2}}, x::StridedArray{Num, 1}; simple::Bool = true) where {T<:Union{Float32, Float64}} = _MTK_MatVecProd(A, x; simple=simple)
+=#
 
 Base.deleteat!(a::VectorOfArray, i::Integer) = (Base._deleteat!(a.u, i, 1); a.u)