Cleanup of script testing (#62)

* some more elegant tests for the scripts

* missed a little

* with bsk

Author: landinjm

Project.toml

@@ -1,7 +1,7 @@
 name = "PoissonBoltzmannIPAM2025"
 uuid = "480e0da0-c3fa-11f0-becb-6b90fb73a51e"
-version = "0.1.0"
 authors = ["Jürgen Fuhrmann <juergen-fuhrmann@web.de>"]
+version = "0.1.0"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
@@ -13,6 +13,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 HypertextLiteral = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JuliaMPBSolver = "d8b18f01-5396-498d-b34d-247825c18ff0"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LessUnitful = "f29f6376-6e90-4d80-80c9-fb8ec61203d5"
@@ -42,6 +43,7 @@ ForwardDiff = "1"
 HypertextLiteral = "0.9"
 InteractiveUtils = "1.11.0"
 Interpolations = "0.16"
+JLD2 = "0.6.3"
 LaTeXStrings = "1.0.1"
 LessUnitful = "1"
 Markdown = "1.11.0"

draftnotebooks/demo-notebook.jl

@@ -7,6 +7,7 @@ using Interpolations
 using VoronoiFVM
 using PythonPlot
 using JuliaMPBSolver
+using JLD2
 
 const nel = 20.0 * JuliaMPBSolver.Units.el_surface_density # number of electrons/nm^2 at interfaces
 const M_bulk = 1 # (bulk) molarity at center of domain
@@ -51,120 +52,4 @@ solution, X, nv, ε_r =
     computed_parameters,
 )
 
-function bee!(y, ϕ)
-    N = length(user_parameters.charge_numbers)
-    for i in 1:N
-        y[i] =
-            JuliaMPBSolver.Units.thermal_energy(user_parameters.temperature) *
-            log(c_bulk[i] / computed_parameters.bulk_solvent_concentration) /
-            JuliaMPBSolver.Units.F - user_parameters.charge_numbers[i] * ϕ
-    end
-    return nothing
-end
-
-function bee(sol)
-    n = size(sol, 2)
-    N = length(user_parameters.charge_numbers)
-    e = zeros(N, n)
-    y = zeros(N)
-    for i in 1:n
-        bee!(y, sol[1, i])
-        for j in 1:N
-            e[j, i] = y[j]
-        end
-    end
-    return e
-end
-
-function plotsol(sol; size = (600, 400))
-    PythonPlot.clf()
-    fig, ax = pyplot.subplots(2, 1)
-    ax1 = ax[0]
-    ax2 = ax[1]
-    ax1.grid()
-    ax1r = ax1.twinx()
-
-    c = JuliaMPBSolver.Postprocess.compute_concentrations(
-        sol[1, :],
-        user_parameters,
-        computed_parameters,
-    )
-    cm = c[1, :] ⋅ nv / (JuliaMPBSolver.Units.M) / grid_parameters.domain_size
-    cp = c[2, :] ⋅ nv / (JuliaMPBSolver.Units.M) / grid_parameters.domain_size
-    c0 = -(sum(c, dims = 1) .- c̄)
-    e = bee(sol)
-    ax1.set_title(
-        "ϕ∈$(round.(Float64.(extrema(sol[1, :])), sigdigits = 3)), ε_r ∈$(round.(Float64.(extrema(ε_r)), sigdigits = 3))",
-    )
-    ax1.plot(
-        X / JuliaMPBSolver.Units.nm,
-        sol[1, :],
-        color = "green",
-        linewidth = 2,
-        label = "ϕ",
-    )
-    ax1.plot(
-        X / JuliaMPBSolver.Units.nm,
-        e[1, :],
-        color = "blue",
-        linewidth = 2,
-        label = L"ψ^+",
-        linestyle = "dotted",
-    )
-    ax1.plot(
-        X / JuliaMPBSolver.Units.nm,
-        e[2, :],
-        color = "red",
-        linewidth = 2,
-        linestyle = "dotted",
-        label = L"ψ^+",
-    )
-    ax1r.plot(
-        X[1:(end - 1)] / JuliaMPBSolver.Units.nm,
-        ε_r,
-        color = "pink",
-        linewidth = 3,
-        label = L"ε_r",
-    )
-    ax1.set_ylim(-10, 10)
-    ax1.set_xlabel("z/nm")
-    ax1.set_ylabel("ϕ/V")
-    ax1.legend(loc = (0.1, 0.1))
-    ax1r.legend(loc = (0.8, 0.1))
-    ax1r.set_ylim(0, 80)
-
-    ax2.grid()
-    ax2.set_title("M_avg=$(round.((cm, cp), sigdigits = 3))")
-    ax2.set_xlabel("z/nm")
-    ax2.set_ylabel("c/(mol/L)")
-    ax2.set_ylim(0, 60)
-
-    ax2.plot(
-        X / JuliaMPBSolver.Units.nm,
-        c[1, :] / (JuliaMPBSolver.Units.M),
-        color = "blue",
-        linewidth = 2,
-        label = L"c^-",
-    )
-    ax2.plot(
-        X / JuliaMPBSolver.Units.nm,
-        c[2, :] / (JuliaMPBSolver.Units.M),
-        color = "red",
-        linewidth = 2,
-        label = L"c^+",
-    )
-    ax2.plot(
-        X / JuliaMPBSolver.Units.nm,
-        c0[1, :] / (JuliaMPBSolver.Units.M),
-        color = "green",
-        linewidth = 2,
-        label = L"c_{solvent}",
-    )
-    ax2.legend(loc = (0.4, 0.1))
-
-    tight_layout()
-    savefig("simplecell-bsk.jpg", dpi = 300)
-    return PythonPlot.gcf()
-end
-
-plotsol(solution)
+save("simplecell-bsk.jld2", "solution", solution, "X", X, "nv", nv, "ε_r", ε_r)

scripts/simplecell-bsk.jl

@@ -7,6 +7,7 @@ using Interpolations
 using VoronoiFVM
 using PythonPlot
 using JuliaMPBSolver
+using JLD2
 
 const nel = 20.0 * JuliaMPBSolver.Units.el_surface_density # number of electrons/nm^2 at interfaces
 const M_bulk = 1 # (bulk) molarity at center of domain
@@ -52,120 +53,4 @@ solution, X, nv, ε_r =
     computed_parameters,
 )
 
-function bee!(y, ϕ)
-    N = length(user_parameters.charge_numbers)
-    for i in 1:N
-        y[i] =
-            JuliaMPBSolver.Units.thermal_energy(user_parameters.temperature) *
-            log(c_bulk[i] / computed_parameters.bulk_solvent_concentration) /
-            JuliaMPBSolver.Units.F - user_parameters.charge_numbers[i] * ϕ
-    end
-    return nothing
-end
-
-function bee(sol)
-    n = size(sol, 2)
-    N = length(user_parameters.charge_numbers)
-    e = zeros(N, n)
-    y = zeros(N)
-    for i in 1:n
-        bee!(y, sol[1, i])
-        for j in 1:N
-            e[j, i] = y[j]
-        end
-    end
-    return e
-end
-
-function plotsol(sol; size = (600, 400))
-    PythonPlot.clf()
-    fig, ax = pyplot.subplots(2, 1)
-    ax1 = ax[0]
-    ax2 = ax[1]
-    ax1.grid()
-    ax1r = ax1.twinx()
-
-    c = JuliaMPBSolver.Postprocess.compute_concentrations(
-        sol[1, :],
-        user_parameters,
-        computed_parameters,
-    )
-    cm = c[1, :] ⋅ nv / (JuliaMPBSolver.Units.M) / grid_parameters.domain_size
-    cp = c[2, :] ⋅ nv / (JuliaMPBSolver.Units.M) / grid_parameters.domain_size
-    c0 = -(sum(c, dims = 1) .- c̄)
-    e = bee(sol)
-    ax1.set_title(
-        "ϕ∈$(round.(Float64.(extrema(sol[1, :])), sigdigits = 3)), ε_r ∈$(round.(Float64.(extrema(ε_r)), sigdigits = 3))",
-    )
-    ax1.plot(
-        X / JuliaMPBSolver.Units.nm,
-        sol[1, :],
-        color = "green",
-        linewidth = 2,
-        label = "ϕ",
-    )
-    ax1.plot(
-        X / JuliaMPBSolver.Units.nm,
-        e[1, :],
-        color = "blue",
-        linewidth = 2,
-        label = L"ψ^+",
-        linestyle = "dotted",
-    )
-    ax1.plot(
-        X / JuliaMPBSolver.Units.nm,
-        e[2, :],
-        color = "red",
-        linewidth = 2,
-        linestyle = "dotted",
-        label = L"ψ^+",
-    )
-    ax1r.plot(
-        X[1:(end - 1)] / JuliaMPBSolver.Units.nm,
-        ε_r,
-        color = "pink",
-        linewidth = 3,
-        label = L"ε_r",
-    )
-    ax1.set_ylim(-10, 10)
-    ax1.set_xlabel("z/nm")
-    ax1.set_ylabel("ϕ/V")
-    ax1.legend(loc = (0.1, 0.1))
-    ax1r.legend(loc = (0.8, 0.1))
-    ax1r.set_ylim(0, 80)
-
-    ax2.grid()
-    ax2.set_title("M_avg=$(round.((cm, cp), sigdigits = 3))")
-    ax2.set_xlabel("z/nm")
-    ax2.set_ylabel("c/(mol/L)")
-    ax2.set_ylim(0, 60)
-
-    ax2.plot(
-        X / JuliaMPBSolver.Units.nm,
-        c[1, :] / (JuliaMPBSolver.Units.M),
-        color = "blue",
-        linewidth = 2,
-        label = L"c^-",
-    )
-    ax2.plot(
-        X / JuliaMPBSolver.Units.nm,
-        c[2, :] / (JuliaMPBSolver.Units.M),
-        color = "red",
-        linewidth = 2,
-        label = L"c^+",
-    )
-    ax2.plot(
-        X / JuliaMPBSolver.Units.nm,
-        c0[1, :] / (JuliaMPBSolver.Units.M),
-        color = "green",
-        linewidth = 2,
-        label = L"c_{solvent}",
-    )
-    ax2.legend(loc = (0.4, 0.1))
-
-    tight_layout()
-    savefig("simplecell.jpg", dpi = 300)
-    return PythonPlot.gcf()
-end
-
-plotsol(solution)
+save("simplecell.jld2", "solution", solution, "X", X, "nv", nv, "ε_r", ε_r)

scripts/simplecell.jl

@@ -1,4 +1,9 @@
-using Pkg, JuliaMPBSolver, ExampleJuggler, Test, Markdown
+using Pkg
+using JuliaMPBSolver
+using ExampleJuggler
+using Test
+using Markdown
+using JLD2
 
 ExampleJuggler.verbose!(true)
 
@@ -11,8 +16,31 @@ notebooks = [
     "SymmetricCellSurfaceCharge.jl",
 ]
 
+scripts = [
+    "simplecell.jl",
+    "simplecell-bsk.jl",
+]
+
 @testset "Notebooks" begin
     @testscripts(joinpath(@__DIR__, "..", "notebooks"), notebooks)
 end
 
+@testset "Scripts" begin
+    @testscripts(joinpath(@__DIR__, "..", "scripts"), scripts)
+
+    for script in scripts
+        script_name = script[1:(end - 3)]
+        solution, X, nv, ε_r = load(script_name * ".jld2", "solution", "X", "nv", "ε_r")
+        solution_reference, X_reference, nv_reference, ε_r_reference = load(script_name * "-reference.jld2", "solution", "X", "nv", "ε_r")
+
+        @test solution ≈ solution_reference
+        @test X ≈ X_reference
+        @test nv ≈ nv_reference
+        @test ε_r ≈ ε_r_reference
+
+        rm(script_name * ".jld2")
+    end
+
+end
+
 Pkg.test("JuliaMPBSolver")