Skip to content

Commit

Permalink
lint
Browse files Browse the repository at this point in the history
  • Loading branch information
@tlroy
tlroy committed Dec 5, 2023
1 parent 35412a5 commit 4f9d843
Show file tree
Hide file tree
Showing 34 changed files with 618 additions and 584 deletions.
7 changes: 5 additions & 2 deletions echemfem/cylindricalmeasure.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,10 +8,12 @@
from ufl.constantvalue import as_ufl
from ufl.domain import extract_domains


class CylindricalMeasure(Measure):
def __init__(self, radius, *args, **kwargs):
self.radius = radius
super().__init__(*args, **kwargs)

def __rmul__(self, integrand):
"""Multiply a scalar expression with measure to construct a form with
a single integral.
Expand Down Expand Up @@ -75,6 +77,7 @@ def __rmul__(self, integrand):
metadata=self.metadata(),
subdomain_data=self.subdomain_data())
return Form([integral])

def reconstruct(self,
integral_type=None,
subdomain_id=None,
Expand All @@ -90,5 +93,5 @@ def reconstruct(self,
if subdomain_data is None:
subdomain_data = self.subdomain_data()
return CylindricalMeasure(self.radius, self.integral_type(),
domain=domain, subdomain_id=subdomain_id,
metadata=metadata, subdomain_data=subdomain_data)
domain=domain, subdomain_id=subdomain_id,
metadata=metadata, subdomain_data=subdomain_data)
200 changes: 98 additions & 102 deletions echemfem/solver.py
View file

Large diffs are not rendered by default.

5 changes: 2 additions & 3 deletions echemfem/utility_meshes.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@


def IntervalBoundaryLayerMesh(ncells, length_or_left, ncells_bdlayer, length_bdlayer,
boundary = (2,), right=None, distribution_parameters=None, comm=COMM_WORLD):
boundary=(2,), right=None, distribution_parameters=None, comm=COMM_WORLD):
"""
Generate a boundary layer mesh of an interval.
Expand Down Expand Up @@ -73,7 +73,7 @@ def IntervalBoundaryLayerMesh(ncells, length_or_left, ncells_bdlayer, length_bdl


def RectangleBoundaryLayerMesh(nx, ny, Lx, Ly, n_bdlayer, L_bdlayer, Ly_bdlayer=None, ny_bdlayer=None,
boundary = (1,), reorder=None, distribution_parameters=None, comm=COMM_WORLD):
boundary=(1,), reorder=None, distribution_parameters=None, comm=COMM_WORLD):
"""Generate a boundary layer rectangular mesh using quadrilateral elements
:arg nx: The number of cells in the x direction outside the boundary layers
Expand Down Expand Up @@ -193,4 +193,3 @@ def RectangleBoundaryLayerMesh(nx, ny, Lx, Ly, n_bdlayer, L_bdlayer, Ly_bdlayer=
plex.setLabelValue(dmcommon.FACE_SETS_LABEL, face, 4)

return mesh.Mesh(plex, reorder=reorder, distribution_parameters=distribution_parameters, comm=comm)

58 changes: 26 additions & 32 deletions examples/BMCSL.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,73 +15,68 @@
import matplotlib.gridspec as gridspec
from matplotlib.ticker import FormatStrFormatter
import numpy as np
import math
import math

#Initialize bulk concentrations (1 mM CsCl + 9 mM LiCl)
# Initialize bulk concentrations (1 mM CsCl + 9 mM LiCl)

C_Cs_bulk = 1

C_Li_bulk = 9

C_Cl_bulk = 10


class BMCSLSolver(EchemSolver):

def __init__(self):

delta = 0.00008 #Boundary layer thickness (m)
delta = 0.00008 # Boundary layer thickness (m)

mesh = IntervalBoundaryLayerMesh(200, delta, 800, 1e-8,boundary = (2,))
mesh = IntervalBoundaryLayerMesh(200, delta, 800, 1e-8, boundary=(2,))

x = SpatialCoordinate(mesh)[0]


conc_params = []


conc_params.append({"name": "Cl",
"diffusion coefficient": 20.6e-10, # m^2/s
"bulk": C_Cl_bulk, # mol/m3
"z": -1,
"solvated diameter": 0.0 # co-ion size not relevant at sufficiently high negative surface charge density
})


conc_params.append({"name": "Li",
"diffusion coefficient": 10.6e-10, # m^2/s
"bulk": C_Li_bulk, # mol/m3
"z": 1,
"solvated diameter": 7.6e-10 # m (a_Li = 3.8 A)
"solvated diameter": 7.6e-10 # m (a_Li = 3.8 A)
})


conc_params.append({"name": "Cs",
"diffusion coefficient": 20.6e-10, # m^2/s
"bulk": C_Cs_bulk, # mol/m3
"z": 1,
"solvated diameter": 6.6e-10 # m (a_Cs = 3.3 A)
"solvated diameter": 6.6e-10 # m (a_Cs = 3.3 A)
})

self.U_solid = Constant(0)

physical_params = {"flow": ["migration", "poisson", "diffusion finite size_BMCSL"], #Poisson with BMCSL finite size correction
physical_params = {"flow": ["migration", "poisson", "diffusion finite size_BMCSL"], # Poisson with BMCSL finite size correction
"F": 96485., # C/mol
"R": 8.3144598, # J/K/mol
"T": 273.15 + 25., # K
"U_app": conditional(gt(x,1e-6), 0.0, 0.0),
"U_app": conditional(gt(x, 1e-6), 0.0, 0.0),
"vacuum permittivity": 8.8541878128e-12, # F/m
"relative permittivity": 78.4,
"Avogadro constant": 6.02214076e23, #1/mol
"surface charge density": Constant(0), # C/m^2
"Avogadro constant": 6.02214076e23, # 1/mol
"surface charge density": Constant(0), # C/m^2
}


super().__init__(conc_params, physical_params, mesh, family="CG", p=2)


def set_boundary_markers(self):
self.boundary_markers = {"bulk dirichlet": (1,), #C = C_0
"applied": (1,), # U_liquid = 0
self.boundary_markers = {"bulk dirichlet": (1,), # C = C_0
"applied": (1,), # U_liquid = 0
"poisson neumann": (2,),
}

Expand All @@ -98,33 +93,32 @@ def set_boundary_markers(self):

xmesh = Function(solver.V).interpolate(solver.mesh.coordinates[0])

np.savetxt('xmesh.tsv',xmesh.dat.data)
np.savetxt('xmesh.tsv', xmesh.dat.data)


# Surface charge densities in C/m2

sigmas = np.arange(0, -0.625, -0.025)



store_sigmas = [-0.3, -0.4, -0.5, -0.6] # surface charge densities reported in Figure 1
store_sigmas = [-0.3, -0.4, -0.5, -0.6] # surface charge densities reported in Figure 1

for sigma in sigmas:

sigma = round(sigma,2)

solver.physical_params["surface charge density"].assign(sigma)
sigma = round(sigma, 2)

solver.physical_params["surface charge density"].assign(sigma)

solver.solve()

solver.solve()
C_Cl, C_Li, C_Cs, Phi = solver.u.subfunctions

C_Cl, C_Li , C_Cs, Phi = solver.u.subfunctions
if (sigma in store_sigmas):

if(sigma in store_sigmas):

print(sigma)
print(sigma)

np.savetxt('Cs_{0:.1g}.tsv'.format(sigma),C_Cs.dat.data)
np.savetxt('Cs_{0:.1g}.tsv'.format(sigma), C_Cs.dat.data)

np.savetxt('Li_{0:.1g}.tsv'.format(sigma),C_Li.dat.data)
np.savetxt('Li_{0:.1g}.tsv'.format(sigma), C_Li.dat.data)

np.savetxt('phi_{0:.1g}.tsv'.format(sigma),Phi.dat.data)
np.savetxt('phi_{0:.1g}.tsv'.format(sigma), Phi.dat.data)
102 changes: 50 additions & 52 deletions examples/bicarb.py
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
from firedrake import *
from echemfem import EchemSolver


class CarbonateSolver(EchemSolver):
def __init__(self):
"""
Expand All @@ -10,9 +11,9 @@ def __init__(self):
Reactors. Industrial & Engineering Chemistry Research, 60(31),
pp.11824-11833.
"""
Ly = 1e-3#.6e-3 # m
Lx = 1e-2 # m

Ly = 1e-3 # .6e-3 # m
Lx = 1e-2 # m

mesh = RectangleMesh(200, 100, Lx, Ly, quadrilateral=True)

Expand All @@ -37,48 +38,47 @@ def __init__(self):
C_5_inf = (k5r/k5f)/C_2_inf

def bulk_reaction(y):
yCO2=y[0];
yOH=y[1];
yHCO3=y[2];
yCO3=y[3];
yH=y[4];


dCO2 = -(k1f) *yCO2*yOH \
+(k1r) *yHCO3 \
-(k2f) *yCO2 \
+(k2r) *yHCO3*yH

dOH = -(k1f) *yCO2*yOH \
+(k1r) *yHCO3 \
+(k3f) *yCO3 \
-(k3r) *yOH*yHCO3\
-(k5f) *yOH*yH\
+(k5r)

dHCO3 = (k1f) *yCO2*yOH\
-(k1r) *yHCO3\
+(k3f) *yCO3\
-(k3r) *yOH*yHCO3\
+(k2f) *yCO2 \
-(k2r) *yHCO3*yH \
+(k4f) *yCO3*yH\
-(k4r) *yHCO3

dCO3 = -(k3f) *yCO3 \
+(k3r) *yOH*yHCO3\
-(k4f)*yCO3*yH\
+(k4r) *yHCO3

dH = (k2f) *yCO2 \
-(k2r) *yHCO3*yH \
-(k4f) *yCO3*yH\
+(k4r) *yHCO3\
-(k5f) *yOH*yH\
+(k5r)

return [dCO2, dOH, dHCO3, dCO3, dH]#, 0.]

yCO2 = y[0]
yOH = y[1]
yHCO3 = y[2]
yCO3 = y[3]
yH = y[4]

dCO2 = -(k1f) * yCO2*yOH \
+ (k1r) * yHCO3 \
- (k2f) * yCO2 \
+ (k2r) * yHCO3*yH

dOH = -(k1f) * yCO2*yOH \
+ (k1r) * yHCO3 \
+ (k3f) * yCO3 \
- (k3r) * yOH*yHCO3\
- (k5f) * yOH*yH\
+ (k5r)

dHCO3 = (k1f) * yCO2*yOH\
- (k1r) * yHCO3\
+ (k3f) * yCO3\
- (k3r) * yOH*yHCO3\
+ (k2f) * yCO2 \
- (k2r) * yHCO3*yH \
+ (k4f) * yCO3*yH\
- (k4r) * yHCO3

dCO3 = -(k3f) * yCO3 \
+ (k3r) * yOH*yHCO3\
- (k4f)*yCO3*yH\
+ (k4r) * yHCO3

dH = (k2f) * yCO2 \
- (k2r) * yHCO3*yH \
- (k4f) * yCO3*yH\
+ (k4r) * yHCO3\
- (k5f) * yOH*yH\
+ (k5r)

return [dCO2, dOH, dHCO3, dCO3, dH] # , 0.]

C_CO2_bulk = C_1_inf
C_OH_bulk = C_2_inf
C_HCO3_bulk = C_3_inf
Expand Down Expand Up @@ -108,19 +108,17 @@ def bulk_reaction(y):
"bulk": C_CO32_bulk, # mol/m3
})


conc_params.append({"name": "H",
"diffusion coefficient": 9.311E-9, # m^2/s
"bulk": C_H_bulk, # mol/m3
})

#conc_params.append({"name": "K",
# conc_params.append({"name": "K",
# "diffusion coefficient": 1.96E-9, # m^2/s
# "bulk": C_K_bulk, # mol/m3
# })


physical_params = {"flow": ["advection","diffusion"],
physical_params = {"flow": ["advection", "diffusion"],
"bulk reaction": bulk_reaction,
}

Expand All @@ -136,17 +134,17 @@ def neumann(self, C, conc_params, u):
if name == "OH":
return 2*(1.91E-1)*u[0]


def set_boundary_markers(self):
self.boundary_markers = {"inlet": (1),
"bulk dirichlet": (4),
"outlet": (2,),
"neumann": (3,),
"neumann": (3,),
}

def set_velocity(self):
_, y = SpatialCoordinate(self.mesh)
self.vel = as_vector([1.91*y,Constant(0)]) # m/s
self.vel = as_vector([1.91*y, Constant(0)]) # m/s


solver = CarbonateSolver()
solver.setup_solver()
Expand Down
4 changes: 2 additions & 2 deletions examples/bortels_threeion.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,8 @@ def reaction(u, V):
R = physical_params["R"]
T = physical_params["T"]
eta = V - U
return J0 * (exp(F / R / T * (eta)) -
(C / C_b) * exp(-F / R / T * (eta)))
return J0 * (exp(F / R / T * (eta))
- (C / C_b) * exp(-F / R / T * (eta)))

def reaction_cathode(u):
return reaction(u, Constant(0))
Expand Down
Loading

0 comments on commit 4f9d843

Please sign in to comment.