diff --git a/Project.toml b/Project.toml index 39aaf29..360287e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "BattMo" uuid = "6f0c0536-3c2c-4762-a987-c605a8a6f898" authors = ["Olav Møyner "] -version = "0.1.1" +version = "0.1.2" [deps] JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" @@ -10,6 +10,7 @@ Jutul = "2b460a1a-8a2b-45b2-b125-b5c536396eb9" MAT = "23992714-dd62-5051-b70f-ba57cb901cac" Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc" @@ -20,5 +21,6 @@ Jutul = "0.2.10" MAT = "0.10.4" Polynomials = "3.2.11" PrecompileTools = "1.1.1" +StaticArrays = "1.5.17" Tullio = "0.3.5" julia = "1.8" diff --git a/src/BattMo.jl b/src/BattMo.jl index ef9822d..c2434e9 100644 --- a/src/BattMo.jl +++ b/src/BattMo.jl @@ -1,5 +1,7 @@ module BattMo using PrecompileTools +using StaticArrays + import JSON import Jutul: number_of_cells, number_of_faces, diff --git a/src/models/activematerial.jl b/src/models/activematerial.jl index fd4716c..acd08fc 100644 --- a/src/models/activematerial.jl +++ b/src/models/activematerial.jl @@ -5,8 +5,8 @@ const ActiveMaterialParameters = JutulStorage abstract type SolidDiffusionDiscretization end -struct P2Ddiscretization <: SolidDiffusionDiscretization - data::Dict{Symbol, Any} +struct P2Ddiscretization{T} <: SolidDiffusionDiscretization + data::T # At the moment the following keys are included : # N::Integer # Discretization size for solid diffusion # rp::Real # Particle radius @@ -18,8 +18,8 @@ end struct NoParticleDiffusion <: SolidDiffusionDiscretization end -struct ActiveMaterial{D} <: ElectroChemicalComponent where {D<:SolidDiffusionDiscretization} - params::ActiveMaterialParameters +struct ActiveMaterial{D, T} <: ElectroChemicalComponent where {D<:SolidDiffusionDiscretization, T<:ActiveMaterialParameters} + params::T # At the moment the following keys are include # - ocp_func::F where {F <: Function} # - n_charge_carriers @@ -30,6 +30,9 @@ struct ActiveMaterial{D} <: ElectroChemicalComponent where {D<:SolidDiffusionDis discretization::D end +const ActiveMaterialP2D{D, T} = ActiveMaterial{D, T} +const ActiveMaterialNoParticleDiffusion{T} = ActiveMaterial{NoParticleDiffusion, T} + struct Ocp <: ScalarVariable end struct DiffusionCoef <: ScalarVariable end struct ReactionRateConst <: ScalarVariable end @@ -46,16 +49,18 @@ Jutul.local_discretization(::SolidDiffusionBc, i) = nothing const ActiveMaterialModel = SimulationModel{O, S} where {O<:JutulDomain, S<:ActiveMaterial} ## Create ActiveMaterial with full p2d solid diffusion -function ActiveMaterial{P2Ddiscretization}(params::ActiveMaterialParameters, rp, N, D) +function ActiveMaterialP2D(params::ActiveMaterialParameters, rp, N, D) data = setupSolidDiffusionDiscretization(rp, N, D) discretization = P2Ddiscretization(data) - return ActiveMaterial{P2Ddiscretization}(params, discretization) + params = Jutul.convert_to_immutable_storage(params) + return ActiveMaterial{typeof(discretization), typeof(params)}(params, discretization) end ## Create ActiveMaterial with no solid diffusion -function ActiveMaterial{NoParticleDiffusion}(params::ActiveMaterialParameters) +function ActiveMaterialNoParticleDiffusion(params::ActiveMaterialParameters) discretization = NoParticleDiffusion() - return ActiveMaterial{NoParticleDiffusion}(params, discretization) + params = Jutul.convert_to_immutable_storage(params) + return ActiveMaterial{NoParticleDiffusion, typeof(params)}(params, discretization) end @@ -68,11 +73,11 @@ function Base.getindex(disc::P2Ddiscretization, key::Symbol) return disc.data[key] end -function discretisation_type(system::ActiveMaterial{P2Ddiscretization}) +function discretisation_type(system::ActiveMaterialP2D) return :P2Ddiscretization end -function discretisation_type(system::ActiveMaterial{NoParticleDiffusion}) +function discretisation_type(system::ActiveMaterialNoParticleDiffusion) return :NoParticleDiffusion end @@ -80,7 +85,7 @@ function discretisation_type(model::ActiveMaterialModel) discretisation_type(model.system) end -function solid_diffusion_discretization_number(system::ActiveMaterial{P2Ddiscretization}) +function solid_diffusion_discretization_number(system::ActiveMaterialP2D) return system.discretization[:N] end @@ -108,7 +113,7 @@ function setupSolidDiffusionDiscretization(rp, N, D) for i = 1 : N + 1 hT[i] = (4*pi*rf[i]^2/(dr/2)) end - + div = Vector{Tuple{Int64, Int64, Float64}}(undef, 2*(N - 1)) k = 1 @@ -118,7 +123,10 @@ function setupSolidDiffusionDiscretization(rp, N, D) div[k] = (j + 1, j, -1) k += 1 end - + hT = SVector(hT...) + div = SVector(div...) + vols = SVector(vols...) + data = Dict(:N => N , :D => D , :rp => rp , @@ -126,9 +134,8 @@ function setupSolidDiffusionDiscretization(rp, N, D) :vols => vols, :div => div , ) - - return data - + # Convert to concrete type + return NamedTuple(pairs(data)) end ##################################################### @@ -136,7 +143,7 @@ end #################################################### function select_primary_variables!(S, - system::ActiveMaterial{P2Ddiscretization}, + system::ActiveMaterialP2D, model::SimulationModel ) S[:Phi] = Phi() @@ -156,7 +163,7 @@ function degrees_of_freedom_per_entity(model::ActiveMaterialModel, end function select_parameters!(S, - system::ActiveMaterial{P2Ddiscretization}, + system::ActiveMaterialP2D, model::SimulationModel) S[:Temperature] = Temperature() @@ -166,7 +173,7 @@ function select_parameters!(S, end function select_secondary_variables!(S, - system::ActiveMaterial{P2Ddiscretization}, + system::ActiveMaterialP2D, model::SimulationModel ) S[:Charge] = Charge() @@ -177,7 +184,7 @@ function select_secondary_variables!(S, end function select_equations!(eqs, - system::ActiveMaterial{P2Ddiscretization}, + system::ActiveMaterialP2D, model::SimulationModel ) @@ -197,7 +204,7 @@ function Jutul.number_of_equations_per_entity(model::ActiveMaterialModel, ::Soli end function select_minimum_output_variables!(out , - system::ActiveMaterial{P2Ddiscretization}, + system::ActiveMaterialP2D, model::SimulationModel ) push!(out, :Charge) @@ -209,10 +216,10 @@ end @jutul_secondary( function update_vocp!(Ocp, tv::Ocp, - model:: SimulationModel{<:Any, ActiveMaterial{P2Ddiscretization}, <:Any, <:Any}, + model:: SimulationModel{<:Any, ActiveMaterialP2D{D, T}, <:Any, <:Any}, Cs, ix - ) + ) where {D, T} ocp_func = model.system.params[:ocp_func] cmax = model.system.params[:maximum_concentration] refT = 298.15 @@ -226,10 +233,10 @@ end @jutul_secondary( function update_reaction_rate!(ReactionRateConst , tv::ReactionRateConst , - model::SimulationModel{<:Any, ActiveMaterial{P2Ddiscretization}, <:Any, <:Any}, + model::SimulationModel{<:Any, ActiveMaterialP2D{D, T}, <:Any, <:Any}, Cs , ix - ) + ) where {D, T} rate_func = model.system.params[:reaction_rate_constant_func] refT = 298.15 for cell in ix @@ -241,9 +248,9 @@ end @jutul_secondary( function update_solid_diffusion_flux!(SolidDiffFlux, tv::SolidDiffFlux, - model::SimulationModel{<:Any, ActiveMaterial{P2Ddiscretization}, <:Any, <:Any}, + model::SimulationModel{<:Any, ActiveMaterialP2D{D, T}, <:Any, <:Any}, Cp, - ix) + ix) where {D, T} s = model.system for cell in ix @inbounds @views update_solid_flux!(SolidDiffFlux[:, cell], Cp[:, cell], s) @@ -252,7 +259,7 @@ end ) -function update_solid_flux!(flux, Cp, system::ActiveMaterial{P2Ddiscretization}) +function update_solid_flux!(flux, Cp, system::ActiveMaterialP2D) # compute lithium flux in particle, using harmonic averaging. At the moment D has a constant value within particle # but this is going to change. @@ -286,15 +293,15 @@ function Jutul.update_equation_in_entity!(eq_buf , ldisc = Nothing) disc = model.system.discretization - N = disc[:N] + N = length(eq_buf) hT = disc[:hT] vols = disc[:vols] div = disc[:div] D = disc[:D] - Cp = state.Cp[:, self_cell] - Cp0 = state0.Cp[:, self_cell] - flux = state.SolidDiffFlux[:, self_cell] + Cp = @views state.Cp[:, self_cell] + Cp0 = @views state0.Cp[:, self_cell] + flux = @views state.SolidDiffFlux[:, self_cell] Cs = state.Cs[self_cell] for i = 1 : N @@ -339,7 +346,7 @@ end function select_primary_variables!(S, - system::ActiveMaterial{NoParticleDiffusion}, + system::ActiveMaterialNoParticleDiffusion, model::SimulationModel ) S[:Phi] = Phi() @@ -348,7 +355,7 @@ function select_primary_variables!(S, end function select_secondary_variables!(S, - system::ActiveMaterial{NoParticleDiffusion}, + system::ActiveMaterialNoParticleDiffusion, model::SimulationModel ) @@ -360,7 +367,7 @@ function select_secondary_variables!(S, end function select_parameters!(S, - system::ActiveMaterial{NoParticleDiffusion}, + system::ActiveMaterialNoParticleDiffusion, model::SimulationModel) S[:Temperature] = Temperature() @@ -370,7 +377,7 @@ function select_parameters!(S, end function select_equations!(eqs, - system::ActiveMaterial{NoParticleDiffusion}, + system::ActiveMaterialNoParticleDiffusion, model::SimulationModel ) @@ -382,7 +389,7 @@ end function select_minimum_output_variables!(out , - system::ActiveMaterial{NoParticleDiffusion}, + system::ActiveMaterialNoParticleDiffusion, model::SimulationModel ) push!(out, :Charge) @@ -397,10 +404,10 @@ end @jutul_secondary( function update_vocp!(Ocp , tv::Ocp , - model::M, + model::SimulationModel{<:Any, ActiveMaterialNoParticleDiffusion{T}, <:Any, <:Any}, C , ix - ) where {M <: SimulationModel{<:Any, ActiveMaterial{NoParticleDiffusion}, <:Any, <:Any}} + ) where T ocp_func = model.system.params[:ocp_func] cmax = model.system.params[:maximum_concentration] @@ -417,10 +424,10 @@ end @jutul_secondary( function update_reaction_rate!(ReactionRateConst , tv::ReactionRateConst, - model::SimulationModel{<:Any, ActiveMaterial{NoParticleDiffusion}, <:Any, <:Any}, + model::SimulationModel{<:Any, ActiveMaterialNoParticleDiffusion{T}, <:Any, <:Any}, C , ix - ) + ) where T rate_func = model.system.params[:reaction_rate_constant_func] refT = 298.15 for i in ix diff --git a/src/models/battery_cross_terms.jl b/src/models/battery_cross_terms.jl index 16addfd..2152c2c 100644 --- a/src/models/battery_cross_terms.jl +++ b/src/models/battery_cross_terms.jl @@ -87,13 +87,13 @@ function reaction_rate(phi_a , n = activematerial.params[:n_charge_carriers] cmax = activematerial.params[:maximum_concentration] vsa = activematerial.params[:volumetric_surface_area] - + eta = (phi_a - phi_e - ocp) th = 1e-3*cmax j0 = R0*regularized_sqrt(c_e*(cmax - c_a)*c_a, th)*n*FARADAY_CONST R = vsa*butler_volmer_equation(j0, 0.5, n, eta, T) - return R./(n*FARADAY_CONST) + return R/(n*FARADAY_CONST) end @@ -159,9 +159,7 @@ function Jutul.update_cross_term_in_entity!(out , @assert cs == :Charge v = 1.0*vols*R*n*FARADAY_CONST end - out[] = -v - end diff --git a/src/models/elyte.jl b/src/models/elyte.jl index f4fc0b4..a817dbc 100644 --- a/src/models/elyte.jl +++ b/src/models/elyte.jl @@ -154,7 +154,7 @@ end t = transference(sys) F = FARADAY_CONST for i in ix - @inbounds chemCoef[i] = 1/F*(1 - t)*Conductivity[i]*2*DmuDc[i] + @inbounds chemCoef[i] = 1.0/F*(1.0 - t)*Conductivity[i]*2.0*DmuDc[i] end end @@ -172,10 +172,10 @@ function Jutul.face_flux!(::T, c, other, face, face_sign, eq::ConservationLaw{:C end -function Jutul.face_flux!(::T, c, other, face, face_sign, eq::ConservationLaw{:Mass, <:Any}, state, model::ElectrolyteModel, dt, flow_disc) where T +function Jutul.face_flux!(q::T, c, other, face, face_sign, eq::ConservationLaw{:Mass, <:Any}, state, model::ElectrolyteModel, dt, flow_disc) where T t = transference(model.system) - z = 1 + z = 1.0 F = FARADAY_CONST @inbounds trans = state.ECTransmissibilities[face] @@ -187,7 +187,5 @@ function Jutul.face_flux!(::T, c, other, face, face_sign, eq::ConservationLaw{:M j = j - jchem*(1.0) massFlux = diffFlux + t/(z*F)*j - - return T(massFlux) - + return setindex(q, massFlux, 1)::T end diff --git a/src/mrst_test_utils.jl b/src/mrst_test_utils.jl index 5fa2e7c..074d663 100644 --- a/src/mrst_test_utils.jl +++ b/src/mrst_test_utils.jl @@ -123,7 +123,7 @@ function setup_model_1d(jsondict; use_groups = false, kwarg...) include_cc = true - model = setup_battery_model_1d(jsondict, include_cc = include_cc) + model = setup_battery_model_1d(jsondict, include_cc = include_cc; kwarg...) parameters = setup_battery_parameters_1d(jsondict, model) initState = setup_battery_initial_state_1d(jsondict, model) @@ -156,11 +156,11 @@ function setup_geomparams(jsondict) end -function setup_battery_model_1d(jsondict; include_cc = true, use_groups = false) +function setup_battery_model_1d(jsondict; include_cc = true, use_groups = false, general_ad = false) geomparams = setup_geomparams(jsondict) - function setup_component(geomparam::Dict, sys; addDirichlet = false) + function setup_component(geomparam::Dict, sys; addDirichlet = false, general_ad = false) facearea = geomparam[:facearea] @@ -193,7 +193,11 @@ function setup_battery_model_1d(jsondict; include_cc = true, use_groups = false) end - flow = TwoPointPotentialFlowHardCoded(g) + if general_ad + flow = PotentialFlow(g) + else + flow = TwoPointPotentialFlowHardCoded(g) + end disc = (charge_flow = flow,) domain = DiscretizedDomain(domain, disc) @@ -202,7 +206,7 @@ function setup_battery_model_1d(jsondict; include_cc = true, use_groups = false) end - function setup_component(geomparams::Dict, sys::Electrolyte, bcfaces = nothing) + function setup_component(geomparams::Dict, sys::Electrolyte, bcfaces = nothing; general_ad = false) # specific implementation for electrolyte # requires geometric parameters for :NAM, :SEP, :PAM @@ -234,7 +238,11 @@ function setup_battery_model_1d(jsondict; include_cc = true, use_groups = false) domain[:halfTrans, HalfFaces()] = facearea*T_hf domain[:bcTrans, BoundaryFaces()] = facearea*T_b - flow = TwoPointPotentialFlowHardCoded(g) + if general_ad + flow = PotentialFlow(g) + else + flow = TwoPointPotentialFlowHardCoded(g) + end disc = (charge_flow = flow,) domain = DiscretizedDomain(domain, disc) @@ -276,13 +284,13 @@ function setup_battery_model_1d(jsondict; include_cc = true, use_groups = false) rp = inputparams_am["SolidDiffusion"]["rp"] N = Int64(inputparams_am["SolidDiffusion"]["N"]) D = inputparams_am["SolidDiffusion"]["D0"] - sys_am = ActiveMaterial{P2Ddiscretization}(am_params, rp, N, D) + sys_am = ActiveMaterialP2D(am_params, rp, N, D) else - sys_am = ActiveMaterial{NoParticleDiffusion}(am_params) + sys_am = ActiveMaterialNoParticleDiffusion(am_params) end geomparam = geomparams[name] - model_am = setup_component(geomparam, sys_am) + model_am = setup_component(geomparam, sys_am, general_ad = general_ad) return model_am @@ -292,7 +300,7 @@ function setup_battery_model_1d(jsondict; include_cc = true, use_groups = false) if include_cc sys_cc = CurrentCollector() - model_cc = setup_component(geomparams[:CC], sys_cc, addDirichlet = true) + model_cc = setup_component(geomparams[:CC], sys_cc, addDirichlet = true, general_ad = general_ad) end # Setup NAM @@ -319,7 +327,7 @@ function setup_battery_model_1d(jsondict; include_cc = true, use_groups = false) params[:conductivity] = func elyte = Electrolyte(params) - model_elyte = setup_component(geomparams, elyte) + model_elyte = setup_component(geomparams, elyte, general_ad = general_ad) # Setup PAM model_pam = setup_active_material(:PAM, geomparams) @@ -327,7 +335,7 @@ function setup_battery_model_1d(jsondict; include_cc = true, use_groups = false) # Setup negative current collector if any if include_cc sys_pp = CurrentCollector() - model_pp = setup_component(geomparams[:PP], sys_pp) + model_pp = setup_component(geomparams[:PP], sys_pp, general_ad = general_ad) end # Setup control model @@ -454,14 +462,14 @@ function setup_battery_model(exported; include_cc = true, use_p2d = true, use_gr else error("not recongized") end - + T_prm = typeof(am_params) if use_p2d rp = inputparams_am["SolidDiffusion"]["rp"] N = Int64(inputparams_am["SolidDiffusion"]["N"]) D = inputparams_am["SolidDiffusion"]["D0"] - sys_am = ActiveMaterial{P2Ddiscretization}(am_params, rp, N, D) + sys_am = ActiveMaterialP2D(am_params, rp, N, D) else - sys_am = ActiveMaterial{NoParticleDiffusion}(am_params) + sys_am = ActiveMaterialNoParticleDiffusion(am_params) end if include_cc