Skip to content

Commit

Permalink
Expanded options
Browse files Browse the repository at this point in the history 
+ Speed increases
+ Custom stop conditions
+ Revamped the file saving/loading system to default to the current working directory
+ Added an option to switch between LU and KLU decomposition during initialization (defaults to KLU)
+ Bug fixes
  • Loading branch information
@MarcBerliner
MarcBerliner committed Dec 27, 2021
1 parent f860371 commit 1030104
Show file tree
Hide file tree
Showing 11 changed files with 132 additions and 69 deletions.
7 changes: 4 additions & 3 deletions Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,15 +1,15 @@
name = "PETLION"
uuid = "5e0a28e4-193c-47fa-bbb8-9c901cc1ac2c"
authors = ["Marc D. Berliner", "Richard D. Braatz"]
version = "0.2.0"
version = "0.2.1"

[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"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
Expand All @@ -18,6 +18,7 @@ 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"

Expand All @@ -35,7 +36,7 @@ SpecialFunctions = "1"
StatsBase = "0.3 - 0.33"
Sundials = "4"
Symbolics = "0.1, 1, 2, 3"
julia = "1, 1.5, 1.6, 1.7"
julia = "1"

[extras]
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
Expand Down
1 change: 1 addition & 0 deletions src/PETLION.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ using RecipesBase
using SpecialFunctions: erf
using SHA: sha1

import SuiteSparse
import Sundials
import LinearAlgebra

Expand Down
19 changes: 11 additions & 8 deletions src/checks.jl
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
@inline function check_simulation_stop!(model, t::Float64, Y, YP, run::AbstractRun, p, bounds, opts::options_model;
@inline function check_simulation_stop!(model, t::Float64, Y, YP, run::AbstractRun, p, bounds, opts::options_model_immutable{T};
ϵ::Float64 = t < 1.0 ? opts.reltol : 0.0,
)
) where T<:Function

if t ≥ run.tf
run.info.flag = 0
Expand All @@ -21,6 +21,7 @@
check_stop_c_e( p, run, model, Y, YP, bounds, ϵ, I)
check_stop_η_plating( p, run, model, Y, YP, bounds, ϵ, I)
check_stop_dfilm( p, run, model, Y, YP, bounds, ϵ, I)
opts.stop_function( p, run, model, Y, YP, bounds, ϵ, I)

return nothing
end
Expand Down Expand Up @@ -220,7 +221,7 @@ end



@inline function check_solve(run::run_constant, model::R1, int::R2, p, bounds, opts::R5, funcs, keep_Y::Bool, iter::Int64, Y::Vector{Float64}, t::Float64) where {R1<:model_output,R2<:Sundials.IDAIntegrator,R5<:options_model}
@inline function check_solve(run::run_constant, model::R1, int::R2, p, bounds, opts::R5, funcs, keep_Y::Bool, iter::Int64, Y::Vector{Float64}, t::Float64) where {R1<:model_output,R2<:Sundials.IDAIntegrator,R5<:AbstractOptionsModel}
if t === int.tprev
# Sometimes the initial step at t = 0 can be too large. This reduces the step size
if t === 0.0
Expand All @@ -245,7 +246,7 @@ end
return true
end

@inline function check_solve(run::run_function, model::R1, int::R2, p::param, bounds::boundary_stop_conditions, opts::R5, funcs, keep_Y::Bool, iter::Int64, Y::Vector{Float64}, t::Float64) where {R1<:model_output,R2<:Sundials.IDAIntegrator,R5<:options_model}
@inline function check_solve(run::run_function, model::R1, int::R2, p::param, bounds::boundary_stop_conditions, opts::R5, funcs, keep_Y::Bool, iter::Int64, Y::Vector{Float64}, t::Float64) where {R1<:model_output,R2<:Sundials.IDAIntegrator,R5<:AbstractOptionsModel}
if iter === opts.maxiters
error("Reached max iterations of $(opts.maxiters) at t = $(int.t)")
elseif within_bounds(run)
Expand All @@ -264,7 +265,7 @@ end
end
end

@inline function check_reinitialization!(model::R1, int::R2, run::R3, p::R4, bounds::R5, opts::R6, funcs) where {R1<:model_output, R2<:Sundials.IDAIntegrator, R3<:AbstractRun,R4<:param,R5<:boundary_stop_conditions,R6<:options_model}
@inline function check_reinitialization!(model::R1, int::R2, run::R3, p::R4, bounds::R5, opts::R6, funcs) where {R1<:model_output, R2<:Sundials.IDAIntegrator, R3<:AbstractRun,R4<:param,R5<:boundary_stop_conditions,R6<:AbstractOptionsModel}
"""
Checking the current function for discontinuities.
If there is a significant change in current after a step size of dt = reltol,
Expand All @@ -273,21 +274,23 @@ end

Y = int.u.v
YP = int.du.v
# take a step of Δt = the relative tolerance
t_new = int.t + opts.reltol

value_old = value(run)
value_new = run.func(t_new,Y,YP,p)

# if the function values at t vs. t + Δt are very different (i.e., there is a discontinuity)
# then reinitialize the DAE at t + Δt
if !≈(value_old, value_new, atol=opts.abstol, rtol=opts.reltol)
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
return nothing
end

@inline function check_errors_parameters_runtime(p::R1,opts::R2) where {R1<:param,R2<:options_model}
@inline function check_errors_parameters_runtime(p::R1,opts::R2) where {R1<:param,R2<:AbstractOptionsModel}
ϵ_sp, ϵ_sn = active_material(p)

if ( ϵ_sp > 1 ) error("ϵ_p + ϵ_fp must be ∈ [0, 1)") end
Expand All @@ -307,4 +310,4 @@ function check_errors_initial(θ, numerics, N)
end

check_is_hold(x::Symbol,model::model_output) = (x===:hold) && (!isempty(model) ? true : error("Cannot use `:hold` without a previous model."))
check_is_hold(x,model) = false
check_is_hold(::Any,::model_output) = false
6 changes: 5 additions & 1 deletion src/external.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -413,7 +413,11 @@ end

function strings_directory_func(N::discretizations_per_section, numerics::T; create_dir=false) where T<:options_numerical

dir_saved_models = joinpath(options[:FILE_DIRECTORY], "saved_models")
dir_saved_models = "saved_models"
# If the file directory is not specified, use the current working directory
if !isnothing(options[:FILE_DIRECTORY])
dir_saved_models = joinpath(options[:FILE_DIRECTORY], dir_saved_models)
end

if create_dir && !isdir(dir_saved_models)
mkdir(dir_saved_models)
Expand Down
5 changes: 3 additions & 2 deletions src/generate_functions.jl
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,8 @@
const PETLION_VERSION = (0,2,0)
const PETLION_VERSION = (0,2,1)
const options = Dict{Symbol,Any}(
:SAVE_SYMBOLIC_FUNCTIONS => true,
:FILE_DIRECTORY => pwd(),
:FILE_DIRECTORY => nothing,
:FACTORIZATION_METHOD => :KLU, # :KLU or :LU
)

function load_functions(p::AbstractParam)
Expand Down
70 changes: 40 additions & 30 deletions src/model_evaluation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,21 +5,23 @@
I = nothing, # constant current C-rate. also accepts :rest
V = nothing, # constant voltage. also accepts :hold if continuing simulation
P = nothing, # constant power. also accepts :hold if continuing simulation
SOC::Number = p.opts.SOC, # initial SOC of the simulation. only valid if not continuing simulation
outputs::R2 = p.opts.outputs, # model output states
abstol = p.opts.abstol, # absolute tolerance in DAE solve
reltol = p.opts.reltol, # relative tolerance in DAE solve
abstol_init = abstol, # absolute tolerance in initialization
reltol_init = reltol, # relative tolerance in initialization
maxiters = p.opts.maxiters, # maximum solver iterations
check_bounds = p.opts.check_bounds, # check if the boundaries (V_min, SOC_max, etc.) are satisfied
reinit = p.opts.reinit, # reinitialize the initial guess
verbose = p.opts.verbose, # print information about the run
interp_final = p.opts.interp_final, # interpolate the final points if a boundary is hit
tstops = p.opts.tstops, # times the solver explicitly solves for
tdiscon = p.opts.tdiscon, # times of known discontinuities in the current function
interp_bc = p.opts.interp_bc, # :interpolate or :extrapolate
save_start = p.opts.save_start, # warm-start for the initial guess
SOC::Number = p.opts.SOC, # initial SOC of the simulation. only valid if not continuing simulation
outputs::R2 = p.opts.outputs, # model output states
abstol = p.opts.abstol, # absolute tolerance in DAE solve
reltol = p.opts.reltol, # relative tolerance in DAE solve
abstol_init = abstol, # absolute tolerance in initialization
reltol_init = reltol, # relative tolerance in initialization
maxiters = p.opts.maxiters, # maximum solver iterations
check_bounds = p.opts.check_bounds, # check if the boundaries (V_min, SOC_max, etc.) are satisfied
reinit = p.opts.reinit, # reinitialize the initial guess
verbose = p.opts.verbose, # print information about the run
interp_final = p.opts.interp_final, # interpolate the final points if a boundary is hit
tstops = p.opts.tstops, # times the solver explicitly solves for
tdiscon = p.opts.tdiscon, # times of known discontinuities in the current function
interp_bc = p.opts.interp_bc, # :interpolate or :extrapolate
save_start = p.opts.save_start, # warm-start for the initial guess
stop_function = p.opts.stop_function,
calc_integrator = p.opts.calc_integrator,
V_max = p.bounds.V_max,
V_min = p.bounds.V_min,
SOC_max = p.bounds.SOC_max,
Expand Down Expand Up @@ -59,7 +61,7 @@
run = get_run(I,V,P,t0,tf,p,model)

# putting opts and bounds into a structure.
opts = options_model(outputs, Float64(SOC), abstol, reltol, abstol_init, reltol_init, maxiters, check_bounds, reinit, verbose, interp_final, tstops, tdiscon, interp_bc, save_start, var_keep)
opts = options_model_immutable(outputs, Float64(SOC), abstol, reltol, abstol_init, reltol_init, maxiters, check_bounds, reinit, verbose, interp_final, tstops, tdiscon, interp_bc, save_start, var_keep, stop_function, calc_integrator)
bounds = boundary_stop_conditions(V_max, V_min, SOC_max, SOC_min, T_max, c_s_n_max, I_max, I_min, η_plating_min, c_e_min, dfilm_max)

# getting the initial conditions and run setup
Expand Down Expand Up @@ -146,7 +148,7 @@ end

@inline within_bounds(run::AbstractRun) = run.info.flag === -1

@inline function initialize_model!(model::model_struct, p::param{jac}, run::T, bounds::boundary_stop_conditions, opts::options_model, res_I_guess=1.0) where {jac<:AbstractJacobian,method<:AbstractMethod,T<:AbstractRun{method,<:Any},model_struct<:Union{model_output,Vector{Float64}}}
@inline function initialize_model!(model::model_struct, p::param{jac}, run::T, bounds::boundary_stop_conditions, opts::AbstractOptionsModel, res_I_guess=nothing) where {jac<:AbstractJacobian,method<:AbstractMethod,T<:AbstractRun{method,<:Any},model_struct<:Union{model_output,Vector{Float64}}}
if !haskey(p.funcs,run)
get_method_funcs!(p,run)
funcs = p.funcs(run)
Expand Down Expand Up @@ -198,13 +200,13 @@ end
return int, funcs, model
end

@inline function retrieve_integrator(run::T, p::param, funcs::Jac_and_res{<:Sundials.IDAIntegrator}, Y0, YP0, opts::options_model, new_run::Bool) where {method<:AbstractMethod,T<:AbstractRun{method,<:Any}}
@inline function retrieve_integrator(run::T, p::param, funcs::Jac_and_res{<:Sundials.IDAIntegrator}, Y0, YP0, opts::AbstractOptionsModel, new_run::Bool) where {method<:AbstractMethod,T<:AbstractRun{method,<:Any}}
"""
If the model has previously been evaluated for a constant run simulation, you can reuse
the integrator function with its cache instead of creating a new one
"""

if isempty(funcs.int)
if isempty(funcs.int) || opts.calc_integrator
int = create_integrator(run, p, funcs, Y0, YP0, opts)
else
# reuse the integrator cache
Expand All @@ -223,7 +225,7 @@ end
return int
end

@inline function create_integrator(run::T, p::param, funcs::Q, Y0, YP0, opts::options_model) where {T<:AbstractRun,Q<:Jac_and_res}
@inline function create_integrator(run::T, p::param, funcs::Q, Y0, YP0, opts::AbstractOptionsModel) where {T<:AbstractRun,Q<:Jac_and_res}
R_full = funcs.R_full
J_full = funcs.J_full
DAEfunc = DAEFunction(
Expand All @@ -242,7 +244,7 @@ end

return int
end
@inline function postfix_integrator!(int::Sundials.IDAIntegrator, run::AbstractRun, opts::options_model, new_run::Bool)
@inline function postfix_integrator!(int::Sundials.IDAIntegrator, run::AbstractRun, opts::AbstractOptionsModel, new_run::Bool)
tstops = int.opts.tstops.valtree
empty!(tstops)

Expand All @@ -258,7 +260,7 @@ end
return nothing
end

@inline function solve!(model,int::R1,run::R2,p,bounds,opts::R3,funcs) where {R1<:Sundials.IDAIntegrator,R2<:AbstractRun,R3<:options_model}
@inline function solve!(model,int::R1,run::R2,p,bounds,opts::R3,funcs) where {R1<:Sundials.IDAIntegrator,R2<:AbstractRun,R3<:AbstractOptionsModel}
keep_Y = opts.var_keep.Y
Y = int.u.v
YP = int.du.v
Expand All @@ -281,10 +283,10 @@ end
return nothing
end

@inline function exit_simulation!(p::R1, model::R2, run::R3, bounds::R4, int::R5, opts::R6; cancel_interp::Bool=false) where {R1<:param,R2<:model_output,R3<:AbstractRun,R4<:boundary_stop_conditions,R5<:Sundials.IDAIntegrator,R6<:options_model}
@inline function exit_simulation!(p::R1, model::R2, run::R3, bounds::R4, int::R5, opts::R6; cancel_interp::Bool=false) where {R1<:param,R2<:model_output,R3<:AbstractRun,R4<:boundary_stop_conditions,R5<:Sundials.IDAIntegrator,R6<:AbstractOptionsModel}
# if a stop condition (besides t = tf) was reached
if !cancel_interp
if opts.interp_final && !(run.info.flag === 0)
if opts.interp_final && !(run.info.flag === 0) && int.t > 1
interp_final_points!(p, model, run, bounds, int, opts)
else
set_var!(model.Y, opts.var_keep.Y ? copy(int.u.v) : int.u.v, opts.var_keep.Y)
Expand Down Expand Up @@ -315,8 +317,13 @@ end
return nothing
end

@views @inbounds @inline function interp_final_points!(p::R1, model::R2, run::R3, bounds::R4, int::R5, opts::R6) where {R1<:param,R2<:model_output,R3<:AbstractRun,R4<:boundary_stop_conditions,R5<:Sundials.IDAIntegrator,R6<:options_model}
YP = opts.var_keep.YP ? model.YP[end] : Float64[]
@views @inbounds @inline function interp_final_points!(p::R1, model::R2, run::R3, bounds::R4, int::R5, opts::R6) where {R1<:param,R2<:model_output,R3<:AbstractRun,R4<:boundary_stop_conditions,R5<:Sundials.IDAIntegrator,R6<:AbstractOptionsModel}
if opts.var_keep.YP
YP = length(model.YP) > 1 ? bounds.t_final_interp_frac.*(model.YP[end] .- model.YP[end-1]) .+ model.YP[end-1] : model.YP[end]
else
YP = Float64[]
end

t = bounds.t_final_interp_frac*(int.t - int.tprev) + int.tprev

set_var!(model.Y, bounds.t_final_interp_frac.*(int.u.v .- model.Y[end]) .+ model.Y[end], opts.var_keep.Y)
Expand All @@ -341,7 +348,7 @@ end

return key, key_exists
end
@inline function initialize_states!(p::param, Y0::T, YP0::T, run::AbstractRun, opts::options_model, funcs::Jac_and_res, SOC::Float64;kw...) where {T<:Vector{Float64}}
@inline function initialize_states!(p::param, Y0::T, YP0::T, run::AbstractRun, opts::AbstractOptionsModel, funcs::Jac_and_res, SOC::Float64;kw...) where {T<:Vector{Float64}}
if opts.save_start
key, key_exists = save_start_init!(Y0, run, p, SOC)

Expand All @@ -357,7 +364,10 @@ 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;

@inline factorization!(L::SuiteSparse.UMFPACK.UmfpackLU{Float64, Int64},A::SparseMatrixCSC{Float64, Int64}) = LinearAlgebra.lu!(L,A)
@inline factorization!(L::KLUFactorization{Float64, Int64},A::SparseMatrixCSC{Float64, Int64}) = klu!(L,A)
@inline function newtons_method!(p::param,Y::R1,YP::R1,run,opts::AbstractOptionsModel,R_alg::T1,R_diff::T2,J_alg::T3;
itermax::Int64=100, t::Float64=0.0
) where {R1<:Vector{Float64},T1<:residual_combined,T2<:residual_combined,T3<:jacobian_combined}

Expand All @@ -367,14 +377,14 @@ end
YP .= 0.0
J = J_alg.sp
γ = 0.0
L = J_alg.L # KLU factorization
L = J_alg.L # 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)
klu!(L, J)
factorization!(L, J)

Y_old .= Y_new
Y_new .-= L\res
Expand Down
23 changes: 12 additions & 11 deletions src/outputs.jl
View file
Original file line number Diff line number Diff line change
@@ -1,20 +1,21 @@
const empty_vector_of_array = VectorOfArray(Array{Vector{Float64}}(Float64[]))
Base.@kwdef struct model_states{vec1D,vec2D,R1}
"""
If you want to add anything to this struct, you must also check/modify set_vars!`.
Otherwise, it may not work as intended
"""
# Matrices (vectors in space and time)
Y::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
YP::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
c_e::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
c_s_avg::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
T::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
film::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
Q::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
j::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
j_s::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
Φ_e::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
Φ_s::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? VectorOfArray(Array{Vector{Float64}}([])) : nothing
Y::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
YP::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
c_e::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
c_s_avg::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
T::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
film::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
Q::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
j::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
j_s::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
Φ_e::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
Φ_s::vec2D = ( vec2D === VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
# Vectors (vectors in time, not space)
I::vec1D = ( vec1D === Array{Float64,1} ) ? Float64[] : nothing
t::vec1D = ( vec1D === Array{Float64,1} ) ? Float64[] : nothing
Expand Down
13 changes: 8 additions & 5 deletions src/physics_equations/auxiliary_states_and_coefficients.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -489,10 +489,13 @@ function limiting_electrode(p::AbstractParam)
θ = p.θ
ϵ_sp, ϵ_sn = active_material(p)

if ϵ_sp*θ[:l_p]*θ[:c_max_p]*(θ[:θ_min_p] - θ[:θ_max_p]) > ϵ_sn*θ[:l_n]*θ[:c_max_n]*(θ[:θ_max_n] - θ[:θ_min_n])
return :p
Q_p = ϵ_sp*θ[:l_p]*θ[:c_max_p]*(θ[:θ_min_p] - θ[:θ_max_p])
Q_n = ϵ_sn*θ[:l_n]*θ[:c_max_n]*(θ[:θ_max_n] - θ[:θ_min_n])

if Q_p > Q_n
return :p, Q_p
else
return :n
return :n, Q_n
end
end

Expand All @@ -501,8 +504,8 @@ end
"""
Calculate the 1C current density (A⋅hr/m²) based on the limiting electrode
"""
F = 96485.3321233

F = const_Faradays
ϵ_sp = 1.0 - (θ[:ϵ_fp] + θ[:ϵ_p])
ϵ_sn = 1.0 - (θ[:ϵ_fn] + θ[:ϵ_n])

Expand Down
Loading

2 comments on commit 1030104

@MarcBerliner
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register()

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/51280

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.2.1 -m "<description of version>" 10301046a2dcc9ac849d4ae8b650cae4e9181741
git push origin v0.2.1

Please sign in to comment.