diff --git a/.gitignore b/.gitignore index 65e5528..f3b6f7d 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,14 @@ Manifest*.toml __pycache__ condaenv **/venv +draftresults +draft.pdf +*.nav +*.snm +*.toc +*.aux +*.log +*.bbl +*.blg +*.out +*.tdo diff --git a/Project.toml b/Project.toml index 626cb05..fe81d8b 100644 --- a/Project.toml +++ b/Project.toml @@ -11,11 +11,13 @@ DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821" ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" HypertextLiteral = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2" +InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" JuliaMPBSolver = "d8b18f01-5396-498d-b34d-247825c18ff0" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LessUnitful = "f29f6376-6e90-4d80-80c9-fb8ec61203d5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" Pluto = "c3e4b0f8-55cb-11ea-2926-15256bba5781" PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8" @@ -23,6 +25,7 @@ PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46" PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" PythonPlot = "274fc56d-3b97-40fa-a1cd-1b4a50311bf9" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" VoronoiFVM = "82b139dc-5afc-11e9-35da-9b9bdfd336f3" @@ -33,12 +36,15 @@ JuliaMPBSolver = {path = "packages/JuliaMPBSolver"} [compat] DoubleFloats = "1" DynamicQuantities = "1.10.0" +ExampleJuggler = "2.4.0" ExtendableGrids = "1" ForwardDiff = "1" HypertextLiteral = "0.9" +InteractiveUtils = "1.11.0" Interpolations = "0.16" LaTeXStrings = "1.0.1" LessUnitful = "1" +Markdown = "1.11.0" PlutoUI = "0.7" PreallocationTools = "0.4" PythonCall = "0.9" @@ -48,9 +54,10 @@ VoronoiFVM = "2.14" julia = "1.11" [extras] +ExampleJuggler = "3bbe58f8-ed81-4c4e-a134-03e85fcf4a1a" JuliaMPBSolver = "d8b18f01-5396-498d-b34d-247825c18ff0" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Pkg", "Test", "JuliaMPBSolver"] +test = ["Pkg", "Test", "JuliaMPBSolver", "ExampleJuggler"] diff --git a/notebooks/HalfCellAppliedPotential.jl b/notebooks/HalfCellAppliedPotential.jl new file mode 100644 index 0000000..744e562 --- /dev/null +++ b/notebooks/HalfCellAppliedPotential.jl @@ -0,0 +1,261 @@ +### A Pluto.jl notebook ### +# v0.20.19 + +using Markdown +using InteractiveUtils + +# ╔═╡ a70cef7d-2a2f-4155-bdf3-fec9df94c63f +begin + using Pkg + Pkg.activate(joinpath(@__DIR__, "..")) + using Revise + using PlutoUI, HypertextLiteral, UUIDs + using VoronoiFVM + using ExtendableGrids + using LinearAlgebra + using LessUnitful + using Test + using DoubleFloats + using PythonPlot + using PythonCall + using PythonPlot: pyplot + using LaTeXStrings + using Colors + using JuliaMPBSolver.ICMPBP: ICMPBP, ICMPBData, AppliedPotentialHalfCell, set_molarity!, calc_cmol, calc_c0mol, calc_χ, get_E, get_φ, get_p, get_c0, + set_κ!, set_φ! + using DrWatson, PoissonBoltzmannIPAM2025 + +end + +# ╔═╡ af6ae00d-f032-4743-878b-e575466b6e84 +md""" +# Half cell with applied potential +""" + +# ╔═╡ 87d66a6f-b40d-4d76-afc6-4b4f086e80a4 +begin + const nm = ufac"nm" + const V = ufac"V" + const cm = ufac"cm" + const dm = ufac"dm" + const μF = ufac"μF" + const N_A = ph"N_A" + const mol = ufac"mol" +end + +# ╔═╡ 087614cc-a4f6-4867-a6bd-201d1f2a7fc2 +begin + const L = 10nm + const n = 101 +end + +# ╔═╡ c7a51cb9-790a-4818-95d7-c32f2252e8b3 +Xhalf = (10 .^ range(-3, 1, length = 101)) * nm; + +# ╔═╡ 4ef3d13c-2f57-4575-bc88-c8092ae6910f +gridhalf = simplexgrid(Xhalf) + +# ╔═╡ 3a99a9bb-7862-41f8-a535-d3c580da6909 +md""" +All values are given with respect to SI basic units (m, kg, s, V, A) +""" + +# ╔═╡ b24b7e23-61ea-41fc-a345-286e904c042b +datavhalf = ICMPBData(χvar = true) + +# ╔═╡ 1bb47749-edde-4bee-be9f-059a7652b354 +begin + datavhalf.qscale = sqrt(datavhalf.qscale) + halfcell = AppliedPotentialHalfCell(gridhalf, datavhalf, dielectric_decrement = false, valuetype = Float64) + + halfcelldd = AppliedPotentialHalfCell(gridhalf, datavhalf, dielectric_decrement = true, valuetype = Float64) +end; + + +# ╔═╡ dd3c4807-3972-4e9c-a44a-3347b065d01c +p3 = plotcells(halfcell, halfcelldd); p3 + +# ╔═╡ 935f8897-4d68-447b-8316-1a0fe0285d54 +@test isa(p3, Figure) + +# ╔═╡ d07ac411-7985-4b5f-a88b-8aa4037b7d65 +function makecolors(V) + h = 0.5 / length(V) + return colors = [ (i * h, 0, 1 - i * h) for i in 1:length(V) ] +end + +# ╔═╡ 2cf54b71-d99b-40bd-b9a6-0a7cf919614b +p1 = let + molarities = [0.001, 0.01, 0.1, 1] + φ_max = 0.5 + colors = makecolors(molarities) + pyplot.close() + clf() + fig, ax = pyplot.subplots(1, 1) + fig.set_size_inches(600 / 100, 300 / 100) + κ = 10 + ax.set_title("Double layer capacitance, κ=$(κ)\n Dashed: witghout dielectric decrement") + set_κ!(halfcell, κ) + set_κ!(halfcelldd, κ) + for i in 1:length(molarities) + M = molarities[i] + set_molarity!(halfcell, M) + + volts, dlcaps = ICMPBP.dlcapsweep(halfcell; φ_max) + plot(volts, dlcaps / (μF / cm^2), color = colors[i], linestyle = "--") + + set_molarity!(halfcelldd, M) + volts, dlcaps = ICMPBP.dlcapsweep(halfcelldd; φ_max, damp_initial = 0.5) + plot(volts, dlcaps / (μF / cm^2), label = "M=$(M)", color = colors[i]) + end + ax.set_xlabel("U/V") + ax.set_ylabel(L"C_{dl}/(μF/cm^2)") + ax.legend() + ax.grid() + tight_layout() + + !haskey(ENV, "CI") && savefig(draftresultsdir("halfcell_dlcap_M"), dpi = 600) + + gcf() +end; p1 + +# ╔═╡ 52d5fb3c-072d-4a6e-a392-7e236b4ec933 +@test isa(p1, Figure) + +# ╔═╡ cc8e8ba6-4cc6-4e8a-bbd6-07b9f61c2ef5 +p2 = let + kappas = [1, 5, 10, 20] + φ_max = 0.5 + M = 0.1 + colors = makecolors(kappas) + pyplot.close() + clf() + fig, ax = pyplot.subplots(1, 1) + fig.set_size_inches(600 / 100, 300 / 100) + ax.set_title("Double layer capacitance, M=$(M)\n Dashed: without dielectric decrement") + set_molarity!(halfcell, M) + set_molarity!(halfcelldd, M) + for i in 1:length(kappas) + κ = kappas[i] + set_κ!(halfcell, κ) + set_κ!(halfcelldd, κ) + + volts, dlcaps = ICMPBP.dlcapsweep(halfcelldd; φ_max, damp_initial = 0.5) + plot(volts, dlcaps / (μF / cm^2), label = "κ=$(κ)", color = colors[i]) + + volts, dlcaps = ICMPBP.dlcapsweep(halfcell; φ_max) + plot(volts, dlcaps / (μF / cm^2), color = colors[i], linestyle = "--") + + end + ax.set_xlabel("U/V") + ax.set_ylabel(L"C_{dl}/(μF/cm^2)") + ax.legend() + ax.grid() + tight_layout() + !haskey(ENV, "CI") && savefig(draftresultsdir("halfcell_dlcap_kappa"), dpi = 600) + + gcf() +end; p2 + +# ╔═╡ 48896662-d87f-4606-9118-6471184b4dc7 +@test isa(p2, Figure) + +# ╔═╡ 8af12f1c-d35b-4cc9-8185-1bb5adbb69e8 +html"""
""" + +# ╔═╡ 784b4c3e-bb2a-4940-a83a-ed5e5898dfd4 +html"""""" + +# ╔═╡ afe4745f-f9f1-4e23-8735-cbec6fb79c41 +begin + function floataside(text::Markdown.MD; top = 1) + uuid = uuid1() + return @htl( + """ + + + + + """ + ) + end + floataside(stuff; kwargs...) = floataside(md"""$(stuff)"""; kwargs...) +end; + + +# ╔═╡ b8fd36a7-d8d1-45f7-b66e-df9132168bfc +# https://discourse.julialang.org/t/adding-a-restart-process-button-in-pluto/76812/5 +restart_button() = html""" + +"""; + +# ╔═╡ Cell order: +# ╟─af6ae00d-f032-4743-878b-e575466b6e84 +# ╠═a70cef7d-2a2f-4155-bdf3-fec9df94c63f +# ╠═87d66a6f-b40d-4d76-afc6-4b4f086e80a4 +# ╠═087614cc-a4f6-4867-a6bd-201d1f2a7fc2 +# ╠═c7a51cb9-790a-4818-95d7-c32f2252e8b3 +# ╠═4ef3d13c-2f57-4575-bc88-c8092ae6910f +# ╟─3a99a9bb-7862-41f8-a535-d3c580da6909 +# ╠═b24b7e23-61ea-41fc-a345-286e904c042b +# ╠═1bb47749-edde-4bee-be9f-059a7652b354 +# ╠═2cf54b71-d99b-40bd-b9a6-0a7cf919614b +# ╠═52d5fb3c-072d-4a6e-a392-7e236b4ec933 +# ╠═cc8e8ba6-4cc6-4e8a-bbd6-07b9f61c2ef5 +# ╠═48896662-d87f-4606-9118-6471184b4dc7 +# ╠═dd3c4807-3972-4e9c-a44a-3347b065d01c +# ╠═935f8897-4d68-447b-8316-1a0fe0285d54 +# ╠═d07ac411-7985-4b5f-a88b-8aa4037b7d65 +# ╟─8af12f1c-d35b-4cc9-8185-1bb5adbb69e8 +# ╟─784b4c3e-bb2a-4940-a83a-ed5e5898dfd4 +# ╟─afe4745f-f9f1-4e23-8735-cbec6fb79c41 +# ╟─b8fd36a7-d8d1-45f7-b66e-df9132168bfc diff --git a/notebooks/ICMPBP-DD-Draft.jl b/notebooks/ICMPBP-DD-Draft.jl index d2c1c9c..f7be684 100644 --- a/notebooks/ICMPBP-DD-Draft.jl +++ b/notebooks/ICMPBP-DD-Draft.jl @@ -28,8 +28,7 @@ begin using LessUnitful using Test using PythonPlot - using Colors - using JuliaMPBSolver.ICMPBP: ICMPBData, ICMPBSystem, L_Debye, set_molarity!, dlcap0, DerivedData, apply_charge!, ysum, qsweep, capscalc, calc_cmol, calc_c0mol, calc_χ + using JuliaMPBSolver.ICMPBP: ICMPBData, ICMPBSystem, L_Debye, set_molarity!, dlcap0, DerivedData, apply_charge!, ysum, qsweep, capscalc, calc_cmol, calc_c0mol, calc_χ, W end # ╔═╡ ef660f6f-9de3-4896-a65e-13c60df5de1e @@ -145,20 +144,6 @@ md""" ## Simulation result """ -# ╔═╡ 7caea4c4-9783-403c-931c-681f86166a25 -md""" -## Next steps -""" - -# ╔═╡ 65c0b7fa-de4f-454e-ac63-2723d5acc26c -md""" -- Finalize pyiron integration -- Discuss comparison with molecular simulation results -- Add dielectric decrement model -- Discuss additional molecular interactio terms -- Write paper -""" - # ╔═╡ 760e5861-7a6f-41bb-8aec-5e7466c6ec9f md""" ## Simulation setup and run @@ -173,7 +158,7 @@ begin const μF = ufac"μF" const N_A = ph"N_A" const mol = ufac"mol" -end +end; # ╔═╡ ae11bded-9f67-4004-8786-ed54e1ccb932 surfcharge(n) = n * ph"e" / ufac"nm^2" @@ -261,7 +246,7 @@ end; # ╔═╡ 2ef6b8d7-e5f3-4700-8a40-8feffab3569f floataside( md""" - - ``L/nm``: $(@bind L0 PlutoUI.Slider(2:2:10, default=10, show_value=true)) + - ``L/nm``: $(@bind L0 PlutoUI.Slider(2:1:10, default=10, show_value=true)) - ``M_{avg}/(mol/dm^3)``: $(@bind M1_avg PlutoUI.Slider(0.1:0.1:2, default=1, show_value=true)) - ``n_e/(e/nm^2)``: $(@bind n1_e PlutoUI.Slider(0:0.5:5, default=1, show_value=true)) - ``κ``: $(@bind kappa1 PlutoUI.Slider(0:1:20, default=10, show_value=true)) @@ -270,6 +255,24 @@ floataside( ) +# ╔═╡ eacdd772-1869-406a-b601-64cdd6453ec1 +begin + data1 = ICMPBData( + κ = [kappa1, kappa1] + ) + set_molarity!(data1, M1_avg) + data1.conserveions = true + data1.χvar = true +end + +# ╔═╡ f8c1c2bd-7466-491e-9132-4f15edcfa4c7 +let + clf() + X = -10.0e8:1.0e6:10.0e8 + plot(X, W.(data1.δ0 * X / data1.kT)) + gcf() +end + # ╔═╡ a629e8a1-b1d7-42d8-8c17-43475785218e begin L = L0 * nm @@ -280,6 +283,12 @@ begin i3 = grid[BFaceNodes][3][1] end +# ╔═╡ 8433319f-2f78-494c-9b2e-a5390cf93b00 +sys1 = ICMPBSystem(grid, data1, valuetype = Float64); + +# ╔═╡ 70910bd5-b8ca-4021-9b40-233b50ea5601 +inival1 = unknowns(sys1, data1); + # ╔═╡ f579cf2d-9511-48a8-bf11-7400ef76ee3d function plotsol( sol, sys; data = data(sys), @@ -314,9 +323,11 @@ function plotsol( nv = nodevolumes(sys) cm, cp = c[1, :] ⋅ nv / L, c[2, :] ⋅ nv / L - M_bulk = sol[6:7, i3] * data.cscale / (ph"N_A" * ufac"mol/dm^3") - crm, crp = M_bulk[1], M_bulk[2] - ax2.set_title("M_bulk=$(myround.((crm, crp))), M_avg=$(myround.((cm, cp)))") + if data1.conserveions + M_bulk = sol[7:8, i3] * data.cscale / (ph"N_A" * ufac"mol/dm^3") + crm, crp = M_bulk[1], M_bulk[2] + # ax2.set_title("M_bulk=$(myround.((crm, crp))), M_avg=$(myround.((cm, cp)))") + end # ax2.set_title("M_avg=$(myround.((cm, cp)))") ax2.set_xlabel("z/nm") ax2.set_ylabel("c/(mol/L)") @@ -334,31 +345,22 @@ function plotsol( ax3r.plot(X / nm, calc_χ(sol, sys), color = "orange", label = "χ") ax3.legend(loc = (0.1, 0.1)) ax3r.legend(loc = (0.8, 0.1)) + ax3r.set_ylim(0, 100) tight_layout() return PythonPlot.gcf() end -# ╔═╡ eacdd772-1869-406a-b601-64cdd6453ec1 +# ╔═╡ d2df6ed0-e6f5-4677-b790-bfc40de7fd6a begin - data1 = ICMPBData( - q = [surfcharge(n1_e), -surfcharge(n1_e)], - κ = [kappa1, kappa1] - ) - set_molarity!(data1, M1_avg) - data1.conserveions = true - data1.χvar = false + Q = surfcharge(n1_e) + sol1 = inival1 + for q in range(0, Q, length = 11) + data1.q .= [-q, q] + global sol1 = solve(sys1; inival = sol1, verbose = "") + end end -# ╔═╡ 8433319f-2f78-494c-9b2e-a5390cf93b00 -sys1 = ICMPBSystem(grid, data1); - -# ╔═╡ 70910bd5-b8ca-4021-9b40-233b50ea5601 -inival1 = unknowns(sys1, data1); - -# ╔═╡ d2df6ed0-e6f5-4677-b790-bfc40de7fd6a -sol1 = solve(sys1; inival = inival1, damp_initial = 0.1, verbose = "n", log = true) - # ╔═╡ 5f153fe4-4476-401e-8674-b6902213c19e md""" Newton steps: $(length(history(sol1))) @@ -367,6 +369,12 @@ Newton steps: $(length(history(sol1))) # ╔═╡ c7a08779-53f3-4fab-8bd4-3dffe6135c3b plotsol(sol1, sys1; Mscale) +# ╔═╡ 0e4ec7f0-0aa8-4a32-96a3-40f63f32a12d +sol1 + +# ╔═╡ 7d7ebb45-2fb3-40ea-83f2-62d0a240b2db +@test isa(sol1, AbstractMatrix) + # ╔═╡ Cell order: # ╠═60941eaa-1aea-11eb-1277-97b991548781 # ╟─ef660f6f-9de3-4896-a65e-13c60df5de1e @@ -380,20 +388,21 @@ plotsol(sol1, sys1; Mscale) # ╟─bab7b779-339d-4702-a88a-e6cbf5a72cd0 # ╟─16a79117-e89b-4478-a571-8c011b5784c1 # ╟─5f153fe4-4476-401e-8674-b6902213c19e -# ╟─c7a08779-53f3-4fab-8bd4-3dffe6135c3b -# ╟─7caea4c4-9783-403c-931c-681f86166a25 -# ╟─65c0b7fa-de4f-454e-ac63-2723d5acc26c +# ╠═c7a08779-53f3-4fab-8bd4-3dffe6135c3b +# ╠═0e4ec7f0-0aa8-4a32-96a3-40f63f32a12d +# ╠═eacdd772-1869-406a-b601-64cdd6453ec1 # ╟─760e5861-7a6f-41bb-8aec-5e7466c6ec9f # ╠═f4facb34-1f4a-432d-8a1e-30299e542bcd -# ╠═2ef6b8d7-e5f3-4700-8a40-8feffab3569f +# ╟─2ef6b8d7-e5f3-4700-8a40-8feffab3569f # ╠═a629e8a1-b1d7-42d8-8c17-43475785218e # ╠═ae11bded-9f67-4004-8786-ed54e1ccb932 -# ╠═eacdd772-1869-406a-b601-64cdd6453ec1 # ╠═8433319f-2f78-494c-9b2e-a5390cf93b00 # ╠═70910bd5-b8ca-4021-9b40-233b50ea5601 # ╠═d2df6ed0-e6f5-4677-b790-bfc40de7fd6a +# ╠═7d7ebb45-2fb3-40ea-83f2-62d0a240b2db +# ╠═f8c1c2bd-7466-491e-9132-4f15edcfa4c7 # ╟─f75f1d3a-47e5-475b-97b1-bb275a510783 # ╠═dc05f31c-a28e-4470-8916-72dda567b149 # ╠═f579cf2d-9511-48a8-bf11-7400ef76ee3d # ╠═000e87b5-cd1b-4b23-99be-6b7006502312 -# ╠═0a824568-997a-4a48-9fe3-71cc5f9b7471 +# ╟─0a824568-997a-4a48-9fe3-71cc5f9b7471 diff --git a/notebooks/ICMPBP-Draft.jl b/notebooks/ICMPBP-Draft.jl index b9bfcb5..74e4e02 100644 --- a/notebooks/ICMPBP-Draft.jl +++ b/notebooks/ICMPBP-Draft.jl @@ -1,5 +1,5 @@ ### A Pluto.jl notebook ### -# v0.20.21 +# v0.20.19 using Markdown using InteractiveUtils @@ -140,21 +140,9 @@ begin sol0 = solve!(state; inival, verbose = "n", damp_initial = 0.1) end -# ╔═╡ d6f1acbc-d2b0-4d26-b9f6-ad46b9ddfb05 -sol0[:, i3] - # ╔═╡ 1394cfd7-aecc-4559-a460-98082fd43a9c state.matrix -# ╔═╡ a4c2f75e-6f20-4e28-8252-90d58fbf1b24 -data.conserveions - -# ╔═╡ efb12e12-825b-4dfd-aa10-c6afb304b6bf -ph"e" / ufac"nm^2" - -# ╔═╡ 9108a7e6-176e-481b-8ca6-3c8445051e1c -data.n_avg / ph"N_A" - # ╔═╡ 6f037b32-e2a8-4693-b46c-952d6b140e8e begin data1 = apply_charge!(deepcopy(data), 1 * ph"e" / ufac"nm^2") @@ -163,6 +151,14 @@ begin sol1[:, i3] end +# ╔═╡ 388c4956-211c-40c5-b5ff-6e2059e47a38 +begin + Y = ysum(sys, sol1) +end + +# ╔═╡ 53b7a161-87d5-4702-87ad-82e691571c2c +@test all(x -> ≈(x, 1.0, atol = 2.0e-2), Y) + # ╔═╡ 1c0145d5-76b1-48c1-8852-de1a2668285a molarities = [0.1, 1] @@ -172,11 +168,6 @@ qsweep(sys, verbose = "n", qmax = 2) # ╔═╡ a7f2692e-a15f-47b7-8486-8948ce7ab3f7 result_pp = capscalc(sys, molarities; qmax = 1) -# ╔═╡ e114ec0d-13d3-4455-b1c9-d1c5d76671d9 -md""" -#### Pressure poisson problem -""" - # ╔═╡ 0b6f33b9-41d4-48fd-8026-8a3bddcc1989 md""" #### Result plot @@ -199,7 +190,7 @@ function capsplot(ax, result, title) c = (1 - imol * hmol, 0, imol * hmol) ax.plot( - result[imol].voltages, result[imol].dlcaps / (μF / cm^2), + result[imol].voltages, -result[imol].dlcaps / (μF / cm^2), color = c, label = "$(result[imol].molarity)M" ) ax.scatter( @@ -216,14 +207,14 @@ let res = (600, 200) fig = figure(1, figsize = (6, 3)) ax1 = fig.add_subplot(1, 1, 1) - capsplot(ax1, result_pp, "Pressure poisson") + capsplot(ax1, result_pp, "Double layer capacitance") tight_layout() matplotlib.pyplot.close() fig end # ╔═╡ Cell order: -# ╠═ef660f6f-9de3-4896-a65e-13c60df5de1e +# ╟─ef660f6f-9de3-4896-a65e-13c60df5de1e # ╠═60941eaa-1aea-11eb-1277-97b991548781 # ╟─4082c3d3-b728-4bcc-b480-cdee41d9ab99 # ╟─920b7d84-56c6-4958-aed9-fc67ba0c43f6 @@ -239,16 +230,13 @@ end # ╠═31a1f686-f0b6-430a-83af-187df411b293 # ╠═684aa24b-046f-426f-9b99-f0c45c70f654 # ╠═14ac1c80-cae5-42f1-b0d3-33aa5bba4de6 -# ╠═d6f1acbc-d2b0-4d26-b9f6-ad46b9ddfb05 # ╠═1394cfd7-aecc-4559-a460-98082fd43a9c -# ╠═a4c2f75e-6f20-4e28-8252-90d58fbf1b24 -# ╠═efb12e12-825b-4dfd-aa10-c6afb304b6bf -# ╠═9108a7e6-176e-481b-8ca6-3c8445051e1c # ╠═6f037b32-e2a8-4693-b46c-952d6b140e8e +# ╠═388c4956-211c-40c5-b5ff-6e2059e47a38 +# ╠═53b7a161-87d5-4702-87ad-82e691571c2c # ╠═1c0145d5-76b1-48c1-8852-de1a2668285a # ╠═f1c33101-00e6-4af9-9e68-6cdf5fe92b59 # ╠═a7f2692e-a15f-47b7-8486-8948ce7ab3f7 -# ╟─e114ec0d-13d3-4455-b1c9-d1c5d76671d9 # ╟─0b6f33b9-41d4-48fd-8026-8a3bddcc1989 # ╟─1859db0c-1c1e-4a78-9d4a-e2ec45c3cffc # ╠═000e87b5-cd1b-4b23-99be-6b7006502312 diff --git a/notebooks/ICMPBP-EndOfHackathon.jl b/notebooks/ICMPBP-EndOfHackathon.jl index 313ccf8..0c0d562 100644 --- a/notebooks/ICMPBP-EndOfHackathon.jl +++ b/notebooks/ICMPBP-EndOfHackathon.jl @@ -1,5 +1,5 @@ ### A Pluto.jl notebook ### -# v0.20.21 +# v0.20.19 using Markdown using InteractiveUtils @@ -358,6 +358,9 @@ Newton steps: $(length(history(sol1))) # ╔═╡ c7a08779-53f3-4fab-8bd4-3dffe6135c3b plotsol(sol1, sys1; Mscale) +# ╔═╡ 7b34987e-93cc-4882-ab76-44cf761fbc78 +@test isa(sol1, AbstractMatrix) + # ╔═╡ Cell order: # ╠═60941eaa-1aea-11eb-1277-97b991548781 # ╟─ef660f6f-9de3-4896-a65e-13c60df5de1e @@ -383,6 +386,7 @@ plotsol(sol1, sys1; Mscale) # ╠═8433319f-2f78-494c-9b2e-a5390cf93b00 # ╠═70910bd5-b8ca-4021-9b40-233b50ea5601 # ╠═d2df6ed0-e6f5-4677-b790-bfc40de7fd6a +# ╠═7b34987e-93cc-4882-ab76-44cf761fbc78 # ╟─f75f1d3a-47e5-475b-97b1-bb275a510783 # ╠═dc05f31c-a28e-4470-8916-72dda567b149 # ╠═f579cf2d-9511-48a8-bf11-7400ef76ee3d diff --git a/notebooks/MPBP-Draft.jl b/notebooks/MPBP-Draft.jl index 4655b19..d25fdf9 100644 --- a/notebooks/MPBP-Draft.jl +++ b/notebooks/MPBP-Draft.jl @@ -1,5 +1,5 @@ ### A Pluto.jl notebook ### -# v0.20.21 +# v0.20.19 using Markdown using InteractiveUtils @@ -143,11 +143,11 @@ inival = unknowns(sys, data) # ╔═╡ 14ac1c80-cae5-42f1-b0d3-33aa5bba4de6 begin sol0 = solve(sys; inival, verbose = "n", damp_initial = 0.1) - ysum(sys, sol0) + Y = ysum(sys, sol0) end -# ╔═╡ efb12e12-825b-4dfd-aa10-c6afb304b6bf -ph"e" / ufac"nm^2" +# ╔═╡ f9e954e5-d45a-4229-a2f8-e3f407d49faf +@test all(x -> ≈(x, 1.0, atol = 5.0e-3), Y) # ╔═╡ 6f037b32-e2a8-4693-b46c-952d6b140e8e begin @@ -159,12 +159,12 @@ end # ╔═╡ 1c0145d5-76b1-48c1-8852-de1a2668285a molarities = [0.001, 0.01, 0.1, 1] -# ╔═╡ f1c33101-00e6-4af9-9e68-6cdf5fe92b59 -qsweep(sys) - # ╔═╡ a7f2692e-a15f-47b7-8486-8948ce7ab3f7 result_pp = capscalc(sys, molarities) +# ╔═╡ b155bc1f-e645-4b75-aad7-4b6bb0dc3d27 +@test isa(result_pp, Vector) + # ╔═╡ e114ec0d-13d3-4455-b1c9-d1c5d76671d9 md""" #### Pressure poisson problem @@ -192,7 +192,7 @@ function capsplot(ax, result, title) c = (1 - imol * hmol, 0, imol * hmol) ax.plot( - result[imol].voltages, result[imol].dlcaps / (μF / cm^2), + result[imol].voltages, -result[imol].dlcaps / (μF / cm^2), color = c, label = "$(result[imol].molarity)M" ) ax.scatter( @@ -209,7 +209,7 @@ let res = (600, 200) fig = figure(1, figsize = (6, 3)) ax1 = fig.add_subplot(1, 1, 1) - capsplot(ax1, result_pp, "Pressure poisson") + capsplot(ax1, result_pp, "Double Layer capacitance") tight_layout() matplotlib.pyplot.close() fig @@ -234,11 +234,11 @@ end # ╠═31a1f686-f0b6-430a-83af-187df411b293 # ╠═9b16c019-f2ce-4b42-97e1-6d13a463a232 # ╠═14ac1c80-cae5-42f1-b0d3-33aa5bba4de6 -# ╠═efb12e12-825b-4dfd-aa10-c6afb304b6bf +# ╠═f9e954e5-d45a-4229-a2f8-e3f407d49faf # ╠═6f037b32-e2a8-4693-b46c-952d6b140e8e # ╠═1c0145d5-76b1-48c1-8852-de1a2668285a -# ╠═f1c33101-00e6-4af9-9e68-6cdf5fe92b59 # ╠═a7f2692e-a15f-47b7-8486-8948ce7ab3f7 +# ╟─b155bc1f-e645-4b75-aad7-4b6bb0dc3d27 # ╟─e114ec0d-13d3-4455-b1c9-d1c5d76671d9 # ╟─0b6f33b9-41d4-48fd-8026-8a3bddcc1989 # ╟─1859db0c-1c1e-4a78-9d4a-e2ec45c3cffc diff --git a/notebooks/SymmetricCellSurfaceCharge.jl b/notebooks/SymmetricCellSurfaceCharge.jl new file mode 100644 index 0000000..c6a1410 --- /dev/null +++ b/notebooks/SymmetricCellSurfaceCharge.jl @@ -0,0 +1,199 @@ +### A Pluto.jl notebook ### +# v0.20.19 + +using Markdown +using InteractiveUtils + +# ╔═╡ a70cef7d-2a2f-4155-bdf3-fec9df94c63f +begin + using Pkg + Pkg.activate(joinpath(@__DIR__, "..")) + using Revise + using PlutoUI, HypertextLiteral, UUIDs + using VoronoiFVM + using ExtendableGrids + using LinearAlgebra + using LessUnitful + using Test + using DoubleFloats + using PythonPlot + using PythonCall + using PythonPlot: pyplot + using LaTeXStrings + using Colors + using DrWatson, PoissonBoltzmannIPAM2025 + using JuliaMPBSolver.ICMPBP: ICMPBP, ICMPBData, SurfaceChargedSymmetricCell, AbstractSymmetricCell, set_molarity!, calc_cmol, calc_c0mol, calc_χ, get_E, get_φ, get_p, get_c0, + set_κ!, set_q! +end + +# ╔═╡ af6ae00d-f032-4743-878b-e575466b6e84 +md""" +# Symmetric cell with ion conservation and applied charge +""" + +# ╔═╡ 87d66a6f-b40d-4d76-afc6-4b4f086e80a4 +begin + const nm = ufac"nm" + const V = ufac"V" + const cm = ufac"cm" + const dm = ufac"dm" + const μF = ufac"μF" + const N_A = ph"N_A" + const mol = ufac"mol" +end + +# ╔═╡ 087614cc-a4f6-4867-a6bd-201d1f2a7fc2 +begin + const Ll = 10nm + const Ls = 2nm + const n = 101 +end + +# ╔═╡ c7a51cb9-790a-4818-95d7-c32f2252e8b3 +Xl = range(0, Ll, length = n) + +# ╔═╡ d87781e8-af89-41d0-a56d-94744406042d +Xs = range(0, Ls, length = n) + +# ╔═╡ 4ef3d13c-2f57-4575-bc88-c8092ae6910f +begin + gridl = simplexgrid(Xl) + bfacemask!(gridl, [Ll / 2], [Ll / 2], 3, tol = 1.0e-10 * nm) + +end + +# ╔═╡ 4d78d81c-f1a3-4386-adce-2a7fa1e3eaef +begin + grids = simplexgrid(Xs) + bfacemask!(grids, [Ls / 2], [Ls / 2], 3, tol = 1.0e-10 * nm) + +end + +# ╔═╡ 3a99a9bb-7862-41f8-a535-d3c580da6909 +md""" +All values are given with respect to SI basic units (m, kg, s, V, A) +""" + +# ╔═╡ b24b7e23-61ea-41fc-a345-286e904c042b +data = ICMPBData(χvar = true, conserveions = true) + +# ╔═╡ 1bb47749-edde-4bee-be9f-059a7652b354 +begin + smallcell = SurfaceChargedSymmetricCell(grids, data, dielectric_decrement = false) + largecell = SurfaceChargedSymmetricCell(gridl, data, dielectric_decrement = false) + + smallcelldd = SurfaceChargedSymmetricCell(grids, data, dielectric_decrement = true) + largecelldd = SurfaceChargedSymmetricCell(gridl, data, dielectric_decrement = true) + +end; + +# ╔═╡ 6aae7cca-01f1-45be-9d4b-a3942bd042f3 +p = plotcells(largecell, largecelldd; figname = "largecell");p + +# ╔═╡ 27c9e166-df00-4a30-a6ba-32e9a071136b +@test isa(p, Figure) + +# ╔═╡ 773e6c0a-08ea-47b2-8924-42759f944c6b +p2 = plotcells(smallcell, smallcelldd, figname = "smallcell");p2 + +# ╔═╡ ceaff98b-1ba6-43b3-bd2f-0e7fe585f428 +@test isa(p2, Figure) + +# ╔═╡ 8af12f1c-d35b-4cc9-8185-1bb5adbb69e8 +html"""
""" + +# ╔═╡ 784b4c3e-bb2a-4940-a83a-ed5e5898dfd4 +html"""""" + +# ╔═╡ afe4745f-f9f1-4e23-8735-cbec6fb79c41 +begin + function floataside(text::Markdown.MD; top = 1) + uuid = uuid1() + return @htl( + """ + + + + + """ + ) + end + floataside(stuff; kwargs...) = floataside(md"""$(stuff)"""; kwargs...) +end; + + +# ╔═╡ b8fd36a7-d8d1-45f7-b66e-df9132168bfc +# https://discourse.julialang.org/t/adding-a-restart-process-button-in-pluto/76812/5 +restart_button() = html""" + +"""; + +# ╔═╡ Cell order: +# ╟─af6ae00d-f032-4743-878b-e575466b6e84 +# ╠═a70cef7d-2a2f-4155-bdf3-fec9df94c63f +# ╠═87d66a6f-b40d-4d76-afc6-4b4f086e80a4 +# ╠═087614cc-a4f6-4867-a6bd-201d1f2a7fc2 +# ╠═c7a51cb9-790a-4818-95d7-c32f2252e8b3 +# ╠═d87781e8-af89-41d0-a56d-94744406042d +# ╠═4ef3d13c-2f57-4575-bc88-c8092ae6910f +# ╠═4d78d81c-f1a3-4386-adce-2a7fa1e3eaef +# ╟─3a99a9bb-7862-41f8-a535-d3c580da6909 +# ╠═b24b7e23-61ea-41fc-a345-286e904c042b +# ╠═1bb47749-edde-4bee-be9f-059a7652b354 +# ╠═6aae7cca-01f1-45be-9d4b-a3942bd042f3 +# ╠═27c9e166-df00-4a30-a6ba-32e9a071136b +# ╠═773e6c0a-08ea-47b2-8924-42759f944c6b +# ╠═ceaff98b-1ba6-43b3-bd2f-0e7fe585f428 +# ╟─8af12f1c-d35b-4cc9-8185-1bb5adbb69e8 +# ╟─784b4c3e-bb2a-4940-a83a-ed5e5898dfd4 +# ╟─afe4745f-f9f1-4e23-8735-cbec6fb79c41 +# ╟─b8fd36a7-d8d1-45f7-b66e-df9132168bfc diff --git a/packages/JuliaMPBSolver/Project.toml b/packages/JuliaMPBSolver/Project.toml index 1dcb87d..4d9a9cd 100644 --- a/packages/JuliaMPBSolver/Project.toml +++ b/packages/JuliaMPBSolver/Project.toml @@ -1,6 +1,6 @@ name = "JuliaMPBSolver" uuid = "d8b18f01-5396-498d-b34d-247825c18ff0" -version = "0.3.1" +version = "0.3.2" authors = ["Jürgen Fuhrmann "] [deps] @@ -11,6 +11,7 @@ LessUnitful = "f29f6376-6e90-4d80-80c9-fb8ec61203d5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46" +SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" VoronoiFVM = "82b139dc-5afc-11e9-35da-9b9bdfd336f3" @@ -21,6 +22,7 @@ LessUnitful = "1.2.1" LinearAlgebra = "1.10.0" Markdown = "1.10.0" PreallocationTools = "0.4.34" +SciMLBase = "2.127.0" Unitful = "1.25.1" VoronoiFVM = "2.13.0" julia = "1.10" diff --git a/packages/JuliaMPBSolver/src/JuliaMPBSolver.jl b/packages/JuliaMPBSolver/src/JuliaMPBSolver.jl index f4abe43..8658b7f 100644 --- a/packages/JuliaMPBSolver/src/JuliaMPBSolver.jl +++ b/packages/JuliaMPBSolver/src/JuliaMPBSolver.jl @@ -24,7 +24,9 @@ module ICMPBP using ExtendableGrids using VoronoiFVM using LinearAlgebra + using SciMLBase include("icmbp-p.jl") + include("cells.jl") end include("api.jl") diff --git a/packages/JuliaMPBSolver/src/cells.jl b/packages/JuliaMPBSolver/src/cells.jl new file mode 100644 index 0000000..0983211 --- /dev/null +++ b/packages/JuliaMPBSolver/src/cells.jl @@ -0,0 +1,195 @@ +abstract type AbstractMPBCell end + +abstract type AbstractHalfCell <: AbstractMPBCell end +abstract type AbstractSymmetricCell <: AbstractMPBCell end + +struct AppliedPotentialHalfCell <: AbstractHalfCell + sys::VoronoiFVM.System +end + +struct AppliedPotentialSymmetricCell <: AbstractMPBCell + sys::VoronoiFVM.System +end + +struct SurfaceChargedHalfCell <: AbstractHalfCell + sys::VoronoiFVM.System +end + +struct SurfaceChargedSymmetricCell <: AbstractSymmetricCell + sys::VoronoiFVM.System +end + +function VoronoiFVM.unknowns(cell::AbstractMPBCell) + sys = cell.sys + data = sys.physics.data + (; i0, iφ, ip, iE, coffset, N) = data + u = unknowns(sys; inival = 0) + i3 = 0 + if data.conserveions + i3 = sys.grid[BFaceNodes][3][1] + end + for α in 1:N + u[α, :] .= 0.1 + if data.conserveions + u[coffset + α, i3] = data.n_E[α] / data.cscale + end + end + + u[i0, :] .= 1 - N * 0.1 + return u +end + +mpbdata(cell::AbstractMPBCell) = cell.sys.physics.data + +calc_cmol(sol, cell::AbstractMPBCell) = calc_cmol(sol, cell.sys) +calc_c0mol(sol, cell::AbstractMPBCell) = calc_c0mol(sol, cell.sys) +calc_χ(sol, cell::AbstractMPBCell) = calc_χ(sol, cell.sys) +get_E(sol, cell::AbstractMPBCell) = sol[mpbdata(cell).iE, :] * mpbdata(cell).Escale +get_φ(sol, cell::AbstractMPBCell) = sol[mpbdata(cell).iφ, :] +get_p(sol, cell::AbstractMPBCell) = sol[mpbdata(cell).ip, :] * mpbdata(cell).pscale +get_c0(sol, cell::AbstractMPBCell) = sol[mpbdata(cell).i0, :] +set_κ!(cell::AbstractMPBCell, κ::Number) = mpbdata(cell).κ = [κ, κ] +set_molarity!(cell::AbstractMPBCell, M) = set_molarity!(mpbdata(cell), M) +set_φ!(cell::AbstractMPBCell, φ::Number) = mpbdata(cell).φ = φ +set_q!(cell::AbstractMPBCell, q::Number) = mpbdata(cell).q .= [-q, q] + +function SciMLBase.solve(cell::AbstractMPBCell; inival = unknowns(cell), verbose = "", damp_initial = 0.1, kwargs...) + sys = cell.sys + return solve(sys; inival, damp_initial, verbose, kwargs...) +end + + +function halfcell_applied_potential_bcondition!(y, u, bnode, data) + (; iφ, ip) = data + boundary_dirichlet!(y, u, bnode, species = iφ, region = 1, value = data.φ) + boundary_dirichlet!(y, u, bnode, species = iφ, region = 2, value = 0) + boundary_dirichlet!(y, u, bnode, species = ip, region = 2, value = 0) + return nothing +end + +function AppliedPotentialHalfCell(grid, data; dielectric_decrement = false, valuetype = Float64) + data = deepcopy(data) + data.nv = ones(num_nodes(grid)) # help to satisfy sparsity detector + data.conserveions = false + data.χvar = dielectric_decrement + + sys = VoronoiFVM.System( + grid; + data = data, + flux = poisson_and_p_flux!, + reaction = reaction!, + bcondition = halfcell_applied_potential_bcondition!, + generic = ionconservation!, + unknown_storage = :dense, + valuetype + ) + + # Enable species for all fields + for i in 1:data.N + enable_species!(sys, i, [1]) + end + + enable_species!(sys, data.i0, [1]) + enable_species!(sys, data.iφ, [1]) + enable_species!(sys, data.ip, [1]) + enable_species!(sys, data.iE, [1]) + data.nv = nodevolumes(sys) + return AppliedPotentialHalfCell(sys) +end + +function dlcapsweep( + cell::AppliedPotentialHalfCell; φ_max = 1.0, δφ = 1.0e-5, + steps = 51, hmin = 1.0e-5, damp_initial = 1, kwargs... + ) + set_φ!(cell, 0) + sol0 = solve(cell; damp_initial) + volts = zeros(0) + dlcaps = zeros(0) + hmax = φ_max / steps + for dir in [-1, 1] + sol = sol0 + h = 0.001 + φ = 0.0 + first = true + while φ < φ_max + try + if first + _φ = φ + else + _φ = min(φ + h, φ_max) + end + set_φ!(cell, dir * _φ) + sol = solve(cell; inival = sol, damp_initial, kwargs...) + Q = calc_spacecharge(cell.sys, sol) + + set_φ!(cell, dir * (_φ + δφ)) + sol = solve(cell; inival = sol, damp_initial, kwargs...) + Qδ = calc_spacecharge(cell.sys, sol) + + cdl = (Qδ - Q) / (dir * δφ) + push!(volts, dir * φ) + push!(dlcaps, cdl) + φ = _φ + h = min(h * 1.2, hmax) + first = false + catch e + h = h / 2 + if h < hmin || first + rethrow(e) + end + end + end + if dir == -1 + volts = reverse(volts)[1:(end - 1)] + dlcaps = reverse(dlcaps)[1:(end - 1)] + end + end + return volts, dlcaps +end + + +function symmcell_surfacecharge_bcondition!(y, u, bnode, data) + (; iφ, ip) = data + (; iφ, ip) = data + boundary_neumann!(y, u, bnode, species = iφ, region = 2, value = data.q[2] * data.qscale) + boundary_neumann!(y, u, bnode, species = iφ, region = 1, value = data.q[1] * data.qscale) + boundary_dirichlet!(y, u, bnode, species = ip, region = 3, value = 0) + return nothing +end + + +function SurfaceChargedSymmetricCell(grid, data; dielectric_decrement = false, valuetype = Float64) + data = deepcopy(data) + data.nv = ones(num_nodes(grid)) # help to satisfy sparsity detector + data.conserveions = true + data.χvar = dielectric_decrement + + sys = VoronoiFVM.System( + grid; + data = data, + flux = poisson_and_p_flux!, + reaction = reaction!, + bcondition = symmcell_surfacecharge_bcondition!, + generic = ionconservation!, + unknown_storage = :sparse, + valuetype + ) + + # Enable species for all fields + for i in 1:data.N + enable_species!(sys, i, [1]) + end + + enable_species!(sys, data.i0, [1]) + enable_species!(sys, data.iφ, [1]) + enable_species!(sys, data.ip, [1]) + enable_species!(sys, data.iE, [1]) + + # Enable species in mid of the domain for ion conservation and + # electroneutrality constraint + for α in 1:data.N + enable_boundary_species!(sys, data.coffset + α, [3]) + end + data.nv = nodevolumes(sys) + return SurfaceChargedSymmetricCell(sys) +end diff --git a/packages/JuliaMPBSolver/src/icmbp-p.jl b/packages/JuliaMPBSolver/src/icmbp-p.jl index 0c99cfe..4c02ea6 100644 --- a/packages/JuliaMPBSolver/src/icmbp-p.jl +++ b/packages/JuliaMPBSolver/src/icmbp-p.jl @@ -94,7 +94,10 @@ begin n_avg::Vector{Float64} = fill(molarity, N) "Surface charges" - q::Vector{Float64} = [0, 0] * ufac"C/m^2" + q::Vector{Float64} = fill(0 * ufac"C/m^2", 2) + + "Applied potential" + φ::Float64 = 0 * ufac"C/m^2" "Solvent molarity" n0_ref::Float64 = 55.508 * ph"N_A" / ufac"dm^3" @@ -138,6 +141,12 @@ begin "Concentration scaling parameter" cscale::Float64 = ph"N_A" + "Charge scaling parameter" + qscale::Float64 = 1.0 / ph"e" + + "Electric field scaling parameter" + Escale::Float64 = ufac"V/nm" + "Reference voltage" E_ref::Float64 = 0.0 * ufac"V" @@ -168,14 +177,33 @@ begin end end +# ╔═╡ 0089bc79-cc71-45b6-985b-932db659df98 +ICMPBData() + # ╔═╡ 858ed8e1-84b1-4105-8ea0-45209aea40c6 md""" #### `apply_charge!(data,q)` """ # ╔═╡ 4929c105-4c01-4c83-ad2f-2056a8c51d29 -function apply_charge!(data::ICMPBData, q) - data.q .= [q, - q] +function apply_charge!(data::ICMPBData, q::Vector) + data.q = q + return data +end + +function apply_charge!(data::ICMPBData, q::Number) + data.q = [-q, q] + return data +end + +# ╔═╡ 308d7e40-1e15-475e-90a2-394a64e0e8d0 +md""" +#### `apply_voltage!(data, q)` +""" + +# ╔═╡ acc60604-bcc5-4ecf-9772-e2fc851a9232 +function apply_voltage!(data::ICMPBData, φ) + data.φ = φ return data end @@ -210,7 +238,7 @@ L_{Debye}=\sqrt{ \frac{(1+χ)ε_0k_BT}{e^2n_E}} # ╔═╡ a41c6c1f-ceb5-4590-a421-cae5078d167b function L_Debye(data) return sqrt( - (1 + data.χ) * data.ε_0 * ph"k_B" * data.T / (data.e^2 * data.n_E[1]), + (1 + data.χ0) * data.ε_0 * ph"k_B" * data.T / (data.e^2 * data.n_E[1]), ) end; @@ -270,7 +298,7 @@ function W(x) result = 1.0 / 3.0 + u * result ysmall = 3 * result ylarge = 3 * ((coth(x) - 1.0 / x) / x) - return ifelse(abs(x) < 0.5, ysmall, ylarge) + return ifelse(abs(x) < 0.1, ysmall, ylarge) end # ╔═╡ 9aa654e2-4cd9-4b31-af84-8060729cd3a3 @@ -285,7 +313,7 @@ function Λ(x) # thx 深度求索 u = x * x ysmall = log(1.0 + u * (1.0 / 6.0 + u * (1.0 / 120.0 + u * (1.0 / 5040.0 + u * (1.0 / 362880.0))))) ylarge = log(sinh(x) / x) - return ifelse(abs(x) < 0.5, ysmall, ylarge) + return ifelse(abs(x) < 0.1, ysmall, ylarge) end # ╔═╡ a26cf11b-0ce1-4c1d-a64d-1917178ff676 @@ -307,7 +335,7 @@ Ion molar fractions # ╔═╡ 188f67d8-2ae8-474c-8e58-68b8b4fde02e function y_α(φ, p, E, α, data, ddata) η_φ = data.z[α] * data.e * (φ - data.E_ref) - η_p = ddata.v[α] * (p * data.pscale - data.p_ref) + η_p = ddata.v[α] * (p - data.p_ref) return ddata.y_E[α] * exp(-(η_φ + η_p) / (data.kT) + data.χvar * Λ(data.δ[α] * E / data.kT)) end; @@ -320,7 +348,7 @@ Solvent molar fraction # ╔═╡ d7531d5f-fc2d-42b2-9cf9-6a737b0f0f8d function y0(p, E, data, ddata) - return ddata.y0_E * exp(-data.v0 * (p * data.pscale - data.p_ref) / (data.kT) + data.χvar * Λ(data.δ0 * E / data.kT)) + return ddata.y0_E * exp(-data.v0 * (p - data.p_ref) / (data.kT) + data.χvar * Λ(data.δ0 * E / data.kT)) end; # ╔═╡ f6f004a6-d71b-4813-a363-9f51dc37e42a @@ -383,18 +411,22 @@ function spacecharge(u, data) return data.e * sumyz / sumyv end + # ╔═╡ bf7f8bab-4807-440b-8796-fcf75ad313d7 function susceptibility(u, data) - (; iE, i0) = data + (; iE, i0, Escale) = data + E = u[iE] * Escale y = u[i0] - χ = data.v0 * u[i0] * data.χ0 * W(data.δ0 * u[iE] / data.kT) + χ = data.v0 * y * data.χ0 * W(data.δ0 * E / data.kT) sumyv = data.v0 * y for α in 1:(data.N) + y = u[α] v = data.vu[α] + data.κ[α] * data.v0 - χ += v * u[α] * data.χ[α] * W(data.δ[α] * u[iE] / data.kT) - sumyv += v * u[α] + χ += v * y * data.χ[α] * W(data.δ[α] * E / data.kT) + sumyv += v * y end - return χ / sumyv + χ = χ / sumyv + return χ end # ╔═╡ b41838bb-3d5b-499c-9eb5-137c252ae366 @@ -529,8 +561,8 @@ Obtain ion number densities from system # ╔═╡ 800dfed8-9f29-4138-96f8-e8bf1f2f00e6 function calc_cnum(sol, sys) data = sys.physics.data - i3 = sys.grid[BFaceNodes][3][1] if data.conserveions + i3 = sys.grid[BFaceNodes][3][1] ddata = DerivedData(data, sol[(data.coffset + 1):end, i3]) else ddata = DerivedData(data) @@ -549,8 +581,8 @@ end; function calc_c0num(sol, sys) data = sys.physics.data grid = sys.grid - i3 = sys.grid[BFaceNodes][3][1] if data.conserveions + i3 = sys.grid[BFaceNodes][3][1] ddata = DerivedData(data, sol[(data.coffset + 1):end, i3]) else ddata = DerivedData(data) @@ -564,6 +596,11 @@ function calc_c0num(sol, sys) return c0 end; +# ╔═╡ 4baae81f-e68a-4259-a095-4194fde390d2 +md""" +#### calc_χ(sol,sys) +""" + # ╔═╡ 3b7a90cd-8f58-4abc-805a-2891ad6ceb9a function calc_χ(sol, sys) data = sys.physics.data @@ -592,6 +629,12 @@ calc_cmol(sol, sys) = calc_cnum(sol, sys) / (ph"N_A" * ufac"mol/dm^3"); # ╔═╡ 79cc671b-ef6e-42da-8641-61e43f221cb1 calc_c0mol(sol, sys) = calc_c0num(sol, sys) / (ph"N_A" * ufac"mol/dm^3"); +# ╔═╡ dad1836b-17f0-42ca-b527-799a6f5c3d07 +function calc_spacecharge(sys, sol) + data = sys.physics.data + return VoronoiFVM.integrate(sys, (y, u, node, data) -> y[data.iφ] = -spacecharge(u, data), sol)[data.iφ, 1] +end + # ╔═╡ 02d69f1c-4525-4f69-9938-cb0495171c3a md""" #### ysum(sys,sol) @@ -645,7 +688,7 @@ md""" # ╔═╡ 05695820-fa21-49b7-b52f-8a94cf2fa0fa md""" -We use N+3 fields of unknowns in the following sequence: +We use N+4 fields of unknowns in the following sequence: ``y_1 \dots y_N, y_0, \varphi, p, E``. Pressures are scaled by `pscale` (default: 1GPa). """ @@ -660,10 +703,8 @@ Callback which runs in every grid point. # ╔═╡ e1c13f1e-5b67-464b-967b-25e3a93e33d9 function reaction!(f, u, node, data) - (; i0, iφ, ip, N) = data - φ = u[iφ] - p = u[ip] - f[iφ] = -spacecharge(u, data) + # Everything dependent on u is on the left side of the equation, therefore the minus sign + f[data.iφ] = -spacecharge(u, data) * data.qscale return end; @@ -672,37 +713,46 @@ md""" #### `poisson_and_p_flux!(f, u, edge, data)` Runs on every grid edge. Calculate fluxes for the Poisson and the pressure equations. +```math + -\nabla \cdot (\varepsilon_0(1+χ) \nabla \varphi) - q =0 +``` """ # ╔═╡ 64e47917-9c61-4d64-a6a1-c6e8c7b28c59 function poisson_and_p_flux!(f, u, edge, data) - (; iφ, ip, iE, N) = data + (; iφ, ip, iE, Escale, N) = data nspec = size(u, 1) - E1 = u[iE, 1] - E2 = u[iE, 2] - uu = zeros(eltype(u), nspec) + T = eltype(u) + uu1 = zeros(eltype(u), nspec) + uu2 = zeros(eltype(u), nspec) for i in 1:nspec - uu[i] = u[i, 1] + uu1[i] = u[i, 1] + uu2[i] = u[i, 2] end - χ1 = susceptibility(uu, data) - q1 = spacecharge(uu, data) - for i in 1:nspec - uu[i] = u[i, 2] - end - χ2 = susceptibility(uu, data) - q2 = spacecharge(uu, data) + + q1 = spacecharge(uu1, data) + q2 = spacecharge(uu2, data) + + χ = T(data.χ0) + E = zero(T) + χ1 = zero(T) + χ2 = zero(T) if data.χvar - χ = (χ1 + χ2) / 2 - else - χ = data.χ0 + + χ1 = susceptibility(uu1, data) + χ2 = susceptibility(uu2, data) + χ = 2 / (1 / χ1 + 1 / χ2) + E1 = u[iE, 1] * Escale + E2 = u[iE, 2] * Escale + E = (E1 + E2) / 2 end - χ = data.χ0 - f[iφ] = (1.0 + χ) * data.ε_0 * (u[iφ, 1] - u[iφ, 2]) - f[ip] = - (u[ip, 1] - u[ip, 2]) + (u[iφ, 1] - u[iφ, 2]) * (q1 + q2) / (2 * data.pscale) - - data.χvar * (data.ε_0 / 2) * ((E1 + E2) / 2)^2 * (χ1 - χ2) + # Poisson equation is scaled with qscale + f[iφ] = (1.0 + χ) * data.ε_0 * (u[iφ, 1] - u[iφ, 2]) * data.qscale + + # pressure equation is scaled with 1/pscale + f[ip] = (u[ip, 1] - u[ip, 2]) + (u[iφ, 1] - u[iφ, 2]) * (q1 + q2) / (2 * data.pscale) + data.ε_0 * E^2 * (χ1 - χ2) / (2 * data.pscale) return end; @@ -716,8 +766,8 @@ Boundary condition callback. The Dirichlet condition for the pressure in the mid # ╔═╡ 743b9a7a-d6ac-4da0-8538-2045d965b547 function bcondition!(y, u, bnode, data) (; iφ, ip) = data - boundary_neumann!(y, u, bnode, species = iφ, region = 2, value = data.q[2]) - boundary_neumann!(y, u, bnode, species = iφ, region = 1, value = data.q[1]) + boundary_neumann!(y, u, bnode, species = iφ, region = 2, value = data.q[2] * data.qscale) + boundary_neumann!(y, u, bnode, species = iφ, region = 1, value = data.q[1] * data.qscale) boundary_dirichlet!(y, u, bnode, species = ip, region = 3, value = 0) return nothing end @@ -759,11 +809,12 @@ end # ╔═╡ dbccaa88-65d9-47ab-be78-83df64a6db24 function ionconservation!(f, u, sys, data) - (; coffset, i0, iφ, iE, ip, N, z, nv, n_avg) = data + (; coffset, i0, iφ, ip, iE, Escale, pscale, N, z, nv, n_avg) = data # Set the result to zero f .= 0 - # Find mid-of-the-domain node number from boundary region 3 - i3 = sys.grid[BFaceNodes][3][1] + + i3 = 0 + X = sys.grid[Coordinates][1, :] # Parameters u and f come as vectors, `idx` allows to access their contents with @@ -772,6 +823,8 @@ function ionconservation!(f, u, sys, data) # Obtain values of the bulk molecular densities if data.conserveions + # Find mid-of-the-domain node number from boundary region 3 + i3 = sys.grid[BFaceNodes][3][1] n_E = [u[idx[coffset + i, i3]] * data.cscale for i in 1:N] # Calculate derived data ddata = DerivedData(data, n_E) @@ -779,18 +832,21 @@ function ionconservation!(f, u, sys, data) ddata = DerivedData(data) end - # Calculate molar fractions for all nodes of the grid for i in 1:num_nodes(sys.grid) - f[idx[i0, i]] = u[idx[i0, i]] - y0(u[idx[ip, i]], u[idx[iE, i]], data, ddata) + + # Calculate molar fractions + f[idx[i0, i]] = u[idx[i0, i]] - y0(u[idx[ip, i]] * pscale, u[idx[iE, i]] * Escale, data, ddata) for α in 1:N - f[idx[α, i]] = u[idx[α, i]] - y_α(u[idx[iφ, i]], u[idx[ip, i]], u[idx[iE, i]], α, data, ddata) + f[idx[α, i]] = u[idx[α, i]] - y_α(u[idx[iφ, i]], u[idx[ip, i]] * pscale, u[idx[iE, i]] * Escale, α, data, ddata) end + + # Calculate electric field strength if i == 1 - f[idx[iE, i]] = u[idx[iE, i]] - abs((u[idx[iφ, i + 1]] - u[idx[iφ, i]]) / (X[i + 1] - X[i])) + f[idx[iE, i]] = u[idx[iE, i]] - (u[idx[iφ, i + 1]] - u[idx[iφ, i]]) / (X[i + 1] - X[i]) / Escale elseif i == num_nodes(sys.grid) - f[idx[iE, i]] = u[idx[iE, i]] - abs((u[idx[iφ, i - 1]] - u[idx[iφ, i]]) / (X[i - 1] - X[i])) + f[idx[iE, i]] = u[idx[iE, i]] - (u[idx[iφ, i - 1]] - u[idx[iφ, i]]) / (X[i - 1] - X[i]) / Escale else - f[idx[iE, i]] = u[idx[iE, i]] - abs((u[idx[iφ, i + 1]] - u[idx[iφ, i - 1]]) / (X[i + 1] - X[i - 1])) + f[idx[iE, i]] = u[idx[iE, i]] - (u[idx[iφ, i + 1]] - u[idx[iφ, i - 1]]) / (X[i + 1] - X[i - 1]) / Escale end end @@ -807,9 +863,9 @@ function ionconservation!(f, u, sys, data) f[idx[coffset + N, i3]] += z[α] * u[idx[coffset + α, i3]] / z[N] end + # Calculate number density integrals y = zeros(eltype(u), N + 1) uu = zeros(eltype(u), N + 1) - # Calculate number density integrals for iv in 1:length(nv) for α in 1:(N + 1) uu[α] = u[idx[α, iv]] @@ -826,7 +882,7 @@ function ionconservation!(f, u, sys, data) end # ╔═╡ 7bf3a130-3b47-428e-916f-4a0ec1237844 -function ICMPBSystem(grid, data) +function ICMPBSystem(grid, data; valuetype = Float64) data.nv = ones(num_nodes(grid)) # trigger sparsity detector sys = VoronoiFVM.System( @@ -836,7 +892,8 @@ function ICMPBSystem(grid, data) reaction = reaction!, bcondition = bcondition!, generic = ionconservation!, - unknown_storage = :sparse + unknown_storage = :sparse, + valuetype ) # Enable species for all fields @@ -855,8 +912,8 @@ function ICMPBSystem(grid, data) for α in 1:data.N enable_boundary_species!(sys, data.coffset + α, [3]) end - data.nv = nodevolumes(sys) end + data.nv = nodevolumes(sys) return sys end; @@ -874,9 +931,6 @@ Sweep over series of surface charges and calculate resulting potential difference. """ -# ╔═╡ 3b389ecf-4c63-4eb8-b8be-76442eacef80 - - # ╔═╡ 178b947f-3fef-44ed-9eca-fdb9916bc2b6 function qsweep(sys; qmax = 10, nsteps = 100, verbose = "", kwargs...) data = deepcopy(sys.physics.data) @@ -943,8 +997,11 @@ end # ╠═55b2ee36-c4f9-4ba3-84ed-faeb556aa026 # ╟─6e4aaa60-29c5-4f75-a3f1-24e340c25e6c # ╠═0d825f88-cd67-4368-90b3-29f316b72e6e +# ╠═0089bc79-cc71-45b6-985b-932db659df98 # ╟─858ed8e1-84b1-4105-8ea0-45209aea40c6 # ╠═4929c105-4c01-4c83-ad2f-2056a8c51d29 +# ╟─308d7e40-1e15-475e-90a2-394a64e0e8d0 +# ╠═acc60604-bcc5-4ecf-9772-e2fc851a9232 # ╟─f3049938-2637-401d-9411-4d7be07c19ca # ╟─e69e10cc-e21a-418d-90b2-ae218dca0c73 # ╠═5d6340c4-2ddd-429b-a60b-3de5570a7398 @@ -982,10 +1039,12 @@ end # ╟─0c54efd0-f279-4dc6-8b00-ba092dd13f44 # ╠═800dfed8-9f29-4138-96f8-e8bf1f2f00e6 # ╠═24910762-7d56-446b-a758-d8e830fe9a09 +# ╟─4baae81f-e68a-4259-a095-4194fde390d2 # ╠═3b7a90cd-8f58-4abc-805a-2891ad6ceb9a # ╟─9fe3ca93-c051-426e-8b9a-cc59f59319ad # ╠═2ee34d76-7238-46c2-94d1-a40d8b017af6 # ╠═79cc671b-ef6e-42da-8641-61e43f221cb1 +# ╠═dad1836b-17f0-42ca-b527-799a6f5c3d07 # ╟─02d69f1c-4525-4f69-9938-cb0495171c3a # ╠═48670f54-d303-4c3a-a191-06e6592a2e0a # ╟─7a607454-7b75-4313-920a-2dbdad258015 @@ -995,7 +1054,7 @@ end # ╟─05695820-fa21-49b7-b52f-8a94cf2fa0fa # ╟─b9b0cb4f-cf72-418e-a65e-0f4c8a10e34c # ╠═e1c13f1e-5b67-464b-967b-25e3a93e33d9 -# ╟─c1168002-a716-4568-9a52-ac00f32f0134 +# ╠═c1168002-a716-4568-9a52-ac00f32f0134 # ╠═64e47917-9c61-4d64-a6a1-c6e8c7b28c59 # ╟─05acd04f-74af-42f9-b039-7ee5b2ba63ff # ╠═743b9a7a-d6ac-4da0-8538-2045d965b547 @@ -1007,7 +1066,6 @@ end # ╠═b0a45e53-8b98-4e18-8b41-7f6d0bc1f76e # ╟─3dc0d408-ab00-4da0-9999-2ddf6e4fbf60 # ╟─a67c0d46-456d-4b0c-8519-961b37043350 -# ╠═3b389ecf-4c63-4eb8-b8be-76442eacef80 # ╠═178b947f-3fef-44ed-9eca-fdb9916bc2b6 # ╟─4032bc46-7820-45d4-bc11-9350ecf1797a # ╠═fc84996b-02c0-4c16-8632-79f4e1900f78 diff --git a/results/halfcell.png b/results/halfcell.png new file mode 100644 index 0000000..0286dc9 Binary files /dev/null and b/results/halfcell.png differ diff --git a/results/halfcell_dlcap_M.png b/results/halfcell_dlcap_M.png new file mode 100644 index 0000000..3237105 Binary files /dev/null and b/results/halfcell_dlcap_M.png differ diff --git a/results/halfcell_dlcap_kappa.png b/results/halfcell_dlcap_kappa.png new file mode 100644 index 0000000..43cf3a2 Binary files /dev/null and b/results/halfcell_dlcap_kappa.png differ diff --git a/results/largecell.png b/results/largecell.png new file mode 100644 index 0000000..5f7a903 Binary files /dev/null and b/results/largecell.png differ diff --git a/results/smallcell.png b/results/smallcell.png new file mode 100644 index 0000000..996bb7e Binary files /dev/null and b/results/smallcell.png differ diff --git a/src/PoissonBoltzmannIPAM2025.jl b/src/PoissonBoltzmannIPAM2025.jl index 530e875..994476b 100644 --- a/src/PoissonBoltzmannIPAM2025.jl +++ b/src/PoissonBoltzmannIPAM2025.jl @@ -1,7 +1,18 @@ module PoissonBoltzmannIPAM2025 -using DrWatson +using DrWatson, PoissonBoltzmannIPAM2025 +using PythonPlot +using PythonPlot: pyplot +using LessUnitful +using ExtendableGrids +using VoronoiFVM +using Test +using JuliaMPBSolver.ICMPBP: ICMPBP, SurfaceChargedSymmetricCell, AbstractHalfCell, AbstractSymmetricCell, set_molarity!, calc_cmol, calc_c0mol, calc_χ, get_E, get_φ, get_p, get_c0, + set_κ!, set_q!, set_φ! resultsdir(args...) = projectdir("results", args...) +draftresultsdir(args...) = projectdir("draftresults", args...) -export resultsdir +include("plotcells.jl") + +export resultsdir, draftresultsdir, plotcells end diff --git a/src/plotcells.jl b/src/plotcells.jl new file mode 100644 index 0000000..24c6256 --- /dev/null +++ b/src/plotcells.jl @@ -0,0 +1,251 @@ +surfcharge(n) = n * ph"e" / ufac"nm^2" + +function makecolors(V) + h = 0.5 / length(V) + return colors = [ (i * h, 0, 1 - i * h) for i in 1:length(V) ] +end + +function plotcells( + cell::AbstractSymmetricCell, celldd::AbstractSymmetricCell; + figname = "symmcell", + κ = 10, M = 1, q = 1.0, + title = latexstring(""" \\text{Solutions for } \\kappa=$(κ), M_{avg}=$(M)mol/L, q=$(q)e/nm^2""") * + "\n" * + latexstring("""\\text{Dashed: without dielectric decrement} """) + + ) + pyplot.close() + clf() + fig, ax = pyplot.subplots(3, 1) + fig.set_size_inches(600 / 100, 600 / 100) + ax1 = ax[0] + ax2 = ax[1] + ax3 = ax[2] + set_κ!(cell, κ) + set_κ!(celldd, κ) + set_molarity!(cell, M) + set_molarity!(celldd, M) + fig.suptitle(title) + grid = cell.sys.grid + X = grid[Coordinates][1, :] / ufac"nm" + solc = unknowns(cell) + for _q in range(0, q, length = 5) + set_q!(cell, surfcharge(_q)) + solc = solve(cell, inival = solc) + end + + cc = calc_cmol(solc, cell) + c0c = calc_c0mol(solc, cell) + φc = get_φ(solc, cell) + pc = get_p(solc, cell) / ufac"GPa" + Ec = get_E(solc, cell) .|> abs + χc = calc_χ(solc, cell) + + + solv = unknowns(celldd) + for _q in range(0, q, length = 10) + set_q!(celldd, surfcharge(_q)) + solv = solve(celldd, inival = solv, damp_initial = 0.1) + end + cv = calc_cmol(solv, celldd) + c0v = calc_c0mol(solv, celldd) + φv = get_φ(solv, celldd) + pv = get_p(solv, celldd) / ufac"GPa" + Ev = get_E(solv, celldd) .|> abs + χv = calc_χ(solv, celldd) + + ax1r = ax1.twinx() + ax1.plot(X, φc, color = "olive", linestyle = "--") + ax1r.plot(X, pc, color = "darkorange", linestyle = "--") + ax1.plot(X, φv, color = "olive", label = "ϕ") + ax1r.plot(X, pv, color = "darkorange", label = "p") + ax1.set_xlabel("z/nm") + ax1.legend(loc = (0.2, 0.1)) + ax1.set_ylabel("ϕ/V") + φ_max = ceil(maximum(φv) * 10) / 10 + ax1r.set_ylabel("p/GPa") + ax1r.legend(loc = (0.7, 0.1)) + ax1.set_ylim((-φ_max, φ_max)) + ax1.set_yticks(range(-φ_max, φ_max, length = 5)) + + pmax = ceil(20 * maximum(pc)) / 10 + ax1r.set_ylim((0, pmax)) + ax1r.set_yticks(range(0, pmax, length = 5)) + ax1.grid() + + + ax2.grid() + + ax2.set_xlabel("z/nm") + ax2.set_ylabel("c/(mol/L)") + ax2r = ax2.twinx() + ax2r.set_ylabel(L"c_{solvent}/(mol/L)") + + + ax2.plot(X, cv[1, :], color = "blue", label = L"c^-") + ax2.plot(X, cv[2, :], color = "red", label = L"c^+") + ax2r.plot(X, c0v, color = "green", label = L"c_{solvent}") + + ax2.plot(X, cc[1, :], color = "blue", linestyle = "--") + ax2.plot(X, cc[2, :], color = "red", linestyle = "--") + ax2r.plot(X, c0c, color = "green", linestyle = "--") + ax2.legend(loc = (0.15, 0.6)) + ax2r.legend(loc = (0.7, 0.7)) + ax2.set_ylim((0, 5)) + ax2.set_yticks(range(0, 5, length = 5)) + ax2r.set_ylim((0, 60)) + ax2r.set_yticks(range(0, 60, length = 5)) + + ax3.grid() + ax3r = ax3.twinx() + ax3.plot(X, Ev / ufac"V / nm", color = "darkviolet", label = "|E|") + ax3r.plot(X, χv, color = "steelblue", label = "χ") + ax3.plot(X, Ec / ufac"V / nm", color = "darkviolet", linestyle = "--") + ax3.legend(loc = (0.2, 0.2)) + ax3r.plot(X, χc, color = "steelblue", linestyle = "--") + ax3r.legend(loc = (0.65, 0.2)) + ax3.set_xlabel("z/nm") + Emax = ceil(maximum(Ev / ufac"V / nm")) + ax3.set_ylabel("|E|/(V/nm)") + ax3.set_ylim((0, Emax)) + ax3.set_yticks(range(0, Emax, length = 5)) + ax3r.set_ylabel("χ") + ax3r.set_ylim((0, 100)) + ax3r.set_yticks(range(0, 100, length = 5)) + + tight_layout() + !haskey(ENV, "CI") && savefig(draftresultsdir(figname), dpi = 600) + return gcf() + +end + + +function plotcells( + cell::AbstractHalfCell, celldd::AbstractHalfCell; + figname = "halfcell", + κ = 10, M = 1, Δφ = 0.5, + title = latexstring(""" \\text{Solutions for } κ=$(κ), M=$(M)mol/L, \\Delta \\varphi=$(Δφ)V""") * + "\n" * + latexstring("""\\text{Dashed: without dielectric decrement} """) + ) + pyplot.close() + clf() + fig, ax = pyplot.subplots(3, 1) + fig.set_size_inches(600 / 100, 600 / 100) + ax1 = ax[0] + ax2 = ax[1] + ax3 = ax[2] + set_κ!(cell, κ) + set_κ!(celldd, κ) + set_molarity!(cell, M) + set_molarity!(celldd, M) + fig.suptitle(title) + grid = cell.sys.grid + X = grid[Coordinates][1, :] / ufac"nm" + solc = unknowns(cell) + for _Δφ in range(0, Δφ, length = 5) + set_φ!(cell, _Δφ) + solc = solve(cell, inival = solc) + end + + cc = calc_cmol(solc, cell) + c0c = calc_c0mol(solc, cell) + φc = get_φ(solc, cell) + pc = get_p(solc, cell) / ufac"GPa" + Ec = get_E(solc, cell) .|> abs + χc = calc_χ(solc, cell) + + + solv = unknowns(celldd) + set_φ!(celldd, 0) + solv = solve(celldd, inival = solv, damp_initial = 0.5) + h = 0.001 + v = 0.0 + while v < Δφ + try + vv = min(v + h, Δφ) + set_φ!(celldd, vv) + solv = solve(celldd, inival = solv, damp_initial = 0.5) + v = vv + h = h * 1.2 + catch e + h = h / 2 + if h < 1.0e-5 + rethrow(e) + end + end + end + #set_φ!(celldd, Δφ) + #solv = solve(celldd, inival = solc, damp_initial = 0.1) + + cv = calc_cmol(solv, celldd) + c0v = calc_c0mol(solv, celldd) + φv = get_φ(solv, celldd) + pv = get_p(solv, celldd) / ufac"GPa" + Ev = get_E(solv, celldd) .|> abs + χv = calc_χ(solv, celldd) + + ax1r = ax1.twinx() + ax1.plot(X, φc, color = "olive", linestyle = "--") + ax1r.plot(X, pc, color = "darkorange", linestyle = "--") + ax1.plot(X, φv, color = "olive", label = "ϕ") + ax1r.plot(X, pv, color = "darkorange", label = L"p") + ax1.set_xlabel("z/nm") + ax1.legend(loc = (0.2, 0.7)) + ax1.set_ylabel("ϕ/V") + φ_max = ceil(maximum(φv) * 10) / 10 + ax1r.set_ylabel("p/GPa") + ax1r.legend(loc = (0.7, 0.7)) + ax1.set_ylim((0, φ_max)) + ax1.set_yticks(range(0, φ_max, length = 5)) + + pmax = ceil(20 * maximum(pc)) / 10 + ax1r.set_ylim((0, pmax)) + ax1r.set_yticks(range(0, pmax, length = 5)) + ax1.grid() + + + ax2.grid() + + ax2.set_xlabel("z/nm") + ax2.set_ylabel("c/(mol/L)") + ax2r = ax2.twinx() + ax2r.set_ylabel(L"c_{solvent}/(mol/L)") + + + ax2.plot(X, cv[1, :], color = "blue", label = L"c^-") + ax2.plot(X, cv[2, :], color = "red", label = L"c^+") + ax2r.plot(X, c0v, color = "green", label = L"c_{solvent}") + + ax2.plot(X, cc[1, :], color = "blue", linestyle = "--") + ax2.plot(X, cc[2, :], color = "red", linestyle = "--") + ax2r.plot(X, c0c, color = "green", linestyle = "--") + ax2.legend(loc = (0.2, 0.6)) + ax2r.legend() + ax2.set_ylim((0, 5)) + ax2.set_yticks(range(0, 5, length = 5)) + ax2r.set_ylim((0, 60)) + ax2r.set_yticks(range(0, 60, length = 5)) + + ax3.grid() + ax3r = ax3.twinx() + ax3.plot(X, Ev / ufac"V / nm", color = "darkviolet", label = "E") + ax3r.plot(X, χv, color = "steelblue", label = "χ") + ax3.plot(X, Ec / ufac"V / nm", color = "darkviolet", linestyle = "--") + ax3.legend(loc = (0.2, 0.1)) + ax3r.plot(X, χc, color = "steelblue", linestyle = "--") + ax3r.legend(loc = (0.8, 0.1)) + ax3.set_xlabel("z/nm") + Emax = ceil(maximum(Ev / ufac"V / nm")) + ax3.set_ylabel("|E|/(V/nm)") + ax3.set_ylim((0, Emax)) + ax3.set_yticks(range(0, Emax, length = 5)) + ax3r.set_ylabel("χ") + ax3r.set_ylim((0, 100)) + ax3r.set_yticks(range(0, 100, length = 5)) + + tight_layout() + !haskey(ENV, "CI") && savefig(draftresultsdir(figname), dpi = 600) + return gcf() + +end diff --git a/test/runtests.jl b/test/runtests.jl index 9f86fad..0a49493 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,3 +1,18 @@ -using Pkg, JuliaMPBSolver +using Pkg, JuliaMPBSolver, ExampleJuggler, Test, Markdown + +ExampleJuggler.verbose!(true) + +notebooks = [ + "MPBP-Draft.jl", + "ICMPBP-Draft.jl", + "ICMPBP-EndOfHackathon.jl", + "ICMPBP-DD-Draft.jl", + "HalfCellAppliedPotential.jl", + "SymmetricCellSurfaceCharge.jl", +] + +@testset "Notebooks" begin + @testscripts(joinpath(@__DIR__, "..", "notebooks"), notebooks) +end Pkg.test("JuliaMPBSolver") diff --git a/tex/draft.tex b/tex/draft.tex index 5fdc18b..b2f53fc 100644 --- a/tex/draft.tex +++ b/tex/draft.tex @@ -76,7 +76,7 @@ \subsection{Poisson equation} We set the polarization free energy contribution as $$ -\rho\psi^{Pol} = -k_BT\sum_{\alpha} n_\alpha\Lambda\left(\frac{\delta_\alpha |E|}{k_BT}\right) + \rho\psi^{Pol} = -k_BT\sum_{\alpha} n_\alpha\Lambda\left(\frac{\delta_\alpha |E|}{k_BT}\right) $$ where $E=\nabla \phi$ is the electric field. @@ -153,7 +153,7 @@ \subsubsection{Second order pressure PDE (``Pressure Poisson'')} First, one can use the momentum balance in mechanical equilibrium \begin{align} \label{eq:momentumbalance} - \nabla p = q\nabla \phi - \frac{\varepsilon_0}{2} |E|^2\nabla \chi + \nabla p = - q\nabla \phi - \frac{\varepsilon_0}{2} |E|^2\nabla \chi \end{align} for the pressure derived as part of the generalized Nernst-Planck-Poisson @@ -161,7 +161,7 @@ \subsubsection{Second order pressure PDE (``Pressure Poisson'')} a second order PDE \begin{equation}\label{eq:pressurePDE} \begin{aligned} - - \Delta p & = q\nabla\phi \; \text{in} \Omega. + - \Delta p & = \nabla\cdot (q\nabla\phi) + \nabla\cdot ( \frac{\varepsilon_0}{2} |E|^2\nabla \chi) \; \text{in} \Omega. \end{aligned} \end{equation} The boundary conditions are: