From 9290d2263bbc88a2ddff9caf67723251fcbdf634 Mon Sep 17 00:00:00 2001 From: Erasdna Date: Wed, 23 Aug 2023 13:11:19 +0200 Subject: [PATCH 1/2] mrst_test_utils cleanup from matlab-bridge --- src/BattMo.jl | 7 +- src/mrst_test_utils.jl | 2111 +++++++++++++++++++++------------------- src/types.jl | 34 + 3 files changed, 1165 insertions(+), 987 deletions(-) create mode 100644 src/types.jl diff --git a/src/BattMo.jl b/src/BattMo.jl index c2434e9..5fca732 100644 --- a/src/BattMo.jl +++ b/src/BattMo.jl @@ -71,14 +71,15 @@ include("models/current_and_voltage_boundary.jl") include("models/battery_cross_terms.jl") # Works now include("models/convergence.jl") include("physics.jl") +include("types.jl") include("mrst_test_utils.jl") include("linsolve.jl") # Precompilation of solver. Run a small battery simulation to precompile everything. @compile_workload begin - for use_general_ad in [false, true] - run_battery("p2d_40", info_level = -1, general_ad = use_general_ad); - end + for use_general_ad in [false, true] + run_battery("p2d_40", general_ad = use_general_ad,info_level = -1); + end end end # module diff --git a/src/mrst_test_utils.jl b/src/mrst_test_utils.jl index 074d663..0c5f03d 100644 --- a/src/mrst_test_utils.jl +++ b/src/mrst_test_utils.jl @@ -1,285 +1,694 @@ -struct Constants - F - R - hour - function Constants() - new(96485.3329, - 8.31446261815324, - 3600) - end +############################################################################################ +#Exported functions: +############################################################################################ +export run_battery +export inputRefToStates + +############################################################################################ +#Run battery +############################################################################################ + +function run_battery(init::InputFile; + use_p2d::Bool=true, + extra_timing::Bool=false, + max_step::Union{Integer, Nothing}=nothing, + linear_solver::Symbol=:direct, + general_ad::Bool=false, + use_groups::Bool=false, + kwarg...) + """ + Run battery wrapper method. Can use inputs from either Matlab or Json files and performs + simulation using a simple discharge CV policy + """ + + #Setup simulation + sim, forces, state0, parameters, init, model = setup_sim(init, use_p2d=use_p2d, use_groups=use_groups, general_ad=general_ad; kwarg...) + + #Set up config and timesteps + timesteps=setup_timesteps(init;max_step=max_step) + cfg=setup_config(sim,model,linear_solver,extra_timing; kwarg...) + + #Perform simulation + states, reports = simulate(state0,sim, timesteps, forces=forces, config=cfg; kwarg ...) + + extra = Dict(:model => model, + :state0 => state0, + :parameters => parameters, + :init => init, + :timesteps => timesteps, + :config => cfg, + :forces => forces, + :simulator => sim) + + return (states=states, reports=reports, extra=extra,exported=init) end -struct SourceAtCell - cell - src - function SourceAtCell(cell, src) - new(cell, src) - end +#Allows running run_battery with simple option for loading mat files +function run_battery(init::String; + use_p2d::Bool=true, + extra_timing::Bool=false, + max_step=nothing, + linear_solver::Symbol=:direct, + general_ad::Bool=false, + use_groups::Bool=false, + kwarg...) + """ + Simplifies reading pre-generated .mat files from the data repository inside BattMo.jl + """ + + #Path to mat files + path = string(dirname(pathof(BattMo)), "/../test/battery/data/") + suffix = ".mat" + + return run_battery(MatlabFile(path * init * suffix); + use_p2d=use_p2d, + extra_timing=extra_timing, + max_step=max_step, + linear_solver=linear_solver, + general_ad=general_ad, + use_groups=use_groups, + kwarg...) end -function getTrans(model1, model2, faces, cells, quantity) - """ setup transmissibility for coupling between models at boundaries""" +##################################################################################### +#Setup config +##################################################################################### + +function setup_config(sim::Jutul.JutulSimulator, + model::MultiModel, + linear_solver::Symbol, + extra_timing::Bool; + kwarg...) + """ + Sets up the config object used during simulation. In this current version this + setup is the same for json and mat files. The specific setup values should + probably be given as inputs in future versions of BattMo.jl + """ - T_all1 = model1["operators"]["T_all"][faces[:, 1]] - T_all2 = model2["operators"]["T_all"][faces[:, 2]] + cfg = simulator_config(sim; kwarg...) + cfg[:linear_solver] = battery_linsolve(model, linear_solver) + cfg[:debug_level] = 0 + #cfg[:max_timestep_cuts] = 0 + cfg[:max_residual] = 1e20 + cfg[:min_nonlinear_iterations] = 1 + cfg[:extra_timing] = extra_timing + # cfg[:max_nonlinear_iterations] = 5 + cfg[:safe_mode] = false + cfg[:error_on_incomplete] = false + #Original matlab steps will be too large! + cfg[:failure_cuts_timestep]=true - s1 = model1[quantity][cells[:, 1]] - s2 = model2[quantity][cells[:, 2]] - - T = 1.0./((1.0./(T_all1.*s1))+(1.0./(T_all2.*s2))) + # if false + # cfg[:info_level] = 5 + # cfg[:max_nonlinear_iterations] = 1 + # cfg[:max_timestep_cuts] = 0 + # end - return T + cfg[:tolerances][:PP][:default] = 1e-1 + cfg[:tolerances][:BPP][:default] = 1e-1 + + return cfg end -function getTrans_1d(model1, model2, bcfaces, bccells, parameters1, parameters2, quantity) - """ setup transmissibility for coupling between models at boundaries. Intermediate 1d version""" +##################################################################################### +#Setup timestepping +##################################################################################### - d1 = physical_representation(model1) - d2 = physical_representation(model2) +function setup_timesteps(init::JSONFile; + kwarg ...) + """ + Method setting up the timesteps from a json file object. + """ - bcTrans1 = d1[:bcTrans][bcfaces[:, 1]] - bcTrans2 = d2[:bcTrans][bcfaces[:, 2]] - cells1 = bccells[:, 1] - cells2 = bccells[:, 2] + total = init.object["TimeStepping"]["totalTime"] + n = init.object["TimeStepping"]["N"] - s1 = parameters1[quantity][cells1] - s2 = parameters2[quantity][cells2] - - T = 1.0./((1.0./(bcTrans1.*s1))+(1.0./(bcTrans2.*s2))) + dt = total / n + timesteps = rampupTimesteps(total, dt, 5) - return T - + return timesteps end -function getHalfTrans_1d(model, bcfaces, bccells, parameters, quantity) - """ recover half transmissibilities for boundary faces and weight them by the coefficient sent as quantity for the corresponding given cells. Intermediate 1d version. Note the indexing in BoundaryFaces is used""" +function setup_timesteps(init::MatlabFile; + max_step::Union{Integer,Nothing}=nothing, + kwarg...) + """ + Method setting up the timesteps from a mat file object. If use_state_ref is true + the simulation will use the same timesteps as the pre-run matlab simulation. + """ + + if init.use_state_ref + steps = size(init.object["states"], 1) + alltimesteps = Vector{Float64}(undef, steps) + time = 0 + end_step = 0 + + #Alternative to minE=3.2 + minE = init.object["model"]["Control"]["lowerCutoffVoltage"] + + for i = 1:steps + alltimesteps[i] = init.object["states"][i]["time"] - time + time = init.object["states"][i]["time"] + E = init.object["states"][i]["Control"]["E"] + if (E > minE + 0.001) + end_step = i + end + end + if !isnothing(max_step) + end_step = min(max_step, end_step) + end + timesteps = alltimesteps[1:end_step] + else + timesteps=init.object["schedule"]["step"]["val"][:] - d = physical_representation(model) - bcTrans = d[:bcTrans][bcfaces] - s = parameters[quantity][bccells] - - T = bcTrans.*s + end - return T + return timesteps end -function getHalfTrans(model, faces, cells, quantity) - """ recover half transmissibilities for boundary faces and weight them by the coefficient sent as quantity for the given cells. -Here, the faces should belong the corresponding cells at the same index""" +#################################################################################### +#Setup simulation +#################################################################################### - T_all = model["operators"]["T_all"] - s = model[quantity][cells] - T = T_all[faces].*s +#Replaces setup_sim_1d +function setup_sim(init::JSONFile; + use_groups::Bool=false, + general_ad::Bool=false, + kwarg ... ) - return T - -end + model, state0, parameters = setup_model(init, use_groups=use_groups, general_ad=general_ad; kwarg...) -function getHalfTrans(model, faces) - """ recover the half transmissibilities for boundary faces""" - - T_all = model["operators"]["T_all"] - T = T_all[faces] - - return T - -end + setup_coupling!(init,model,parameters) -function my_number_of_cells(model::MultiModel) - - cells = 0 - for smodel in model.models - cells += number_of_cells(smodel.domain) - end + minE = init.object["Control"]["lowerCutoffVoltage"] - return cells - -end + CRate = init.object["Control"]["CRate"] + cap = computeCellCapacity(model) + con = Constants() -function convert_to_int_vector(x::Float64) - vec = Int64.(Vector{Float64}([x])) - return vec -end + inputI = (cap / con.hour) * CRate -function convert_to_int_vector(x::Matrix{Float64}) - vec = Int64.(Vector{Float64}(x[:,1])) - return vec -end + @. state0[:BPP][:Phi] = minE * 1.5 -function setup_model(exported; use_p2d = true, use_groups = false, kwarg...) + tup = Float64(init.object["TimeStepping"]["rampupTime"]) + cFun(time) = currentFun(time, inputI, tup) - include_cc = true + currents = setup_forces(model[:BPP], policy=SimpleCVPolicy(cFun, minE)) + forces = setup_forces(model, BPP=currents) - model = setup_battery_model(exported, use_p2d = use_p2d, include_cc = include_cc) - parameters = setup_battery_parameters(exported, model) - initState = setup_battery_initial_state(exported, model) - - return model, initState, parameters - -end + sim = Simulator(model, state0=state0, parameters=parameters, copy_state=true) -function setup_model_1d(jsondict; use_groups = false, kwarg...) + return sim, forces, state0, parameters, init, model - include_cc = true - - 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) - - return model, initState, parameters - end -function setup_geomparams(jsondict) - - names = (:CC, :NAM, :SEP, :PAM, :PP) - geomparams = Dict(name => Dict() for name in names) - - geomparams[:CC][:N] = jsondict["NegativeElectrode"]["CurrentCollector"]["N"] - geomparams[:CC][:thickness] = jsondict["NegativeElectrode"]["CurrentCollector"]["thickness"] - geomparams[:NAM][:N] = jsondict["NegativeElectrode"]["ActiveMaterial"]["N"] - geomparams[:NAM][:thickness] = jsondict["NegativeElectrode"]["ActiveMaterial"]["thickness"] - geomparams[:SEP][:N] = jsondict["Electrolyte"]["Separator"]["N"] - geomparams[:SEP][:thickness] = jsondict["Electrolyte"]["Separator"]["thickness"] - geomparams[:PAM][:N] = jsondict["PositiveElectrode"]["ActiveMaterial"]["N"] - geomparams[:PAM][:thickness] = jsondict["PositiveElectrode"]["ActiveMaterial"]["thickness"] - geomparams[:PP][:N] = jsondict["PositiveElectrode"]["CurrentCollector"]["N"] - geomparams[:PP][:thickness] = jsondict["PositiveElectrode"]["CurrentCollector"]["thickness"] +function setup_sim(init::MatlabFile; + use_p2d::Bool=true, + use_groups::Bool=false, + general_ad::Bool=false) - for name in names - geomparams[name][:facearea] = jsondict["Geometry"]["faceArea"] - end - - return geomparams - -end + model, state0, parameters = setup_model(init, use_p2d=use_p2d, use_groups=use_groups, general_ad=general_ad) + setup_coupling!(init,model) + #quantities from matlab + minE=init.object["model"]["Control"]["lowerCutoffVoltage"] + inputI=init.object["model"]["Control"]["Imax"] -function setup_battery_model_1d(jsondict; include_cc = true, use_groups = false, general_ad = false) - - geomparams = setup_geomparams(jsondict) + #@. state0[:BPP][:Phi] = state0[:PAM][:Phi][end] #minE * 1.5 - function setup_component(geomparam::Dict, sys; addDirichlet = false, general_ad = false) + # if isempty(init.object["model"]["Control"]["tup"]) + # cFun(time) = currentFun(time, inputI) + # else + # cFun(time) = currentFun(time, inputI,init.object["model"]["Control"]["tup"]) + # end + cFun(time) = currentFun(time,inputI) + forces_pp = nothing - facearea = geomparam[:facearea] - - g = CartesianMesh(Tuple(geomparam[:N]), Tuple(geomparam[:thickness])) - domain = DataDomain(g) + currents = setup_forces(model[:BPP], policy=SimpleCVPolicy(cFun, minE)) - domain[:face_weighted_volumes] = facearea*domain[:volumes] - - k = ones(geomparam[:N]) + forces = Dict( + :CC => nothing, + :NAM => nothing, + :ELYTE => nothing, + :PAM => nothing, + :PP => forces_pp, + :BPP => currents + ) - T = compute_face_trans(domain, k) - T_hf = compute_half_face_trans(domain, k) - T_b = compute_boundary_trans(domain, k) + sim = Simulator(model, state0=state0, parameters=parameters, copy_state=true) - domain[:trans, Faces()] = facearea*T - domain[:halfTrans, HalfFaces()] = facearea*T_hf - domain[:bcTrans, BoundaryFaces()] = facearea*T_b + return sim, forces, state0, parameters, init, model - # We add Dirichlet on negative current collector. This is hacky - if addDirichlet +end - domain.entities[BoundaryDirichletFaces()] = 1 +####################################################################### +#Setup coupling +####################################################################### - bcDirFace = 1 # in BoundaryFaces indexing - bcDirCell = 1 - bcDirInd = 1 - domain[:bcDirHalfTrans, BoundaryDirichletFaces()] = facearea*domain[:bcTrans][bcDirFace] - domain[:bcDirCells, BoundaryDirichletFaces()] = facearea*bcDirCell # - domain[:bcDirInds, BoundaryDirichletFaces()] = facearea*bcDirInd # - - end - - if general_ad - flow = PotentialFlow(g) - else - flow = TwoPointPotentialFlowHardCoded(g) - end - disc = (charge_flow = flow,) - domain = DiscretizedDomain(domain, disc) - - model = SimulationModel(domain, sys, context = DefaultContext()) - return model - - end +#Replaces setup_coupling_1d! +function setup_coupling!(init::JSONFile, + model::MultiModel, + parameters::Dict{Symbol,<:Any} + ) - function setup_component(geomparams::Dict, sys::Electrolyte, bcfaces = nothing; general_ad = false) - # specific implementation for electrolyte - # requires geometric parameters for :NAM, :SEP, :PAM + jsondict=init.object - facearea = geomparams[:SEP][:facearea] - - names = (:NAM, :SEP, :PAM) + geomparams = setup_geomparams(init) + + # setup coupling CC <-> NAM :charge_conservation + + skip_cc = false # we consider only case with current collector (for simplicity, for the moment) + + if !skip_cc - deltas = Vector{Float64}() - for name in names - dx = geomparams[name][:thickness]/geomparams[name][:N] - dx = dx*ones(geomparams[name][:N]) - deltas = vcat(deltas, dx) - end + Ncc = geomparams[:CC][:N] - N = sum([geomparams[name][:N] for name in names]) - deltas = (deltas,) - g = CartesianMesh((N,), deltas) + srange = Ncc + trange = 1 - domain = DataDomain(g) - - domain[:face_weighted_volumes] = facearea*domain[:volumes] - - k = ones(N) - T = compute_face_trans(domain, k) - T_hf = compute_half_face_trans(domain, k) - T_b = compute_boundary_trans(domain, k) + msource = model[:CC] + mtarget = model[:NAM] + + psource = parameters[:CC] + ptarget = parameters[:NAM] - domain[:trans, Faces()] = facearea*T - domain[:halfTrans, HalfFaces()] = facearea*T_hf - domain[:bcTrans, BoundaryFaces()] = facearea*T_b + # Here, the indexing in BoundaryFaces is used + couplingfaces = Array{Int64}(undef, 1, 2) + couplingfaces[1, 1] = 2 + couplingfaces[1, 2] = 1 - if general_ad - flow = PotentialFlow(g) - else - flow = TwoPointPotentialFlowHardCoded(g) - end - disc = (charge_flow = flow,) - domain = DiscretizedDomain(domain, disc) + couplingcells = Array{Int64}(undef, 1, 2) + couplingcells[1, 1] = Ncc + couplingcells[1, 2] = 1 - model = SimulationModel(domain, sys, context = DefaultContext()) + trans = getTrans(msource, mtarget, + couplingfaces, + couplingcells, + psource, ptarget, + :Conductivity) - return model + ct = TPFAInterfaceFluxCT(trange, srange, trans) + ct_pair = setup_cross_term(ct, target = :NAM, source = :CC, equation = :charge_conservation) + add_cross_term!(model, ct_pair) end - - jsonNames = Dict( - :NAM => "NegativeElectrode", - :PAM => "PositiveElectrode", - ) - - function setup_active_material(name::Symbol, geomparams::Dict) - - jsonName = jsonNames[name] + # setup coupling NAM <-> ELYTE charge - inputparams_am = jsondict[jsonName]["ActiveMaterial"] + Nnam = geomparams[:NAM][:N] + + srange = collect(1 : Nnam) # negative electrode + trange = collect(1 : Nnam) # electrolyte (negative side) - am_params = JutulStorage() - am_params[:volumeFraction] = inputparams_am["Interface"]["volumeFraction"] - am_params[:n_charge_carriers] = inputparams_am["Interface"]["n"] - am_params[:maximum_concentration] = inputparams_am["Interface"]["cmax"] - am_params[:volumetric_surface_area] = inputparams_am["Interface"]["volumetricSurfaceArea"] - am_params[:theta0] = inputparams_am["Interface"]["theta0"] - am_params[:theta100] = inputparams_am["Interface"]["theta100"] + if discretisation_type(model[:NAM]) == :P2Ddiscretization - k0 = inputparams_am["Interface"]["k0"] - Eak = inputparams_am["Interface"]["Eak"] - am_params[:reaction_rate_constant_func] = (c, T) -> compute_reaction_rate_constant(c, T, k0, Eak) - - funcname = inputparams_am["Interface"]["OCP"]["functionname"] - am_params[:ocp_func] = getfield(BattMo, Symbol(funcname)) + ct = ButlerVolmerActmatToElyteCT(trange, srange) + ct_pair = setup_cross_term(ct, target = :ELYTE, source = :NAM, equation = :charge_conservation) + add_cross_term!(model, ct_pair) + ct_pair = setup_cross_term(ct, target = :ELYTE, source = :NAM, equation = :mass_conservation) + add_cross_term!(model, ct_pair) + + ct = ButlerVolmerElyteToActmatCT(srange, trange) + ct_pair = setup_cross_term(ct, target = :NAM, source = :ELYTE, equation = :charge_conservation) + add_cross_term!(model, ct_pair) + ct_pair = setup_cross_term(ct, target = :NAM, source = :ELYTE, equation = :solid_diffusion_bc) + add_cross_term!(model, ct_pair) + + else + + @assert discretisation_type(model[:NAM]) == :NoParticleDiffusion + + ct = ButlerVolmerInterfaceFluxCT(trange, srange) + ct_pair = setup_cross_term(ct, target = :ELYTE, source = :NAM, equation = :charge_conservation) + add_cross_term!(model, ct_pair) + ct_pair = setup_cross_term(ct, target = :ELYTE, source = :NAM, equation = :mass_conservation) + add_cross_term!(model, ct_pair) + + end + + # setup coupling ELYTE <-> PAM charge + + Nnam = geomparams[:NAM][:N] + Nsep = geomparams[:SEP][:N] + Npam = geomparams[:PAM][:N] + + srange = collect(1 : Npam) # positive electrode + trange = collect(Nnam + Nsep .+ (1 : Npam)) # electrolyte (positive side) + + if discretisation_type(model[:PAM]) == :P2Ddiscretization + + ct = ButlerVolmerActmatToElyteCT(trange, srange) + ct_pair = setup_cross_term(ct, target = :ELYTE, source = :PAM, equation = :charge_conservation) + add_cross_term!(model, ct_pair) + ct_pair = setup_cross_term(ct, target = :ELYTE, source = :PAM, equation = :mass_conservation) + add_cross_term!(model, ct_pair) + + ct = ButlerVolmerElyteToActmatCT(srange, trange) + ct_pair = setup_cross_term(ct, target = :PAM, source = :ELYTE, equation = :charge_conservation) + add_cross_term!(model, ct_pair) + ct_pair = setup_cross_term(ct, target = :PAM, source = :ELYTE, equation = :solid_diffusion_bc) + add_cross_term!(model, ct_pair) + + else + + @assert discretisation_type(model[:PAM]) == :NoParticleDiffusion + + ct = ButlerVolmerInterfaceFluxCT(trange, srange) + ct_pair = setup_cross_term(ct, target = :ELYTE, source = :PAM, equation = :charge_conservation) + add_cross_term!(model, ct_pair) + ct_pair = setup_cross_term(ct, target = :ELYTE, source = :PAM, equation = :mass_conservation) + add_cross_term!(model, ct_pair) + + end + + + if !skip_cc + # setup coupling PP <-> PAM charge + + Npam = geomparams[:PAM][:N] + + srange = 1 + trange = Npam + + msource = model[:PP] + mtarget = model[:PAM] + + psource = parameters[:PP] + ptarget = parameters[:PAM] + + # Here, the indexing in BoundaryFaces is used + couplingfaces = Array{Int64}(undef, 1, 2) + couplingfaces[1, 1] = 1 + couplingfaces[1, 2] = 2 + + couplingcells = Array{Int64}(undef, 1, 2) + couplingcells[1, 1] = 1 + couplingcells[1, 2] = Npam + + + trans = getTrans(msource, mtarget, + couplingfaces, + couplingcells, + psource, ptarget, + :Conductivity) + + ct = TPFAInterfaceFluxCT(trange, srange, trans) + ct_pair = setup_cross_term(ct, + target = :PAM, + source = :PP, + equation = :charge_conservation) + + add_cross_term!(model, ct_pair) + + # setup coupling with control + + Npp = geomparams[:PP][:N] + + trange = Npp + srange = Int64.(ones(size(trange))) + + msource = model[:PP] + mparameters = parameters[:PP] + # Here the indexing in BoundaryFaces in used + couplingfaces = 2 + couplingcells = Npp + trans = getHalfTrans(msource, couplingfaces, couplingcells, mparameters, :Conductivity) + + ct = TPFAInterfaceFluxCT(trange, srange, trans, symmetric = false) + ct_pair = setup_cross_term(ct, target = :PP, source = :BPP, equation = :charge_conservation) + add_cross_term!(model, ct_pair) + + # Accmulation of charge + ct = AccumulatorInterfaceFluxCT(1, trange, trans) + ct_pair = setup_cross_term(ct, target = :BPP, source = :PP, equation = :charge_conservation) + add_cross_term!(model, ct_pair) + + end + +end + +function setup_coupling!(init::MatlabFile, + model #MultiModel? + ) + # setup coupling CC <-> NAM :charge_conservation + + exported_all=init.object + skip_pp = size(exported_all["model"]["include_current_collectors"]) == (0,0) #! unused + skip_cc = false + + if !skip_cc + + srange = Int64.( + exported_all["model"]["NegativeElectrode"]["couplingTerm"]["couplingcells"][:, 1] + ) + trange = Int64.( + exported_all["model"]["NegativeElectrode"]["couplingTerm"]["couplingcells"][:, 2] + ) + + msource = exported_all["model"]["NegativeElectrode"]["CurrentCollector"] + mtarget = exported_all["model"]["NegativeElectrode"]["ActiveMaterial"] + couplingfaces = Int64.(exported_all["model"]["NegativeElectrode"]["couplingTerm"]["couplingfaces"]) + couplingcells = Int64.(exported_all["model"]["NegativeElectrode"]["couplingTerm"]["couplingcells"]) + trans = getTrans(msource, mtarget, couplingfaces, couplingcells, "EffectiveElectricalConductivity") + + ct = TPFAInterfaceFluxCT(trange, srange, trans) + ct_pair = setup_cross_term(ct, target = :NAM, source = :CC, equation = :charge_conservation) + add_cross_term!(model, ct_pair) + + end + + # setup coupling NAM <-> ELYTE charge + + srange = Int64.(exported_all["model"]["couplingTerms"][1]["couplingcells"][:, 1]) # negative electrode + trange = Int64.(exported_all["model"]["couplingTerms"][1]["couplingcells"][:, 2]) # electrolyte (negative side) + + if discretisation_type(model[:NAM]) == :P2Ddiscretization + + ct = ButlerVolmerActmatToElyteCT(trange, srange) + ct_pair = setup_cross_term(ct, target = :ELYTE, source = :NAM, equation = :charge_conservation) + add_cross_term!(model, ct_pair) + ct_pair = setup_cross_term(ct, target = :ELYTE, source = :NAM, equation = :mass_conservation) + add_cross_term!(model, ct_pair) + + ct = ButlerVolmerElyteToActmatCT(srange, trange) + ct_pair = setup_cross_term(ct, target = :NAM, source = :ELYTE, equation = :charge_conservation) + add_cross_term!(model, ct_pair) + ct_pair = setup_cross_term(ct, target = :NAM, source = :ELYTE, equation = :solid_diffusion_bc) + add_cross_term!(model, ct_pair) + + else + + @assert discretisation_type(model[:NAM]) == :NoParticleDiffusion + + ct = ButlerVolmerInterfaceFluxCT(trange, srange) + ct_pair = setup_cross_term(ct, target = :ELYTE, source = :NAM, equation = :charge_conservation) + add_cross_term!(model, ct_pair) + ct_pair = setup_cross_term(ct, target = :ELYTE, source = :NAM, equation = :mass_conservation) + add_cross_term!(model, ct_pair) + + end + + # setup coupling ELYTE <-> PAM charge + + srange = Int64.(exported_all["model"]["couplingTerms"][2]["couplingcells"][:,1]) # postive electrode + trange = Int64.(exported_all["model"]["couplingTerms"][2]["couplingcells"][:,2]) # electrolyte (positive side) + + if discretisation_type(model[:PAM]) == :P2Ddiscretization + + ct = ButlerVolmerActmatToElyteCT(trange, srange) + ct_pair = setup_cross_term(ct, target = :ELYTE, source = :PAM, equation = :charge_conservation) + add_cross_term!(model, ct_pair) + ct_pair = setup_cross_term(ct, target = :ELYTE, source = :PAM, equation = :mass_conservation) + add_cross_term!(model, ct_pair) + + ct = ButlerVolmerElyteToActmatCT(srange, trange) + ct_pair = setup_cross_term(ct, target = :PAM, source = :ELYTE, equation = :charge_conservation) + add_cross_term!(model, ct_pair) + ct_pair = setup_cross_term(ct, target = :PAM, source = :ELYTE, equation = :solid_diffusion_bc) + add_cross_term!(model, ct_pair) + + else + + @assert discretisation_type(model[:PAM]) == :NoParticleDiffusion + + ct = ButlerVolmerInterfaceFluxCT(trange, srange) + ct_pair = setup_cross_term(ct, target = :ELYTE, source = :PAM, equation = :charge_conservation) + add_cross_term!(model, ct_pair) + ct_pair = setup_cross_term(ct, target = :ELYTE, source = :PAM, equation = :mass_conservation) + add_cross_term!(model, ct_pair) + + end + + + if !skip_cc + # setup coupling PP <-> PAM charge + target = Dict( + :model => :PAM, + :equation => :charge_conservation + ) + source = Dict( + :model => :PP, + :equation => :charge_conservation + ) + srange = Int64.( + exported_all["model"]["PositiveElectrode"]["couplingTerm"]["couplingcells"][:,1] + ) + trange = Int64.( + exported_all["model"]["PositiveElectrode"]["couplingTerm"]["couplingcells"][:,2] + ) + msource = exported_all["model"]["PositiveElectrode"]["CurrentCollector"] + ct = exported_all["model"]["PositiveElectrode"]["couplingTerm"] + couplingfaces = Int64.(ct["couplingfaces"]) + couplingcells = Int64.(ct["couplingcells"]) + trans = getTrans(msource, mtarget, couplingfaces, couplingcells, "EffectiveElectricalConductivity") + ct = TPFAInterfaceFluxCT(trange, srange, trans) + ct_pair = setup_cross_term(ct, target = :PAM, source = :PP, equation = :charge_conservation) + add_cross_term!(model, ct_pair) + end + + if skip_cc + #setup coupling PP <-> BPP charge + target = Dict( + :model => :PAM, + :equation => :charge_conservation + ) + + trange = convert_to_int_vector( + exported_all["model"]["PositiveElectrode"]["ActiveMaterial"]["externalCouplingTerm"]["couplingcells"] + ) + srange = Int64.(ones(size(trange))) + msource = exported_all["model"]["PositiveElectrode"]["ActiveMaterial"] + couplingfaces = Int64.(msource["externalCouplingTerm"]["couplingfaces"]) + couplingcells = Int64.(msource["externalCouplingTerm"]["couplingcells"]) + else + #setup coupling PP <-> BPP charge + target = Dict( + :model => :PP, + :equation => :charge_conservation + ) + + trange = convert_to_int_vector( + exported_all["model"]["PositiveElectrode"]["CurrentCollector"]["externalCouplingTerm"]["couplingcells"] + ) + srange = Int64.(ones(size(trange))) + msource = exported_all["model"]["PositiveElectrode"]["CurrentCollector"] + couplingfaces = Int64.(msource["externalCouplingTerm"]["couplingfaces"]) + couplingcells = Int64.(msource["externalCouplingTerm"]["couplingcells"]) + end + + source = Dict( + :model => :BPP, + :equation => :charge_conservation + ) + + #effcond = exported_all["model"]["PositiveElectrode"]["CurrentCollector"]["EffectiveElectricalConductivity"] + trans = getHalfTrans(msource, couplingfaces, couplingcells, "EffectiveElectricalConductivity") + + if skip_cc + + ct = TPFAInterfaceFluxCT(trange, srange, trans, symmetric = false) + ct_pair = setup_cross_term(ct, target = :PAM, source = :BPP, equation = :charge_conservation) + add_cross_term!(model, ct_pair) + + # Accmulation of charge + ct = AccumulatorInterfaceFluxCT(1, trange, trans) + ct_pair = setup_cross_term(ct, target = :BPP, source = :PAM, equation = :charge_conservation) + add_cross_term!(model, ct_pair) + + else + + ct = TPFAInterfaceFluxCT(trange, srange, trans, symmetric = false) + ct_pair = setup_cross_term(ct, target = :PP, source = :BPP, equation = :charge_conservation) + add_cross_term!(model, ct_pair) + + # Accmulation of charge + ct = AccumulatorInterfaceFluxCT(1, trange, trans) + ct_pair = setup_cross_term(ct, target = :BPP, source = :PP, equation = :charge_conservation) + add_cross_term!(model, ct_pair) + + end + +end + +######################################################################## +#Setup model +######################################################################## + +function setup_model(init::InputFile; + use_p2d::Bool=true, + use_groups::Bool=false, + kwarg...) + + include_cc = true #! + + model = setup_battery_model(init, include_cc=include_cc,use_groups=use_groups,use_p2d=use_p2d; kwarg ... ) + parameters = setup_battery_parameters(init, model) + initState = setup_battery_initial_state(init, model) + + return model, initState, parameters + +end + +################################################################################### +#Setup battery model +################################################################################## + +function setup_battery_model(init::MatlabFile; + include_cc::Bool = true, + use_groups::Bool = false, + use_p2d::Bool = true, + general_ad::Bool = false, + kwarg ... + ) + + + function setup_component(obj::Dict, + sys, + bcfaces = nothing, + general_ad::Bool = false) - use_p2d = true + domain = exported_model_to_domain(obj, bcfaces = bcfaces, general_ad=general_ad) + G = MRSTWrapMesh(obj["G"]) + data_domain = DataDomain(G) + for (k, v) in domain.entities + data_domain.entities[k] = v + end + model = SimulationModel(domain, sys, context = DefaultContext(), data_domain = data_domain) + return model + + end + + jsonNames = Dict( + :NAM => "NegativeElectrode", + :PAM => "PositiveElectrode", + ) + + inputparams = init.object["model"] + + function setup_active_material(name::Symbol, general_ad::Bool) + jsonName = jsonNames[name] + + inputparams_am = inputparams[jsonName]["ActiveMaterial"] + + am_params = JutulStorage() + am_params[:n_charge_carriers] = inputparams_am["Interface"]["n"] + am_params[:maximum_concentration] = inputparams_am["Interface"]["cmax"] + am_params[:volumetric_surface_area] = inputparams_am["Interface"]["volumetricSurfaceArea"] + am_params[:volume_fraction] = inputparams_am["Interface"]["volumeFraction"] + + k0 = inputparams_am["Interface"]["k0"] + Eak = inputparams_am["Interface"]["Eak"] + am_params[:reaction_rate_constant_func] = (c, T) -> compute_reaction_rate_constant(c, T, k0, Eak) + + if name == :NAM + am_params[:ocp_func] = compute_ocp_graphite + elseif name == :PAM + am_params[:ocp_func] = compute_ocp_nmc111 + else + error("not recongized") + end + T_prm = typeof(am_params) if use_p2d rp = inputparams_am["SolidDiffusion"]["rp"] N = Int64(inputparams_am["SolidDiffusion"]["N"]) @@ -289,53 +698,75 @@ function setup_battery_model_1d(jsondict; include_cc = true, use_groups = false, sys_am = ActiveMaterialNoParticleDiffusion(am_params) end - geomparam = geomparams[name] - model_am = setup_component(geomparam, sys_am, general_ad = general_ad) + if include_cc + model_am = setup_component(inputparams_am, sys_am,nothing,general_ad) + else + bcfaces_am = convert_to_int_vector(inputparams_am["externalCouplingTerm"]["couplingfaces"]) + model_am = setup_component(inputparams_am, sys_am, bcfaces_am,general_ad) + # We add also boundary parameters (if any) + S = model_am.parameters + nbc = count_active_entities(model_am.domain, BoundaryDirichletFaces()) + if nbc > 0 + bcvalue_zeros = zeros(nbc) + # add parameters to the model + S[:BoundaryPhi] = BoundaryPotential(:Phi) + S[:BoundaryC] = BoundaryPotential(:C) + S[:BCCharge] = BoundaryCurrent(srccells, :Charge) + S[:BCMass] = BoundaryCurrent(srccells, :Mass) + end + end return model_am end - # Setup negative current collector + # Setup positive current collector if any if include_cc - sys_cc = CurrentCollector() - model_cc = setup_component(geomparams[:CC], sys_cc, addDirichlet = true, general_ad = general_ad) + + inputparams_cc = inputparams["NegativeElectrode"]["CurrentCollector"] + sys_cc = CurrentCollector() + bcfaces = convert_to_int_vector(inputparams_cc["externalCouplingTerm"]["couplingfaces"]) + + model_cc = setup_component(inputparams_cc, sys_cc, bcfaces, general_ad) + end # Setup NAM - model_nam = setup_active_material(:NAM, geomparams) - ## Setup ELYTE + model_nam = setup_active_material(:NAM,general_ad) + ## Setup ELYTE params = JutulStorage(); - inputparams_elyte = jsondict["Electrolyte"] - - params[:transference] = inputparams_elyte["sp"]["t"] - params[:charge] = inputparams_elyte["sp"]["z"] - params[:separator_porosity] = inputparams_elyte["Separator"]["porosity"] - params[:bruggeman] = inputparams_elyte["BruggemanCoefficient"] + inputparams_elyte = inputparams["Electrolyte"] + params[:transference] = inputparams_elyte["sp"]["t"] + params[:charge] = inputparams_elyte["sp"]["z"] + params[:bruggeman] = inputparams_elyte["BruggemanCoefficient"] - # setup diffusion coefficient function - funcname = inputparams_elyte["DiffusionCoefficient"]["functionname"] + # setup diffusion coefficient function, hard coded for the moment because function name is not passed throught model + # TODO : add general code + funcname = "computeDiffusionCoefficient_default" func = getfield(BattMo, Symbol(funcname)) params[:diffusivity] = func # setup diffusion coefficient function - funcname = inputparams_elyte["Conductivity"]["functionname"] + # TODO : add general code + funcname = "computeElectrolyteConductivity_default" func = getfield(BattMo, Symbol(funcname)) params[:conductivity] = func elyte = Electrolyte(params) - model_elyte = setup_component(geomparams, elyte, general_ad = general_ad) + model_elyte = setup_component(inputparams["Electrolyte"], + elyte,nothing,general_ad) # Setup PAM - model_pam = setup_active_material(:PAM, geomparams) + + model_pam = setup_active_material(:PAM,general_ad) # Setup negative current collector if any if include_cc - sys_pp = CurrentCollector() - model_pp = setup_component(geomparams[:PP], sys_pp, general_ad = general_ad) + model_pp = setup_component(inputparams["PositiveElectrode"]["CurrentCollector"], + CurrentCollector(),nothing,general_ad) end # Setup control model @@ -375,59 +806,116 @@ function setup_battery_model_1d(jsondict; include_cc = true, use_groups = false, model = MultiModel(models, groups = groups, reduction = reduction) end - - setup_volume_fractions_1d!(model, geomparams) return model end -function setup_volume_fractions_1d!(model, geomparams) - - names = (:NAM, :SEP, :PAM) - Nelyte = sum([geomparams[name][:N] for name in names]) - vfelyte = zeros(Nelyte) - - names = (:NAM, :PAM) +function setup_battery_model(init::JSONFile; + include_cc::Bool = true, + use_groups::Bool = false, + general_ad::Bool = false, + kwarg ... ) - for name in names - ammodel = model[name] - vf = ammodel.system[:volumeFraction] - Nam = geomparams[name][:N] - ammodel.domain.representation[:volumeFraction] = vf*ones(Nam) - if name == :NAM - nstart = 1 - nend = Nam - elseif name == :PAM - nstart = geomparams[:NAM][:N] + geomparams[:SEP][:N] + 1 - nend = Nelyte + geomparams = setup_geomparams(init) + + jsondict=init.object + + function setup_component(geomparam::Dict, + sys; + addDirichlet::Bool = false, + general_ad::Bool = false) + + facearea = geomparam[:facearea] + + g = CartesianMesh(Tuple(geomparam[:N]), Tuple(geomparam[:thickness])) + domain = DataDomain(g) + + domain[:face_weighted_volumes] = facearea*domain[:volumes] + + k = ones(geomparam[:N]) + + T = compute_face_trans(domain, k) + T_hf = compute_half_face_trans(domain, k) + T_b = compute_boundary_trans(domain, k) + + domain[:trans, Faces()] = facearea*T + domain[:halfTrans, HalfFaces()] = facearea*T_hf + domain[:bcTrans, BoundaryFaces()] = facearea*T_b + + # We add Dirichlet on negative current collector. This is hacky + if addDirichlet + + domain.entities[BoundaryDirichletFaces()] = 1 + + bcDirFace = 1 # in BoundaryFaces indexing + bcDirCell = 1 + bcDirInd = 1 + domain[:bcDirHalfTrans, BoundaryDirichletFaces()] = facearea*domain[:bcTrans][bcDirFace] + domain[:bcDirCells, BoundaryDirichletFaces()] = facearea*bcDirCell # + domain[:bcDirInds, BoundaryDirichletFaces()] = facearea*bcDirInd # + + end + + if general_ad + flow = PotentialFlow(g) else - error("name not recognized") + flow = TwoPointPotentialFlowHardCoded(g) end - vfelyte[nstart : nend] .= 1 - vf + disc = (charge_flow = flow,) + domain = DiscretizedDomain(domain, disc) + + model = SimulationModel(domain, sys, context = DefaultContext()) + return model + end - nstart = geomparams[:NAM][:N] + 1 - nend = nstart + geomparams[:SEP][:N] - separator_porosity = model[:ELYTE].system[:separator_porosity] - vfelyte[nstart : nend] .= separator_porosity*ones(nend - nstart + 1) - - model[:ELYTE].domain.representation[:volumeFraction] = vfelyte + function setup_component(geomparams::Dict, + sys::Electrolyte, + bcfaces = nothing; + general_ad::Bool = false) -end + # specific implementation for electrolyte + # requires geometric parameters for :NAM, :SEP, :PAM + facearea = geomparams[:SEP][:facearea] + + names = (:NAM, :SEP, :PAM) -function setup_battery_model(exported; include_cc = true, use_p2d = true, use_groups = false) + deltas = Vector{Float64}() + for name in names + dx = geomparams[name][:thickness]/geomparams[name][:N] + dx = dx*ones(geomparams[name][:N]) + deltas = vcat(deltas, dx) + end - function setup_component(exported, sys, bcfaces = nothing) + N = sum([geomparams[name][:N] for name in names]) + deltas = (deltas,) + g = CartesianMesh((N,), deltas) - domain = exported_model_to_domain(exported, bcfaces = bcfaces) - G = MRSTWrapMesh(exported["G"]) - data_domain = DataDomain(G) - for (k, v) in domain.entities - data_domain.entities[k] = v + domain = DataDomain(g) + + domain[:face_weighted_volumes] = facearea*domain[:volumes] + + k = ones(N) + T = compute_face_trans(domain, k) + T_hf = compute_half_face_trans(domain, k) + T_b = compute_boundary_trans(domain, k) + + domain[:trans, Faces()] = facearea*T + domain[:halfTrans, HalfFaces()] = facearea*T_hf + domain[:bcTrans, BoundaryFaces()] = facearea*T_b + + if general_ad + flow = PotentialFlow(g) + else + flow = TwoPointPotentialFlowHardCoded(g) end - model = SimulationModel(domain, sys, context = DefaultContext(), data_domain = data_domain) + disc = (charge_flow = flow,) + domain = DiscretizedDomain(domain, disc) + + model = SimulationModel(domain, sys, context = DefaultContext()) + return model end @@ -437,32 +925,30 @@ function setup_battery_model(exported; include_cc = true, use_p2d = true, use_gr :PAM => "PositiveElectrode", ) - inputparams = exported["model"] - function setup_active_material(name) + function setup_active_material(name::Symbol, + geomparams::Dict{Symbol,<:Any}) jsonName = jsonNames[name] - inputparams_am = inputparams[jsonName]["ActiveMaterial"] + inputparams_am = jsondict[jsonName]["ActiveMaterial"] am_params = JutulStorage() + am_params[:volumeFraction] = inputparams_am["Interface"]["volumeFraction"] am_params[:n_charge_carriers] = inputparams_am["Interface"]["n"] am_params[:maximum_concentration] = inputparams_am["Interface"]["cmax"] am_params[:volumetric_surface_area] = inputparams_am["Interface"]["volumetricSurfaceArea"] - am_params[:volume_fraction] = inputparams_am["Interface"]["volumeFraction"] - + am_params[:theta0] = inputparams_am["Interface"]["theta0"] + am_params[:theta100] = inputparams_am["Interface"]["theta100"] + k0 = inputparams_am["Interface"]["k0"] Eak = inputparams_am["Interface"]["Eak"] am_params[:reaction_rate_constant_func] = (c, T) -> compute_reaction_rate_constant(c, T, k0, Eak) - if name == :NAM - am_params[:ocp_func] = compute_ocp_graphite - elseif name == :PAM - am_params[:ocp_func] = compute_ocp_nmc111 - else - error("not recongized") - end - T_prm = typeof(am_params) + funcname = inputparams_am["Interface"]["OCP"]["functionname"] + am_params[:ocp_func] = getfield(BattMo, Symbol(funcname)) + + use_p2d = true if use_p2d rp = inputparams_am["SolidDiffusion"]["rp"] N = Int64(inputparams_am["SolidDiffusion"]["N"]) @@ -472,75 +958,53 @@ function setup_battery_model(exported; include_cc = true, use_p2d = true, use_gr sys_am = ActiveMaterialNoParticleDiffusion(am_params) end - if include_cc - model_am = setup_component(inputparams_am, sys_am) - else - bcfaces_am = convert_to_int_vector(["externalCouplingTerm"]["couplingfaces"]) - model_am = setup_component(inputparams_am, sys_am, bcfaces_am) - # We add also boundary parameters (if any) - S = model_am.parameters - nbc = count_active_entities(model_am.domain, BoundaryDirichletFaces()) - if nbc > 0 - bcvalue_zeros = zeros(nbc) - # add parameters to the model - S[:BoundaryPhi] = BoundaryPotential(:Phi) - S[:BoundaryC] = BoundaryPotential(:C) - S[:BCCharge] = BoundaryCurrent(srccells, :Charge) - S[:BCMass] = BoundaryCurrent(srccells, :Mass) - end - end + geomparam = geomparams[name] + model_am = setup_component(geomparam, sys_am, general_ad = general_ad) return model_am end - # Setup positive current collector if any + # Setup negative current collector if include_cc - - inputparams_cc = inputparams["NegativeElectrode"]["CurrentCollector"] - sys_cc = CurrentCollector() - bcfaces = convert_to_int_vector(inputparams_cc["externalCouplingTerm"]["couplingfaces"]) - - model_cc = setup_component(inputparams_cc, sys_cc, bcfaces) - + sys_cc = CurrentCollector() + model_cc = setup_component(geomparams[:CC], sys_cc, addDirichlet = true, general_ad = general_ad) end # Setup NAM - - model_nam = setup_active_material(:NAM) + model_nam = setup_active_material(:NAM, geomparams) ## Setup ELYTE + params = JutulStorage(); - inputparams_elyte = inputparams["Electrolyte"] - params[:transference] = inputparams_elyte["sp"]["t"] - params[:charge] = inputparams_elyte["sp"]["z"] - params[:bruggeman] = inputparams_elyte["BruggemanCoefficient"] + inputparams_elyte = jsondict["Electrolyte"] - # setup diffusion coefficient function, hard coded for the moment because function name is not passed throught model - # TODO : add general code - funcname = "computeDiffusionCoefficient_default" + params[:transference] = inputparams_elyte["sp"]["t"] + params[:charge] = inputparams_elyte["sp"]["z"] + params[:separator_porosity] = inputparams_elyte["Separator"]["porosity"] + params[:bruggeman] = inputparams_elyte["BruggemanCoefficient"] + + # setup diffusion coefficient function + funcname = inputparams_elyte["DiffusionCoefficient"]["functionname"] func = getfield(BattMo, Symbol(funcname)) params[:diffusivity] = func # setup diffusion coefficient function - # TODO : add general code - funcname = "computeElectrolyteConductivity_default" + funcname = inputparams_elyte["Conductivity"]["functionname"] func = getfield(BattMo, Symbol(funcname)) params[:conductivity] = func elyte = Electrolyte(params) - model_elyte = setup_component(inputparams["Electrolyte"], - elyte) + model_elyte = setup_component(geomparams, elyte, general_ad = general_ad) # Setup PAM - - model_pam = setup_active_material(:PAM) + model_pam = setup_active_material(:PAM, geomparams) # Setup negative current collector if any if include_cc - model_pp = setup_component(inputparams["PositiveElectrode"]["CurrentCollector"], - CurrentCollector()) + sys_pp = CurrentCollector() + model_pp = setup_component(geomparams[:PP], sys_pp, general_ad = general_ad) end # Setup control model @@ -580,14 +1044,24 @@ function setup_battery_model(exported; include_cc = true, use_p2d = true, use_gr model = MultiModel(models, groups = groups, reduction = reduction) end + + setup_volume_fractions!(model, geomparams) return model end -function setup_battery_parameters(exported, model) +################################################################################### +#Setup battery parameters +################################################################################### - parameters = Dict{Symbol, Any}() +function setup_battery_parameters(init::MatlabFile, + model::MultiModel + ) + + parameters = Dict{Symbol, Any}() #NB + + exported=init.object T0 = exported["model"]["initT"] @@ -667,10 +1141,14 @@ function setup_battery_parameters(exported, model) end -function setup_battery_parameters_1d(jsonstruct, model) +function setup_battery_parameters(init::JSONFile, + model::MultiModel + ) parameters = Dict{Symbol, Any}() + jsonstruct=init.object + T0 = jsonstruct["initT"] # Negative current collector (if any) @@ -754,90 +1232,15 @@ function setup_battery_parameters_1d(jsonstruct, model) end -function setup_battery_initial_state_1d(jsonstruct, model) - - if haskey(model.models, :CC) - use_cc = true - else - use_cc = false - end - - T = jsonstruct["initT"] - SOC = jsonstruct["SOC"] - - function setup_init_am(name, model) - - theta0 = model[name].system[:theta0] - theta100 = model[name].system[:theta100] - cmax = model[name].system[:maximum_concentration] - N = model[name].system.discretization[:N] - - theta = SOC*(theta100 - theta0) + theta0; - c = theta*cmax - nc = count_entities(model[name].data_domain, Cells()) - init = Dict() - init[:Cs] = c*ones(nc) - init[:Cp] = c*ones(nc, N) - - OCP = model[name].system[:ocp_func](c, T, cmax) - return (init, nc, OCP) - - end - - function setup_cc(name, phi, model) - nc = count_entities(model[name].data_domain, Cells()) - init = Dict(); - init[:Phi] = phi*ones(nc) - return init - end - - initState = Dict() - - # Setup initial state in negative active material - - init, nc, negOCP = setup_init_am(:NAM, model) - init[:Phi] = zeros(nc) - initState[:NAM] = init - - # Setup initial state in electrolyte - - nc = count_entities(model[:ELYTE].data_domain, Cells()) - - init = Dict() - init[:C] = 1000*ones(nc) - init[:Phi] = - negOCP*ones(nc) - - initState[:ELYTE] = init - - # Setup initial state in positive active material - - init, nc, posOCP = setup_init_am(:PAM, model) - init[:Phi] = (posOCP - negOCP)*ones(nc) - - initState[:PAM] = init - - # Setup negative current collector - - initState[:CC] = setup_cc(:CC, 0, model) - - # Setup positive current collector - - initState[:PP] = setup_cc(:PP, posOCP - negOCP, model) - - init = Dict() - init[:Phi] = [1.0] - init[:Current] = [1.0] - - initState[:BPP] = init - - initState = setup_state(model, initState) - - return initState - -end +################################################################################### +#Setup initial state +################################################################################### +function setup_battery_initial_state(init::MatlabFile, + model::MultiModel + ) -function setup_battery_initial_state(exported, model) + exported=init.object state0 = exported["state0"] @@ -905,405 +1308,147 @@ function setup_battery_initial_state(exported, model) init[:C] = c end - initState[name] = init - - end - - function initialize_electrolyte!(initState) - - init = Dict() - - init[:Phi] = state0["Electrolyte"]["phi"][1] - init[:C] = state0["Electrolyte"]["c"][1] - - initState[:ELYTE] = init - - end - - function initialize_bpp!(initState) - - init = Dict(:Phi => 1.0, :Current => 1.0) - - initState[:BPP] = init - - end - - initState = Dict() - - initialize_current_collector!(initState, :CC) - initialize_active_material!(initState, :NAM) - initialize_electrolyte!(initState) - initialize_active_material!(initState, :PAM) - initialize_current_collector!(initState, :PP) - initialize_bpp!(initState) - - - initState = setup_state(model, initState) - - return initState - -end - -function setup_coupling_1d!(model, parameters, jsondict) - - geomparams = setup_geomparams(jsondict) - - # setup coupling CC <-> NAM :charge_conservation - - skip_cc = false # we consider only case with current collector (for simplicity, for the moment) - - if !skip_cc - - Ncc = geomparams[:CC][:N] - - srange = Ncc - trange = 1 - - msource = model[:CC] - mtarget = model[:NAM] - - psource = parameters[:CC] - ptarget = parameters[:NAM] - - # Here, the indexing in BoundaryFaces is used - couplingfaces = Array{Int64}(undef, 1, 2) - couplingfaces[1, 1] = 2 - couplingfaces[1, 2] = 1 - - couplingcells = Array{Int64}(undef, 1, 2) - couplingcells[1, 1] = Ncc - couplingcells[1, 2] = 1 - - trans = getTrans_1d(msource, mtarget, - couplingfaces, - couplingcells, - psource, ptarget, - :Conductivity) - - ct = TPFAInterfaceFluxCT(trange, srange, trans) - ct_pair = setup_cross_term(ct, target = :NAM, source = :CC, equation = :charge_conservation) - add_cross_term!(model, ct_pair) - - end - - # setup coupling NAM <-> ELYTE charge - - Nnam = geomparams[:NAM][:N] - - srange = collect(1 : Nnam) # negative electrode - trange = collect(1 : Nnam) # electrolyte (negative side) - - if discretisation_type(model[:NAM]) == :P2Ddiscretization - - ct = ButlerVolmerActmatToElyteCT(trange, srange) - ct_pair = setup_cross_term(ct, target = :ELYTE, source = :NAM, equation = :charge_conservation) - add_cross_term!(model, ct_pair) - ct_pair = setup_cross_term(ct, target = :ELYTE, source = :NAM, equation = :mass_conservation) - add_cross_term!(model, ct_pair) - - ct = ButlerVolmerElyteToActmatCT(srange, trange) - ct_pair = setup_cross_term(ct, target = :NAM, source = :ELYTE, equation = :charge_conservation) - add_cross_term!(model, ct_pair) - ct_pair = setup_cross_term(ct, target = :NAM, source = :ELYTE, equation = :solid_diffusion_bc) - add_cross_term!(model, ct_pair) - - else - - @assert discretisation_type(model[:NAM]) == :NoParticleDiffusion - - ct = ButlerVolmerInterfaceFluxCT(trange, srange) - ct_pair = setup_cross_term(ct, target = :ELYTE, source = :NAM, equation = :charge_conservation) - add_cross_term!(model, ct_pair) - ct_pair = setup_cross_term(ct, target = :ELYTE, source = :NAM, equation = :mass_conservation) - add_cross_term!(model, ct_pair) - - end - - # setup coupling ELYTE <-> PAM charge - - Nnam = geomparams[:NAM][:N] - Nsep = geomparams[:SEP][:N] - Npam = geomparams[:PAM][:N] - - srange = collect(1 : Npam) # positive electrode - trange = collect(Nnam + Nsep .+ (1 : Npam)) # electrolyte (positive side) - - if discretisation_type(model[:PAM]) == :P2Ddiscretization - - ct = ButlerVolmerActmatToElyteCT(trange, srange) - ct_pair = setup_cross_term(ct, target = :ELYTE, source = :PAM, equation = :charge_conservation) - add_cross_term!(model, ct_pair) - ct_pair = setup_cross_term(ct, target = :ELYTE, source = :PAM, equation = :mass_conservation) - add_cross_term!(model, ct_pair) - - ct = ButlerVolmerElyteToActmatCT(srange, trange) - ct_pair = setup_cross_term(ct, target = :PAM, source = :ELYTE, equation = :charge_conservation) - add_cross_term!(model, ct_pair) - ct_pair = setup_cross_term(ct, target = :PAM, source = :ELYTE, equation = :solid_diffusion_bc) - add_cross_term!(model, ct_pair) - - else - - @assert discretisation_type(model[:PAM]) == :NoParticleDiffusion - - ct = ButlerVolmerInterfaceFluxCT(trange, srange) - ct_pair = setup_cross_term(ct, target = :ELYTE, source = :PAM, equation = :charge_conservation) - add_cross_term!(model, ct_pair) - ct_pair = setup_cross_term(ct, target = :ELYTE, source = :PAM, equation = :mass_conservation) - add_cross_term!(model, ct_pair) - + initState[name] = init + end - - if !skip_cc - # setup coupling PP <-> PAM charge - - Npam = geomparams[:PAM][:N] - - srange = 1 - trange = Npam - - msource = model[:PP] - mtarget = model[:PAM] - - psource = parameters[:PP] - ptarget = parameters[:PAM] + function initialize_electrolyte!(initState) - # Here, the indexing in BoundaryFaces is used - couplingfaces = Array{Int64}(undef, 1, 2) - couplingfaces[1, 1] = 1 - couplingfaces[1, 2] = 2 - - couplingcells = Array{Int64}(undef, 1, 2) - couplingcells[1, 1] = 1 - couplingcells[1, 2] = Npam - + init = Dict() - trans = getTrans_1d(msource, mtarget, - couplingfaces, - couplingcells, - psource, ptarget, - :Conductivity) + init[:Phi] = state0["Electrolyte"]["phi"][1] + init[:C] = state0["Electrolyte"]["c"][1] - ct = TPFAInterfaceFluxCT(trange, srange, trans) - ct_pair = setup_cross_term(ct, - target = :PAM, - source = :PP, - equation = :charge_conservation) - - add_cross_term!(model, ct_pair) - - # setup coupling with control - - Npp = geomparams[:PP][:N] - - trange = Npp - srange = Int64.(ones(size(trange))) + initState[:ELYTE] = init - msource = model[:PP] - mparameters = parameters[:PP] - # Here the indexing in BoundaryFaces in used - couplingfaces = 2 - couplingcells = Npp - trans = getHalfTrans_1d(msource, couplingfaces, couplingcells, mparameters, :Conductivity) + end - ct = TPFAInterfaceFluxCT(trange, srange, trans, symmetric = false) - ct_pair = setup_cross_term(ct, target = :PP, source = :BPP, equation = :charge_conservation) - add_cross_term!(model, ct_pair) + function initialize_bpp!(initState) - # Accmulation of charge - ct = AccumulatorInterfaceFluxCT(1, trange, trans) - ct_pair = setup_cross_term(ct, target = :BPP, source = :PP, equation = :charge_conservation) - add_cross_term!(model, ct_pair) - + init = Dict(:Phi => state0["Control"]["E"], :Current => 0*state0["Control"]["I"]) + + initState[:BPP] = init + end -end + initState = Dict() + initialize_current_collector!(initState, :CC) + initialize_active_material!(initState, :NAM) + initialize_electrolyte!(initState) + initialize_active_material!(initState, :PAM) + initialize_current_collector!(initState, :PP) + initialize_bpp!(initState) -function setup_coupling!(model, exported_all) - # setup coupling CC <-> NAM :charge_conservation - - skip_pp = size(exported_all["model"]["include_current_collectors"]) == (0,0) - skip_cc = false - if !skip_cc - - srange = Int64.( - exported_all["model"]["NegativeElectrode"]["couplingTerm"]["couplingcells"][:, 1] - ) - trange = Int64.( - exported_all["model"]["NegativeElectrode"]["couplingTerm"]["couplingcells"][:, 2] - ) - - msource = exported_all["model"]["NegativeElectrode"]["CurrentCollector"] - mtarget = exported_all["model"]["NegativeElectrode"]["ActiveMaterial"] - couplingfaces = Int64.(exported_all["model"]["NegativeElectrode"]["couplingTerm"]["couplingfaces"]) - couplingcells = Int64.(exported_all["model"]["NegativeElectrode"]["couplingTerm"]["couplingcells"]) - trans = getTrans(msource, mtarget, couplingfaces, couplingcells, "EffectiveElectricalConductivity") + initState = setup_state(model, initState) - ct = TPFAInterfaceFluxCT(trange, srange, trans) - ct_pair = setup_cross_term(ct, target = :NAM, source = :CC, equation = :charge_conservation) - add_cross_term!(model, ct_pair) - - end + return initState - # setup coupling NAM <-> ELYTE charge +end - srange = Int64.(exported_all["model"]["couplingTerms"][1]["couplingcells"][:, 1]) # negative electrode - trange = Int64.(exported_all["model"]["couplingTerms"][1]["couplingcells"][:, 2]) # electrolyte (negative side) +function setup_battery_initial_state(init::JSONFile, model::MultiModel) - if discretisation_type(model[:NAM]) == :P2Ddiscretization + jsonstruct=init.object - ct = ButlerVolmerActmatToElyteCT(trange, srange) - ct_pair = setup_cross_term(ct, target = :ELYTE, source = :NAM, equation = :charge_conservation) - add_cross_term!(model, ct_pair) - ct_pair = setup_cross_term(ct, target = :ELYTE, source = :NAM, equation = :mass_conservation) - add_cross_term!(model, ct_pair) - - ct = ButlerVolmerElyteToActmatCT(srange, trange) - ct_pair = setup_cross_term(ct, target = :NAM, source = :ELYTE, equation = :charge_conservation) - add_cross_term!(model, ct_pair) - ct_pair = setup_cross_term(ct, target = :NAM, source = :ELYTE, equation = :solid_diffusion_bc) - add_cross_term!(model, ct_pair) - + if haskey(model.models, :CC) + use_cc = true else - - @assert discretisation_type(model[:NAM]) == :NoParticleDiffusion - - ct = ButlerVolmerInterfaceFluxCT(trange, srange) - ct_pair = setup_cross_term(ct, target = :ELYTE, source = :NAM, equation = :charge_conservation) - add_cross_term!(model, ct_pair) - ct_pair = setup_cross_term(ct, target = :ELYTE, source = :NAM, equation = :mass_conservation) - add_cross_term!(model, ct_pair) - + use_cc = false end - # setup coupling ELYTE <-> PAM charge - - srange = Int64.(exported_all["model"]["couplingTerms"][2]["couplingcells"][:,1]) # postive electrode - trange = Int64.(exported_all["model"]["couplingTerms"][2]["couplingcells"][:,2]) # electrolyte (positive side) - - if discretisation_type(model[:PAM]) == :P2Ddiscretization + T = jsonstruct["initT"] + SOC = jsonstruct["SOC"] - ct = ButlerVolmerActmatToElyteCT(trange, srange) - ct_pair = setup_cross_term(ct, target = :ELYTE, source = :PAM, equation = :charge_conservation) - add_cross_term!(model, ct_pair) - ct_pair = setup_cross_term(ct, target = :ELYTE, source = :PAM, equation = :mass_conservation) - add_cross_term!(model, ct_pair) - - ct = ButlerVolmerElyteToActmatCT(srange, trange) - ct_pair = setup_cross_term(ct, target = :PAM, source = :ELYTE, equation = :charge_conservation) - add_cross_term!(model, ct_pair) - ct_pair = setup_cross_term(ct, target = :PAM, source = :ELYTE, equation = :solid_diffusion_bc) - add_cross_term!(model, ct_pair) + function setup_init_am(name, model) - else + theta0 = model[name].system[:theta0] + theta100 = model[name].system[:theta100] + cmax = model[name].system[:maximum_concentration] + N = model[name].system.discretization[:N] - @assert discretisation_type(model[:PAM]) == :NoParticleDiffusion + theta = SOC*(theta100 - theta0) + theta0; + c = theta*cmax + nc = count_entities(model[name].data_domain, Cells()) + init = Dict() + init[:Cs] = c*ones(nc) + init[:Cp] = c*ones(nc, N) - ct = ButlerVolmerInterfaceFluxCT(trange, srange) - ct_pair = setup_cross_term(ct, target = :ELYTE, source = :PAM, equation = :charge_conservation) - add_cross_term!(model, ct_pair) - ct_pair = setup_cross_term(ct, target = :ELYTE, source = :PAM, equation = :mass_conservation) - add_cross_term!(model, ct_pair) + OCP = model[name].system[:ocp_func](c, T, cmax) + return (init, nc, OCP) end - - if !skip_cc - # setup coupling PP <-> PAM charge - target = Dict( - :model => :PAM, - :equation => :charge_conservation - ) - source = Dict( - :model => :PP, - :equation => :charge_conservation - ) - srange = Int64.( - exported_all["model"]["PositiveElectrode"]["couplingTerm"]["couplingcells"][:,1] - ) - trange = Int64.( - exported_all["model"]["PositiveElectrode"]["couplingTerm"]["couplingcells"][:,2] - ) - msource = exported_all["model"]["PositiveElectrode"]["CurrentCollector"] - ct = exported_all["model"]["PositiveElectrode"]["couplingTerm"] - couplingfaces = Int64.(ct["couplingfaces"]) - couplingcells = Int64.(ct["couplingcells"]) - trans = getTrans(msource, mtarget, couplingfaces, couplingcells, "EffectiveElectricalConductivity") - ct = TPFAInterfaceFluxCT(trange, srange, trans) - ct_pair = setup_cross_term(ct, target = :PAM, source = :PP, equation = :charge_conservation) - add_cross_term!(model, ct_pair) + function setup_cc(name, phi, model) + nc = count_entities(model[name].data_domain, Cells()) + init = Dict(); + init[:Phi] = phi*ones(nc) + return init end - if skip_cc - #setup coupling PP <-> BPP charge - target = Dict( - :model => :PAM, - :equation => :charge_conservation - ) + initState = Dict() - trange = convert_to_int_vector( - exported_all["model"]["PositiveElectrode"]["ActiveMaterial"]["externalCouplingTerm"]["couplingcells"] - ) - srange = Int64.(ones(size(trange))) - msource = exported_all["model"]["PositiveElectrode"]["ActiveMaterial"] - couplingfaces = Int64.(msource["externalCouplingTerm"]["couplingfaces"]) - couplingcells = Int64.(msource["externalCouplingTerm"]["couplingcells"]) - else - #setup coupling PP <-> BPP charge - target = Dict( - :model => :PP, - :equation => :charge_conservation - ) + # Setup initial state in negative active material + + init, nc, negOCP = setup_init_am(:NAM, model) + init[:Phi] = zeros(nc) + initState[:NAM] = init + + # Setup initial state in electrolyte + + nc = count_entities(model[:ELYTE].data_domain, Cells()) + + init = Dict() + init[:C] = 1000*ones(nc) + init[:Phi] = - negOCP*ones(nc) - trange = convert_to_int_vector( - exported_all["model"]["PositiveElectrode"]["CurrentCollector"]["externalCouplingTerm"]["couplingcells"] - ) - srange = Int64.(ones(size(trange))) - msource = exported_all["model"]["PositiveElectrode"]["CurrentCollector"] - couplingfaces = Int64.(msource["externalCouplingTerm"]["couplingfaces"]) - couplingcells = Int64.(msource["externalCouplingTerm"]["couplingcells"]) - end + initState[:ELYTE] = init + + # Setup initial state in positive active material - source = Dict( - :model => :BPP, - :equation => :charge_conservation - ) + init, nc, posOCP = setup_init_am(:PAM, model) + init[:Phi] = (posOCP - negOCP)*ones(nc) - #effcond = exported_all["model"]["PositiveElectrode"]["CurrentCollector"]["EffectiveElectricalConductivity"] - trans = getHalfTrans(msource, couplingfaces, couplingcells, "EffectiveElectricalConductivity") + initState[:PAM] = init - if skip_cc - - ct = TPFAInterfaceFluxCT(trange, srange, trans, symmetric = false) - ct_pair = setup_cross_term(ct, target = :PAM, source = :BPP, equation = :charge_conservation) - add_cross_term!(model, ct_pair) + # Setup negative current collector + + initState[:CC] = setup_cc(:CC, 0, model) - # Accmulation of charge - ct = AccumulatorInterfaceFluxCT(1, trange, trans) - ct_pair = setup_cross_term(ct, target = :BPP, source = :PAM, equation = :charge_conservation) - add_cross_term!(model, ct_pair) - - else - - ct = TPFAInterfaceFluxCT(trange, srange, trans, symmetric = false) - ct_pair = setup_cross_term(ct, target = :PP, source = :BPP, equation = :charge_conservation) - add_cross_term!(model, ct_pair) + # Setup positive current collector - # Accmulation of charge - ct = AccumulatorInterfaceFluxCT(1, trange, trans) - ct_pair = setup_cross_term(ct, target = :BPP, source = :PP, equation = :charge_conservation) - add_cross_term!(model, ct_pair) + initState[:PP] = setup_cc(:PP, posOCP - negOCP, model) + + init = Dict() + init[:Phi] = [1.0] + init[:Current] = [1.0] - end + initState[:BPP] = init + + initState = setup_state(model, initState) + + return initState end -function currentFun(t::T, inputI::T) where T - #inputI = 9.4575 - tup = 0.1 + +################################################################################## +#Current function +################################################################################## +# function currentFun(t::T, inputI::T) where T +# #inputI = 9.4575 +# tup = 0.1 +# val::T = 0.0 +# if t <= tup +# val = sineup(0.0, inputI, 0.0, tup, t) +# else +# val = inputI +# end +# return val +# end + + +function currentFun(t::T, inputI::T, tup::T=0.1) where T val::T = 0.0 if t <= tup val = sineup(0.0, inputI, 0.0, tup, t) @@ -1313,42 +1458,176 @@ function currentFun(t::T, inputI::T) where T return val end +################################################################################## +#Setup volume fraction +################################################################################## +function setup_volume_fractions!(model::MultiModel, geomparams::Dict{Symbol,<:Any}) -function currentFun(t::T, inputI::T, tup::T) where T - val::T = 0.0 - if t <= tup - val = sineup(0.0, inputI, 0.0, tup, t) - else - val = inputI + names = (:NAM, :SEP, :PAM) + Nelyte = sum([geomparams[name][:N] for name in names]) + vfelyte = zeros(Nelyte) + + names = (:NAM, :PAM) + + for name in names + ammodel = model[name] + vf = ammodel.system[:volumeFraction] + Nam = geomparams[name][:N] + ammodel.domain.representation[:volumeFraction] = vf*ones(Nam) + if name == :NAM + nstart = 1 + nend = Nam + elseif name == :PAM + nstart = geomparams[:NAM][:N] + geomparams[:SEP][:N] + 1 + nend = Nelyte + else + error("name not recognized") + end + vfelyte[nstart : nend] .= 1 - vf end - return val + + nstart = geomparams[:NAM][:N] + 1 + nend = nstart + geomparams[:SEP][:N] + separator_porosity = model[:ELYTE].system[:separator_porosity] + vfelyte[nstart : nend] .= separator_porosity*ones(nend - nstart + 1) + + model[:ELYTE].domain.representation[:volumeFraction] = vfelyte + end -function amg_precond(; max_levels = 10, max_coarse = 10, type = :smoothed_aggregation) +################################################################################## +#Transmissibilities +################################################################################## + +function getTrans(model1::Dict{String,<:Any}, + model2::Dict{String, Any}, + faces, + cells, + quantity::String) + """ setup transmissibility for coupling between models at boundaries""" + + T_all1 = model1["operators"]["T_all"][faces[:, 1]] + T_all2 = model2["operators"]["T_all"][faces[:, 2]] + + s1 = model1[quantity][cells[:, 1]] + s2 = model2[quantity][cells[:, 2]] - gs_its = 1 - cyc = AlgebraicMultigrid.V() - if type == :smoothed_aggregation - m = smoothed_aggregation - else - m = ruge_stuben + T = 1.0./((1.0./(T_all1.*s1))+(1.0./(T_all2.*s2))) + + return T + +end + +function getTrans(model1::Jutul.SimulationModel, + model2::Jutul.SimulationModel, + bcfaces, + bccells, + parameters1, + parameters2, + quantity) + """ setup transmissibility for coupling between models at boundaries. Intermediate 1d version""" + + d1 = physical_representation(model1) + d2 = physical_representation(model2) + + bcTrans1 = d1[:bcTrans][bcfaces[:, 1]] + bcTrans2 = d2[:bcTrans][bcfaces[:, 2]] + cells1 = bccells[:, 1] + cells2 = bccells[:, 2] + + s1 = parameters1[quantity][cells1] + s2 = parameters2[quantity][cells2] + + T = 1.0./((1.0./(bcTrans1.*s1))+(1.0./(bcTrans2.*s2))) + + return T + +end + +function getHalfTrans(model::Jutul.SimulationModel, + bcfaces, + bccells, + parameters, + quantity) + """ recover half transmissibilities for boundary faces and weight them by the coefficient sent as quantity for the corresponding given cells. Intermediate 1d version. Note the indexing in BoundaryFaces is used""" + + d = physical_representation(model) + bcTrans = d[:bcTrans][bcfaces] + s = parameters[quantity][bccells] + + T = bcTrans.*s + + return T +end + +function getHalfTrans(model::Dict{String, Any}, + faces, + cells, + quantity::String) + """ recover half transmissibilities for boundary faces and weight them by the coefficient sent as quantity for the given cells. + Here, the faces should belong the corresponding cells at the same index""" + + T_all = model["operators"]["T_all"] + s = model[quantity][cells] + T = T_all[faces].*s + + return T + +end + +function getHalfTrans(model::Dict{String,<:Any}, + faces) + """ recover the half transmissibilities for boundary faces""" + + T_all = model["operators"]["T_all"] + T = T_all[faces] + + return T + +end + +################################################################################## +#Setup geomparams +################################################################################## + +function setup_geomparams(init::JSONFile) + + jsondict=init.object + + names = (:CC, :NAM, :SEP, :PAM, :PP) + geomparams = Dict(name => Dict() for name in names) + + geomparams[:CC][:N] = jsondict["NegativeElectrode"]["CurrentCollector"]["N"] + geomparams[:CC][:thickness] = jsondict["NegativeElectrode"]["CurrentCollector"]["thickness"] + geomparams[:NAM][:N] = jsondict["NegativeElectrode"]["ActiveMaterial"]["N"] + geomparams[:NAM][:thickness] = jsondict["NegativeElectrode"]["ActiveMaterial"]["thickness"] + geomparams[:SEP][:N] = jsondict["Electrolyte"]["Separator"]["N"] + geomparams[:SEP][:thickness] = jsondict["Electrolyte"]["Separator"]["thickness"] + geomparams[:PAM][:N] = jsondict["PositiveElectrode"]["ActiveMaterial"]["N"] + geomparams[:PAM][:thickness] = jsondict["PositiveElectrode"]["ActiveMaterial"]["thickness"] + geomparams[:PP][:N] = jsondict["PositiveElectrode"]["CurrentCollector"]["N"] + geomparams[:PP][:thickness] = jsondict["PositiveElectrode"]["CurrentCollector"]["thickness"] + + for name in names + geomparams[name][:facearea] = jsondict["Geometry"]["faceArea"] end - gs = GaussSeidel(ForwardSweep(), gs_its) - return AMGPreconditioner(m, max_levels = max_levels, max_coarse = max_coarse, presmoother = gs, postsmoother = gs, cycle = cyc) + return geomparams end +################################################################################## +#Compute cell capactity +################################################################################## -function computeCellCapacity(model) +function computeCellCapacity(model::MultiModel) con = Constants() - - function computeHalfCellCapacity(name) + + function computeHalfCellCapacity(name::Symbol) ammodel = model[name] - sys = ammodel.system - + sys = ammodel.system F = con.F n = sys[:n_charge_carriers] cMax = sys[:maximum_concentration] @@ -1379,226 +1658,75 @@ function computeCellCapacity(model) end -function setup_sim_1d(jsondict; use_groups = false, general_ad = false) - - model, state0, parameters = setup_model_1d(jsondict, use_groups = use_groups, general_ad = general_ad) - - setup_coupling_1d!(model, parameters, jsondict) - - minE = jsondict["Control"]["lowerCutoffVoltage"] - - CRate = jsondict["Control"]["CRate"] - cap = computeCellCapacity(model) - con = Constants() - - inputI = (cap/con.hour)*CRate - - @. state0[:BPP][:Phi] = minE*1.5 - - tup = Float64(jsondict["TimeStepping"]["rampupTime"]) - cFun(time) = currentFun(time, inputI, tup) - - currents = setup_forces(model[:BPP], policy = SimpleCVPolicy(cFun, minE)) - forces = setup_forces(model, BPP = currents) - - sim = Simulator(model, state0 = state0, parameters = parameters, copy_state = true) - - return sim, forces, state0, parameters, model - -end - -function setup_sim(name; use_p2d = true, use_groups = false, general_ad = false) +################################################################################### +#Other utils +################################################################################### - fn = string(dirname(pathof(BattMo)), "/../test/battery/data/", name, ".mat") - exported = MAT.matread(fn) - - model, state0, parameters = setup_model(exported, use_p2d = use_p2d, use_groups = use_groups, general_ad = general_ad) - setup_coupling!(model, exported) - - inputI = 0; - minE = 10 - steps = size(exported["states"],1) - - for i = 1:steps - - inputI = max(inputI, exported["states"][i]["Control"]["I"]) - minE = min(minE, exported["states"][i]["Control"]["E"]) - +struct Constants + F + R + hour + function Constants() + new(96485.3329, + 8.31446261815324, + 3600) end - - @. state0[:BPP][:Phi] = minE*1.5 - cFun(time) = currentFun(time, inputI) - forces_pp = nothing - - currents = setup_forces(model[:BPP], policy = SimpleCVPolicy(cFun, minE)) - - forces = Dict( - :CC => nothing, - :NAM => nothing, - :ELYTE => nothing, - :PAM => nothing, - :PP => forces_pp, - :BPP => currents - ) - - sim = Simulator(model, state0 = state0, parameters = parameters, copy_state = true) - - return sim, forces, state0, parameters, exported, model - end -export run_battery - -function run_battery(name; - use_p2d = true, - extra_timing = false, - max_step = nothing, - linear_solver = :direct, - general_ad = false, - use_groups = false, - kwarg...) - - sim, forces, state0, parameters, exported, model = setup_sim(name, use_p2d = use_p2d, use_groups = use_groups, general_ad = general_ad) - - steps = size(exported["states"], 1) - alltimesteps = Vector{Float64}(undef, steps) - time = 0; - end_step = 0 - minE = 3.2 - - for i = 1 : steps - alltimesteps[i] = exported["states"][i]["time"] - time - time = exported["states"][i]["time"] - E = exported["states"][i]["Control"]["E"] - if (E > minE + 0.001) - end_step = i - end - end - if !isnothing(max_step) - end_step = min(max_step, end_step) - end - - timesteps = alltimesteps[1 : end_step] - - cfg = simulator_config(sim; kwarg...) - cfg[:linear_solver] = battery_linsolve(model, linear_solver) - cfg[:debug_level] = 0 - #cfg[:max_timestep_cuts] = 0 - cfg[:max_residual] = 1e20 - cfg[:min_nonlinear_iterations] = 1 - cfg[:extra_timing] = extra_timing - # cfg[:max_nonlinear_iterations] = 5 - cfg[:safe_mode] = false - cfg[:error_on_incomplete] = true - if false - cfg[:info_level] = 5 - cfg[:max_nonlinear_iterations] = 1 - cfg[:max_timestep_cuts] = 0 +struct SourceAtCell + cell + src + function SourceAtCell(cell, src) + new(cell, src) end - - cfg[:tolerances][:PP][:default] = 1e-1 - cfg[:tolerances][:BPP][:default] = 1e-1 - # Run simulation - - states, report = simulate(sim, timesteps, forces = forces, config = cfg) - stateref = exported["states"] - - extra = Dict(:model => model, - :state0 => state0, - :states_ref => stateref, - :parameters => parameters, - :exported => exported, - :timesteps => timesteps, - :config => cfg, - :forces => forces, - :simulator => sim) - - return (states = states, reports = report, extra = extra, exported = exported) end - -function rampupTimesteps(time, dt, n = 8) +function rampupTimesteps(time::Real, dt::Real, n::Integer=8) ind = [8; collect(range(n, 1, step=-1))] - dt_init = [dt/2^k for k in ind] + dt_init = [dt / 2^k for k in ind] cs_time = cumsum(dt_init) if any(cs_time .> time) - dt_init = dt_init[cs_time .< time]; + dt_init = dt_init[cs_time. model, - :state0 => state0, - :parameters => parameters, - :timesteps => timesteps, - :config => cfg, - :forces => forces, - :simulator => sim) +end - return (states = states, reports = reports, extra = extra) +function convert_to_int_vector(x::Float64) + vec = Int64.(Vector{Float64}([x])) + return vec end +function convert_to_int_vector(x::Matrix{Float64}) + vec = Int64.(Vector{Float64}(x[:, 1])) + return vec +end -export inputRefToStates function inputRefToStates(states, stateref) statesref = deepcopy(states); for i in 1:size(states,1) @@ -1656,3 +1784,18 @@ function test_mrst_battery(name) end + +function amg_precond(; max_levels = 10, max_coarse = 10, type = :smoothed_aggregation) + + gs_its = 1 + cyc = AlgebraicMultigrid.V() + if type == :smoothed_aggregation + m = smoothed_aggregation + else + m = ruge_stuben + end + gs = GaussSeidel(ForwardSweep(), gs_its) + + return AMGPreconditioner(m, max_levels = max_levels, max_coarse = max_coarse, presmoother = gs, postsmoother = gs, cycle = cyc) + +end \ No newline at end of file diff --git a/src/types.jl b/src/types.jl new file mode 100644 index 0000000..42ceba0 --- /dev/null +++ b/src/types.jl @@ -0,0 +1,34 @@ +import JSON +import MAT + +export InputFile, MatlabFile, JSONFile + +######################################################## +#Types to simplify reading data from different formats +######################################################## +abstract type InputFile end + +struct MatlabFile <: InputFile + source::String + object::Dict{String,Any} + use_state_ref::Bool + function MatlabFile(file::String; + state_ref::Bool=true) #Read from stored mat file + + return new(file, MAT.matread(file),state_ref) + end + function MatlabFile(file::String, + object::Dict{String,Any}; + state_ref::Bool=false) #Used when we launch julia from matlab + + return new(file,object,state_ref) + end +end + +struct JSONFile <: InputFile + source::String + object::Dict{String,Any} + function JSONFile(file::String) + return new(file,JSON.parsefile(file)) + end +end From d4960e336a5192c60b0629b93631cd1b5d31b9bd Mon Sep 17 00:00:00 2001 From: Erasdna Date: Wed, 23 Aug 2023 13:51:37 +0200 Subject: [PATCH 2/2] test compatibility --- examples/battery/example_battery.jl | 11 ++++++----- examples/battery/example_cycle.jl | 21 +++++++++++---------- src/BattMo.jl | 2 +- src/mrst_test_utils.jl | 9 +++++---- test/battery/test_battery_boundary.jl | 2 +- 5 files changed, 24 insertions(+), 21 deletions(-) diff --git a/examples/battery/example_battery.jl b/examples/battery/example_battery.jl index a3b6ee3..e8bfb3c 100644 --- a/examples/battery/example_battery.jl +++ b/examples/battery/example_battery.jl @@ -43,11 +43,12 @@ end # sim, forces, state0, parameters, exported, model = BattMo.setup_sim(name, use_p2d = use_p2d) -do_json = true +do_json = false # true if do_json - - states, reports, extra = run_battery_1d(info_level = 0, extra_timing = false); + fn = string(dirname(pathof(BattMo)), "/../test/battery/data/jsonfiles/", name, "_jl.json") + init = JSONFile(fn) + states, reports, extra = run_battery(init;use_p2d=use_p2d,info_level = 0, extra_timing = false); nsteps = size(states) @@ -66,7 +67,7 @@ if do_json else - states, reports, extra, exported = run_battery(name, use_p2d = use_p2d, info_level = 5, max_step = nothing); + states, reports, extra, exported = run_battery(name; use_p2d = use_p2d, info_level = 0, max_step = nothing); nsteps = size(states, 1) @@ -75,7 +76,7 @@ else time = cumsum(timesteps[1 : nsteps]) E = [state[:BPP][:Phi][1] for state in states] - statesref = extra[:states_ref] + statesref = extra[:init].object["states"] timeref = time Eref = [state["Control"]["E"] for state in statesref[1 : nsteps]] diff --git a/examples/battery/example_cycle.jl b/examples/battery/example_cycle.jl index 8f4dd39..b3cde4a 100644 --- a/examples/battery/example_cycle.jl +++ b/examples/battery/example_cycle.jl @@ -1,12 +1,13 @@ using BattMo -name="model1D_50" +#name="model1D_50" +name = "p2d_40" # Run base case and plot the results against BattMo-MRST reference states, reports, extra = run_battery(name, info_level = 1, max_step = nothing); prm = extra[:parameters] model = extra[:model] timesteps = extra[:timesteps] steps = size(states, 1) -stateref = extra[:states_ref] +stateref = extra[:init].object["states"] E = Matrix{Float64}(undef,steps,2) for step in 1:steps phi = states[step][:BPP][:Phi][1] @@ -22,15 +23,15 @@ Plots.plot!(T, E[:, 2], label = "BattMo") closeall() display(plot1) ## Compute and plot the charge flux -f = :ELYTE -ix = 1 -get_j = ix -> output_flux(model[f], states[ix][f], prm[f], :charge_conservation) +# f = :ELYTE +# ix = 1 +# get_j = ix -> output_flux(model[f], states[ix][f], prm[f], :charge_conservation) -plot2 = Plots.plot(get_j(1)) -for ix in 5:5:length(states) - Plots.plot!(get_j(ix)) -end -display(plot2) +# plot2 = Plots.plot(get_j(1)) +# for ix in 5:5:length(states) +# Plots.plot!(get_j(ix)) +# end +# display(plot2) ## Set up a cycle policy and simulate using Jutul diff --git a/src/BattMo.jl b/src/BattMo.jl index 5fca732..fc03dc7 100644 --- a/src/BattMo.jl +++ b/src/BattMo.jl @@ -78,7 +78,7 @@ include("linsolve.jl") # Precompilation of solver. Run a small battery simulation to precompile everything. @compile_workload begin for use_general_ad in [false, true] - run_battery("p2d_40", general_ad = use_general_ad,info_level = -1); + run_battery("p2d_40"; general_ad = use_general_ad,info_level = -1); end end diff --git a/src/mrst_test_utils.jl b/src/mrst_test_utils.jl index 0c5f03d..8b95b64 100644 --- a/src/mrst_test_utils.jl +++ b/src/mrst_test_utils.jl @@ -22,7 +22,7 @@ function run_battery(init::InputFile; """ #Setup simulation - sim, forces, state0, parameters, init, model = setup_sim(init, use_p2d=use_p2d, use_groups=use_groups, general_ad=general_ad; kwarg...) + sim, forces, state0, parameters, init, model = setup_sim(init, use_p2d=use_p2d, use_groups=use_groups, general_ad=general_ad) #Set up config and timesteps timesteps=setup_timesteps(init;max_step=max_step) @@ -197,7 +197,7 @@ function setup_sim(init::JSONFile; currents = setup_forces(model[:BPP], policy=SimpleCVPolicy(cFun, minE)) forces = setup_forces(model, BPP=currents) - sim = Simulator(model, state0=state0, parameters=parameters, copy_state=true) + sim = Simulator(model; state0=state0, parameters=parameters, copy_state=true) return sim, forces, state0, parameters, init, model @@ -206,7 +206,8 @@ end function setup_sim(init::MatlabFile; use_p2d::Bool=true, use_groups::Bool=false, - general_ad::Bool=false) + general_ad::Bool=false, + kwarg ... ) model, state0, parameters = setup_model(init, use_p2d=use_p2d, use_groups=use_groups, general_ad=general_ad) setup_coupling!(init,model) @@ -236,7 +237,7 @@ function setup_sim(init::MatlabFile; :BPP => currents ) - sim = Simulator(model, state0=state0, parameters=parameters, copy_state=true) + sim = Simulator(model; state0=state0, parameters=parameters, copy_state=true) return sim, forces, state0, parameters, init, model diff --git a/test/battery/test_battery_boundary.jl b/test/battery/test_battery_boundary.jl index 9234b31..88a5436 100644 --- a/test/battery/test_battery_boundary.jl +++ b/test/battery/test_battery_boundary.jl @@ -11,7 +11,7 @@ testcases =[ @testset "$modelname" begin for use_general_ad in [false, true] states, report, extra = run_battery(modelname, info_level = -1, general_ad = use_general_ad); - stateref = extra[:states_ref] + stateref = extra[:init].object["states"] steps = size(states, 1) E = Matrix{Float64}(undef,steps,2) for step in 1:steps