import numpy as np
"""
Module that contains the implementation of the diffusion term diffusion tensor.
"""
[docs]
class DiffusionTensor:
r"""
Class containing the anisotropic and inhomogeneous diffusion tensor :math:`\mathbf{K}(x,y,t)`.
The diffusion tensor is given by
a diffusion coefficient :math:`K_u` in the wind direction and
a diffusion coefficient :math:`K_{u^\perp}` in the crosswind direction.
Hence, the anistropic diffusion tensor is given by
:math:`\mathbf{K} = R(\vec{u})diag(K_u,K_{u^\perp})R(\vec{u})^T`
with :math:`R(\vec{u})` the rotation matrice
of angle :math:`\theta` between the wind field and the cartesian frame.
Attributes
----------
U : ~pheromone_dispersion.velocity.Velocity
The wind field :math:`\vec{u}(x,y,t)`.
K_u : float
The diffusion coefficient in the wind direction :math:`K_u`.
K_u_t : float
The diffusion coefficient in the crosswind direction :math:`K_{u^\perp}`.
at_vertical_interface: ~numpy.ndarray
The diffusion tensor :math:`\mathbf{K}(x,y)` at the vertical interfaces of each cells of the mesh at the current time :math:`t`.
at_horizontal_interface: ~numpy.ndarray
The diffusion tensor :math:`\mathbf{K}(x,y)` at the horizontal interfaces of each cells of the mesh at the current time :math:`t`.
"""
[docs]
def __init__(self, U, K_u, K_u_t):
r"""
Constructor method
Parameters
----------
U : ~pheromone_dispersion.velocity.Velocity
The wind field :math:`\vec{u}(x,y,t)`.
K_u : float
The diffusion coefficient in the wind direction :math:`K_u`.
K_u_t : float
The diffusion coefficient in the crosswind direction :math:`K_{u^\perp}`.
"""
# To do
# add exceptions to make sure that all the inputs have the right type
self.U = U
self.K_u = K_u
self.K_u_t = K_u_t
self.diffusion_tensor_from_velocity_field()
[docs]
def diffusion_tensor_from_velocity_field(self):
r"""
Compute the attributes :attr:`at_vertical_interface`
and :attr:`at_horizontal_interface`
from the attributes :attr:`U`, :attr:`K_u` and :attr:`K_u_t`
using the rotation formula given above.
"""
# On the vertical interfaces of the mesh
norm_U_at_vertical_interface = np.linalg.norm(self.U.at_vertical_interface, axis=2)
U_at_vertical_interface = self.U.at_vertical_interface / norm_U_at_vertical_interface[:, :, None]
self.at_vertical_interface = np.zeros((U_at_vertical_interface.shape[0], U_at_vertical_interface.shape[1], 2, 2))
self.at_vertical_interface[:, :, 0, 0] = self.K_u_t
self.at_vertical_interface[:, :, 1, 1] = self.K_u_t
self.at_vertical_interface[:, :, 0, 0] += (
-(self.K_u_t - self.K_u) * U_at_vertical_interface[:, :, 0] * U_at_vertical_interface[:, :, 0]
)
self.at_vertical_interface[:, :, 0, 1] += (
-(self.K_u_t - self.K_u) * U_at_vertical_interface[:, :, 0] * U_at_vertical_interface[:, :, 1]
)
self.at_vertical_interface[:, :, 1, 0] += (
-(self.K_u_t - self.K_u) * U_at_vertical_interface[:, :, 0] * U_at_vertical_interface[:, :, 1]
)
self.at_vertical_interface[:, :, 1, 1] += (
-(self.K_u_t - self.K_u) * U_at_vertical_interface[:, :, 1] * U_at_vertical_interface[:, :, 1]
)
# On the horizontal interfaces of the mesh
norm_U_at_horizontal_interface = np.linalg.norm(self.U.at_horizontal_interface, axis=2)
U_at_horizontal_interface = self.U.at_horizontal_interface / norm_U_at_horizontal_interface[:, :, None]
self.at_horizontal_interface = np.zeros((U_at_horizontal_interface.shape[0], U_at_horizontal_interface.shape[1], 2, 2))
self.at_horizontal_interface[:, :, 0, 0] = self.K_u_t
self.at_horizontal_interface[:, :, 1, 1] = self.K_u_t
self.at_horizontal_interface[:, :, 0, 0] += (
-(self.K_u_t - self.K_u) * U_at_horizontal_interface[:, :, 0] * U_at_horizontal_interface[:, :, 0]
)
self.at_horizontal_interface[:, :, 0, 1] += (
-(self.K_u_t - self.K_u) * U_at_horizontal_interface[:, :, 0] * U_at_horizontal_interface[:, :, 1]
)
self.at_horizontal_interface[:, :, 1, 0] += (
-(self.K_u_t - self.K_u) * U_at_horizontal_interface[:, :, 0] * U_at_horizontal_interface[:, :, 1]
)
self.at_horizontal_interface[:, :, 1, 1] += (
-(self.K_u_t - self.K_u) * U_at_horizontal_interface[:, :, 1] * U_at_horizontal_interface[:, :, 1]
)
[docs]
def at_current_time(self, tc):
"""
Update all the attributes at a given time.
Parameters
----------
tc : float or integer
The current time.
Notes
-----
Updates the velocity field :attr:`U` and its own attributes
using the method :meth:`~pheromone_dispersion.velocity.Velocity.at_current_time`
of the class :class:`~pheromone_dispersion.velocity.Velocity`.
Then, recomputes the attributes :attr:`at_vertical_interface`
and :attr:`at_horizontal_interface` from the updated velocity field
using the method :meth:`diffusion_tensor_from_velocity_field`.
"""
self.U.at_current_time(tc)
self.diffusion_tensor_from_velocity_field()