diff --git a/examples/Data/ParameterAnalysis/FractionsHeightLeft.npz b/examples/Data/ParameterAnalysis/FractionsHeightLeft.npz
index d61941b95f8becc100880ad255c7d1fa80f09ede..2153287dd6e8839bbdd97f01f77e3222d4ef5707 100644
Binary files a/examples/Data/ParameterAnalysis/FractionsHeightLeft.npz and b/examples/Data/ParameterAnalysis/FractionsHeightLeft.npz differ
diff --git a/examples/Figures/ParameterAnalysis/FractionsHeightLeft.svg b/examples/Figures/ParameterAnalysis/FractionsHeightLeft.svg
index b4999bd80d4c822961bf8d919deabcaee15e3ba8..1e48149d9fe265fa9d5c6a784f1e85baa9d07635 100644
--- a/examples/Figures/ParameterAnalysis/FractionsHeightLeft.svg
+++ b/examples/Figures/ParameterAnalysis/FractionsHeightLeft.svg
@@ -6,7 +6,7 @@
<rdf:RDF xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:cc="http://creativecommons.org/ns#" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
<cc:Work>
<dc:type rdf:resource="http://purl.org/dc/dcmitype/StillImage"/>
- <dc:date>2024-08-19T12:51:22.628209</dc:date>
+ <dc:date>2024-08-20T15:15:16.805709</dc:date>
<dc:format>image/svg+xml</dc:format>
<dc:creator>
<cc:Agent>
@@ -42,16 +42,16 @@ z
<g id="line2d_1">
<path d="M 139.043011 701.92
L 139.043011 75
-" clip-path="url(#p1871188695)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
+" clip-path="url(#pb9693e132a)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
</g>
<g id="line2d_2">
<defs>
- <path id="mc368a271b7" d="M 0 0
+ <path id="m547502dcbc" d="M 0 0
L 0 3.5
" style="stroke: #000000; stroke-width: 0.8"/>
</defs>
<g>
- <use xlink:href="#mc368a271b7" x="139.043011" y="701.92" style="stroke: #000000; stroke-width: 0.8"/>
+ <use xlink:href="#m547502dcbc" x="139.043011" y="701.92" style="stroke: #000000; stroke-width: 0.8"/>
</g>
</g>
<g id="text_1">
@@ -120,11 +120,11 @@ z
<g id="line2d_3">
<path d="M 249.127102 701.92
L 249.127102 75
-" clip-path="url(#p1871188695)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
+" clip-path="url(#pb9693e132a)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
</g>
<g id="line2d_4">
<g>
- <use xlink:href="#mc368a271b7" x="249.127102" y="701.92" style="stroke: #000000; stroke-width: 0.8"/>
+ <use xlink:href="#m547502dcbc" x="249.127102" y="701.92" style="stroke: #000000; stroke-width: 0.8"/>
</g>
</g>
<g id="text_2">
@@ -178,11 +178,11 @@ z
<g id="line2d_5">
<path d="M 359.211193 701.92
L 359.211193 75
-" clip-path="url(#p1871188695)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
+" clip-path="url(#pb9693e132a)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
</g>
<g id="line2d_6">
<g>
- <use xlink:href="#mc368a271b7" x="359.211193" y="701.92" style="stroke: #000000; stroke-width: 0.8"/>
+ <use xlink:href="#m547502dcbc" x="359.211193" y="701.92" style="stroke: #000000; stroke-width: 0.8"/>
</g>
</g>
<g id="text_3">
@@ -199,11 +199,11 @@ L 359.211193 75
<g id="line2d_7">
<path d="M 469.295284 701.92
L 469.295284 75
-" clip-path="url(#p1871188695)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
+" clip-path="url(#pb9693e132a)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
</g>
<g id="line2d_8">
<g>
- <use xlink:href="#mc368a271b7" x="469.295284" y="701.92" style="stroke: #000000; stroke-width: 0.8"/>
+ <use xlink:href="#m547502dcbc" x="469.295284" y="701.92" style="stroke: #000000; stroke-width: 0.8"/>
</g>
</g>
<g id="text_4">
@@ -246,11 +246,11 @@ z
<g id="line2d_9">
<path d="M 579.379375 701.92
L 579.379375 75
-" clip-path="url(#p1871188695)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
+" clip-path="url(#pb9693e132a)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
</g>
<g id="line2d_10">
<g>
- <use xlink:href="#mc368a271b7" x="579.379375" y="701.92" style="stroke: #000000; stroke-width: 0.8"/>
+ <use xlink:href="#m547502dcbc" x="579.379375" y="701.92" style="stroke: #000000; stroke-width: 0.8"/>
</g>
</g>
<g id="text_5">
@@ -266,11 +266,11 @@ L 579.379375 75
<g id="line2d_11">
<path d="M 689.463466 701.92
L 689.463466 75
-" clip-path="url(#p1871188695)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
+" clip-path="url(#pb9693e132a)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
</g>
<g id="line2d_12">
<g>
- <use xlink:href="#mc368a271b7" x="689.463466" y="701.92" style="stroke: #000000; stroke-width: 0.8"/>
+ <use xlink:href="#m547502dcbc" x="689.463466" y="701.92" style="stroke: #000000; stroke-width: 0.8"/>
</g>
</g>
<g id="text_6">
@@ -286,11 +286,11 @@ L 689.463466 75
<g id="line2d_13">
<path d="M 799.547557 701.92
L 799.547557 75
-" clip-path="url(#p1871188695)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
+" clip-path="url(#pb9693e132a)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
</g>
<g id="line2d_14">
<g>
- <use xlink:href="#mc368a271b7" x="799.547557" y="701.92" style="stroke: #000000; stroke-width: 0.8"/>
+ <use xlink:href="#m547502dcbc" x="799.547557" y="701.92" style="stroke: #000000; stroke-width: 0.8"/>
</g>
</g>
<g id="text_7">
@@ -306,11 +306,11 @@ L 799.547557 75
<g id="line2d_15">
<path d="M 909.631648 701.92
L 909.631648 75
-" clip-path="url(#p1871188695)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
+" clip-path="url(#pb9693e132a)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
</g>
<g id="line2d_16">
<g>
- <use xlink:href="#mc368a271b7" x="909.631648" y="701.92" style="stroke: #000000; stroke-width: 0.8"/>
+ <use xlink:href="#m547502dcbc" x="909.631648" y="701.92" style="stroke: #000000; stroke-width: 0.8"/>
</g>
</g>
<g id="text_8">
@@ -326,11 +326,11 @@ L 909.631648 75
<g id="line2d_17">
<path d="M 1019.715739 701.92
L 1019.715739 75
-" clip-path="url(#p1871188695)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
+" clip-path="url(#pb9693e132a)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
</g>
<g id="line2d_18">
<g>
- <use xlink:href="#mc368a271b7" x="1019.715739" y="701.92" style="stroke: #000000; stroke-width: 0.8"/>
+ <use xlink:href="#m547502dcbc" x="1019.715739" y="701.92" style="stroke: #000000; stroke-width: 0.8"/>
</g>
</g>
<g id="text_9">
@@ -450,16 +450,16 @@ z
<g id="line2d_19">
<path d="M 95.009375 673.423638
L 1063.749375 673.423638
-" clip-path="url(#p1871188695)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
+" clip-path="url(#pb9693e132a)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
</g>
<g id="line2d_20">
<defs>
- <path id="m294d782db7" d="M 0 0
+ <path id="m2c2a3b9dc2" d="M 0 0
L -3.5 0
" style="stroke: #000000; stroke-width: 0.8"/>
</defs>
<g>
- <use xlink:href="#m294d782db7" x="95.009375" y="673.423638" style="stroke: #000000; stroke-width: 0.8"/>
+ <use xlink:href="#m2c2a3b9dc2" x="95.009375" y="673.423638" style="stroke: #000000; stroke-width: 0.8"/>
</g>
</g>
<g id="text_11">
@@ -473,18 +473,18 @@ L -3.5 0
</g>
<g id="ytick_2">
<g id="line2d_21">
- <path d="M 95.009375 559.433042
-L 1063.749375 559.433042
-" clip-path="url(#p1871188695)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
+ <path d="M 95.009375 559.433008
+L 1063.749375 559.433008
+" clip-path="url(#pb9693e132a)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
</g>
<g id="line2d_22">
<g>
- <use xlink:href="#m294d782db7" x="95.009375" y="559.433042" style="stroke: #000000; stroke-width: 0.8"/>
+ <use xlink:href="#m2c2a3b9dc2" x="95.009375" y="559.433008" style="stroke: #000000; stroke-width: 0.8"/>
</g>
</g>
<g id="text_12">
<!-- 0.2 -->
- <g transform="translate(40.3 570.830698) scale(0.3 -0.3)">
+ <g transform="translate(40.3 570.830664) scale(0.3 -0.3)">
<use xlink:href="#DejaVuSans-30"/>
<use xlink:href="#DejaVuSans-2e" x="63.623047"/>
<use xlink:href="#DejaVuSans-32" x="95.410156"/>
@@ -493,18 +493,18 @@ L 1063.749375 559.433042
</g>
<g id="ytick_3">
<g id="line2d_23">
- <path d="M 95.009375 445.442446
-L 1063.749375 445.442446
-" clip-path="url(#p1871188695)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
+ <path d="M 95.009375 445.442378
+L 1063.749375 445.442378
+" clip-path="url(#pb9693e132a)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
</g>
<g id="line2d_24">
<g>
- <use xlink:href="#m294d782db7" x="95.009375" y="445.442446" style="stroke: #000000; stroke-width: 0.8"/>
+ <use xlink:href="#m2c2a3b9dc2" x="95.009375" y="445.442378" style="stroke: #000000; stroke-width: 0.8"/>
</g>
</g>
<g id="text_13">
<!-- 0.4 -->
- <g transform="translate(40.3 456.840103) scale(0.3 -0.3)">
+ <g transform="translate(40.3 456.840035) scale(0.3 -0.3)">
<defs>
<path id="DejaVuSans-34" d="M 2419 4116
L 825 1625
@@ -534,18 +534,18 @@ z
</g>
<g id="ytick_4">
<g id="line2d_25">
- <path d="M 95.009375 331.451851
-L 1063.749375 331.451851
-" clip-path="url(#p1871188695)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
+ <path d="M 95.009375 331.451749
+L 1063.749375 331.451749
+" clip-path="url(#pb9693e132a)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
</g>
<g id="line2d_26">
<g>
- <use xlink:href="#m294d782db7" x="95.009375" y="331.451851" style="stroke: #000000; stroke-width: 0.8"/>
+ <use xlink:href="#m2c2a3b9dc2" x="95.009375" y="331.451749" style="stroke: #000000; stroke-width: 0.8"/>
</g>
</g>
<g id="text_14">
<!-- 0.6 -->
- <g transform="translate(40.3 342.849507) scale(0.3 -0.3)">
+ <g transform="translate(40.3 342.849405) scale(0.3 -0.3)">
<defs>
<path id="DejaVuSans-36" d="M 2113 2584
Q 1688 2584 1439 2293
@@ -586,18 +586,18 @@ z
</g>
<g id="ytick_5">
<g id="line2d_27">
- <path d="M 95.009375 217.461255
-L 1063.749375 217.461255
-" clip-path="url(#p1871188695)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
+ <path d="M 95.009375 217.461119
+L 1063.749375 217.461119
+" clip-path="url(#pb9693e132a)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
</g>
<g id="line2d_28">
<g>
- <use xlink:href="#m294d782db7" x="95.009375" y="217.461255" style="stroke: #000000; stroke-width: 0.8"/>
+ <use xlink:href="#m2c2a3b9dc2" x="95.009375" y="217.461119" style="stroke: #000000; stroke-width: 0.8"/>
</g>
</g>
<g id="text_15">
<!-- 0.8 -->
- <g transform="translate(40.3 228.858911) scale(0.3 -0.3)">
+ <g transform="translate(40.3 228.858775) scale(0.3 -0.3)">
<defs>
<path id="DejaVuSans-38" d="M 2034 2216
Q 1584 2216 1326 1975
@@ -647,18 +647,18 @@ z
</g>
<g id="ytick_6">
<g id="line2d_29">
- <path d="M 95.009375 103.47066
-L 1063.749375 103.47066
-" clip-path="url(#p1871188695)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
+ <path d="M 95.009375 103.470489
+L 1063.749375 103.470489
+" clip-path="url(#pb9693e132a)" style="fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square"/>
</g>
<g id="line2d_30">
<g>
- <use xlink:href="#m294d782db7" x="95.009375" y="103.47066" style="stroke: #000000; stroke-width: 0.8"/>
+ <use xlink:href="#m2c2a3b9dc2" x="95.009375" y="103.470489" style="stroke: #000000; stroke-width: 0.8"/>
</g>
</g>
<g id="text_16">
<!-- 1.0 -->
- <g transform="translate(40.3 114.868316) scale(0.3 -0.3)">
+ <g transform="translate(40.3 114.868146) scale(0.3 -0.3)">
<use xlink:href="#DejaVuSans-31"/>
<use xlink:href="#DejaVuSans-2e" x="63.623047"/>
<use xlink:href="#DejaVuSans-30" x="95.410156"/>
@@ -729,21 +729,279 @@ z
</g>
<g id="line2d_31">
<path d="M 139.043011 673.423636
-L 579.379375 483.439312
+L 390.034739 673.320129
+L 416.45492 673.08386
+L 434.068375 672.676851
+L 447.278466 672.081037
+L 456.085193 671.443747
+L 464.89192 670.511791
+L 469.295284 669.896604
+L 473.698648 669.155399
+L 478.102011 668.263752
+L 482.505375 667.192998
+L 486.908739 665.90967
+L 491.312102 664.374924
+L 495.715466 662.543985
+L 500.11883 660.365657
+L 504.522193 657.781932
+L 508.925557 654.727798
+L 513.32892 651.131302
+L 517.732284 646.913998
+L 522.135648 641.991879
+L 526.539011 636.276908
+L 530.942375 629.679274
+L 535.345739 622.110428
+L 539.749102 613.486931
+L 544.152466 603.735059
+L 548.55583 592.795972
+L 552.959193 580.631179
+L 557.362557 567.227845
+L 561.76592 552.603418
+L 566.169284 536.809014
+L 570.572648 519.930981
+L 574.976011 502.090247
+L 579.379375 483.439255
+L 588.186102 444.439501
+L 601.396193 384.753544
+L 605.799557 365.341827
+L 610.20292 346.462331
+L 614.606284 328.253846
+L 619.009648 310.827695
+L 623.413011 294.267451
+L 627.816375 278.630166
+L 632.219739 263.948683
+L 636.623102 250.234669
+L 641.026466 237.481983
+L 645.42983 225.670104
+L 649.833193 214.767385
+L 654.236557 204.734016
+L 658.63992 195.524595
+L 663.043284 187.090289
+L 667.446648 179.380594
+L 671.850011 172.344723
+L 676.253375 165.93266
+L 680.656739 160.095942
+L 685.060102 154.788197
+L 689.463466 149.965504
+L 693.86683 145.586601
+L 698.270193 141.612982
+L 702.673557 138.008902
+L 707.07692 134.741333
+L 711.480284 131.779864
+L 715.883648 129.096582
+L 720.287011 126.665928
+L 724.690375 124.464547
+L 729.093739 122.471136
+L 733.497102 120.666289
+L 737.900466 119.03235
+L 742.30383 117.553271
+L 746.707193 116.214475
+L 755.51392 113.906044
+L 764.320648 112.015313
+L 773.127375 110.466901
+L 781.934102 109.198941
+L 790.74083 108.160698
+L 799.547557 107.310588
+L 812.757648 106.315357
+L 825.967739 105.578039
+L 843.581193 104.883229
+L 861.194648 104.417477
+L 883.211466 104.044865
+L 914.035011 103.755714
+L 962.472011 103.565431
L 1019.715739 103.496364
-" clip-path="url(#p1871188695)" style="fill: none; stroke: #1f77b4; stroke-width: 4; stroke-linecap: square"/>
+L 1019.715739 103.496364
+" clip-path="url(#pb9693e132a)" style="fill: none; stroke: #1f77b4; stroke-width: 4; stroke-linecap: square"/>
</g>
<g id="line2d_32">
<path d="M 139.043011 103.496364
-L 579.379375 483.439312
+L 227.110284 103.661681
+L 266.740557 103.940748
+L 293.160739 104.327359
+L 315.177557 104.883229
+L 332.791011 105.578039
+L 346.001102 106.315357
+L 359.211193 107.310588
+L 368.01792 108.160698
+L 376.824648 109.198941
+L 385.631375 110.466901
+L 394.438102 112.015313
+L 403.24483 113.906044
+L 407.648193 115.002733
+L 412.051557 116.214475
+L 416.45492 117.553271
+L 420.858284 119.03235
+L 425.261648 120.666289
+L 429.665011 122.471136
+L 434.068375 124.464547
+L 438.471739 126.665928
+L 442.875102 129.096582
+L 447.278466 131.779864
+L 451.68183 134.741333
+L 456.085193 138.008902
+L 460.488557 141.612982
+L 464.89192 145.586601
+L 469.295284 149.965504
+L 473.698648 154.788197
+L 478.102011 160.095942
+L 482.505375 165.93266
+L 486.908739 172.344723
+L 491.312102 179.380594
+L 495.715466 187.090289
+L 500.11883 195.524595
+L 504.522193 204.734016
+L 508.925557 214.767385
+L 513.32892 225.670104
+L 517.732284 237.481983
+L 522.135648 250.234669
+L 526.539011 263.948683
+L 530.942375 278.630166
+L 535.345739 294.267451
+L 539.749102 310.827695
+L 544.152466 328.253846
+L 548.55583 346.462331
+L 552.959193 365.341827
+L 561.76592 404.533295
+L 574.976011 464.156561
+L 579.379375 483.439255
+L 583.782739 502.090247
+L 588.186102 519.930981
+L 592.589466 536.809014
+L 596.99283 552.603418
+L 601.396193 567.227845
+L 605.799557 580.631179
+L 610.20292 592.795972
+L 614.606284 603.735059
+L 619.009648 613.486931
+L 623.413011 622.110428
+L 627.816375 629.679274
+L 632.219739 636.276908
+L 636.623102 641.991879
+L 641.026466 646.913998
+L 645.42983 651.131302
+L 649.833193 654.727798
+L 654.236557 657.781932
+L 658.63992 660.365657
+L 663.043284 662.543985
+L 667.446648 664.374924
+L 671.850011 665.90967
+L 676.253375 667.192998
+L 680.656739 668.263752
+L 685.060102 669.155399
+L 689.463466 669.896604
+L 698.270193 671.021672
+L 707.07692 671.792747
+L 715.883648 672.318964
+L 729.093739 672.81
+L 746.707193 673.144781
+L 773.127375 673.338781
+L 825.967739 673.415873
+L 1019.715739 673.423636
L 1019.715739 673.423636
-" clip-path="url(#p1871188695)" style="fill: none; stroke: #ff7f0e; stroke-width: 4; stroke-linecap: square"/>
+" clip-path="url(#pb9693e132a)" style="fill: none; stroke: #ff7f0e; stroke-width: 4; stroke-linecap: square"/>
</g>
<g id="line2d_33">
- <path d="M 139.043011 673.397935
-L 579.379375 483.439312
-L 1019.715739 673.397935
-" clip-path="url(#p1871188695)" style="fill: none; stroke: #2ca02c; stroke-width: 4; stroke-linecap: square"/>
+ <path d="M 139.043011 673.397764
+L 227.110284 673.23251
+L 266.740557 672.953767
+L 293.160739 672.568054
+L 315.177557 672.014391
+L 332.791011 671.323853
+L 346.001102 670.5929
+L 359.211193 669.60924
+L 368.01792 668.771713
+L 376.824648 667.752195
+L 385.631375 666.512082
+L 394.438102 665.005056
+L 403.24483 663.175781
+L 412.051557 660.958508
+L 416.45492 659.680634
+L 420.858284 658.275678
+L 425.261648 656.731887
+L 429.665011 655.036628
+L 434.068375 653.176366
+L 438.471739 651.136668
+L 442.875102 648.902219
+L 447.278466 646.456863
+L 451.68183 643.783684
+L 456.085193 640.865116
+L 460.488557 637.683111
+L 464.89192 634.219372
+L 469.295284 630.455656
+L 473.698648 626.374169
+L 478.102011 621.958071
+L 482.505375 617.192106
+L 486.908739 612.063372
+L 491.312102 606.562247
+L 495.715466 600.68349
+L 500.11883 594.427512
+L 504.522193 587.801816
+L 508.925557 580.822582
+L 513.32892 573.516359
+L 517.732284 565.921783
+L 526.539011 550.092174
+L 539.749102 526.003138
+L 544.152466 518.328859
+L 548.55583 511.059462
+L 552.959193 504.344758
+L 557.362557 498.336376
+L 561.76592 493.181051
+L 566.169284 489.013186
+L 570.572648 485.947283
+L 574.976011 484.070957
+L 579.379375 483.439255
+L 583.782739 484.070957
+L 588.186102 485.947283
+L 592.589466 489.013186
+L 596.99283 493.181051
+L 601.396193 498.336376
+L 605.799557 504.344758
+L 610.20292 511.059462
+L 614.606284 518.328859
+L 619.009648 526.003138
+L 627.816375 542.008325
+L 636.623102 558.091217
+L 641.026466 565.921783
+L 645.42983 573.516359
+L 649.833193 580.822582
+L 654.236557 587.801816
+L 658.63992 594.427512
+L 663.043284 600.68349
+L 667.446648 606.562247
+L 671.850011 612.063372
+L 676.253375 617.192106
+L 680.656739 621.958071
+L 685.060102 626.374169
+L 689.463466 630.455656
+L 693.86683 634.219372
+L 698.270193 637.683111
+L 702.673557 640.865116
+L 707.07692 643.783684
+L 711.480284 646.456863
+L 715.883648 648.902219
+L 720.287011 651.136668
+L 724.690375 653.176366
+L 729.093739 655.036628
+L 733.497102 656.731887
+L 737.900466 658.275678
+L 742.30383 659.680634
+L 746.707193 660.958508
+L 755.51392 663.175781
+L 764.320648 665.005056
+L 773.127375 666.512082
+L 781.934102 667.752195
+L 790.74083 668.771713
+L 799.547557 669.60924
+L 812.757648 670.5929
+L 825.967739 671.323853
+L 843.581193 672.014391
+L 861.194648 672.47822
+L 883.211466 672.84984
+L 914.035011 673.138555
+L 962.472011 673.328711
+L 1019.715739 673.397764
+L 1019.715739 673.397764
+" clip-path="url(#pb9693e132a)" style="fill: none; stroke: #2ca02c; stroke-width: 4; stroke-linecap: square"/>
</g>
<g id="patch_3">
<path d="M 95.009375 701.92
@@ -900,7 +1158,7 @@ z
</g>
</g>
<defs>
- <clipPath id="p1871188695">
+ <clipPath id="pb9693e132a">
<rect x="95.009375" y="75" width="968.74" height="626.92"/>
</clipPath>
</defs>
diff --git a/examples/ReproducableCode/DoubleLayerCapacity/CompressibleValidationAnalyticalMethod.py b/examples/ReproducableCode/DoubleLayerCapacity/CompressibleValidationAnalyticalMethod.py
index cf5abb92597d4273da9257e3de8a446285e43617..db4b3e5726b0c28bfde61bd083054754a35c3809 100644
--- a/examples/ReproducableCode/DoubleLayerCapacity/CompressibleValidationAnalyticalMethod.py
+++ b/examples/ReproducableCode/DoubleLayerCapacity/CompressibleValidationAnalyticalMethod.py
@@ -15,6 +15,7 @@ import os
src_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '../../../', 'src')
sys.path.insert(0, src_path)
+from Eq04 import solve_System_4eq
from Eq02 import solve_System_2eq
from Helper_DoubleLayerCapacity import Phi_pot_center, dx, C_dl, n, Q_num_, Q_DL_dimless_ana, Q_DL_dim_ana
@@ -78,7 +79,7 @@ chi = 80 # [-]
# Parameter and bcs for the electrolyte
Lambda2 = (k*T*epsilon0*(1+chi))/(e0**2 * nR_m * (LR)**2)
a2 = (pR)/(nR_m * k * T)
-K = 20_000#'incompressible'
+K = 1_000#'incompressible'
kappa = 0
Molarity = 0.01
y_R = Molarity / nR_mol
@@ -102,26 +103,26 @@ phi_left_dimless = np.linspace(Vol_start, Volt_end, n_Volts) * e0/(k*T)
# Numerical calculations
# Solution vectors
-y_A_num, y_C_num, y_S_num, phi_num, p_num, x_num = [], [], [], [], [], []
-for i, phi_bcs in enumerate(phi_left_dimless):
- y_A_, y_C_, phi_, p_, x_ = solve_System_2eq(phi_bcs, phi_right, p_right, z_A, z_C, y_R, y_R, K, Lambda2, a2, number_cells, solvation=kappa, refinement_style='hard_log', rtol=1e-4, max_iter=2500, return_type='Vector', relax_param=relax_param)
- y_S_ = 1 - y_A_ - y_C_
- y_A_num.append(y_A_)
- y_C_num.append(y_C_)
- y_S_num.append(y_S_)
- phi_num.append(phi_)
- p_num.append(p_)
- x_num.append(x_)
+# y_A_num, y_C_num, y_S_num, phi_num, p_num, x_num = [], [], [], [], [], []
+# for i, phi_bcs in enumerate(phi_left_dimless):
+# y_A_, y_C_, phi_, p_, x_ = solve_System_4eq(phi_bcs, phi_right, p_right, z_A, z_C, y_R, y_R, K, Lambda2, a2, number_cells, solvation=kappa, refinement_style='hard_log', rtol=1e-4, max_iter=2500, return_type='Vector', relax_param=relax_param)
+# y_S_ = 1 - y_A_ - y_C_
+# y_A_num.append(y_A_)
+# y_C_num.append(y_C_)
+# y_S_num.append(y_S_)
+# phi_num.append(phi_)
+# p_num.append(p_)
+# x_num.append(x_)
-Q_num = []
-for j in range(len(phi_left_dimless)):
- Q_num.append(Q_num_(y_A_num[j], y_C_num[j], n(p_num[j], K), x_num[j]))
-Q_num = np.array(Q_num)
+# Q_num = []
+# for j in range(len(phi_left_dimless)):
+# Q_num.append(Q_num_(y_A_num[j], y_C_num[j], n(p_num[j], K), x_num[j]))
+# Q_num = np.array(Q_num)
-dx_ = phi_left_dimless[1] - phi_left_dimless[0] # [1/V], Assumption: phi^L is uniformly distributed
-C_DL_num = (Q_num[1:] - Q_num[:-1])/dx_ # [µAs/cm³]
-C_DL_num = np.array(C_DL_num)
-C_dl_num = C_dl(Q_num, phi_left_dimless)
+# dx_ = phi_left_dimless[1] - phi_left_dimless[0] # [1/V], Assumption: phi^L is uniformly distributed
+# C_DL_num = (Q_num[1:] - Q_num[:-1])/dx_ # [µAs/cm³]
+# C_DL_num = np.array(C_DL_num)
+# C_dl_num = C_dl(Q_num, phi_left_dimless)
# Analytical calculations
# K = K_vec[0]
@@ -136,7 +137,7 @@ plt.figure()
# plt.plot(phi_left, Q_num[1] - Q_ana, label='Difference')
# plt.plot(phi_left_dimless, Q_num - Q_ana, label='Difference')
# plt.plot(phi_left_dimless, Q_num[1], label='Numerical')
-plt.plot(phi_left_dimless, Q_num, label='Numerical')
+# plt.plot(phi_left_dimless, Q_num, label='Numerical')
plt.plot(phi_left_dimless, Q_ana, label='Analytical')
plt.grid()
plt.legend()
@@ -149,7 +150,7 @@ plt.figure()
# plt.plot(Phi_pot_center(phi_left), C_DL_num[1] - C_DL_ana, label='Difference')
# plt.plot(Phi_pot_center(phi_left_dimless), C_DL_num - C_DL_ana, label='Difference')
# plt.plot(Phi_pot_center(phi_left_dimless), C_DL_num[1], label='Numerical')
-plt.plot(Phi_pot_center(phi_left_dimless), C_DL_num, label='Numerical')
+# plt.plot(Phi_pot_center(phi_left_dimless), C_DL_num, label='Numerical')
plt.plot(Phi_pot_center(phi_left_dimless), C_DL_ana, label='Analytical')
plt.grid()
plt.legend()
diff --git a/examples/ReproducableCode/DoubleLayerCapacity/IncompressibleValidationAnalyticalMethod.py b/examples/ReproducableCode/DoubleLayerCapacity/IncompressibleValidationAnalyticalMethod.py
index c17780a490e4e75201db5ab285f5e42d77e07d05..f23debe52dda6274fbd5bbe9e9b560982c320b66 100644
--- a/examples/ReproducableCode/DoubleLayerCapacity/IncompressibleValidationAnalyticalMethod.py
+++ b/examples/ReproducableCode/DoubleLayerCapacity/IncompressibleValidationAnalyticalMethod.py
@@ -15,6 +15,7 @@ import os
src_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '../../../', 'src')
sys.path.insert(0, src_path)
+from Eq04 import solve_System_4eq
from Eq02 import solve_System_2eq
from Helper_DoubleLayerCapacity import Phi_pot_center, dx, C_dl, n, Q_num_, Q_DL_dimless_ana, Q_DL_dim_ana
@@ -43,7 +44,7 @@ chi = 80 # [-]
Lambda2 = (k*T*epsilon0*(1+chi))/(e0**2 * nR_m * (LR)**2)
a2 = (pR)/(nR_m * k * T)
K = 'incompressible'
-kappa = 5
+kappa = 0
Molarity = 0.01
y_R = Molarity / nR_mol
z_A, z_C = -1.0, 1.0
@@ -60,7 +61,7 @@ rtol = 1e-4 # ! Change back to 1e-8
# phi^L domain
Vol_start = 0.1 # ! Change back to 0
Volt_end = 0.75
-n_Volts = 25#0
+n_Volts = 5#0
phi_left = np.linspace(Vol_start, Volt_end, n_Volts) * e0/(k*T)
@@ -69,7 +70,7 @@ phi_left = np.linspace(Vol_start, Volt_end, n_Volts) * e0/(k*T)
# Solution vectors
y_A_num, y_C_num, y_S_num, phi_num, p_num, x_num = [], [], [], [], [], []
for i, phi_bcs in enumerate(phi_left):
- y_A_, y_C_, phi_, p_, x_ = solve_System_2eq(phi_bcs, phi_right, p_right, z_A, z_C, y_R, y_R, K, Lambda2, a2, number_cells, solvation=kappa, refinement_style='hard_log', rtol=rtol, max_iter=2500, return_type='Vector', relax_param=relax_param)
+ y_A_, y_C_, phi_, p_, x_ = solve_System_4eq(phi_bcs, phi_right, p_right, z_A, z_C, y_R, y_R, K, Lambda2, a2, number_cells, solvation=kappa, refinement_style='hard_log', rtol=rtol, max_iter=2500, return_type='Vector', relax_param=relax_param)
y_S_ = 1 - y_A_ - y_C_
y_A_num.append(y_A_)
y_C_num.append(y_C_)
@@ -98,7 +99,7 @@ C_DL_ana = C_dl(Q_ana, phi_left)
# Plotting
plt.figure()
# plt.plot(phi_left, Q_num - Q_ana, label='Difference')
-plt.plot(phi_left, Q_num, label='Numerical')
+# plt.plot(phi_left, Q_num, label='Numerical')
plt.plot(phi_left, Q_ana, label='Analytical')
plt.grid()
plt.legend()
@@ -109,7 +110,7 @@ plt.show()
plt.figure()
# plt.plot(Phi_pot_center(phi_left), C_DL_num - C_DL_ana, label='Difference')
-plt.plot(Phi_pot_center(phi_left), C_dl_num, label='Numerical')
+# plt.plot(Phi_pot_center(phi_left), C_dl_num, label='Numerical')
plt.plot(Phi_pot_center(phi_left), C_DL_ana, label='Analytical')
plt.grid()
plt.legend()
diff --git a/examples/ReproducableCode/ParameterAnalysis/FractionsHeightLeft.py b/examples/ReproducableCode/ParameterAnalysis/FractionsHeightLeft.py
index cf11f02196958e587386270dcb323501039ff042..9e7d8e12f52f9c6b5ff9b067213523a424607f1d 100644
--- a/examples/ReproducableCode/ParameterAnalysis/FractionsHeightLeft.py
+++ b/examples/ReproducableCode/ParameterAnalysis/FractionsHeightLeft.py
@@ -33,11 +33,12 @@ number_cells = 1024*4
refinement_style = 'hard_log'
rtol = 1e-8
max_iter = 10_000
+relax_param = 0.05
-phi_left_vec = np.linspace(-10, 10, 101)
+phi_left_vec = np.linspace(-10, 10, 201)
y_A, y_C, y_S, phi, p, x = [], [], [], [], [], []
for phi_left_ in phi_left_vec:
- y_A_, y_C_, phi_, p_, x_ = solve_System_2eq(phi_left_, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, refinement_style=refinement_style, return_type='Vector', max_iter=max_iter, rtol=rtol)
+ y_A_, y_C_, phi_, p_, x_ = solve_System_2eq(phi_left_, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, relax_param=relax_param, refinement_style=refinement_style, return_type='Vector', max_iter=max_iter, rtol=rtol)
y_S_ = 1 - y_A_ - y_C_
y_A.append(y_A_)
y_C.append(y_C_)
diff --git a/examples/ReproducableCode/ValidationSimplifiedSystem.py b/examples/ReproducableCode/ValidationSimplifiedSystem.py
index d8a2d835df8ee8ef65c2d3c00197c22d9b4c2a19..9d5c47de097d75288059fcb021f945aa1fed9fb4 100644
--- a/examples/ReproducableCode/ValidationSimplifiedSystem.py
+++ b/examples/ReproducableCode/ValidationSimplifiedSystem.py
@@ -32,7 +32,6 @@ y_A_R = 0.01
y_C_R = 0.01
z_A = -1.0
z_C = 1.0
-K = 'incompressible'
Lambda2 = 8.553e-6
a2 = 7.5412e-4
solvation = 15
@@ -40,6 +39,57 @@ number_cells = 1024#*4
refinement_style = 'log'
rtol = 1e-8
+# # Incompressible
+# K = 'incompressible'
+
+# # solve the complete system
+# y_A_4eq, y_C_4eq, phi_4eq, p_4eq, x_4eq = solve_System_4eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation, relax_param=0.05, refinement_style='hard_log', return_type='Vector', max_iter=1_000, rtol=rtol)
+
+# # solve the simplified system
+# y_A_2eq, y_C_2eq, phi_2eq, p_2eq, x_2eq = solve_System_2eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation, relax_param=0.05, refinement_style='hard_log', return_type='Vector', max_iter=1_000, rtol=rtol)
+
+# # Evaluate the difference
+# plt.figure()
+# # plt.plot(x_4eq, y_A_4eq, label='y_A_4eq')
+# # plt.plot(x_2eq, y_A_2eq, label='y_A_2eq')
+# plt.plot(x_4eq, y_A_4eq - y_A_2eq, label='y_A_4eq - y_A_2eq')
+# plt.grid()
+# plt.xlim(0, 0.05)
+# plt.legend()
+# plt.show()
+
+# plt.figure()
+# # plt.plot(x_4eq, y_C_4eq, label='y_C_4eq')
+# # plt.plot(x_2eq, y_C_2eq, label='y_C_2eq')
+# plt.plot(x_4eq, y_C_4eq - y_C_2eq, label='y_C_4eq - y_C_2eq')
+# plt.grid()
+# plt.xlim(0, 0.05)
+# plt.legend()
+# plt.show()
+
+# plt.figure()
+# # plt.plot(x_4eq, phi_4eq, label='phi_4eq')
+# # plt.plot(x_2eq, phi_2eq, label='phi_2eq')
+# plt.plot(x_4eq, phi_4eq - phi_2eq, label='phi_4eq - phi_2eq')
+# plt.grid()
+# plt.xlim(0, 0.05)
+# plt.legend()
+# plt.show()
+
+# plt.figure()
+# # plt.plot(x_4eq, p_4eq, label='p_4eq')
+# # plt.plot(x_2eq, p_2eq, label='p_2eq')
+# plt.plot(x_4eq, p_4eq - p_2eq, label='p_4eq - p_2eq')
+# plt.grid()
+# plt.xlim(0, 0.05)
+# plt.legend()
+# plt.show()
+
+
+
+# Compressible
+K = 10_000
+
# solve the complete system
y_A_4eq, y_C_4eq, phi_4eq, p_4eq, x_4eq = solve_System_4eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation, relax_param=0.05, refinement_style='hard_log', return_type='Vector', max_iter=1_000, rtol=rtol)
@@ -48,20 +98,20 @@ y_A_2eq, y_C_2eq, phi_2eq, p_2eq, x_2eq = solve_System_2eq(phi_left, phi_right,
# Evaluate the difference
plt.figure()
-# plt.plot(x_4eq, y_A_4eq, label='y_A_4eq')
-# plt.plot(x_2eq, y_A_2eq, label='y_A_2eq')
-plt.plot(x_4eq, y_A_4eq - y_A_2eq, label='y_A_4eq - y_A_2eq')
+plt.plot(x_4eq, y_A_4eq, label='y_A_4eq')
+plt.plot(x_2eq, y_A_2eq, label='y_A_2eq')
+# plt.plot(x_4eq, y_A_4eq - y_A_2eq, label='y_A_4eq - y_A_2eq')
plt.grid()
-plt.xlim(0, 0.05)
+plt.xlim(0, 0.1)
plt.legend()
plt.show()
plt.figure()
-# plt.plot(x_4eq, y_C_4eq, label='y_C_4eq')
-# plt.plot(x_2eq, y_C_2eq, label='y_C_2eq')
-plt.plot(x_4eq, y_C_4eq - y_C_2eq, label='y_C_4eq - y_C_2eq')
+plt.plot(x_4eq, y_C_4eq, label='y_C_4eq')
+plt.plot(x_2eq, y_C_2eq, label='y_C_2eq')
+# plt.plot(x_4eq, y_C_4eq - y_C_2eq, label='y_C_4eq - y_C_2eq')
plt.grid()
-plt.xlim(0, 0.05)
+plt.xlim(0, 0.1)
plt.legend()
plt.show()
@@ -70,7 +120,7 @@ plt.figure()
# plt.plot(x_2eq, phi_2eq, label='phi_2eq')
plt.plot(x_4eq, phi_4eq - phi_2eq, label='phi_4eq - phi_2eq')
plt.grid()
-plt.xlim(0, 0.05)
+plt.xlim(0, 0.1)
plt.legend()
plt.show()
@@ -79,6 +129,6 @@ plt.figure()
# plt.plot(x_2eq, p_2eq, label='p_2eq')
plt.plot(x_4eq, p_4eq - p_2eq, label='p_4eq - p_2eq')
plt.grid()
-plt.xlim(0, 0.05)
+plt.xlim(0, 0.1)
plt.legend()
plt.show()
\ No newline at end of file
diff --git a/src/ElectrolyticDiode.py b/src/ElectrolyticDiode.py
new file mode 100644
index 0000000000000000000000000000000000000000..647f6dd821d3ec2ce6002147d864c89c68b3d4a7
--- /dev/null
+++ b/src/ElectrolyticDiode.py
@@ -0,0 +1,335 @@
+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+
+This module solves the dimensionless system of equations for the example of an electric diode
+
+# ! Disclaimer: Comparing the solution with: Entropy and convergence analysis for two finite volume schemes for a Nernst–Planck–Poisson system with ion volume constraints (Gaudeu, Fuhrmann) DOI: 10.1007/s00211-022-01279-y yields different results. The electric potential seems to fit quite good, but shows some slight differences - the distribution fits and the rough values also. The atomic fractions do not fit with the solution in the paper. The rectification effect, observed in the paper, cannot be observed in our solution. The ion concentrations are same high for all three biases and show almost now difference.
+The paper solves without the equation for the pressure and uses different molar volumina for the anions and cations compared to the solvent. These different molar voluma might be the reason for the differences in the atomic fractions. The used dimensionless form of the model does not allow to set different molar volumina. This needs further investigation.
+'''
+
+import numpy as np
+from mpi4py import MPI
+from dolfinx import mesh, fem, log, io
+from dolfinx.fem.petsc import NonlinearProblem
+from dolfinx.nls.petsc import NewtonSolver
+from ufl import TestFunctions, split, dot, grad, dx, inner, ln
+from basix.ufl import element, mixed_element
+import matplotlib.pyplot as plt
+from ufl import Measure
+from dolfinx.mesh import locate_entities, meshtags
+
+def ElectrolyticDiode(Bias_type:str, phi_bias:float, z_A:float, z_C:float, y_A_bath:float, y_C_bath:float, K:float|str, Lambda2:float, a2:float, number_cells:list, solvation:float = 0, PoissonBoltzmann:bool=False, relax_param:float=None, Lx:float=2, Ly:float=10, rtol:float=1e-8, max_iter:float=500, return_type:str='Vector'):
+ '''
+ Solves the system of equations for the example of an electric diode
+
+ Solve the dimensionless system of equations presented in: Numerical Treatment of a Thermodynamically Consistent Electrolyte Model, B.Sc. Thesis Habscheid 2024
+
+ ! If the Newton solver diverges, you may try to reduce the relaxation parameter.
+
+ Parameters
+ ----------
+ Bias_type : str
+ ForwardBias, NoBias, BackwardBias
+ phi_bias : float
+ Bias in φ
+ z_A : float
+ Charge number of species A
+ z_C : float
+ Charge number of species C
+ y_A_bath : float
+ Bath concentration of anions
+ y_C_bath : float
+ Bath concentration of cations
+ K : float | str
+ Dimensioness bulk modulus of the electrolyte. If 'incompressible', the system is solved for an incompressible electrolyte
+ ! Only implemented for incompressible mixtures, yet
+ Lambda2 : float
+ Dimensionless parameter
+ a2 : float
+ Dimensionless parameter
+ number_cells : int
+ Number of cells in the mesh
+ solvation : float, optional
+ solvation number, by default 0
+ PoissonBoltzmann : bool, optional
+ Solve classical Nernst-Planck model with the use of the Poisson-Boltzmann formulation if True, else solve the presented model by Dreyer, Guhlke, Müller, by default False
+ relax_param : float, optional
+ Relaxation parameter for the Newton solver
+ xₙ₊₁ = γ xₙ f(xₙ)/f'(xₙ) with γ the relaxation parameter
+ , by default None -> Determined automatically
+ Lx : float, optional
+ Length of domain in x-direction, by default 2
+ Ly : float, optional
+ Length of domain in y-direction, by default 10
+ rtol : float, optional
+ Relative tolerance for Newton solver, by default 1e-8
+ max_iter : float, optional
+ Maximum number of Newton iterations, by default 500
+ return_type : str, optional
+ Type of return value, by default 'Vector'
+ 'Vector' -> Returns the solution as a numpy array for triangle elements
+ 'Extended' -> Returns the solution as both, a numpy array for triangle elements and the fenicsx function u
+
+
+ Returns
+ -------
+ y_A_vector, y_C_vector, phi_vector, p_vector, x_vector
+ Returns atomic fractions for species A and C, electric potential, pressure, and the mesh as numpy arrays for triangle elements. x_vector is a list of both dimensions [x, y]
+ If return_type is 'Extended', also returns the fenicsx function u
+ '''
+ if K != 'incompressible':
+ raise NotImplementedError('Only incompressible electrolytes are implemented yet')
+ x0 = np.array([0, 0])
+ x1 = np.array([Lx, Ly])
+ match Bias_type:
+ case 'ForwardBias': phi_bias = phi_bias
+ case 'NoBias': phi_bias = 0
+ case 'BackwardBias': phi_bias = -phi_bias
+
+ # Define boundaries
+ geom_tol = 1E-4 # ! Geometric tolerance, may need to be adjusted
+ def Left(x):
+ return np.isclose(x[0], x0[0], geom_tol, geom_tol)
+
+ def Right(x):
+ return np.isclose(x[0], x1[0], geom_tol, geom_tol)
+
+ def Top(x):
+ return np.isclose(x[1], x1[1], geom_tol, geom_tol)
+
+ def Bottom(x):
+ return np.isclose(x[1], x0[1], geom_tol, geom_tol)
+
+ def Right_Top(x):
+ NeumannPart = np.logical_and(1/2*Ly < x[1], x[1] < 3/4*Ly)
+ RightPart = np.isclose(x[0], x1[0], geom_tol, geom_tol)
+ return np.logical_and(RightPart, NeumannPart)
+
+ def Right_Bottom(x):
+ NeumannPart = np.logical_and(1/4*Ly < x[1], x[1] < 1/2*Ly)
+ RightPart = np.isclose(x[0], x1[0], geom_tol, geom_tol)
+ return np.logical_and(RightPart, NeumannPart)
+
+ def Left_Top(x):
+ NeumannPart = np.logical_and(1/2*Ly < x[1], x[1] < 3/4*Ly)
+ LeftPart = np.isclose(x[0], x0[0], geom_tol, geom_tol)
+ return np.logical_and(LeftPart, NeumannPart)
+
+ def Left_Bottom(x):
+ NeumannPart = np.logical_and(1/4*Ly < x[1], x[1] < 1/2*Ly)
+ LeftPart = np.isclose(x[0], x0[0], geom_tol, geom_tol)
+ return np.logical_and(LeftPart, NeumannPart)
+
+ # Create the mesh
+ msh = mesh.create_rectangle(MPI.COMM_WORLD, [x0, x1], number_cells)
+
+ # Define Finite Elements
+ CG1_elem = element('Lagrange', msh.basix_cell(), 1)
+
+ # Define Mixed Function Space
+ W_elem = mixed_element([CG1_elem, CG1_elem, CG1_elem, CG1_elem])
+ W = fem.functionspace(msh, W_elem)
+
+ # Define Trial- and Testfunctions
+ u = fem.Function(W)
+ y_A, y_C, phi, p = split(u)
+ (v_A, v_C, v_1, v_2) = TestFunctions(W)
+
+ # Define boundary conditions
+ W0, _ = W.sub(0).collapse()
+ W1, _ = W.sub(1).collapse()
+ W2, _ = W.sub(2).collapse()
+ W3, _ = W.sub(3).collapse()
+
+ def phi_bottom_(x):
+ return np.full_like(x[1], 0)
+ def phi_top_(x):
+ return np.full_like(x[1], phi_bias)
+ def y_A_bottom_(x):
+ return np.full_like(x[1], y_A_bath)
+ def y_A_top_(x):
+ return np.full_like(x[1], y_A_bath)
+ def y_C_bottom_(x):
+ return np.full_like(x[1], y_C_bath)
+ def y_C_top_(x):
+ return np.full_like(x[1], y_C_bath)
+
+ phi_bottom_bcs = fem.Function(W2)
+ phi_bottom_bcs.interpolate(phi_bottom_)
+ phi_top_bcs = fem.Function(W2)
+ phi_top_bcs.interpolate(phi_top_)
+ y_A_bottom_bcs = fem.Function(W0)
+ y_A_bottom_bcs.interpolate(y_A_bottom_)
+ y_A_top_bcs = fem.Function(W0)
+ y_A_top_bcs.interpolate(y_A_top_)
+ y_C_bottom_bcs = fem.Function(W1)
+ y_C_bottom_bcs.interpolate(y_C_bottom_)
+ y_C_top_bcs = fem.Function(W1)
+ y_C_top_bcs.interpolate(y_C_top_)
+
+ # Identify face tags and create dirichlet conditions
+ facet_bottom_dofs = fem.locate_dofs_geometrical((W.sub(2), W.sub(2).collapse()[0]), Bottom)
+ facet_top_dofs = fem.locate_dofs_geometrical((W.sub(2), W.sub(2).collapse()[0]), Top)
+ bc_bottom_phi = fem.dirichletbc(phi_bottom_bcs, facet_bottom_dofs, W.sub(2))
+ bc_top_phi = fem.dirichletbc(phi_top_bcs, facet_top_dofs, W.sub(2))
+
+ facet_bottom_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Bottom)
+ bc_bottom_y_A = fem.dirichletbc(y_A_bottom_bcs, facet_bottom_dofs, W.sub(0))
+ facet_top_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Top)
+ bc_top_y_A = fem.dirichletbc(y_A_top_bcs, facet_top_dofs, W.sub(0))
+
+ facet_bottom_dofs = fem.locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]), Bottom)
+ bc_bottom_y_C = fem.dirichletbc(y_C_bottom_bcs, facet_bottom_dofs, W.sub(1))
+ facet_top_dofs = fem.locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]), Top)
+ bc_top_y_C = fem.dirichletbc(y_C_top_bcs, facet_top_dofs, W.sub(1))
+
+ # Collect boundary conditions
+ bcs = [bc_bottom_phi, bc_top_phi, bc_bottom_y_A, bc_top_y_A, bc_bottom_y_C, bc_top_y_C]
+
+ # Variational formulation
+ if K == 'incompressible':
+ def nF(y_A, y_C):
+ return (z_C * y_C + z_A * y_A)
+
+ A = (
+ inner(grad(phi), grad(v_1)) * dx
+ - Lambda2 * nF(y_A, y_C) * v_1 * dx
+ ) + (
+ inner(grad(p), grad(v_2)) * dx
+ + 1 / a2 * nF(y_A, y_C) * dot(grad(phi), grad(v_2)) * dx
+ ) + (
+ inner(grad(ln(y_A) + a2 * solvation * p - ln(1-y_A-y_C) + z_A * phi), grad(v_A)) * dx
+ + inner(grad(ln(y_C) + a2 * solvation * p- ln(1-y_A-y_C) + z_C * phi), grad(v_C)) * dx
+ )
+
+ # Define Neumann boundaries
+ boundaries = [(0, Right_Top), (1, Right_Bottom), (2, Left_Top), (3, Left_Bottom)]
+ facet_indices, facet_markers = [], []
+ fdim = msh.topology.dim - 1
+ for (marker, locator) in boundaries:
+ facets = locate_entities(msh, fdim, locator)
+ facet_indices.append(facets)
+ facet_markers.append(np.full_like(facets, marker))
+ facet_indices = np.hstack(facet_indices).astype(np.int32)
+ facet_markers = np.hstack(facet_markers).astype(np.int32)
+ sorted_facets = np.argsort(facet_indices)
+ facet_tag = meshtags(msh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
+ ds = Measure("ds", domain=msh, subdomain_data=facet_tag)
+ L = 0
+ L += (
+ inner(g_phi, v_1) * ds(0)
+ + inner(-g_phi, v_1) * ds(1)
+ # Uncomment, if you want to apply the Neumann bcs on both sides/ left side
+ # + inner(g_phi, v_1) * ds(2)
+ # + inner(-g_phi, v_1) * ds(3)
+ )
+
+ F = A + L
+
+ # Initialize initial guess for u
+ y_C_init = fem.Function(W1)
+ y_A_init = fem.Function(W0)
+ y_C_init.interpolate(lambda x: np.full_like(x[0], y_C_bath))
+ y_A_init.interpolate(lambda x: np.full_like(x[0], y_A_bath))
+
+ with u.vector.localForm() as u_loc:
+ u_loc.set(0)
+ u.sub(0).interpolate(y_A_init)
+ u.sub(1).interpolate(y_C_init)
+
+ # Define Nonlinear Problem
+ problem = NonlinearProblem(F, u, bcs=bcs)
+
+ solver = NewtonSolver(MPI.COMM_WORLD, problem)
+ solver.convergence_criterion = 'residual' #"incremental"
+ solver.rtol = rtol
+ solver.relaxation_parameter = relax_param
+ solver.max_it = max_iter
+ solver.report = True
+
+ # Solve the problem
+ log.set_log_level(log.LogLevel.INFO)
+ n, converged = solver.solve(u)
+ assert (converged)
+ print(f"Number of interations: {n:d}")
+
+ # Extract solution
+ y_A_vector = u.sub(0).collapse().x.array
+ y_C_vector = u.sub(1).collapse().x.array
+ phi_vector = u.sub(2).collapse().x.array
+ p_vector = u.sub(3).collapse().x.array
+
+ X = W.sub(0).collapse()[0].tabulate_dof_coordinates()
+ x_vector = [X[:, 0], X[:, 1]]
+
+ if return_type == 'Vector':
+ return y_A_vector, y_C_vector, phi_vector, p_vector, x_vector
+ elif return_type == 'Extended':
+ return y_A_vector, y_C_vector, phi_vector, p_vector, x_vector, u
+ else:
+ raise ValueError('Invalid return_type')
+
+if __name__ == '__main__':
+ phi_bias = 10
+ Bias_type = 'BackwardBias' # 'ForwardBias', 'NoBias', 'BackwardBias'
+ g_phi = 5
+ y_fixed = 0.01
+ z_A = -1.0
+ z_C = 1.0
+ K = 'incompressible'
+ # Lambda2 = 1e-7#7#1e-3
+ # a2 = 1e-4
+ # Lambda2 = 1e-8#8.33e-8
+ # a2 = 7.3656e-4
+ Lambda2 = 8.553e-6
+ a2 = 7.5412e-4
+ number_cells = [20, 200]
+ Lx = 2
+ Ly = 10
+ x0 = np.array([0, 0])
+ x1 = np.array([Lx, Ly])
+ refinement_style = 'uniform'
+ solvation = 0
+ PoissonBoltzmann = False
+ rtol = 1e-3 # ToDo: Change back to 1e-8, currently just for testing
+ relax_param = 0.05
+ max_iter = 15_000
+
+ # Solve the system
+ y_A, y_C, phi, p, X = ElectrolyticDiode(Bias_type, phi_bias, z_A, z_C, y_fixed, y_fixed, K, Lambda2, a2, number_cells, solvation, PoissonBoltzmann, relax_param, Lx, Ly, rtol, max_iter, return_type='Vector')
+ x, y = X[0], X[1]
+ y_S = 1 - y_A - y_C
+
+ # Plot results
+ levelsf = 10
+ levels = 10
+ fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(8,10))
+
+ c = axs[0,0].tricontourf(x, y, y_A, levelsf)
+ fig.colorbar(c)
+ axs[0,0].tricontour(x, y, y_A, levels, colors='black')
+ axs[0,0].set_title('$y_A$')
+
+ c = axs[0,1].tricontourf(x, y, y_C, levelsf)
+ fig.colorbar(c)
+ axs[0,1].tricontour(x, y, y_C, levels, colors='black')
+ axs[0,1].set_title('$y_C$')
+
+ c = axs[0,2].tricontourf(x, y, y_S, levelsf)
+ fig.colorbar(c)
+ axs[0,2].tricontour(x, y, y_S, levels, colors='black')
+ axs[0,2].set_title('$y_S$')
+
+ c = axs[1,0].tricontourf(x, y, phi, levelsf)
+ fig.colorbar(c)
+ axs[1,0].tricontour(x, y, phi, levels, colors='black')
+ axs[1,0].set_title('$\\varphi$')
+
+ c = axs[1,1].tricontourf(x, y, p, levelsf)
+ fig.colorbar(c)
+ axs[1,1].tricontour(x, y, p, levels, colors='black')
+ axs[1,1].set_title('$p$')
+
+ fig.tight_layout()
+ fig.show()
\ No newline at end of file
diff --git a/src/Eq02.py b/src/Eq02.py
index 7ecd7d0e72574f5207edd068d100190aeb56e7c9..7c2c9fc4e929d0ee99ea24844a695f53df48750c 100644
--- a/src/Eq02.py
+++ b/src/Eq02.py
@@ -197,7 +197,7 @@ def solve_System_2eq(phi_left:float, phi_right:float, p_right:float, z_A:float,
# Define variational problem
if K == 'incompressible':
# total free charge density
- def nF(y_A, y_C):
+ def nF(y_A, y_C, p):
return (z_C * y_C + z_A * y_A)
else:
# total number density
@@ -207,13 +207,13 @@ def solve_System_2eq(phi_left:float, phi_right:float, p_right:float, z_A:float,
# total free charge density
def nF(y_A, y_C, p):
return (z_C * y_C + z_A * y_A) * n(p)
- # Variational Form
+ # Variational Form
A = (
inner(grad(phi), grad(v_1)) * dx
- - 1 / Lambda2 * nF(y_A(phi, p), y_C(phi, p)) * v_1 * dx
+ - 1 / Lambda2 * nF(y_A(phi, p), y_C(phi, p), p) * v_1 * dx
) + (
inner(grad(p), grad(v_2)) * dx
- + 1 / a2 * nF(y_A(phi, p), y_C(phi, p)) * dot(grad(phi), grad(v_2)) * dx
+ + 1 / a2 * nF(y_A(phi, p), y_C(phi, p), p) * dot(grad(phi), grad(v_2)) * dx
)
F = A
diff --git a/src/Helper_DoubleLayerCapacity.py b/src/Helper_DoubleLayerCapacity.py
index 4da01ebab88716abe73041b31c0f973a3470e9c4..77149453f1646d83950c3add8b5c16eb56ab010e 100644
--- a/src/Helper_DoubleLayerCapacity.py
+++ b/src/Helper_DoubleLayerCapacity.py
@@ -152,34 +152,69 @@ def Q_DL_dimless_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float
Charge of the system in dimensionless units
'''
z_N = 0
+ E_p = p_R
if K == 'incompressible':
+ # Assume E_p = p_R = 0
D_A = y_A_R / (np.exp(-(solvation+1)*a2*p_R-z_A*phi_R))
D_C = y_C_R / (np.exp(-(solvation+1)*a2*p_R-z_C*phi_R))
D_N = y_N_R / (np.exp(-(solvation+1)*a2*p_R-z_N*phi_R))
- E_p = p_R
- CLambda_L = np.log(D_A * np.exp(-z_A * phi_L) + D_C * np.exp(-z_C * phi_L) + D_N * np.exp(-z_N * phi_L))
- CLambda_R = np.log(D_A * np.exp(-z_A * phi_R) + D_C * np.exp(-z_C * phi_R) + D_N * np.exp(-z_N * phi_R))
+ D_tilde_A = np.log(D_A)
+ D_tilde_C = np.log(D_C)
+ D_tilde_N = np.log(D_N)
+
+ CLambda_L = - z_A * phi_L + D_tilde_A - z_C * phi_L + D_tilde_C - z_N * phi_L + D_tilde_N
+ CLambda_R = - z_A * phi_R + D_tilde_A - z_C * phi_R + D_tilde_C - z_N * phi_R + D_tilde_N
- # dx_Phi_L = np.sqrt((CLambda_L - (solvation+1) * a2 * E_p) * 2 / (solvation + 1))
- # dx_Phi_R = np.sqrt((CLambda_R - (solvation+1) * a2 * E_p) * 2 / (solvation + 1))
- # Q_DL = Lambda2 * (dx_Phi_L - dx_Phi_R)
Lambda = np.sqrt(Lambda2)
- Q_DL = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * np.sqrt(2/(solvation+1)) * (np.sqrt(CLambda_L + E_p * a2) - np.sqrt(CLambda_R + E_p * a2))
+ a = np.sqrt(a2)
+
+ N = 3
+
+ Q_DL = Lambda * a * (np.sqrt(-CLambda_L/(-(N-1)*(solvation+1)-1)) - np.sqrt(-CLambda_R/(-(N-1)*(solvation+1)-1)))
+
+ # CLambda_L = np.log(D_A * np.exp(-z_A * phi_L) + D_C * np.exp(-z_C * phi_L) + D_N * np.exp(-z_N * phi_L))
+ # CLambda_R = np.log(D_A * np.exp(-z_A * phi_R) + D_C * np.exp(-z_C * phi_R) + D_N * np.exp(-z_N * phi_R))
+
+ # Lambda = np.sqrt(Lambda2)
+ # Q_DL = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * np.sqrt(2) * (np.sqrt(CLambda_L/(solvation+1) - E_p * a2) - np.sqrt(CLambda_R/(solvation+1) - E_p * a2))
else:
C_A = y_A_R / ((K + p_R - 1)**(-(solvation+1)*a2*K)*np.exp(-z_A*phi_R))
C_C = y_C_R / ((K + p_R - 1)**(-(solvation+1)*a2*K)*np.exp(-z_C*phi_R))
C_N = y_N_R / ((K + p_R - 1)**(-(solvation+1)*a2*K)*np.exp(-z_N*phi_R))
- E_p = p_R
- Lambda_tilda_L = C_A * np.exp(-z_A*phi_L) + C_C * np.exp(-z_C*phi_L) + C_N * np.exp(-z_N*phi_L)
- Lambda_tilda_R = C_A * np.exp(-z_A*phi_R) + C_C * np.exp(-z_C*phi_R) + C_N * np.exp(-z_N*phi_R)
+ # # Lambda_tilda_L = C_A * np.exp(-z_A*phi_L) + C_C * np.exp(-z_C*phi_L) + C_N * np.exp(-z_N*phi_L)
+ # # Lambda_tilda_R = C_A * np.exp(-z_A*phi_R) + C_C * np.exp(-z_C*phi_R) + C_N * np.exp(-z_N*phi_R)
+
+ # # # Lambda_hat_L = np.exp(np.log(Lambda_tilda_L) / (a2 * K))
+ # # # Lambda_hat_R = np.exp(np.log(Lambda_tilda_R) / (a2 * K))
+
+ # # dx_Phi_L = np.sqrt(2 * a2/Lambda2 * (1 - K - E_p + Lambda_hat_L))
+ # # dx_Phi_R = np.sqrt(2 * a2/Lambda2 * (1 - K - E_p + Lambda_hat_R))
+
+ # # Q_DL = Lambda2 * (dx_Phi_L - dx_Phi_R)
+
+ # Lambda_tilde_L = C_A * np.exp(-z_A*phi_L) + C_C * np.exp(-z_C*phi_L) + C_N * np.exp(-z_N*phi_L)
+ # Lambda_tilde_R = C_A * np.exp(-z_A*phi_R) + C_C * np.exp(-z_C*phi_R) + C_N * np.exp(-z_N*phi_R)
+
+ # Lambda_hat_L = (1/Lambda_tilde_L)**(-1/((solvation+1)*a2*K))
+ # Lambda_hat_R = (1/Lambda_tilde_R)**(-1/((solvation+1)*a2*K))
+
+ # a = np.sqrt(a2)
+ # Lambda = np.sqrt(Lambda2)
+ # K_tilde = 1 - K + E_p
+
+ # Q_DL = Lambda * a * np.sqrt(2) * (np.sqrt(Lambda_hat_L - K_tilde) - np.sqrt(Lambda_hat_R - K_tilde))
+
+ Lambda_L = C_A * np.exp(-z_A*phi_L) + C_C * np.exp(-z_C*phi_L) + C_N * np.exp(-z_N*phi_L)
+ Lambda_R = C_A * np.exp(-z_A*phi_R) + C_C * np.exp(-z_C*phi_R) + C_N * np.exp(-z_N*phi_R)
- Lambda_hat_L = np.exp(np.log(Lambda_tilda_L) / (a2 * K))
- Lambda_hat_R = np.exp(np.log(Lambda_tilda_R) / (a2 * K))
+ B = -a2 * K
+ A = 1/2 * Lambda2/a2
+ K_tilde = K + E_p - 1
- dx_Phi_L = np.sqrt(2 * a2/Lambda2 * (1 - K - E_p + Lambda_hat_L))
- dx_Phi_R = np.sqrt(2 * a2/Lambda2 * (1 - K - E_p + Lambda_hat_R))
+ dx_Phi_L = (np.power(1/Lambda_L, B) - K_tilde) / A
+ dx_Phi_R = (np.power(1/Lambda_R, B) - K_tilde) / A
Q_DL = Lambda2 * (dx_Phi_L - dx_Phi_R)
return Q_DL
diff --git a/src/__pycache__/Eq02.cpython-312.pyc b/src/__pycache__/Eq02.cpython-312.pyc
index 99acdd38417e0f045265675a551c6f42a3df0a2d..972d71c81431a18eff375df5353a304de3d68c3f 100644
Binary files a/src/__pycache__/Eq02.cpython-312.pyc and b/src/__pycache__/Eq02.cpython-312.pyc differ
diff --git a/src/__pycache__/Helper_DoubleLayerCapacity.cpython-312.pyc b/src/__pycache__/Helper_DoubleLayerCapacity.cpython-312.pyc
index ec8a8993fff2cb48986394c354c67b2ed42ec555..82d88e6411714334b2d2d3d89e81ae6b20e02a4c 100644
Binary files a/src/__pycache__/Helper_DoubleLayerCapacity.cpython-312.pyc and b/src/__pycache__/Helper_DoubleLayerCapacity.cpython-312.pyc differ