begin
   using Pkg
    using Revise
        using GridVisualize

   	import CairoMakie	
        using CairoMakie: lines!
  	 	default_plotter!(CairoMakie)
 		CairoMakie.activate!(type="png")
  
end

begin
    using LiquidElectrolytes
    using Printf
    using ExtendableGrids
    using LinearAlgebra
    using PlutoUI
    using ExtendableFEM
    using VoronoiFVM
    using LessUnitful
    using Test

end

begin
    @phconstants e N_A
    const F = e * N_A
    @unitfactors mol m dm cm nm μA s mV V mPa

end

0.001

begin
    const L = 50 * nm
    const W = 5 * nm
    const nref = 0
    const Δp = 1.0
    const μ = 0.89 * mPa * s
    const cyl = false
    nr = 8 * 2^nref + 1
    nl = Int((nr - 1) * L / W) + 1
    Y = range(0, L, length = nl)
    flowX = range(0, W, length = nr)
    pnpX = geomspace(0, W, 2 * W / nr, 0.25 * W / nr)
    pnpgrid = simplexgrid(pnpX, Y)
    flowgrid = simplexgrid(flowX, Y)
    flowgrid = pnpgrid
    xgrid = simplexgrid(1.0e9 * 2.5 * pnpX, 1.0e9 * Y)
    pnpgrid[CellEdges]
    pnpgrid
end

ExtendableGrids.ExtendableGrid{Float64, Int32}
      dim =       2
   nnodes =     972
   ncells =    1760
  nbfaces =     182
   nedges =    2731

let
    if isdefined(Main,:PlutoRunner)
        vis = GridVisualizer(layout = (1, 2), size = (650, 300), linewidth = 0.1)
        gridplot!(vis[1, 1], flowgrid, title = "flowgrid", aspect = 0.5)
        gridplot!(vis[1, 2], pnpgrid, title = "pnpgrid", aspect = 0.5)
        reveal(vis)
    end
end

begin
    const cbulk = 1.0 * mol / dm^3
    const Δϕ = 0.5 * V
    const σ = -15 * μA * s / cm^2
    const σ_embed = [σ]
    embed(data, λ) = σ_embed[1] = σ * λ
end

embed (generic function with 1 method)

function pnpbcond(y, u, bnode, data)
    (; iϕ, nc) = data
    boundary_neumann!(y, u, bnode, species = iϕ, region = 2, value = σ_embed[1])
    boundary_dirichlet!(y, u, bnode, species = iϕ, region = 1, value = -Δϕ / 2)
    boundary_dirichlet!(y, u, bnode, species = iϕ, region = 3, value = Δϕ / 2)
    for i in 1:nc
        boundary_dirichlet!(y, u, bnode, species = i, region = 1, value = cbulk)
        boundary_dirichlet!(y, u, bnode, species = i, region = 3, value = cbulk)
    end
    return
end

pnpbcond (generic function with 1 method)

function flowbcond(flowsolver)
    assign_operator!(
        flowsolver,
        HomogeneousBoundaryData(
            LiquidElectrolytes.velocity_unknown(flowsolver);
            regions = [4], mask = (1, 0, 1)
        )
    )
    return assign_operator!(
        flowsolver,
        HomogeneousBoundaryData(LiquidElectrolytes.velocity_unknown(flowsolver); regions = [2])
    )
end

flowbcond (generic function with 1 method)

default_pnpdata = ElectrolyteData();

function PNPData(; molar_volume = default_pnpdata.v0, scheme = :μex, κ = 20)
    pnpdata = ElectrolyteData()
    pnpdata.edgevelocity = zeros(num_edges(pnpgrid))
    pnpdata.solvepressure = false
    pnpdata.scheme = scheme
    pnpdata.κ .= κ
    pnpdata.D .= 1.0e-9 * m^2 / s
    pnpdata.pscale = 1
    pnpdata.v .= molar_volume
    return pnpdata
end

PNPData (generic function with 1 method)

gcdata = PNPData(molar_volume = 0, κ = 0)

ElectrolyteData(2, 0, 3, 4, [1.0e-9, 1.0e-9], [1, -1], 0.0180153, [0.0180153, 0.0180153], 1.805054151624549e-5, [0.0, 0.0], [0.0, 0.0], [100.0, 100.0], 0.0, 0.0, 2, 0.0, 1, 298.15, 2478.957029602388, 96485.33212331001, 78.49, 8.8541878128e-12, 1.0, false, :μex, false, [1.805054151624549e-5, 1.805054151624549e-5, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], :DGL)

dgldata = PNPData()

ElectrolyteData(2, 0, 3, 4, [1.0e-9, 1.0e-9], [1, -1], 0.0180153, [0.0180153, 0.0180153], 1.805054151624549e-5, [1.805054151624549e-5, 1.805054151624549e-5], [20.0, 20.0], [100.0, 100.0], 0.0, 0.0, 2, 0.0, 1, 298.15, 2478.957029602388, 96485.33212331001, 78.49, 8.8541878128e-12, 1.0, false, :μex, false, [1.805054151624549e-5, 1.805054151624549e-5, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], :DGL)

begin
    gcsolver = LiquidElectrolytes.PNPStokesSolver(;
        flowgrid,
        pnpgrid,
        μ,
        velospace = H1BR,
        pnpdata = gcdata,
        pnpbcond,
        pnpreaction = (y, u, bnode, data) -> nothing,
        flowbcond
    )
end;

begin
    dglsolver = LiquidElectrolytes.PNPStokesSolver(;
        flowgrid,
        pnpgrid,
        μ,
        velospace = H1BR,
        pnpdata = dgldata,
        pnpbcond,
        pnpreaction = (y, u, bnode, data) -> nothing,
        flowbcond
    )
end;

gcpnpsol, gcflowsol = solve(
    gcsolver;
    verbose = "n",
    niter = 20,
    damp0 = 0.1,
    damp_initial = 1,
    tol_round = 1.0e-9,
    max_round = 3,
    embed,
    nembed = 10
)

([999.9999999999998 999.9999999999998 … 999.9999999999998 999.9999999999998; 999.9999999999998 999.9999999999998 … 999.9999999999998 999.9999999999998; -0.25 -0.25 … 0.25 0.25; -13825.929476277674 -14492.056839464882 … -510270.55560461764 -1.736053810490319e6], 
FEVector information
====================
   block  |  ndofs 	|     min  /  max    	| FEType 		 (name/tag)
 [    1]  |    4675	| -1.14e-01/4.72e-01  	| H1BR{2}  	 (u (velocity))
 [    2]  |    1760	| -1.09e+07/8.74e+04  	| L2P0{1} (broken)  	 (p (pressure))
 [    3]  |     972	| -2.91e-01/2.50e-01  	| H1P1{1}  	 (ϕ (voltage))
 [    4]  |     972	| -1.73e+04/7.04e+08  	| H1P1{1}  	 (q (charge)))

dglpnpsol, dglflowsol = solve(
    dglsolver;
    verbose = "n",
    niter = 20,
    damp0 = 0.1,
    damp_initial = 1,
    tol_round = 1.0e-9,
    max_round = 3,
    embed,
    nembed = 10
)

([999.9999999999998 999.9999999999998 … 999.9999999999998 999.9999999999998; 999.9999999999998 999.9999999999998 … 999.9999999999998 999.9999999999998; -0.25 -0.25 … 0.25 0.25; -18847.096619302196 -19481.623308918366 … -1.2480954335217392e6 -2.3756804244562755e6], 
FEVector information
====================
   block  |  ndofs 	|     min  /  max    	| FEType 		 (name/tag)
 [    1]  |    4675	| -1.07e-01/8.46e-01  	| H1BR{2}  	 (u (velocity))
 [    2]  |    1760	| -1.42e+07/1.62e+05  	| L2P0{1} (broken)  	 (p (pressure))
 [    3]  |     972	| -3.29e-01/2.50e-01  	| H1P1{1}  	 (ϕ (voltage))
 [    4]  |     972	| -8.94e+01/2.19e+08  	| H1P1{1}  	 (q (charge)))

if isdefined(Main,:PlutoRunner)
    LiquidElectrolytes.flowplot(gcflowsol, gcsolver.flowsolver;
                                Plotter = CairoMakie,
                                vscale = 0.5,
                                aspect = 0.25,
                                rasterpoints = 30)
end

if isdefined(Main,:PlutoRunner)
    LiquidElectrolytes.flowplot(dglflowsol, dglsolver.flowsolver;
                                Plotter = CairoMakie,
                                vscale = 0.5,
                                aspect = 0.25,
                                rasterpoints = 30)
end

let
    if isdefined(Main,:PlutoRunner)
        vis = GridVisualizer(; layout = (2, 2), size = (650, 600))
        scalarplot!(vis[1, 1], xgrid, gcpnpsol[1, :] / cbulk)
        scalarplot!(vis[1, 2], xgrid, gcpnpsol[2, :] / cbulk)
        scalarplot!(vis[2, 1], xgrid, gcpnpsol[3, :], colormap = :bwr, limits = (-1, 1))
        scalarplot!(vis[2, 2], xgrid, gcpnpsol[4, :])
        reveal(vis)
    end
end

let
    if isdefined(Main,:PlutoRunner)
        vis = GridVisualizer(; layout = (2, 2), size = (650, 600))
        scalarplot!(vis[1, 1], xgrid, dglpnpsol[1, :] / cbulk)
        scalarplot!(vis[1, 2], xgrid, dglpnpsol[2, :] / cbulk)
        scalarplot!(vis[2, 1], xgrid, dglpnpsol[3, :], colormap = :bwr, limits = (-1, 1))
        scalarplot!(vis[2, 2], xgrid, dglpnpsol[4, :])
        reveal(vis)
    end
end

function midvelo(flowsol, pnpsol, solver)
    nodevelo = LiquidElectrolytes.node_velocity(flowsol, solver.flowsolver)
    nx = length(pnpX)
    iymid = (length(Y) - 1) ÷ 2 + 1
    midrange = ((iymid - 1) * nx + 1):(iymid * nx)
    c_p = pnpsol[1, midrange]
    c_m = pnpsol[2, midrange]
    ϕ = -pnpsol[3, midrange]
    p = pnpsol[4, midrange]
    u = nodevelo[2, midrange]
    return c_p, c_m, ϕ, p, u
end

midvelo (generic function with 1 method)

if isdefined(Main, :PlutoRunner)
    fig=CairoMakie.Figure(size=(650,300))
    c_p, c_m, ϕ, p, u = midvelo(gcflowsol,gcpnpsol, gcsolver)
    ax1=CairoMakie.Axis(fig[1,1],
        yaxisposition=:left,
        xticks=[0,1,2,3,4,5],
        ylabel="u/(m/s)",
        xlabel="x/nm",
        title="Zero ion size"
        )
    CairoMakie.ylims!(ax1,0,0.8)
    CairoMakie.xlims!(ax1,0,5)
    ax2=CairoMakie.Axis(fig[1,1],
        yaxisposition=:right,
        ylabel="c/(mol/L); ϕ/(10mV)"
        )
    CairoMakie.ylims!(ax2,0,10)
    CairoMakie.xlims!(ax2,0,5)

    data=[
lines!(ax1, pnpX / nm, u, color=:green,linewidth = 2)
 lines!(ax2, pnpX / nm, c_p / (mol / dm^3), color=:red, linewidth = 2)
 lines!(ax2, pnpX / nm, c_m / (mol / dm^3), color=:blue, linewidth = 2)
lines!(ax2, pnpX / nm, ϕ / (10mV), color=:black, linewidth = 2)
    ]
CairoMakie.axislegend(ax1, data, ["u_z","c^+","c^-","Φ"], position=:ct, backgroundcolor = :transparent)

    
        c_p, c_m, ϕ, p, u = midvelo(dglflowsol,dglpnpsol, dglsolver)
    ax1=CairoMakie.Axis(fig[1,2],
        yaxisposition=:left,
        xticks=[0,1,2,3,4,5],
        ylabel="u/(m/s)",
        xlabel="x/nm",
        title="Finite ion size\n Solvation"
        )
    CairoMakie.ylims!(ax1,0,0.8)
    CairoMakie.xlims!(ax1,0,5)
    ax2=CairoMakie.Axis(fig[1,2],
        yaxisposition=:right,
        ylabel="c/(mol/L); ϕ/(10mV)"
        )
    CairoMakie.ylims!(ax2,0,10)
    CairoMakie.xlims!(ax2,0,5)

    data=[
lines!(ax1, pnpX / nm, u, color=:green, linewidth = 2),
 lines!(ax2, pnpX / nm, c_p / (mol / dm^3), color=:red,linewidth = 2),
 lines!(ax2, pnpX / nm, c_m / (mol / dm^3), color=:blue, linewidth = 2),
lines!(ax2, pnpX / nm, ϕ / (10mV), color=:black, linewidth = 2)
    ]
CairoMakie.axislegend(ax1, data, ["u_z","c^+","c^-","Φ"], position=:ct, backgroundcolor = :transparent)
    fig
    
end

begin
    gc_c_p, gc_c_m, gc_ϕ, gc_p, gc_u = midvelo(gcflowsol, gcpnpsol, gcsolver)
    dgl_c_p, dgl_c_m, dgl_ϕ, dgl_p, dgl_u = midvelo(dglflowsol, dglpnpsol, dglsolver)
    @test gc_u[1] < dgl_u[1]
    @test gc_c_p[end] > dgl_c_p[end]
end

Test Passed