Commit for v0.1.2

Update README.md

Fixed problem with discontinuous initial states and variable voltage functions

Added examples for variable input functions (I, V, P) and
 updating parameters

Gave info on all parameters

Project.toml

@@ -1,7 +1,7 @@
 name = "PETLION"
 uuid = "5e0a28e4-193c-47fa-bbb8-9c901cc1ac2c"
 authors = ["Marc D. Berliner", "Richard D. Braatz"]
-version = "0.1.1"
+version = "0.1.2"
 
 [deps]
 BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"

README.md

@@ -6,7 +6,7 @@
 PETLION is an open-source, high-performance computing implementation of the porous electrode theory (PET) model in Julia. A typical runtime for a dynamic simulation of full charge or discharge with 301 DAEs is 3 ms on a laptop while allocating about 1 MB of total memory. PETLION is built for efficient parameter estimation, controls, and other complex battery simulations using the rigorous PET model.
 
 # Installation
-To add the PETLION package, run the following command in Julia (package will be added to the Julia repo shortly)
+To add the PETLION package, run the following command in Julia
 ```julia
 julia> import Pkg; Pkg.add("PETLION")
 ```

src/PETLION.jl

@@ -20,21 +20,24 @@ import LinearAlgebra
 using BSON: @load, @save
 
 export run_model, run_model!
-export Params
 export model_output
 
+export Params
+export boundary_stop_conditions, options_model, discretizations_per_section, options_numerical
+
 export D_s_eff_isothermal, D_s_eff
 export rxn_rate_isothermal, rxn_rate
 export D_eff_linear, D_eff
 export K_eff, K_eff_isothermal
 export thermodynamic_factor_linear, thermodynamic_factor
 
-export OCV_LCO
-export OCV_LiC6
+export OCV_LCO,  OCV_NMC
+export OCV_LiC6, OCV_LiC6_with_NMC
 
 export rxn_BV, rxn_BV_γMod_01
 export rxn_MHC
 
+
 include("outputs.jl")
 include("structures.jl")
 include("params.jl")

src/checks.jl

@@ -223,10 +223,10 @@ end
 
     value_old = value(run)
     value_new = run.func(t_new,Y,YP,p)
-    run.value .= value_new
 
     if !≈(value_old, value_new, atol=opts.abstol, rtol=opts.reltol)
-        initialize_states!(p,Y,YP,run,opts,funcs,(@inbounds model.SOC[end]))
+        initialize_states!(p,Y,YP,run,opts,funcs,(@inbounds model.SOC[end]); t=t_new)
+        #run.value .= value_new
 
         Sundials.IDAReInit(int.mem, t_new, Y, YP)
     end

src/external.jl

@@ -1,4 +1,33 @@
 ## Functions to complete the param structure
+function Params(cathode::Function;kwargs...)
+    θ = Dict{Symbol,Float64}()
+    funcs = _funcs_numerical()
+
+    anode, system = cathode(θ,funcs)
+    if anode isa Function && system isa Function
+        anode(θ,funcs)
+        θ, bounds, opts, N, numerics = system(θ, funcs, cathode, anode; kwargs...)
+
+        return initialize_param(θ, bounds, opts, N, numerics)
+    else
+        error("Cathode function must return both anode and system functions.")
+    end
+end
+function Params(;cathode=cathode,anode=anode,system=system, # Input chemistry - can be modified
+    kwargs... # keyword arguments for system
+    )
+    θ = Dict{Symbol,Float64}()
+    funcs = _funcs_numerical()
+
+    cathode(θ, funcs)
+    anode(θ, funcs)
+    θ, bounds, opts, N, numerics = system(θ, funcs; kwargs...)
+
+    p = initialize_param(θ, bounds, opts, N, numerics)
+
+    return p
+end
+
 function initialize_param(θ, bounds, opts, _N, numerics)
 
     ind, N_diff, N_alg, N_tot = state_indices(_N, numerics)

src/model_evaluation.jl

@@ -337,31 +337,32 @@ end
 
     return key, key_exists
 end
-@inline function initialize_states!(p::param{T}, Y0::R1, YP0::R1, run::AbstractRun, opts::options_model, funcs::Jac_and_res, SOC::Float64) where {T<:AbstractJacobian,R1<:Vector{Float64}}
+@inline function initialize_states!(p::param{T}, Y0::R1, YP0::R1, run::AbstractRun, opts::options_model, funcs::Jac_and_res, SOC::Float64;kw...) where {T<:AbstractJacobian,R1<:Vector{Float64}}
     if opts.save_start
         key, key_exists = save_start_init!(Y0, run, p, SOC)
 
-        newtons_method!(p,Y0,YP0,run,opts,funcs.R_alg,funcs.R_diff,funcs.J_alg)
+        newtons_method!(p,Y0,YP0,run,opts,funcs.R_alg,funcs.R_diff,funcs.J_alg;kw...)
 
         if !key_exists
             p.cache.save_start_dict[key] = @inbounds Y0[(1:p.N.alg) .+ p.N.diff]
         end
     else
-        newtons_method!(p,Y0,YP0,run,opts,funcs.R_alg,funcs.R_diff,funcs.J_alg)
+        newtons_method!(p,Y0,YP0,run,opts,funcs.R_alg,funcs.R_diff,funcs.J_alg;kw...)
     end
 
     return nothing
 end
 
 @inline function newtons_method!(p::param,Y::R1,YP::R1,run,opts::options_model,R_alg::T1,R_diff::T2,J_alg::T3;
-    itermax::Int64=100, t::Float64=0.0, γ::Float64=0.0
+    itermax::Int64=100, t::Float64=0.0
     ) where {R1<:Vector{Float64},T1<:residual_combined,T2<:residual_combined,T3<:jacobian_combined}
 
     res   = p.cache.res
     Y_old = p.cache.Y_alg
     Y_new = @views @inbounds Y[p.N.diff+1:end]
     YP   .= 0.0
     J     = J_alg.sp
+    γ     = 0.0
 
     # starting loop for Newton's method
     @inbounds for iter in 1:itermax
@@ -455,3 +456,14 @@ end
     end
     return nothing
 end
+@inline function initial_current!(Y0::Vector{Float64},YP0,p,run::run_function{method,func},model::model_output, res_I_guess) where {method<:method_V,func<:Function}
+    @inbounds run.value .= run.func(0.0,Y0,YP0,p)
+    if !isempty(model)
+        @inbounds Y0[end] = model.I[end]
+    else
+        OCV = calc_V(Y0,p)
+        # Arbitrary guess for the initial current. 
+        @inbounds Y0[end] = value(run) > OCV ? +1.0 : -1.0
+    end
+    return nothing
+end