Skip to content

Commit 34e727f

Browse files
j-fuj-fu
authored
Add documentation/comments (#56)
* Add documentation/comments
1 parent b4ace56 commit 34e727f

3 files changed

Lines changed: 350 additions & 112 deletions

File tree

‎notebooks/ICMPBP-EndOfHackathon.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -348,7 +348,7 @@ sys1 = ICMPBSystem(grid, data1);
348348
inival1 = unknowns(sys1, data1);
349349

350350
# ╔═╡ d2df6ed0-e6f5-4677-b790-bfc40de7fd6a
351-
sol1 = solve(sys1; inival = inival1, damp_initial = 0.1, verbose = "n", log=true)
351+
sol1 = solve(sys1; inival = inival1, damp_initial = 0.1, verbose = "n", log = true)
352352

353353
# ╔═╡ 5f153fe4-4476-401e-8674-b6902213c19e
354354
md"""

‎packages/JuliaMPBSolver/src/api.jl‎

Lines changed: 128 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1,17 +1,83 @@
11
using ExtendableGrids
22
using LessUnitful
33
using VoronoiFVM
4+
5+
"""
6+
prepare_grid(n, domain)
7+
8+
Prepare eqidistant simulation grid. Add interface region in the mid of the
9+
domain in order to fix pressure and provide location for ion conservation
10+
constraints.
11+
"""
12+
13+
function prepare_grid(n, domain)
14+
if iseven(n)
15+
n = n + 1
16+
end
17+
X = range(domain..., length = n)
18+
grid = ExtendableGrids.simplexgrid(X)
19+
xmid = (domain[1] + domain[2]) / 2
20+
bfacemask!(grid, [xmid], [xmid], 3, tol = 1.0e-10 * ufac"nm")
21+
return X, grid
22+
end
23+
24+
"""
25+
mpbpsolve(; n=21,
26+
chargenumbers=[-1, 1],
27+
domain=[0, 1.0e-10],
28+
surfacecharge=[0.16, -0.16],
29+
bulkmolarity=0.1)
30+
31+
Solve the modified Poisson-Boltzmann problem with specified bulk molarity.
32+
33+
This function solves the modified Poisson-Boltzmann equation for a 1D domain with specified
34+
surface charges and bulk ion concentrations. The solver uses a finite volume method
35+
with ExtendableGrids and VoronoiFVM.
36+
37+
# Arguments
38+
- `n::Int=21`: Number of grid points along the domain. Default is 21.
39+
- `chargenumbers::Vector{Int}=[-1, 1]`: Charge numbers of the ionic species (e.g., [-1, 1] for monovalent anions and cations).
40+
- `domain::Vector{Float64}=[0, 1.0e-10]`: Domain boundaries in meters [start, end]. Default is [0, 1.0e-10] (0 to 0.1 nm).
41+
- `surfacecharge::Vector{Float64}=[0.16, -0.16]`: Surface charge densities at the domain boundaries in C/m².
42+
- `bulkmolarity::Float64=0.1`: Bulk molarity of the electrolyte solution in mol/L.
43+
44+
# Returns
45+
- `X`: Grid coordinates in m (Vector{Float64})
46+
- `c0`: Solvent concentration in mol/L (Vector{Float64})
47+
- `c1`: Concentration of first ionic species along the grid in mol/L (Vector{Float64})
48+
- `c2`: Concentration of second ionic species along the grid in mol/L(Vector{Float64})
49+
50+
# Examples
51+
```julia
52+
# Solve with default parameters
53+
X, c0, c_anion, c_cation = mpbpsolve()
54+
55+
# Solve with custom parameters
56+
X, c0, c_anion, c_cation = mpbpsolve(
57+
n=51,
58+
chargenumbers=[-1, 1],
59+
domain=[0, 2.0e-10],
60+
surfacecharge=[0.2, -0.2],
61+
bulkmolarity=0.15
62+
)
63+
```
64+
65+
# Notes
66+
- The function automatically creates a uniform grid over the specified domain
67+
- A boundary face mask is applied at the domain midpoint in order fix the pressure value there
68+
- The solver uses damped Newton iteration with initial damping factor of 0.05
69+
- Surface charges should be specified in SI units (C/m²)
70+
- Domain should be specified in SI units (meters)
71+
"""
472
function mpbpsolve(;
573
n = 21,
674
chargenumbers = [-1, 1],
775
domain = [0, 1.0e-10],
876
surfacecharge = [0.16, -0.16],
977
bulkmolarity = 0.1
1078
)
11-
X = range(domain..., length = n)
12-
grid = ExtendableGrids.simplexgrid(X)
13-
xmid = (domain[1] + domain[2]) / 2
14-
bfacemask!(grid, [xmid], [xmid], 3, tol = 1.0e-10 * ufac"nm")
79+
80+
X, grid = prepare_grid(n, domain)
1581
data = ICMPBP.ICMPBData(; z = chargenumbers, q = surfacecharge)
1682
ICMPBP.set_molarity!(data, bulkmolarity)
1783
sys = ICMPBP.ICMPBSystem(grid, data)
@@ -21,21 +87,71 @@ function mpbpsolve(;
2187
return X, c0, c[1, :], c[2, :]
2288
end
2389

90+
"""
91+
icmpbpsolve(; n=21,
92+
chargenumbers=[-1, 1],
93+
domain=[0, 1.0e-10],
94+
surfacecharge=[0.16, -0.16],
95+
averagemolarity=1)
96+
97+
Solve the ion-conserving modified Poisson-Boltzmann problem with specified average molarity.
98+
99+
This function solves the modified Poisson-Boltzmann equation with ion conservation constraints.
100+
Unlike `mpbpsolve`, this solver conserves the total number of ions in the system and uses
101+
an average molarity parameter rather than bulk molarity. This is particularly useful for
102+
systems where ion conservation is physically important.
103+
104+
# Arguments
105+
- `n::Int=21`: Number of grid points along the domain. Must be odd (automatically incremented if even). Default is 21.
106+
- `chargenumbers::Vector{Int}=[-1, 1]`: Charge numbers of the ionic species (e.g., [-1, 1] for monovalent anions and cations).
107+
- `domain::Vector{Float64}=[0, 1.0e-10]`: Domain boundaries in meters [start, end]. Default is [0, 1.0e-10] (0 to 0.1 nm).
108+
- `surfacecharge::Vector{Float64}=[0.16, -0.16]`: Surface charge densities at the domain boundaries in C/m².
109+
- `averagemolarity::Float64=1`: Average molarity of the electrolyte solution in mol/L.
110+
111+
# Returns
112+
- `X`: Grid coordinates in m (Vector{Float64})
113+
- `c0`: Solvent concentration in mol/L (Vector{Float64})
114+
- `c1`: Concentration of first ionic species along the grid in mol/L (Vector{Float64})
115+
- `c2`: Concentration of second ionic species along the grid in mol/L (Vector{Float64})
116+
117+
# Examples
118+
```julia
119+
# Solve with default parameters
120+
X, c0, c_anion, c_cation = icmpbpsolve()
121+
122+
# Solve with custom parameters and higher resolution
123+
X, c0, c_anion, c_cation = icmpbpsolve(
124+
n=101,
125+
chargenumbers=[-2, 1], # divalent anions, monovalent cations
126+
domain=[0, 5.0e-10],
127+
surfacecharge=[0.3, -0.3],
128+
averagemolarity=0.5
129+
)
130+
```
131+
132+
# Notes
133+
- The function enforces an odd number of grid points for numerical stability
134+
- Ion conservation is enforced through the `conserveions=true` parameter in ICMPBData
135+
- The solver employs damped Newton iteration with initial damping factor of 0.05
136+
- Average molarity represents the spatial average concentration, which may differ from bulk concentration due to ion conservation
137+
- Surface charges should be specified in SI units (C/m²)
138+
- Domain should be specified in SI units (meters)
139+
140+
# Differences from mpbpsolve
141+
- Enforces ion conservation constraints
142+
- Uses average rather than bulk molarity
143+
- Requires odd number of grid points
144+
- May show different concentration profiles due to conservation effects
145+
"""
24146
function icmpbpsolve(;
25147
n = 21,
26148
chargenumbers = [-1, 1],
27149
domain = [0, 1.0e-10],
28150
surfacecharge = [0.16, -0.16],
29151
averagemolarity = 1
30152
)
31-
@info averagemolarity
32-
if iseven(n)
33-
n = n + 1
34-
end
35-
X = range(domain..., length = n)
36-
grid = ExtendableGrids.simplexgrid(X)
37-
xmid = (domain[1] + domain[2]) / 2
38-
bfacemask!(grid, [xmid], [xmid], 3, tol = 1.0e-10 * ufac"nm")
153+
154+
X, grid = prepare_grid(n, domain)
39155
data = ICMPBP.ICMPBData(; z = chargenumbers, q = surfacecharge, conserveions = true)
40156
ICMPBP.set_molarity!(data, averagemolarity)
41157
sys = ICMPBP.ICMPBSystem(grid, data)

0 commit comments

Comments
 (0)