Unify parameter embedding

notebooks/HalfCellAppliedPotential.jl

@@ -1,5 +1,5 @@
 ### A Pluto.jl notebook ###
-# v0.20.19
+# v0.20.21
 
 using Markdown
 using InteractiveUtils

notebooks/ICMPBP-DD-Draft.jl

@@ -28,7 +28,7 @@ begin
     using LessUnitful
     using Test
     using PythonPlot
-    using JuliaMPBSolver.ICMPBP: ICMPBData, ICMPBSystem, L_Debye, set_molarity!, dlcap0, DerivedData, apply_charge!, ysum, qsweep, capscalc, calc_cmol, calc_c0mol, calc_χ, W
+    using JuliaMPBSolver.ICMPBP: ICMPBData, ICMPBSystem, L_Debye, set_molarity!, dlcap0, DerivedData, apply_charge!, ysum, qsweep, capscalc, calc_cmol, calc_c0mol, calc_χ, W, pramp
 end
 
 # ╔═╡ ef660f6f-9de3-4896-a65e-13c60df5de1e
@@ -246,12 +246,13 @@ end;
 # ╔═╡ 2ef6b8d7-e5f3-4700-8a40-8feffab3569f
 floataside(
     md"""
-    - ``L/nm``: $(@bind L0 PlutoUI.Slider(2:1: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))
-      ``M_{scale}/(mol/dm^3)``: $(@bind Mscale PlutoUI.Slider(5:5:60, default=60, show_value=true))
-    """, top = 100
+    - dielectric decrement: $(@bind dd PlutoUI.CheckBox())
+    -``M_{scale}/(mol/dm^3)``: $(@bind Mscale PlutoUI.Slider(5:5:60, default=60, show_value=true))
+        """, top = 100
 )
 
 
@@ -262,7 +263,7 @@ begin
     )
     set_molarity!(data1, M1_avg)
     data1.conserveions = true
-    data1.χvar = true
+    data1.χvar = dd
 end
 
 # ╔═╡ f8c1c2bd-7466-491e-9132-4f15edcfa4c7
@@ -351,13 +352,13 @@ function plotsol(
 end
 
 
-# ╔═╡ d2df6ed0-e6f5-4677-b790-bfc40de7fd6a
+# ╔═╡ 0eb0bbce-6f1c-4b7d-9df7-a375d30f5c99
 begin
     Q = surfcharge(n1_e)
     sol1 = inival1
-    for q in range(0, Q, length = 11)
+    pramp(; p = (0, Q), h = Q / 2, verbose = true) do q
         data1.q .= [-q, q]
-        global sol1 = solve(sys1; inival = sol1, verbose = "")
+        global sol1 = solve(sys1; inival = sol1, damp_initial = 0.1, max_round = 4, verbose = "")
     end
 end
 
@@ -372,6 +373,9 @@ plotsol(sol1, sys1; Mscale)
 # ╔═╡ 0e4ec7f0-0aa8-4a32-96a3-40f63f32a12d
 sol1
 
+# ╔═╡ 8d9d737c-0fbe-4664-afac-afa7d5e5fc4e
+Q
+
 # ╔═╡ 7d7ebb45-2fb3-40ea-83f2-62d0a240b2db
 @test isa(sol1, AbstractMatrix)
 
@@ -393,12 +397,13 @@ sol1
 # ╠═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
 # ╠═8433319f-2f78-494c-9b2e-a5390cf93b00
 # ╠═70910bd5-b8ca-4021-9b40-233b50ea5601
-# ╠═d2df6ed0-e6f5-4677-b790-bfc40de7fd6a
+# ╠═8d9d737c-0fbe-4664-afac-afa7d5e5fc4e
+# ╠═0eb0bbce-6f1c-4b7d-9df7-a375d30f5c99
 # ╠═7d7ebb45-2fb3-40ea-83f2-62d0a240b2db
 # ╠═f8c1c2bd-7466-491e-9132-4f15edcfa4c7
 # ╟─f75f1d3a-47e5-475b-97b1-bb275a510783

notebooks/SymmetricCellSurfaceCharge.jl

@@ -1,5 +1,5 @@
 ### A Pluto.jl notebook ###
-# v0.20.19
+# v0.20.21
 
 using Markdown
 using InteractiveUtils

packages/JuliaMPBSolver/src/JuliaMPBSolver.jl

@@ -25,12 +25,15 @@ module ICMPBP
     using VoronoiFVM
     using LinearAlgebra
     using SciMLBase
+
+    include("pramp.jl")
     include("icmbp-p.jl")
     include("cells.jl")
 end
 
 include("api.jl")
 
+export pramp
 export mpbpsolve
 export icmpbpsolve
 end

packages/JuliaMPBSolver/src/cells.jl

@@ -51,7 +51,7 @@ 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]
+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
@@ -99,7 +99,7 @@ end
 
 function dlcapsweep(
         cell::AppliedPotentialHalfCell; φ_max = 1.0, δφ = 1.0e-5,
-        steps = 51, hmin = 1.0e-5, damp_initial = 1, kwargs...
+        steps = 51, damp_initial = 1, kwargs...
     )
     set_φ!(cell, 0)
     sol0 = solve(cell; damp_initial)
@@ -108,36 +108,19 @@ function dlcapsweep(
     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
+        pramp(; p = (0, φ_max), h = hmax, hmax) do φ
+
+            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)
         end
         if dir == -1
             volts = reverse(volts)[1:(end - 1)]

packages/JuliaMPBSolver/src/pramp.jl

@@ -0,0 +1,66 @@
+"""
+     pramp(
+        f;
+        p = [0, 1],
+        hmin = (p[end] - p[begin]) / 1000,
+        hmax = (p[end] - p[begin]) / 10,
+        h = hmax,
+        hgrow = 1.2,
+        hdegrow = 0.5,
+        verbose = false
+    )
+
+Run 'f(p)'  successively with parameter values from the range given by 'p'.
+Stepsize is given by `h`. If `f` throws an error, solution is retried with a lower
+value of `h`, otherwise, `h` is increased unless it exceeded `hmax`.
+"""
+function pramp(
+        f;
+        p = [0, 1],
+        hmin = (p[end] - p[begin]) / 1000,
+        hmax = p[end] - p[begin],
+        h = hmax,
+        hgrow = 1.2,
+        hdegrow = 0.5,
+        verbose = false
+    )
+    pcurrent = p[begin]
+    first = true
+    h = min(h, hmax)
+    while pcurrent < p[end]
+        try
+            if first
+                ptrial = pcurrent
+            else
+                ptrial = min(pcurrent + h, p[end])
+            end
+            f(ptrial)
+            if verbose
+                if first
+                    println("pramp - success: p=$(pcurrent)")
+                else
+                    println("pramp - success: p=$(pcurrent) + $(h)")
+                end
+            end
+            pcurrent = ptrial
+            if first
+                first = false
+            else
+                h = min(h * hgrow, hmax)
+            end
+        catch e
+            if verbose
+                if first
+                    println("pramp - error: p=$(pcurrent)")
+                else
+                    println("pramp - error: p=$(pcurrent) + $(h)")
+                end
+            end
+            h = h * hdegrow
+            if h < hmin || first
+                rethrow(e)
+            end
+        end
+    end
+    return
+end

src/PoissonBoltzmannIPAM2025.jl

@@ -7,7 +7,7 @@ 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_φ!
+    set_κ!, set_q!, set_φ!, pramp
 
 resultsdir(args...) = projectdir("results", args...)
 draftresultsdir(args...) = projectdir("draftresults", args...)

src/plotcells.jl

@@ -142,9 +142,10 @@ function plotcells(
     fig.suptitle(title)
     grid = cell.sys.grid
     X = grid[Coordinates][1, :] / ufac"nm"
+
     solc = unknowns(cell)
-    for _Δφ in range(0, Δφ, length = 5)
-        set_φ!(cell, _Δφ)
+    pramp(p = (0, Δφ)) do φ
+        set_φ!(cell, φ)
         solc = solve(cell, inival = solc)
     end
 
@@ -157,26 +158,10 @@ function plotcells(
 
 
     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
+    pramp(p = (0, Δφ)) do φ
+        set_φ!(celldd, φ)
+        solv = solve(celldd, inival = solv, damp_initial = 0.5)
     end
-    #set_φ!(celldd, Δφ)
-    #solv = solve(celldd, inival = solc, damp_initial = 0.1)
 
     cv = calc_cmol(solv, celldd)
     c0v = calc_c0mol(solv, celldd)