Skip to content
Snippets Groups Projects
Commit 93a49709 authored by MALOU THIBAULT's avatar MALOU THIBAULT
Browse files

update the documentation of the adjoint_convection_diffusion_2D module to the... 

update the documentation of the adjoint_convection_diffusion_2D module to the Docstring NumpyDoc format
parent 9661d1ce
No related branches found
No related tags found
No related merge requests found
Pipeline #221171 passed
Stage: test
    Stage: build
      Hide whitespace changes
      Inline Side-by-side
      src/source_localization/adjoint_convection_diffusion_2D.py
      + 193
      − 116
      View file @ 93a49709
      ......@@ -14,78 +14,119 @@ from pheromone_dispersion.reaction_operator import Reaction
      class AdjointDiffusionConvectionReaction2DEquation:
      """
      Class containing the adjoint model of the 2D diffusion-convection-reaction PDE model and its solvers
      r"""
      Class containing the adjoint model of the 2D diffusion-convection-reaction PDE model :
      .. math::
      \partial_tc^* + \nabla\cdot(\mathbf{K}^T\nabla c^*)+\nabla(\vec{u}c^*)-(\nabla.\vec{u})c^*-\tau_{loss}c^*
      = \left(\frac{dm}{dc}(c(s))\right)^*\cdot2\mathbf{R}^{-1}\left(m(c(s))-m^{obs}\right)
      ~\forall (x,y)\in\Omega~\forall t\in[0;T[
      with the final and boundary conditions:
      - a null final condition :math:`c(x,y,t=T)=0~\forall (x,y)\in\Omega`,
      - a null diffusive flux :math:`\mathbf{K}^T\nabla c\cdot\vec{n}=0~\forall (x,y)\in\partial\Omega`,
      - a null outgoing convective flux
      :math:`\vec{u}c\cdot\vec{n}=0~\forall (x,y)\in\partial\Omega\cap \{(x,y)|\vec{u}(x,y,t)\cdot\vec{n}>0\}~\forall t\in]0;T]`
      with :math:`\vec{n}` the outgoing normal vector,
      and its solvers.
      Attributes
      ----------
      msh : ~pheromone_dispersion.geom.MeshRect2D
      The geometry of the domain.
      A_adjoint : ~pheromone_dispersion.advection_operator.AdvectionAdjoint
      The adjoint of the convection linear operator
      :math:`A^*:c^*\mapsto - \nabla\cdot(\vec{u}c^*)+(\nabla.\vec{u})c^*~\forall (x,y)\in\Omega~\forall t\in]0;T]`
      with the boundary condition
      :math:`\vec{u}c\cdot\vec{n}=0~\forall (x,y)\in\partial\Omega\cap \{(x,y)|\vec{u}(x,y,t)\cdot\vec{n}>0\}~\forall t\in]0;T]`.
      A : ~pheromone_dispersion.advection_operator.Advection
      The advection linear operator with the flux part of its adjoint operator
      :math:`A^*_f:c^*\mapsto - \nabla\cdot(\vec{u}c^*)~\forall (x,y)\in\Omega~\forall t\in]0;T]`
      with :math:`\vec{u}c\cdot\vec{n}=0~\forall (x,y)\in\partial\Omega\cap \{(x,y)|\vec{u}(x,y,t)\cdot\vec{n}>0\}~\forall t\in]0;T]`.
      This operator is useful for semi-implicit numerical schemes when the flux part of adjoint operator
      :math:`- \nabla\cdot(\vec{u}c^*)` is treated explicitly
      and the divergence part of the adjoint operator
      :math:`(\nabla.\vec{u})c^*` is treated semi-implicitly.
      D : ~pheromone_dispersion.diffusion_operator.Diffusion
      The diffusion linear operator with its adjoint operator
      :math:`D^*:c\mapsto \nabla\cdot(\mathbf{K}^T\nabla c)~\forall (x,y)\in\Omega~\forall t\in]0;T]`
      R : ~pheromone_dispersion.reaction_operator.Reaction
      The reaction linear operator with its adjoint operator :math:`R^*:c\mapsto \tau_{loss}c~\forall (x,y)\in\Omega~\forall t\in]0;T]`.
      ID : ~pheromone_dispersion.identity_operator.Id
      The identity operator :math:`Id:c\mapsto c~\forall (x,y)\in\Omega~\forall t\in]0;T]`.
      negative_divU_advection_term : ~pheromone_dispersion.reaction_operator.Reaction
      The negative part of the divergence part of the adjoint operator of the advection operator
      :math:`A^*_{div^-}:c^*\mapsto (\nabla.\vec{u})_-c^*`.
      This operator is useful for semi-implicit numerical schemes
      when this part of the divergence part of the adjoint operator
      :math:`(\nabla.\vec{u})c^*` is treated implicitly.
      positive_divU_advection_term : ~pheromone_dispersion.reaction_operator.Reaction
      The positive part of the divergence part of the adjoint operator of the advection operator
      :math:`A^*_{div^+}:c^*\mapsto (\nabla.\vec{u})_+c^*`.
      This operator is useful for semi-implicit numerical schemes
      when this part of the divergence part of the adjoint operator
      :math:`(\nabla.\vec{u})c^*` is treated explicitly.
      implemented_solver_type: list of str
      The list of keywords of the implemented types of time discretization.
      The implemented types of time discretization are:
      `['implicit', 'semi-implicit', 'implicit with matrix inversion', 'implicit with stationnary matrix inversion']`
      time_discretization: str
      The keyword identifying the type of time discretization.
      The type of time discretization should be the same the one
      used for the pheromone propagation model implemented in the class
      :class:`~pheromone_dispersion.convection_diffusion_2D.DiffusionConvectionReaction2DEquation`
      to ensure that the adjoint of the discretization of the pheromone propagation model
      coincides with the discretization of the adjoint model.
      tol_inversion : float
      Tolerance for the inversion estimation algorithm called at each time step.
      transpose_inv_matrix_semi_implicit_scheme : ~numpy.ndarray or None
      Inverse matrix of the implicit part matrix of the semi-implicit time discretization.
      Initialized to `None`.
      transpose_inv_matrix_implicit_scheme : ~numpy.ndarray or None
      Inverse matrix of the implicit part matrix of the implicit time discretization.
      Initialized to `None`.
      Notes
      -----
      In a future version of the module, the adjoint operator of the convection operator
      :math:`A^*:c^*\mapsto - \nabla\cdot(\vec{u}c^*)+(\nabla.\vec{u})c^*` will be implemented
      in the :class:`~pheromone_dispersion.advection_operator.Advection` class
      (instead of implemented in the :class:`~pheromone_dispersion.advection_operator.AdvectionAdjoint` class).
      The flux part of this adjoint operator :math:`A^*_f:c^*\mapsto - \nabla\cdot(\vec{u}c^*)` will be implemented
      in a separate class
      (and no longer in the :class:`~pheromone_dispersion.advection_operator.Advection` class)
      """
      def __init__(self, U, K, coeff_depot, msh, time_discretization='semi-implicit', tol_inversion=1e-14):
      """
      Instanciation of the class.
      - input:
      * msh:
      | object of the class mesh_rect_2D
      * U:
      | object of the class velocity
      * K:
      | object of the class diffusion_tensor
      * coeff_depot:
      | array of shape (msh.y, msh.x),
      | contains the deposition coefficient
      * obs:
      | object of the class Obs,
      | contains the observation operator and its derivative wrt the state variable
      * msh:
      | object of the class mesh_rect_2D
      * time_discretization:
      | string, 'semi-implicit' by default
      | describes the type of time discretization
      * tol_inversion:
      | float, optional and 1e-14 by default
      | tolerance for the inversion estimation algorithm called at each time step
      - attributes:
      * msh:
      | object of the class mesh_rect_2D
      * A:
      | object of the class Advection,
      | if the time discretization is semi-implicit, contains the convection linear operator
      * A_adjoint:
      | object of the class Advection,
      | if the time discretization is implicit, contains the adjoint of the convection linear operator
      * D:
      | object of the class Diffusion,
      | contains the diffusion linear operator
      * R:
      | object of the class Reaction,
      | contains the reaction linear operator
      * ID:
      | object of the class Id,
      | contains the identity operator
      * obs:
      | object of the class Obs,
      | contains the observation operator and its derivative wrt the state variable
      * negative_divU_advection_term:
      | object of the class Reaction,
      | if the time discretization is semi-implicit, contains negative part of the reaction term in div U
      | coming from the splitting of the adjoint of the advection term into a advection flux part and a reaction part
      * positive_divU_advection_term:
      | object of the class Reaction,
      | if the time discretization is semi-implicit, contains positive part of the reaction term in div U
      | coming from the splitting of the adjoint of the advection term into a advection flux part and a reaction part
      * time_discretization:
      | string,
      | describes the type of time discretization
      * tol_inversion:
      | float, optional and 1e-14 by default
      | tolerance for the inversion estimation algorithm called at each time step
      * transpose_inv_matrix_semi_implicit_scheme:
      | array of shape (msh.y.size * msh.x.size, msh.y.size * msh.x.size) or None (by default)
      | contains the transpose of the inverse of the implicit matrix of the semi-implicit discretization with matrix inversion
      | is set to None if an other type of discretization is used.
      * transpose_inv_matrix_implicit_scheme:
      | array of shape (msh.y.size * msh.x.size, msh.y.size * msh.x.size) or None (by default)
      | contains the transpose of the inverse matrix of the implicit discretization with matrix inversion
      | is set to None if an other type of discretization is used.
      r"""
      Constructor method
      Parameters
      ----------
      msh: ~pheromone_dispersion.geom.MeshRect2D
      The geometry of the domain.
      U: ~pheromone_dispersion.velocity.Velocity
      The wind field :math:`\vec{u}(x,y,t)`.
      K: ~pheromone_dispersion.diffusion_tensor.DiffusionTensor
      The diffusion tensor :math:`\mathbf{K}(x,y,t)`.
      coeff_depot: ~numpy.ndarray
      The deposition coefficient :math:`\tau_{loss}(x,y)`.
      obs : ~source_localization.obs.Obs
      Object containing all the features related to the observations and estimation of the observed variables
      time_discretization: str, default: 'semi-implicit'
      The keyword identifying the type of time discretization.
      tol_inversion: float, optional, default: 1e-14
      Tolerance for the inversion estimation algorithm called at each time step.
      Raises
      ------
      ValueError
      if the type of discretization is not implemented,
      i.e. if :attr:`time_discretization` is not in :attr:`implemented_solver_type`
      """
      self.msh = msh
      ......@@ -93,14 +134,13 @@ class AdjointDiffusionConvectionReaction2DEquation:
      self.R = Reaction(coeff_depot, msh)
      self.Id = Id(msh)
      self.time_discretization = time_discretization
      implemented_solver_type = [
      self.implemented_solver_type = [
      'implicit',
      'semi-implicit',
      'implicit with matrix inversion',
      'semi-implicit with matrix inversion',
      'implicit with stationnary matrix inversion',
      ]
      if self.time_discretization not in implemented_solver_type:
      if self.time_discretization not in self.implemented_solver_type:
      raise ValueError("The given time discretization is not implemented.")
      if time_discretization == 'semi-implicit':
      ......@@ -121,20 +161,37 @@ class AdjointDiffusionConvectionReaction2DEquation:
      def init_inverse_matrix(self, path_to_matrix=None, matrix_file_name=None):
      """
      Initialize the transpose of the inverse matrix of the linear operator corresponding to the implicit part of the time scheme.
      - input:
      * path_to_matrix:
      | string, optional and None by default
      | path where save or load the inverse matrix.
      | if the path is not provided, it will be './data'
      * matrix_file_name:
      | string, optional and None by default
      | Matrix name.
      | If not provided, it will be be 'inv_matrix_**_scheme' with ** depending on the scheme.
      - do:
      * If the inverse matrix has never been computed, it raises an exception.
      * Else, it loads the previously computed inverse matrix and store its transpose.
      Initialize the attribute :attr:`transpose_inv_matrix_implicit_scheme`
      and :attr:`transpose_inv_matrix_semi_implicit_scheme` if needed.
      If the inverse matrix of the implicit part matrix of the time discretization
      of the associated pheromone propagation model has already been computed,
      it loads the previously computed inverse matrix and compute the adjoint.
      Otherwise, raises an error as the matrix should have been computed
      when initializing the pheromone propagation model as object
      of the class :class:`~pheromone_dispersion.convection_diffusion_2D.DiffusionConvectionReaction2DEquation`
      using the method :meth:`~pheromone_dispersion.convection_diffusion_2D.DiffusionConvectionReaction2DEquation.init_inverse_matrix`.
      Parameters
      ----------
      path_to_matrix: str, optional, default: None
      Path where to save or load the inverse matrix.
      If not provided, set to `'./data'`.
      matrix_file_name: str, optional, default: None
      Name of the file where the matrix stored or will be saved.
      If not provided, set to `'inv_matrix_**_scheme'`
      with ** either `'implicit'` or `'semi_implicit'` depending on the time discretization.
      Raises
      ------
      ValueError
      if no files correspond to the provided path and file name.
      Notes
      -----
      This method is usefull only if :attr:`time_discretization` is
      either `'semi-implicit with matrix inversion'` or `'implicit with stationnary matrix inversion'`.
      Otherwise, the linear system to solve at each time steps is solved using conjugate gradient or GMRES algorithm,
      and the attributes :attr:`transpose_inv_matrix_implicit_scheme` and :attr:`transpose_inv_matrix_semi_implicit_scheme` are not used.
      """
      if path_to_matrix is None:
      path_to_matrix = './data'
      ......@@ -162,17 +219,27 @@ class AdjointDiffusionConvectionReaction2DEquation:
      self.transpose_inv_matrix_implicit_scheme = np.transpose(np.load(Path(path_to_matrix) / matrix_file_name))
      def at_current_time(self, tc):
      """
      Update the linear operators of the PDE at a given time,
      i.e. update the attributes D and S of the class
      if the time discretization is semi-implicit, A, positive_divU_advection_term, negative_divU_advection_term,
      if the time discretization is implicit, A_adjoint,
      and the derivative of the observation operator of the attribute obs.
      - input:
      * tc:
      | float
      | the current time
      r"""
      Update the attributes :attr:`A_adjoint`, :attr:`D`,
      :attr:`A`, :attr:`negative_divU_advection_term`
      and :attr:`positive_divU_advection_term` at a given time.
      Parameters
      ----------
      tc : float or integer
      The current time.
      Notes
      -----
      Updates the attributes :attr:`A_adjoint`, :attr:`D` and :attr:`A` and their own attributes
      using the method :meth:`~pheromone_dispersion.advection_operator.Advection.at_current_time` of resp.
      the class :class:`~pheromone_dispersion.advection_operator.AdvectionAdjoint`,
      the class :class:`~pheromone_dispersion.diffusion_operator.Diffusion` and
      the class :class:`~pheromone_dispersion.advection_operator.Advection`,
      and the attributes :attr:`negative_divU_advection_term`
      and :attr:`positive_divU_advection_term`
      using the method :meth:`~pheromone_dispersion.reaction_operator.Reaction.update_reaction_coeff`
      of the class :class:`~pheromone_dispersion.reaction_operator.Reaction`.
      """
      self.D.at_current_time(tc)
      ......@@ -187,25 +254,35 @@ class AdjointDiffusionConvectionReaction2DEquation:
      self.A_adjoint.at_current_time(tc)
      def solver(self, adjoint_derivative_obs_operator, cost, display_flag=True):
      """
      solve the PDE on the whole time window
      - input:
      * adjoint_derivative_obs_operator:
      | callable,
      | method that returns the adjoint of the derivative of the observation operator
      | given a time and an element of the space of the state variable
      * cost:
      | object of the class Cost,
      | contains especially a method that compute the gradient of the objectif
      | wrt the observed variable for the current optimization iteration and estimation of the observed variable
      * display_flag:
      | boolean, True by default,
      | if True, print the evolution in time of the solver
      - output:
      * p:
      | numpy array of shape (msh.t_array.size * msh.y.size * msh.x.size,),
      | contains the adjoint state everywhere and at all time step concatenated in a vector
      r"""
      Compute the adjoint state :math:`c^*(x,y,t)` by solving the adjoint model on the whole time window.
      Parameters
      ----------
      adjoint_derivative_obs_operator: callable
      Function that returns the adjoint of the derivative of the observation operator map
      :math:`\left(\left(\frac{dm}{dc}\right)^*\cdot\delta m\right) (x,y)`
      given the current time :math:`t` and
      an element of the space of the observed variable :math:`\delta m`.
      cost: ~source_localization.cost.Cost
      The cost function,
      containing especially the gap between the observation data
      and the prediction of the observation variable
      :math:`2\mathbf{R}^{-1}\left(m(c(s))-m^{obs}\right)`.
      display_flag: bool, default: True
      If `True`, print the evolution of the solver through the time iterations.
      Returns
      -------
      p_out: ~numpy.ndarray
      The adjoint state :math:`c^*(x,y,t)`.
      Notes
      -----
      The output :math:`c^*(x,y,t)` has been raveled into
      a (:attr:`msh.t_array.size` * :attr:`msh.y.size` * :attr:`msh.x.size`,)-shape array
      to match the format of the attribute :attr:`~source_localization.control.Control.value`
      of the class :class:`~source_localization.control.Control`.
      """
      # initialization of the unknown variable at the final time and of the output
      ......
      0% or .
      You are about to add 0 people to the discussion. Proceed with caution.
      Finish editing this message first!
      Please register or to comment