Scrutiny tests for computational results

.gitignore

@@ -3,3 +3,14 @@ Manifest*.toml
 __pycache__
 condaenv
 **/venv
+draftresults
+draft.pdf
+*.nav
+*.snm
+*.toc
+*.aux
+*.log
+*.bbl
+*.blg
+*.out
+*.tdo

Project.toml

@@ -11,18 +11,21 @@ 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"
 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"]

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"""<hr>"""
+
+# ╔═╡ 784b4c3e-bb2a-4940-a83a-ed5e5898dfd4
+html"""<style>.dont-panic{ display: none }</style>"""
+
+# ╔═╡ afe4745f-f9f1-4e23-8735-cbec6fb79c41
+begin
+    function floataside(text::Markdown.MD; top = 1)
+        uuid = uuid1()
+        return @htl(
+            """
+            		<style>
+
+
+            		@media (min-width: calc(700px + 30px + 300px)) {
+            			aside.plutoui-aside-wrapper-$(uuid) {
+
+            	color: var(--pluto-output-color);
+            	position:fixed;
+            	right: 1rem;
+            	top: $(top)px;
+            	width: 400px;
+            	padding: 10px;
+            	border: 3px solid rgba(0, 0, 0, 0.15);
+            	border-radius: 10px;
+            	box-shadow: 0 0 11px 0px #00000010;
+            	/* That is, viewport minus top minus Live Docs */
+            	max-height: calc(100vh - 5rem - 56px);
+            	overflow: auto;
+            	z-index: 40;
+            	background-color: var(--main-bg-color);
+            	transition: transform 300ms cubic-bezier(0.18, 0.89, 0.45, 1.12);
+
+            			}
+            			aside.plutoui-aside-wrapper > div {
+            #				width: 300px;
+            			}
+            		}
+            		</style>
+
+            		<aside class="plutoui-aside-wrapper-$(uuid)">
+            		<div>
+            		$(text)
+            		</div>
+            		</aside>
+
+            		"""
+        )
+    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"""
+<script>
+	const button = document.createElement("button")
+
+	button.addEventListener("click", () => {
+		editor_state_set(old_state => ({
+			notebook: {
+				...old_state.notebook,
+				process_status: "no_process",
+			},
+		})).then(() => {
+			window.requestAnimationFrame(() => {
+				document.querySelector("#process_status a").click()
+			})
+		})
+	})
+	button.innerText = "Restart notebook"
+
+	return button
+</script>
+""";
+
+# ╔═╡ 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