more phase diagram fixes

Project.toml

@@ -1,7 +1,7 @@
 name = "ElectrochemicalKinetics"
 uuid = "a2c6e634-85ca-418a-9c67-9b5417ce2d04"
 authors = ["Rachel Kurchin <[email protected]>", "Holden Parks <[email protected]>", "Dhairya Gandhi <[email protected]>"]
-version = "0.1.4"
+version = "0.1.5"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"

README.md

@@ -71,6 +71,7 @@ models = [fit_model(exp_data, model_type; dos_file=dosfile) for model_type in mo
 plot_exp_and_models(exp_data, models)
 ```
 The resulting plot looks like:
+
 <img src="img/plot1.png" alt="models_expts">
 
 But we might also want to see the higher-voltage behavior and compare the models absent the experimental data. There's another function for that (by default it plots to +/- 1V)...
@@ -86,7 +87,7 @@ For more on this, see upcoming paper: [arxiv link placeholder]
 ### Butler-Volmer
 Probably the most basic (and largely empirical) kinetic model. The rate constants are given by:
 
-$$k_{\text{ox}}^{\text{BV}}(\eta)=A\exp\left(\frac{\alpha\eta}{k_{\text B}T}\right)\\k_{\text{red}}^{\text{BV}}(\eta)=A\exp\left(\frac{(1-\alpha)\eta}{k_{\text B}T}\right)$$
+$$k_{\text{ox}}^{\text{BV}}(\eta)=A\exp\left(\frac{\alpha\eta}{k_{\text B}T}\right)\\ k_{\text{red}}^{\text{BV}}(\eta)=A\exp\left(\frac{(1-\alpha)\eta}{k_{\text B}T}\right)$$
 
 It is implemented in the package as `ButlerVolmer`, possessing two fields: the prefactor `A` and the transfer coefficient `α`.
 

src/phase_diagrams.jl

@@ -43,6 +43,7 @@ function g_kinetic(I, km::KineticModel; intercalate=true, warn=true, T=room_T, k
     return g
 end
 
+# TODO: figure out T passthrough issue
 # zeros of this function correspond to pairs of x's satisfying the common tangent condition for a given µ function
 # case where we just feed in two points (x should be of length two)
 function common_tangent(x::Vector, I, km::KineticModel; intercalate=true, kwargs...)
@@ -72,18 +73,22 @@ NOTE 1: appropriate values of `I_step` depend strongly on the prefactor of your
 
 NOTE 2: at lower temperatures (<=320K or so), ButlerVolmer models with the default thermodynamic parameters have a two-phase region at every current, so setting a finite value of I_max is necessary for this function to finish running.
 """
-function phase_diagram(km::KineticModel; I_start=0.0, I_step=1.0, I_max=Inf, verbose=false, intercalate=true, start_guess=[0.05, 0.95], kwargs...)
+function phase_diagram(km::KineticModel; I_start=0.0, I_step=1.0, I_max=Inf, verbose=false, intercalate=true, start_guess=[0.05, 0.95], tol=5e-3, kwargs...)
     I = I_start
     pbs_here = find_phase_boundaries(I, km; intercalate=intercalate, guess=start_guess, kwargs...)
     pbs = pbs_here'
     I_vals = [I_start]
-    while abs(pbs_here[2] - pbs_here[1]) > 1e-3 && I < I_max
+    pb_diff = [0.0, 0.0]
+    while abs(pbs_here[2] - pbs_here[1]) > tol && I < I_max
         I = I + I_step
         if verbose
             println("Solving at I=", I, "...")
         end
         try
-            pbs_here = find_phase_boundaries(I, km; intercalate=intercalate, guess=pbs_here, kwargs...)
+            pbs_old = pbs_here
+            pbs_here = find_phase_boundaries(I, km; intercalate=intercalate, guess=pbs_old .+ pb_diff, kwargs...)
+            pb_diff = pbs_here .- pbs_old
+            # TODO: check that they haven't crossed over
             pbs = vcat(pbs, pbs_here')
             push!(I_vals, I)
             if verbose

test/phase_diagram_tests.jl

@@ -117,24 +117,60 @@ end
 end
 
 @testset "Phase Diagram" begin
-    km = bv # TODO: expand this
-
-    # simplest case, just one pair of x values (this function is still pretty slow though)
-    v1 = find_phase_boundaries(100, km)
-
-    @test all(isapprox.(common_tangent(v1, 100, km), Ref(0.0), atol=1e-6))
-    v2 = find_phase_boundaries(100, km, T=350)     
-    @test all(isapprox.(common_tangent(v2, 100, km, T=350), Ref(0.0), atol=1e-5))
-    # they should get "narrower" with temperature
-    @test v2[1] > v1[1]
-    @test v2[2] < v1[2]
-    v3 = find_phase_boundaries(400, km)
-    # ...and also with current
-    @test v3[1] > v1[1]
-    @test v3[2] < v1[2]
-
-    # test actual numerical values too
-    @test all(isapprox.(v1, [0.045547, 0.841501], atol=1e-5))
-    @test all(isapprox.(v2, [0.090478, 0.77212], atol=1e-4))
-    @test all(isapprox.(v3, [0.105747, 0.6809], atol=1e-4))
+    v1_vals = Dict(:ButlerVolmer=>[0.045547, 0.841501],
+                :Marcus=>[0.04901, 0.81884],
+                :AsymptoticMarcusHushChidsey=>[0.04958, 0.81897])
+    v2_vals = Dict(:ButlerVolmer=>[0.090478, 0.77212],
+                :Marcus=>[0.078934, 0.804275],
+                :AsymptoticMarcusHushChidsey=>[0.07609, 0.81652])
+    v3_vals = Dict(:ButlerVolmer=>[0.105747, 0.6809],
+                :Marcus=>[0.154934, 0.54218],
+                :AsymptoticMarcusHushChidsey=>[0.147034,0.566617])
+    I_max_T300 = Dict(:ButlerVolmer=>700.0,
+                :Marcus=>600.0,
+                :AsymptoticMarcusHushChidsey=>650.0)
+    I_max_T400 = Dict(:ButlerVolmer=>250.0,
+                :Marcus=>400.0,
+                :AsymptoticMarcusHushChidsey=>500.0)
+    for km in kms
+        model_type = nameof(typeof(km))
+        @testset "$model_type" begin
+            # simplest case, just one pair of x values (this function is still pretty slow though)
+            v1 = find_phase_boundaries(100, km)
+
+            @test all(isapprox.(common_tangent(v1, 100, km), Ref(0.0), atol=1e-6))
+            v2 = find_phase_boundaries(100, km, T=350)     
+            @test all(isapprox.(common_tangent(v2, 100, km, T=350), Ref(0.0), atol=1e-5))
+            # they should get "narrower" with temperature
+            @test v2[1] > v1[1]
+            @test v2[2] < v1[2]
+            v3 = find_phase_boundaries(400, km, guess=[0.1,0.8])
+            # ...and also with current
+            @test v3[1] > v1[1]
+            @test v3[2] < v1[2]
+
+            # test actual numerical values too
+            @test all(isapprox.(v1, v1_vals[model_type], atol=1e-5))
+            @test all(isapprox.(v2, v2_vals[model_type], atol=1e-4))
+            @test all(isapprox.(v3, v3_vals[model_type], atol=1e-4))
+
+            # and actually build the full map for a couple temps
+            pbs, I = phase_diagram(km, I_max=700, I_step=50)
+            @test maximum(I) == I_max_T300[model_type]
+        end
+    end
+
+    @testset "MHC" begin
+        # do some lighterweight stuff here since this takes awhile
+        mhc = MarcusHushChidsey(amhc.A, amhc.λ)
+        # also give it a good initial guess so it only takes one or two optimizer steps
+        v1 = find_phase_boundaries(100, mhc, guess=v1_vals[:AsymptoticMarcusHushChidsey])
+        @test all(isapprox.(v1, [0.04967204036, 0.81822937543], atol=1e-5))
+
+        v2 = find_phase_boundaries(100, mhc, T=350, guess=v2_vals[:AsymptoticMarcusHushChidsey])
+        @test all(isapprox.(v2, [0.075615, 0.8171636], atol=1e-5))
+
+        v3 = find_phase_boundaries(400, mhc, guess=v3_vals[:AsymptoticMarcusHushChidsey])
+        @test all(isapprox.(v3, [0.1467487, 0.5704125], atol=1e-5))
+    end
 end