#
# Finite Element discretisation class which uses scikit-fem
#
import pybamm
from scipy.sparse import csr_matrix, csc_matrix
from scipy.sparse.linalg import inv
import numpy as np
from pybamm.util import import_optional_dependency
[docs]
class ScikitFiniteElement(pybamm.SpatialMethod):
"""
A class which implements the steps specific to the finite element method during
discretisation. The class uses scikit-fem to discretise the problem to obtain
the mass and stiffness matrices. At present, this class is only used for
solving the Poisson problem -grad^2 u = f in the y-z plane (i.e. not the
through-cell direction).
For broadcast, we follow the default behaviour from SpatialMethod.
"""
def __init__(self, options=None):
super().__init__(options)
pybamm.citations.register("Gustafsson2020")
def build(self, mesh):
super().build(mesh)
# add npts_for_broadcast to mesh domains for this particular discretisation
for dom in mesh.keys():
mesh[dom].npts_for_broadcast_to_nodes = mesh[dom].npts
[docs]
def spatial_variable(self, symbol):
"""
Creates a discretised spatial variable compatible with
the FiniteElement method.
Parameters
----------
symbol : :class:`pybamm.SpatialVariable`
The spatial variable to be discretised.
Returns
-------
:class:`pybamm.Vector`
Contains the discretised spatial variable
"""
symbol_mesh = self.mesh
if symbol.name == "y":
entries = symbol_mesh["current collector"].coordinates[0, :][:, np.newaxis]
elif symbol.name == "z":
entries = symbol_mesh["current collector"].coordinates[1, :][:, np.newaxis]
else:
raise pybamm.GeometryError(
f"Spatial variable must be 'y' or 'z' not {symbol.name}"
)
return pybamm.Vector(entries, domains=symbol.domains)
[docs]
def gradient(self, symbol, discretised_symbol, boundary_conditions):
"""Matrix-vector multiplication to implement the gradient operator. The
gradient w of the function u is approximated by the finite element method
using the same function space as u, i.e. we solve w = grad(u), which
corresponds to the weak form w*v*dx = grad(u)*v*dx, where v is a suitable
test function.
Parameters
----------
symbol: :class:`pybamm.Symbol`
The symbol that we will take the Laplacian of.
discretised_symbol: :class:`pybamm.Symbol`
The discretised symbol of the correct size
boundary_conditions : dict
The boundary conditions of the model
({symbol: {"negative tab": neg. tab bc, "positive tab": pos. tab bc}})
Returns
-------
:class: `pybamm.Concatenation`
A concatenation that contains the result of acting the discretised
gradient on the child discretised_symbol. The first column corresponds
to the y-component of the gradient and the second column corresponds
to the z component of the gradient.
"""
skfem = import_optional_dependency("skfem")
domain = symbol.domain[0]
mesh = self.mesh[domain]
# get gradient matrix
grad_y_matrix, grad_z_matrix = self.gradient_matrix(symbol, boundary_conditions)
# assemble mass matrix (there is no need to zero out entries here, since
# boundary conditions are already accounted for in the governing pde
# for the symbol we are taking the gradient of. we just want to get the
# correct weights)
@skfem.BilinearForm
def mass_form(u, v, w):
return u * v
mass = skfem.asm(mass_form, mesh.basis)
# we need the inverse
mass_inv = pybamm.Matrix(inv(csc_matrix(mass)))
# compute gradient
grad_y = mass_inv @ (grad_y_matrix @ discretised_symbol)
grad_z = mass_inv @ (grad_z_matrix @ discretised_symbol)
# create concatenation
grad = pybamm.Concatenation(
grad_y, grad_z, check_domain=False, concat_fun=np.hstack
)
grad.copy_domains(symbol)
return grad
[docs]
def gradient_squared(self, symbol, discretised_symbol, boundary_conditions):
"""Multiplication to implement the inner product of the gradient operator
with itself. See :meth:`pybamm.SpatialMethod.gradient_squared`
"""
grad = self.gradient(symbol, discretised_symbol, boundary_conditions)
grad_y, grad_z = grad.orphans
return grad_y**2 + grad_z**2
[docs]
def gradient_matrix(self, symbol, boundary_conditions):
"""
Gradient matrix for finite elements in the appropriate domain.
Parameters
----------
symbol: :class:`pybamm.Symbol`
The symbol for which we want to calculate the gradient matrix
boundary_conditions : dict
The boundary conditions of the model
({symbol: {"negative tab": neg. tab bc, "positive tab": pos. tab bc}})
Returns
-------
:class:`pybamm.Matrix`
The (sparse) finite element gradient matrix for the domain
"""
skfem = import_optional_dependency("skfem")
# get primary domain mesh
domain = symbol.domain[0]
mesh = self.mesh[domain]
# make form for the gradient in the y direction
@skfem.BilinearForm
def gradient_dy(u, v, w):
return u.grad[0] * v
# make form for the gradient in the z direction
@skfem.BilinearForm
def gradient_dz(u, v, w):
return u.grad[1] * v
# assemble the matrices
grad_y = skfem.asm(gradient_dy, mesh.basis)
grad_z = skfem.asm(gradient_dz, mesh.basis)
return pybamm.Matrix(grad_y), pybamm.Matrix(grad_z)
[docs]
def divergence(self, symbol, discretised_symbol, boundary_conditions):
"""Matrix-vector multiplication to implement the divergence operator.
See :meth:`pybamm.SpatialMethod.divergence`
"""
raise NotImplementedError
[docs]
def laplacian(self, symbol, discretised_symbol, boundary_conditions):
"""Matrix-vector multiplication to implement the Laplacian operator.
Parameters
----------
symbol: :class:`pybamm.Symbol`
The symbol that we will take the Laplacian of.
discretised_symbol: :class:`pybamm.Symbol`
The discretised symbol of the correct size
boundary_conditions : dict
The boundary conditions of the model
({symbol: {"negative tab": neg. tab bc, "positive tab": pos. tab bc}})
Returns
-------
:class: `pybamm.Array`
Contains the result of acting the discretised gradient on
the child discretised_symbol
"""
skfem = import_optional_dependency("skfem")
domain = symbol.domain[0]
mesh = self.mesh[domain]
stiffness_matrix = self.stiffness_matrix(symbol, boundary_conditions)
# get boundary conditions and type
neg_bc_value, neg_bc_type = boundary_conditions[symbol]["negative tab"]
pos_bc_value, pos_bc_type = boundary_conditions[symbol]["positive tab"]
# boundary load vector is adjusted to account for boundary conditions below
boundary_load = pybamm.Vector(np.zeros(mesh.npts))
# assemble boundary load if Neumann boundary conditions
if "Neumann" in [neg_bc_type, pos_bc_type]:
# make form for unit load over the boundary
@skfem.LinearForm
def unit_bc_load_form(v, w):
return v
if neg_bc_type == "Neumann":
# assemble unit load over tab
neg_bc_load = skfem.asm(unit_bc_load_form, mesh.negative_tab_basis)
# value multiplied by weights
boundary_load = boundary_load + neg_bc_value * pybamm.Vector(neg_bc_load)
elif neg_bc_type == "Dirichlet":
# set Dirichlet value at facets corresponding to tab
neg_bc_load = np.zeros(mesh.npts)
neg_bc_load[mesh.negative_tab_dofs] = 1
boundary_load = boundary_load + neg_bc_value * pybamm.Vector(neg_bc_load)
else:
raise ValueError(
f"boundary condition must be Dirichlet or Neumann, not '{neg_bc_type}'"
)
if pos_bc_type == "Neumann":
# assemble unit load over tab
pos_bc_load = skfem.asm(unit_bc_load_form, mesh.positive_tab_basis)
# value multiplied by weights
boundary_load = boundary_load + pos_bc_value * pybamm.Vector(pos_bc_load)
elif pos_bc_type == "Dirichlet":
# set Dirichlet value at facets corresponding to tab
pos_bc_load = np.zeros(mesh.npts)
pos_bc_load[mesh.positive_tab_dofs] = 1
boundary_load = boundary_load + pos_bc_value * pybamm.Vector(pos_bc_load)
else:
raise ValueError(
f"boundary condition must be Dirichlet or Neumann, not '{pos_bc_type}'"
)
return -stiffness_matrix @ discretised_symbol + boundary_load
[docs]
def stiffness_matrix(self, symbol, boundary_conditions):
"""
Laplacian (stiffness) matrix for finite elements in the appropriate domain.
Parameters
----------
symbol: :class:`pybamm.Symbol`
The symbol for which we want to calculate the Laplacian matrix
boundary_conditions : dict
The boundary conditions of the model
({symbol: {"negative tab": neg. tab bc, "positive tab": pos. tab bc}})
Returns
-------
:class:`pybamm.Matrix`
The (sparse) finite element stiffness matrix for the domain
"""
skfem = import_optional_dependency("skfem")
# get primary domain mesh
domain = symbol.domain[0]
mesh = self.mesh[domain]
# make form for the stiffness
@skfem.BilinearForm
def stiffness_form(u, v, w):
return sum(u.grad * v.grad)
# assemble the stifnness matrix
stiffness = skfem.asm(stiffness_form, mesh.basis)
# get boundary conditions and type
try:
_, neg_bc_type = boundary_conditions[symbol]["negative tab"]
_, pos_bc_type = boundary_conditions[symbol]["positive tab"]
except KeyError as error:
raise pybamm.ModelError(
f"No boundary conditions provided for symbol `{symbol}``"
) from error
# adjust matrix for Dirichlet boundary conditions
if neg_bc_type == "Dirichlet":
self.bc_apply(stiffness, mesh.negative_tab_dofs)
if pos_bc_type == "Dirichlet":
self.bc_apply(stiffness, mesh.positive_tab_dofs)
return pybamm.Matrix(stiffness)
[docs]
def integral(self, child, discretised_child, integration_dimension):
"""Vector-vector dot product to implement the integral operator.
See :meth:`pybamm.SpatialMethod.integral`
"""
# Calculate integration vector
integration_vector = self.definite_integral_matrix(child)
out = integration_vector @ discretised_child
return out
[docs]
def definite_integral_matrix(self, child, vector_type="row"):
"""
Matrix for finite-element implementation of the definite integral over
the entire domain
.. math::
I = \\int_{\\Omega}\\!f(s)\\,dx
for where :math:`\\Omega` is the domain.
Parameters
----------
child : :class:`pybamm.Symbol`
The symbol being integrated
vector_type : str, optional
Whether to return a row or column vector (default is row)
Returns
-------
:class:`pybamm.Matrix`
The finite element integral vector for the domain
"""
skfem = import_optional_dependency("skfem")
# get primary domain mesh
domain = child.domain[0]
mesh = self.mesh[domain]
# make form for the integral
@skfem.LinearForm
def integral_form(v, w):
return v
# assemble
vector = skfem.asm(integral_form, mesh.basis)
if vector_type == "row":
return pybamm.Matrix(vector[np.newaxis, :])
elif vector_type == "column":
return pybamm.Matrix(vector[:, np.newaxis])
[docs]
def indefinite_integral(self, child, discretised_child, direction):
"""Implementation of the indefinite integral operator. The
input discretised child must be defined on the internal mesh edges.
See :meth:`pybamm.SpatialMethod.indefinite_integral`
"""
raise NotImplementedError
[docs]
def boundary_integral(self, child, discretised_child, region):
"""Implementation of the boundary integral operator.
See :meth:`pybamm.SpatialMethod.boundary_integral`
"""
# Calculate integration vector
integration_vector = self.boundary_integral_vector(child.domain, region=region)
out = integration_vector @ discretised_child
out.clear_domains()
return out
[docs]
def boundary_integral_vector(self, domain, region):
"""A node in the expression tree representing an integral operator over the
boundary of a domain
.. math::
I = \\int_{\\partial a}\\!f(u)\\,du,
where :math:`\\partial a` is the boundary of the domain, and
:math:`u\\in\\text{domain boundary}`.
Parameters
----------
domain : list
The domain(s) of the variable in the integrand
region : str
The region of the boundary over which to integrate. If region is
`entire` the integration is carried out over the entire boundary. If
region is `negative tab` or `positive tab` then the integration is only
carried out over the appropriate part of the boundary corresponding to
the tab.
Returns
-------
:class:`pybamm.Matrix`
The finite element integral vector for the domain
"""
skfem = import_optional_dependency("skfem")
# get primary domain mesh
mesh = self.mesh[domain[0]]
# make form for the boundary integral
@skfem.LinearForm
def integral_form(v, w):
return v
if region == "entire":
# assemble over all facets
integration_vector = skfem.asm(integral_form, mesh.facet_basis)
elif region == "negative tab":
# assemble over negative tab facets
integration_vector = skfem.asm(integral_form, mesh.negative_tab_basis)
elif region == "positive tab":
# assemble over positive tab facets
integration_vector = skfem.asm(integral_form, mesh.positive_tab_basis)
return pybamm.Matrix(integration_vector[np.newaxis, :])
[docs]
def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
"""
Returns the average value of the symbol over the negative tab ("negative tab")
or the positive tab ("positive tab") in the Finite Element Method.
Overwrites the default :meth:`pybamm.SpatialMethod.boundary_value`
"""
# Return average value on the negative tab for "negative tab" and positive tab
# for "positive tab"
if isinstance(symbol, pybamm.BoundaryValue):
# get integration_vector
if symbol.side == "negative tab":
region = "negative tab"
elif symbol.side == "positive tab":
region = "positive tab"
domain = symbol.children[0].domain
integration_vector = self.boundary_integral_vector(domain, region=region)
# divide integration weights by (numerical) tab width to give average value
boundary_val_vector = integration_vector / (
integration_vector @ pybamm.Vector(np.ones(integration_vector.shape[1]))
)
elif isinstance(symbol, pybamm.BoundaryGradient):
raise NotImplementedError
# Return boundary value with domain given by symbol
boundary_value = boundary_val_vector @ discretised_child
boundary_value.copy_domains(symbol)
return boundary_value
[docs]
def mass_matrix(self, symbol, boundary_conditions):
"""
Calculates the mass matrix for the finite element method.
Parameters
----------
symbol: :class:`pybamm.Variable`
The variable corresponding to the equation for which we are
calculating the mass matrix.
boundary_conditions : dict
The boundary conditions of the model
({symbol: {"negative tab": neg. tab bc, "positive tab": pos. tab bc}})
Returns
-------
:class:`pybamm.Matrix`
The (sparse) mass matrix for the spatial method.
"""
return self.assemble_mass_form(symbol, boundary_conditions)
[docs]
def boundary_mass_matrix(self, symbol, boundary_conditions):
"""
Calculates the mass matrix for the finite element method assembled
over the boundary.
Parameters
----------
symbol: :class:`pybamm.Variable`
The variable corresponding to the equation for which we are
calculating the mass matrix.
boundary_conditions : dict
The boundary conditions of the model
({symbol: {"negative tab": neg. tab bc, "positive tab": pos. tab bc}})
Returns
-------
:class:`pybamm.Matrix`
The (sparse) mass matrix for the spatial method.
"""
return self.assemble_mass_form(symbol, boundary_conditions, region="boundary")
[docs]
def bc_apply(self, M, boundary, zero=False):
"""
Adjusts the assembled finite element matrices to account for boundary conditions.
Parameters
----------
M: :class:`scipy.sparse.coo_matrix`
The assembled finite element matrix to adjust.
boundary: :class:`numpy.array`
Array of the indices which correspond to the boundary.
zero: bool, optional
If True, the rows of M given by the indicies in boundary are set to zero.
If False, the diagonal element is set to one. default is False.
"""
row = np.arange(0, np.size(boundary))
if zero:
data = np.zeros_like(boundary)
else:
data = np.ones_like(boundary)
bc_matrix = csr_matrix(
(data, (row, boundary)), shape=(np.size(boundary), np.shape(M)[1])
)
M[boundary] = bc_matrix
