begin
    using PlutoUI, HypertextLiteral, UUIDs, Markdown
    using Test
    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 = "png")
    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
   nnodes =     103
   ncells =     102
  nbfaces =       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
       solvepressure = true
             weights = [1.8051e-5, 1.8051e-5, 1.0, 0.0]
        edgevelocity = 0.0
               model = DGL

## 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)
    return 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}}(
  grid = ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=1, nnodes=103, ncells=102,
  nbfaces=2),  
  physics = Physics(data=LiquidElectrolytes.ElectrolyteData, flux=pbflux,
  storage=default_storage, reaction=pbreaction, breaction=bcondition, ),  
  num_species = 2)

begin
    @info "PBSystem"
result = dlcapsweep(
    pbsys;
    inival = unknowns(pbsys, inival=0),
    iϕ = 1, voltages,
    molarity = molarity * ufac"mol/dm^3",
    verbose="ne"
)
end

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.4740886459264626, 0.4771681320581145, 0.4803056367685965, 0.48350295029697676, 0.4867619398696288, 0.4900845540156684  …  0.4900510063343688, 0.48672903834923176, 0.4834706750511675, 0.48027396871597183, 0.47713705292440345, 0.47405813818568454, 0.4710355081039008, 0.468067515477788, 0.46515257891899253, 0.4622891794836548], nothing)

begin
    @info "EquilibriumSystem"
    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)
    @info scheme
    function pnp_bcondition(f, u, bnode, data::ElectrolyteData)
        (; iϕ, Γ_we, ϕ_we) = data
        boundary_dirichlet!(f, u, bnode, species = iϕ, region = Γ_we, value = ϕ_we)
        return bulkbcondition(f, u, bnode, data)
    end
    xcelldata = deepcopy(celldata)
    xcelldata.scheme = scheme
    μcell = PNPSystem(grid; bcondition = pnp_bcondition, celldata = xcelldata)
    return μresult = dlcapsweep(
        μcell;
        voltages,
        molarity = molarity * ufac"mol/dm^3",
        δ = 1.0e-4, damp_initial = 0.5,
        verbose="ne"
    )
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.4623175632068044, 0.46518147068241866, 0.46809693018068366, 0.4710654612349874, 0.4740886459264626, 0.4771681320581145, 0.4803056367685965, 0.483502950298087, 0.4867619398696288, 0.4900845540156684  …  0.49005100633547904, 0.48672903834923176, 0.4834706750511675, 0.4802739687148616, 0.47713705292551367, 0.47405813818679476, 0.4710355081039008, 0.4680675154766778, 0.46515257891899253, 0.4622891794836548], 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.4623175632023635, 0.46518147068241866, 0.46809693018068366, 0.47106546123276694, 0.47408864592424216, 0.4771681320603349, 0.4803056367685965, 0.48350295029697676, 0.4867619398696288, 0.4900845540156684  …  0.49005100633547904, 0.486729038350342, 0.4834706750511675, 0.4802739687137514, 0.4771370529266239, 0.4740581381890152, 0.4710355081016804, 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.46231756320569417, 0.46518147068130844, 0.46809693018068366, 0.4710654612316567, 0.47408864592313193, 0.47716813205700426, 0.4803056367685965, 0.48350295029586654, 0.486761939870739, 0.4900845540156684  …  0.4900510063332586, 0.48672903834812153, 0.4834706750511675, 0.4802739687137514, 0.47713705292440345, 0.474058138187905, 0.47103550810612127, 0.4680675154766778, 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, aresult.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
)

@test isapprox(result.dlcaps, ecaps, rtol = 5.0e-4)

Test Passed

@test isapprox(result.dlcaps, μresult.dlcaps, rtol = 1.0e-11)

Test Passed

@test isapprox(result.dlcaps, aresult.dlcaps, rtol = 1.0e-11)

Test Passed

@test isapprox(result.dlcaps, cresult.dlcaps, rtol = 1.0e-11)

Test Passed

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)