API

EchemSolver classes

class echemfem.EchemSolver(conc_params, physical_params, mesh, echem_params=[], homog_params=[], gas_params=[], stats_file='', overwrite_stats_file=False, p=1, family='DG', p_penalty=None, SUPG=False, cylindrical=False)

Base class for an electrochemical model solver.

This class is used to create solvers for the transport of multiple ion (or non-charged species).

conc_params

List containing one dictionary for each species. Each dictionary contains physical parameters for the species.

Type:

list of dict

physical_params

Dictionary containing physical parameters

Type:

dict

mesh

Mesh object from firedrake

Type:

firedrake.mesh.MeshGeometry

echem_params

List containing one dictionary for each charge-transfer reaction.

Type:

list, optional

gas_params

List containing one dictionary for each gaseous species. Each dictionary contains physical parameters for the gaseous species. Note that this implementation is not yet validated.

Type:

list, optional

stats_file

File name for performance statistics

Type:

str, optional

overwrite_stats_file

Set to True to overwrite new file.

Type:

bool, optional

p

Polynomial degree of finite elements.

Type:

int, optional

family

Finite element family. Choose between “CG” and “DG”.

Type:

str, optional

p_penalty

Polynomial degree of the DG penalization. Only used for P-multigrid.,

Type:

str, optional

SUPG

Streamline upwind diffusion for stablization of the advection-migration term. Only used for CG.

Type:

bool, optional

cylindrical

if True, uses symmetric cylindrical coordinates (r, z).

Type:

bool, optional

charge_conservation_form(u, test_fn, conc_params, W=None, i_bc=None)

Returns weak form of the charge conservation equation for liquid potential.

Using interior penalty for potential gradient

effective_diffusion(D, phase='liquid')

Bruggeman correlation for effective diffusion in a porous medium

electroneutrality_form(u, test_fn, conc_params)

Returns weak form for the electroneutrality condition.

Only useable for CG.

gas_mass_conservation_form(X, test_fn, gas_params, u=None)

Returns weak form of a mass conservation equation for a gaseous species.

Using interior penalty for the diffusion term and upwinding for the advection tterm.

gas_mass_conservation_form_pressure(p, test_fn, gas_params, u=None, W=None, i_bc=None)

Returns weak form of the mass conservation equation for a gaseous species.

DG: Using interior penalty for the pressure diffusion term

gas_pressure_form(p, test_fn, gas_params, u=None)

Returns weak form of the gas Pressure equation, i.e. the total gaseous mass conservation equation.

DG: Using interior penalty for the pressure diffusion term.

liquid_pressure_form(p, test_fn, conc_params, u=None, W=None, i_bc=None)

Returns weak form of Pressure equation, i.e. the water mass conservation equation.

DG: Using interior penalty for pressure gradient

mass_conservation_form(C, test_fn, conc_params, u=None)

Returns weak form of a mass conservation equation for an aqueous species.

DG: Using interior penalty for the diffusion term and upwinding for the advection-diffusion term.

neumann(C, conc_params, u)

Custom Neumann Boundary condition

Parameters:

  • C – concentration of species i

  • conc_params – conc_params[i], i.e. the concentration parameters of species i

  • u – full solution state

Returns:

normal flux for species i

output_state(u, prefix='results/', initiate=True, **kwargs)

Outputs the provided state variables in a pvd file

Parameters:

  • u (firedrake.function.Function) – State to output. Must be same FunctionSpace as self.u

  • prefix (str) – Path to results directory

  • initiate (bool) – if True, create the output file

  • **kwargs – Arbitrary keyword arguments passed to File.write

poisson_neumann(u)

Custom Neumann Boundary condition for Poisson

potential_poisson_form(u, test_fn, conc_params, solid=False, W=None, i_bc=None)

Returns weak form of the Poisson equation for a potential.

DG: Using interior penalty for potential gradient

set_boundary_markers()

Set boundary attributes

This method should set self.boundary_markers, a dict, where the keys are str: representing the type of boundary condition, and the values are tuple containing the boundary indices. For example :

self.boundary_markers = {"bulk": (1,2,)}

This would set the boundary condition "bulk" on boundary indices 1 and 2.

set_function_spaces(layers, family)

Set base function spaces

self.V : FunctionSpace for concentrations

self.Vu : FunctionSpace for potentials (if using Poisson or Electroneutrality)

self.Vp : FunctionSpace for pressures (if using Darcy)

set_velocity()

Set velocity

This method should set self.vel as a Firedrake vector quantity.

setup_bulk_reactions()

Setting up the homogeneous bulk reactions using homog_params

Setup using the law of mass action

setup_forms(us, v)

Setup weak forms

setup_solver(initial_guess=True, initial_solve=True)

Sets up the initial guess and solver

This creates firedrake.variational_solver.NonlinearVariationalProblem and firedrake.variational_solver.NonlinearVariationalSolver objects stored in self.problem and self.echem_solver, respectively.

Parameters:

  • initial_guess (bool) – If True, sets all concentrations to their given “bulk” value.

  • initial_solve (bool) – If True, solves the system assuming constant concentrations, which provides good initial guesses for the potential(s)

solve()

Solves the problem and outputs the solutions

steady_forms(us, v)

Add steady-state portion of equations to weak forms

class echemfem.TransientEchemSolver(conc_params, physical_params, mesh, echem_params=[], homog_params=[], gas_params=[], stats_file='', overwrite_stats_file=False, p=1, family='DG', p_penalty=None, SUPG=False, cylindrical=False)

Transient electrochemical model solver.

Adds accumulation terms and time loop to create a transient model

mass_accumulation_form(C, C_old, test_fn, conc_params, u=None)

Returns weak form of a mass accumulation term for an aqueous species.

setup_forms(us, v)

Setup weak forms for transient case

solve(times)

Solves the transient problem and outputs the solutions

Parameters:

times (list or numpy.ndarray) – array of the discrete temporal grid

transient_forms(us, v, us_old)

Add transient portion of equations to weak forms

FlowSolver classes

class echemfem.FlowSolver(mesh, fluid_params, boundary_markers)

Base class for a flow solver.

This class is used to create solvers for fluid flow decoupled from the chemistry

mesh

Mesh object from firedrake

Type:

firedrake.mesh.MeshGeometry

fluid_params

Dictionary containing physical parameters

Type:

dict

boundary_markers

Dictionary where the keys are str: representing the type of boundary condition, and the values are tuple containing the boundary indices. For example :

boundary_markers = {"inlet velocity": (1,)}

sets the boundary condition "inlet velocity" on boundary 1.

Type:

dict

setup_dirichlet_bcs()

Set up Dirichlet boundary conditions that are the same regardless of model

setup_functions()

Set up FunctionSpaces, solution Function, and TestFunctions

setup_problem()

Set up PDE system - problem specific

setup_solver()

Set up default nonlinear solver

solve()

Solve system and output results

class echemfem.NavierStokesFlowSolver(mesh, fluid_params, boundary_markers)

Incompressible Navier-Stokes solver

For nondimensional Navier-Stokes, pass Reynolds number to fluid_params. For dimensional Navier-Stokes, pass density and kinematic viscosity.

setup_problem()

Set up PDE system - problem specific

setup_solver(ksp_solver='lu')

Optional PCD preconditioner for nondimensional Navier-Stokes

Parameters:

ksp_solver (str) – "lu" or "pcd"

class echemfem.NavierStokesBrinkmanFlowSolver(mesh, fluid_params, boundary_markers)

Incompressible Navier-Stokes-Brinkman solver

For dimensionless form, pass Reynolds and Darcy numbers to fluid_params. For dimensional form, pass density, permeability, and kinematic or dynamic viscosity.

setup_problem()

Set up PDE system - problem specific

Utility meshes

These are boundary layer meshes, which are modified versions of firedrake.utility_meshes.IntervalMesh() and firedrake.utility_meshes.RectangleMesh().

echemfem.utility_meshes.IntervalBoundaryLayerMesh(ncells, length_or_left, ncells_bdlayer, length_bdlayer, boundary=(2, ), right=None, distribution_parameters=None, comm=<mpi4py.MPI.Intracomm object>)

Generate a boundary layer mesh of an interval.

Parameters:

  • ncells – The number of the cells over the interval, outside the boundary layer.

  • length_or_left – The length of the interval (if right is not provided) or else the left hand boundary point.

  • ncells_bdlayer – The number of the cells in each boundary layer.

  • length_bdlayer – The length of the boundary layer(s).

  • boundary – (optional) boundary markers of the boundary layer(s).

  • right – (optional) position of the right boundary point (in which case length_or_left should be the left boundary point).

Keyword Arguments:

comm – Optional communicator to build the mesh on (defaults to COMM_WORLD).

The left hand boundary point has boundary marker 1, while the right hand point has marker 2.

echemfem.utility_meshes.RectangleBoundaryLayerMesh(nx, ny, Lx, Ly, n_bdlayer, L_bdlayer, Ly_bdlayer=None, ny_bdlayer=None, boundary=(1, ), reorder=None, distribution_parameters=None, comm=<mpi4py.MPI.Intracomm object>)

Generate a boundary layer rectangular mesh using quadrilateral elements

Parameters:

  • nx – The number of cells in the x direction outside the boundary layers

  • ny – The number of cells in the y direction outside the boundary layers

  • Lx – The extent in the x direction

  • Ly – The extent in the y direction

  • n_bdlayer – The number of the cells in each boundary layer. If ny_bdlayer is defined, length of the boundary layer in the x-direction.

  • L_bdlayer – The length of the boundary layer(s). If Ly_bdlayer is defined, length of in the x-direction

  • Ly_bdlayer – (optional) Length of the y-direction boundary layer

  • ny_bdlayer – (optional) number of cells of the y-direction boundary layer

  • boundary – (optional) boundary markers of the boundary layer(s).

Keyword Arguments:

  • reorder – (optional), should the mesh be reordered

  • comm – Optional communicator to build the mesh on (defaults to COMM_WORLD).

The boundary edges in this mesh are numbered as follows:

  • 1: plane x == 0

  • 2: plane x == Lx

  • 3: plane y == 0

  • 4: plane y == Ly