begin
    using PlutoUI,HypertextLiteral,UUIDs, Markdown
    using LinearAlgebra
    using LessUnitful
    using VoronoiFVM
    using LiquidElectrolytes
    using ExtendableGrids
    using LessUnitful
    using GridVisualize
    if isdefined(Main,:PlutoRunner)
        using Colors
    using CairoMakie	
    default_plotter!(CairoMakie)
    CairoMakie.activate!(type="svg")
    end
end;

pkgdir(LiquidElectrolytes)

"/home/runner/work/LiquidElectrolytes.jl/LiquidElectrolytes.jl"

begin
    const nm=ufac"nm"
    const V=ufac"V"
    const mol=ufac"mol"
    const dm=ufac"dm"
end;

grid = let
    hmin = 1.0e-1 * nm * 2.0^(-nref)
    hmax = 1.0 * nm * 2.0^(-nref)
    L = 20.0 * nm
    X = geomspace(0, L, hmin, hmax)
    simplexgrid(X)
end

ExtendableGrids.ExtendableGrid{Float64, Int32};
dim: 1 nodes: 103 cells: 102 bfaces: 2

gridplot(grid,size=(600,200))

begin
    const molarity=0.01
    const voltages=-1:0.01:1
    const nref=1
    const solvation=10
    celldata = ElectrolyteData(;
                               nc = 2,
                               Γ_we = 1,
                               Γ_bulk = 2)
    
    celldata.c_bulk.=molarity * mol / dm^3
    celldata.κ.=solvation
    celldata
end

                  nc = 2
                  na = 0
                  iϕ = 3
                  ip = 4
                   D = [2.0e-9, 2.0e-9]
                   z = [1, -1]
                  M0 = 0.018015
                   M = [0.018015, 0.018015]
                  v0 = 1.8051e-5
                   v = [1.8051e-5, 1.8051e-5]
                   κ = [10.0, 10.0]
              c_bulk = [10.0, 10.0]
              ϕ_bulk = 0.0
              p_bulk = 0.0
              Γ_bulk = 2
                ϕ_we = 0.0
                Γ_we = 1
                   T = 298.15
                  RT = 2479.0
                   F = 96485.0
                   ε = 78.49
                 ε_0 = 8.8542e-12
              pscale = 1.0e9
            eneutral = false
              scheme = μex
              epsreg = 1.0e-20
             weights = [1.8051e-5, 1.8051e-5, 1.0, 0.0]

    ## Define boundary conditions 
    function bcondition(f,u,bnode,data)
    (;Γ_we,Γ_bulk,ϕ_we) = data
       	iϕ,ip=1,2 
    ## Dirichlet ϕ=ϕ_we at Γ_we
    boundary_dirichlet!(f,u,bnode,species=iϕ,region=Γ_we,value=ϕ_we)
    boundary_dirichlet!(f,u,bnode,species=iϕ,region=Γ_bulk,value=data.ϕ_bulk)
    boundary_dirichlet!(f,u,bnode,species=ip,region=Γ_bulk,value=data.p_bulk)
    end

bcondition (generic function with 1 method)

pbsys=PBSystem(grid;celldata=deepcopy(celldata),bcondition)

VoronoiFVM.System{Float64, Float64, Int32, Int64, Matrix{Int32}}(num_species=2)

result=dlcapsweep(pbsys;inival=unknowns(pbsys),iϕ=1,voltages,molarity=molarity*ufac"mol/dm^3")

DLCapSweepResult{Vector{Float64}, Vector{Float64}, Nothing}([-1.0, -0.99, -0.98, -0.97, -0.96, -0.95, -0.94, -0.93, -0.92, -0.91  …  0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.0], [0.46231756320458395, 0.46518147068130844, 0.46809693018068366, 0.4710654612349874, 0.4740886459253524, 0.4771681320547838, 0.4803056367685965, 0.48350295029697676, 0.486761939870739, 0.4900845540156684  …  0.4900510063343688, 0.48672903834923176, 0.4834706750511675, 0.4802739687148616, 0.4771370529277341, 0.47405813818568454, 0.4710355081039008, 0.468067515477788, 0.4651525789178823, 0.46228917948254455], nothing)

begin
    ecell=create_equilibrium_system(grid,EquilibriumData(celldata))
    evolts,ecaps=dlcapsweep_equi(ecell;vmax=voltages[end],molarity,δV=1.0e-4,nsteps=101)
end;

function pnpsweep(scheme) 
        
    function pnp_bcondition(f,u,bnode,data::ElectrolyteData)
        (; iϕ,Γ_we,ϕ_we) = data
        boundary_dirichlet!(f,u,bnode,species=iϕ,region=Γ_we,value=ϕ_we)
        bulkbcondition(f,u,bnode,data)
    end
    xcelldata=deepcopy(celldata)
    xcelldata.scheme=scheme
    μcell=PNPSystem(grid;bcondition=pnp_bcondition,celldata=xcelldata)
    μresult=dlcapsweep(μcell;voltages,molarity=molarity*ufac"mol/dm^3",δ=1.0e-4)
end;

μresult=pnpsweep(:μex)

DLCapSweepResult{Vector{Float64}, Vector{Float64}, Nothing}([-1.0, -0.99, -0.98, -0.97, -0.96, -0.95, -0.94, -0.93, -0.92, -0.91  …  0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.0], [0.4623175632034737, 0.46518147068241866, 0.4680969301817939, 0.4710654612349874, 0.4740886459253524, 0.47716813206144515, 0.4803056367697067, 0.48350295029697676, 0.4867619398696288, 0.4900845540156684  …  0.4900510063343688, 0.48672903834812153, 0.4834706750511675, 0.4802739687137514, 0.47713705292551367, 0.47405813818568454, 0.4710355081016804, 0.4680675154766778, 0.46515257891899253, 0.46228917948254455], nothing)

aresult=pnpsweep(:act)

DLCapSweepResult{Vector{Float64}, Vector{Float64}, Nothing}([-1.0, -0.99, -0.98, -0.97, -0.96, -0.95, -0.94, -0.93, -0.92, -0.91  …  0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.0], [0.46231756320458395, 0.46518147068241866, 0.46809693018068366, 0.47106546123276694, 0.47408864592313193, 0.47716813205589403, 0.4803056367697067, 0.48350295029697676, 0.48676193986851857, 0.4900845540156684  …  0.49005100633547904, 0.486729038350342, 0.4834706750511675, 0.48027396871264116, 0.4771370529266239, 0.47405813818568454, 0.4710355081039008, 0.4680675154766778, 0.46515257891899253, 0.462289179484765], nothing)

cresult=pnpsweep(:cent)

DLCapSweepResult{Vector{Float64}, Vector{Float64}, Nothing}([-1.0, -0.99, -0.98, -0.97, -0.96, -0.95, -0.94, -0.93, -0.92, -0.91  …  0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.0], [0.4623175632068044, 0.46518147068241866, 0.4680969301817939, 0.4710654612349874, 0.47408864592424216, 0.4771681320625554, 0.48030563677081695, 0.48350295029697676, 0.48676193986851857, 0.4900845540156684  …  0.4900510063321484, 0.4867290383514522, 0.4834706750511675, 0.4802739687137514, 0.47713705292440345, 0.47405813818679476, 0.47103550810612127, 0.468067515477788, 0.46515257891899253, 0.46228917948254455], nothing)

myround(x)=round(x,sigdigits=5)

myround (generic function with 1 method)

floataside(let
    vis=GridVisualizer(size=(350,200),legend=:rt)
    scalarplot!(vis,result.voltages,result.dlcaps,color=:red,label="pb")
    
    scalarplot!(vis,evolts,ecaps,color=:green,label="equi",clear=false)	
    scalarplot!(vis,μresult.voltages,μresult.dlcaps,color=:blue,label="μpnp",clear=false)
    scalarplot!(vis,μresult.voltages,aresult.dlcaps,color=:magenta,label="apnp",clear=false)
    scalarplot!(vis,cresult.voltages,cresult.dlcaps,color=:cyan,label="cpnp",clear=false)
    
    reveal(vis)
end, top=200)

let
    pb0=findfirst(x->x≈0,result.voltages)
    equi0=findfirst(x->x≈0, evolts)
    pnp0=findfirst(x->x≈0, μresult.voltages)
    
floataside(
    md"""
DLcap0: $(dlcap0(celldata)|>myround), $(dlcap0(EquilibriumData(celldata))|>myround)

DLCapMin(pb,equi,pnp): $(result.dlcaps[pb0]|>myround,ecaps[equi0]|>myround,μresult.dlcaps[pnp0]|>myround))
    
    """,
    
    top=450

)
end

Any

floataside (generic function with 2 methods)