import numpy as np
from scipy.sparse.linalg import LinearOperator as LinOp
"""
Module that contains the implementation of the linear operator of the diffusion term of the pheromone propagation model.
"""
[docs]
class Diffusion(LinOp):
r"""
Class containing the diffusion term linear operator
:math:`D:c\mapsto -\nabla\cdot(\mathbf{K}\nabla c)~\forall (x,y)\in\Omega~\forall t\in]0;T]`
with the related boundary condition
:math:`\mathbf{K}\nabla c\cdot\vec{n}=0~\forall (x,y)\in\partial\Omega`.
The implementation of this linear operator is a subclass of the :py:class:`~scipy.sparse.linalg.LinearOperator` class.
Attributes
----------
msh: ~pheromone_dispersion.geom.MeshRect2D
The geometry of the domain.
K: ~pheromone_dispersion.diffusion_tensor.DiffusionTensor
The diffusion tensor :math:`\mathbf{K}(x,y,t)`.
shape : tuple of int
Shape of the matrix of the linear operator.
The shape is (:attr:`msh.y.size` * :attr:`msh.x.size`, :attr:`msh.y.size` * :attr:`msh.x.size`)
dtype : ~numpy.dtype
Data type of the elements of the matrix of the linear operator.
The type is `float64`
"""
[docs]
def __init__(self, K, msh):
r"""
Constructor method
Parameters
----------
K: ~pheromone_dispersion.diffusion_tensor.DiffusionTensor
The diffusion tensor :math:`\mathbf{K}(x,y,t)`.
msh : ~pheromone_dispersion.geom.MeshRect2D
The geometry of the domain.
"""
self.msh = msh
self.K = K
self.shape = (self.msh.y.size * self.msh.x.size, self.msh.y.size * self.msh.x.size)
self.dtype = np.dtype(float)
[docs]
def at_current_time(self, tc):
r"""
Update the diffusion tensor at a given time.
Parameters
----------
tc : float or integer
The current time.
Notes
-----
Updates the diffusion tensor :attr:`K` and its own attributes
using the method :meth:`~pheromone_dispersion.diffusion_tensor.DiffusionTensor.at_current_time`
of the class :class:`~pheromone_dispersion.diffusion_tensor.DiffusionTensor`.
"""
self.K.at_current_time(tc)
[docs]
def _matvec(self, x_out):
r"""
Compute the image (matrix-vector product) of :math:`D:c\mapsto -\nabla\cdot(\mathbf{K}\nabla c)`
with the related boundary condition
for a given concentration map :math:`c` at the current time.
Parameters
----------
x_out : ~numpy.ndarray
The map of concentration of pheromones :math:`c(x,y)` at a given time :math:`t` raveled into a vector.
Returns
-------
~numpy.ndarray
Array containing the image :math:`-\nabla\cdot(\mathbf{K}\nabla c)`.
Notes
-----
The input :math:`c` has to be raveled into a (:attr:`msh.y.size` * :attr:`msh.x.size`,)-shape array
to match the format of the :class:`~scipy.sparse.linalg.LinearOperator` class.
The same way, the ouput can be reshape into a (:attr:`msh.y.size`, :attr:`msh.x.size`)-shape array
to get the map :math:`Dc(x,y)` at the current time :math:`t`
For now, the numerical scheme is implemented for diagonal diffusion tensor,
i.e. for unidirectional wind field
or for isotropic diffusion tensor (:math:`K_u=K_{u^\perp}`)
"""
# To do:
# implemente for non-diagonal diffusion tensor
# phantom cells should be added that satisfies the boundary conditions
# a _rmatvec method shoud be added for non-diagonal diffusion tensor
# reshape the vector (of (msh.x.size*msh.y.size, )) containg the concentration into a matrix (of shape (msh.y.size, msh.x.size))
x = x_out.reshape((self.msh.y.size, self.msh.x.size))
# Computation of the gradient of concentration on the vertical interfaces between the cells
grad_at_vertical_interface = np.zeros((self.msh.y.size, self.msh.x_vertical_interface.size, 2))
grad_at_vertical_interface[:, 1:-1, 0] = (x[:, 1:] - x[:, :-1]) / self.msh.dx
grad_at_vertical_interface[1:-1, 1:-1, 1] = (x[2:, :-1] + x[2:, 1:] - (x[:-2, :-1] + x[:-2, 1:])) / (4 * self.msh.dy)
# adding phantom cells, TO BE IMPROVED
grad_at_vertical_interface[0, 1:-1, 1] = (x[1, :-1] + x[1, 1:] - (x[0, :-1] + x[0, 1:])) / (4 * self.msh.dy)
grad_at_vertical_interface[-1, 1:-1, 1] = (x[-1, :-1] + x[-1, 1:] - (x[-2, :-1] + x[-2, 1:])) / (4 * self.msh.dy)
# Computation of the gradient concentration on the horizontal interfaces between the cells
grad_at_horizontal_interface = np.zeros((self.msh.y_horizontal_interface.size, self.msh.x.size, 2))
grad_at_horizontal_interface[1:-1, 1:-1, 0] = (x[1:, 2:] + x[:-1, 2:] - (x[1:, :-2] + x[:-1, :-2])) / (4 * self.msh.dx)
grad_at_horizontal_interface[1:-1, :, 1] = (x[1:, :] - x[:-1, :]) / self.msh.dy
# adding phantom cells, TO BE IMPROVED
grad_at_horizontal_interface[1:-1, 0, 0] = (x[1:, 1] + x[:-1, 1] - (x[1:, 0] + x[:-1, 0])) / (4 * self.msh.dx)
grad_at_horizontal_interface[1:-1, -1, 0] = (x[1:, -1] + x[:-1, -1] - (x[1:, -2] + x[:-1, -2])) / (4 * self.msh.dx)
# Computation of the horizontal flux going through the vertical interfaces
flux_horizontal = np.zeros((self.msh.y.size, self.msh.x_vertical_interface.size))
flux_horizontal[:, :] = (
self.K.at_vertical_interface[:, :, 0, 0] * grad_at_vertical_interface[:, :, 0]
+ self.K.at_vertical_interface[:, :, 0, 1] * grad_at_vertical_interface[:, :, 1]
)
# Boundary Condition on the horizontal flux going through the vertical interfaces
flux_horizontal[:, -1] = 0
flux_horizontal[:, 0] = 0
# Computation of the vertical flux going through the horizontal interfaces
flux_vertical = np.zeros((self.msh.y_horizontal_interface.size, self.msh.x.size))
flux_vertical[:, :] = (
self.K.at_horizontal_interface[:, :, 1, 0] * grad_at_horizontal_interface[:, :, 0]
+ self.K.at_horizontal_interface[:, :, 1, 1] * grad_at_horizontal_interface[:, :, 1]
)
# Boundary Condition on the vertical flux going through the horizontal interfaces
flux_vertical[0, :] = 0
flux_vertical[-1, :] = 0
# Computation of the total flux over the border of the cell for a cartesian mesh
sum_faces = (
self.msh.dy * (flux_horizontal[:, 1:] - flux_horizontal[:, :-1]) + self.msh.dx * (flux_vertical[1:, :] - flux_vertical[:-1, :])
) / self.msh.mass_cell
return sum_faces.reshape((self.msh.y.size * self.msh.x.size,))
def _matmat(self, x_out):
# Compute the image (matrix-matrix product)
# using a naive method based on the _matvec method
output = np.zeros((self.shape[0], x_out.shape[1]))
for i_col, col in enumerate(x_out.T):
output[:, i_col] = self.matvec(col)
return output