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

remove the solvers of the adjoint model no longer used, and a bit of cleaning

parent 93a49709
No related branches found
No related tags found
No related merge requests found
Pipeline #221194 passed
Stage: test
    Stage: build
      Hide whitespace changes
      Inline Side-by-side
      features/2D_adjoint_model.feature
      + 0
      − 1
      View file @ 50c7cb26
      ......@@ -12,7 +12,6 @@ Scenario: Generate the AdjointDiffusionConvectionReaction2DEquation object for t
      Then the reaction term of the PDE is a linear operator
      And the diffusion term of the PDE is a linear operator
      And the convection term of the semi implicit PDE is a linear operator
      And the div parts of the convection term of the semi implicit PDE have the expected values
      And the convection term of the implicit PDE is a linear operator
      And the output of the solver has the expected shape
      And the initialization fails if the given time discretization is not implemented
      src/source_localization/adjoint_convection_diffusion_2D.py
      + 43
      − 93
      View file @ 50c7cb26
      ......@@ -6,7 +6,6 @@ import numpy as np
      from scipy.sparse.linalg import cg
      from scipy.sparse.linalg import gmres
      from pheromone_dispersion.advection_operator import Advection
      from pheromone_dispersion.advection_operator import AdvectionAdjoint
      from pheromone_dispersion.diffusion_operator import Diffusion
      from pheromone_dispersion.identity_operator import Id
      ......@@ -41,15 +40,6 @@ class AdjointDiffusionConvectionReaction2DEquation:
      :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]`
      ......@@ -57,22 +47,10 @@ class AdjointDiffusionConvectionReaction2DEquation:
      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']`
      `['implicit', 'semi-implicit', 'semi-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
      ......@@ -82,11 +60,13 @@ class AdjointDiffusionConvectionReaction2DEquation:
      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.
      transpose_inv_matrix_implicit_part: ~numpy.ndarray or None
      Inverse matrix of the implicit part matrix of the time discretization of the adjoint model.
      If the time discretizations are the same, this matrix coincides
      with the transpose of the inverse matrix of the implicit part matrix
      of the time discretization of the pheromone propagation implemented
      in the class
      :class:`~pheromone_dispersion.convection_diffusion_2D.DiffusionConvectionReaction2DEquation`.
      Initialized to `None`.
      ......@@ -94,11 +74,8 @@ class AdjointDiffusionConvectionReaction2DEquation:
      -----
      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)
      in the :class:`~pheromone_dispersion.advection_operator.Advection` class,
      instead of implemented in the :class:`~pheromone_dispersion.advection_operator.AdvectionAdjoint` class.
      """
      def __init__(self, U, K, coeff_depot, msh, time_discretization='semi-implicit', tol_inversion=1e-14):
      ......@@ -137,32 +114,20 @@ class AdjointDiffusionConvectionReaction2DEquation:
      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 self.implemented_solver_type:
      raise ValueError("The given time discretization is not implemented.")
      if time_discretization == 'semi-implicit':
      self.A = Advection(U, msh)
      pos_div_val = np.where(self.A.U.div >= 0, self.A.U.div, 0.0)
      neg_div_val = np.where(self.A.U.div < 0, self.A.U.div, 0.0)
      self.negative_divU_advection_term = Reaction(neg_div_val, msh)
      self.positive_divU_advection_term = Reaction(pos_div_val, msh)
      self.A_adjoint = AdvectionAdjoint(U, msh)
      # elif time_discretization == 'implicit':
      # self.A_adjoint = AdvectionAdjoint(U, msh)
      self.tol_inversion = tol_inversion
      self.transpose_inv_matrix_semi_implicit_scheme = None
      self.transpose_inv_matrix_implicit_scheme = None
      self.transpose_inv_matrix_implicit_part = None
      def init_inverse_matrix(self, path_to_matrix=None, matrix_file_name=None):
      """
      Initialize the attribute :attr:`transpose_inv_matrix_implicit_scheme`
      and :attr:`transpose_inv_matrix_semi_implicit_scheme` if needed.
      Initialize the attribute :attr:`transpose_inv_matrix_implicit_part` 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.
      ......@@ -191,7 +156,7 @@ class AdjointDiffusionConvectionReaction2DEquation:
      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.
      and the attribute :attr:`transpose_inv_matrix_implicit_part` is not used.
      """
      if path_to_matrix is None:
      path_to_matrix = './data'
      ......@@ -200,10 +165,9 @@ class AdjointDiffusionConvectionReaction2DEquation:
      if matrix_file_name is None:
      if self.time_discretization == 'semi-implicit with matrix inversion':
      matrix_file_name = 'inv_matrix_semi_implicit_scheme'
      if self.time_discretization == 'implicit with matrix inversion':
      matrix_file_name = 'inv_matrix_implicit_scheme'
      if self.time_discretization == 'implicit with stationnary matrix inversion':
      matrix_file_name = 'inv_matrix_implicit_scheme'
      if matrix_file_name[-4:] != '.npy':
      matrix_file_name += '.npy'
      if not (Path(path_to_matrix) / matrix_file_name).exists():
      ......@@ -212,17 +176,13 @@ class AdjointDiffusionConvectionReaction2DEquation:
      else:
      print("=== Load of the inverse of the matrix of the implicit part of the " + self.time_discretization + " scheme ===")
      if self.time_discretization == 'semi-implicit with matrix inversion':
      self.transpose_inv_matrix_semi_implicit = np.transpose(np.load(Path(path_to_matrix) / matrix_file_name))
      if self.time_discretization == 'implicit with matrix inversion':
      self.transpose_inv_matrix_implicit_scheme = np.transpose(np.load(Path(path_to_matrix) / matrix_file_name))
      self.transpose_inv_matrix_implicit_part = np.transpose(np.load(Path(path_to_matrix) / matrix_file_name))
      if self.time_discretization == 'implicit with stationnary matrix inversion':
      self.transpose_inv_matrix_implicit_scheme = np.transpose(np.load(Path(path_to_matrix) / matrix_file_name))
      self.transpose_inv_matrix_implicit_part = np.transpose(np.load(Path(path_to_matrix) / matrix_file_name))
      def at_current_time(self, tc):
      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.
      Update the attributes :attr:`A_adjoint` and :attr:`D` at a given time.
      Parameters
      ----------
      ......@@ -231,26 +191,13 @@ class AdjointDiffusionConvectionReaction2DEquation:
      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`.
      Updates the attributes :attr:`A_adjoint` and :attr:`D` their own attributes
      using the method :meth:`~pheromone_dispersion.advection_operator.AdvectionAdjoint.at_current_time` of resp.
      the class :class:`~pheromone_dispersion.advection_operator.AdvectionAdjoint` and
      the class :class:`~pheromone_dispersion.diffusion_operator.Diffusion`.
      """
      self.D.at_current_time(tc)
      if self.time_discretization == 'semi-implicit':
      self.A.at_current_time(tc, self.A.minus_U)
      pos_div_val = np.where(self.A.U.div >= 0, self.A.U.div, 0.0)
      neg_div_val = np.where(self.A.U.div < 0, self.A.U.div, 0.0)
      self.negative_divU_advection_term.update_reaction_coeff(neg_div_val)
      self.positive_divU_advection_term.update_reaction_coeff(pos_div_val)
      # elif self.time_discretization == 'implicit':
      self.A_adjoint.at_current_time(tc)
      def solver(self, adjoint_derivative_obs_operator, cost, display_flag=True):
      ......@@ -305,23 +252,29 @@ class AdjointDiffusionConvectionReaction2DEquation:
      # using a conjugate gradient method for the current time step
      if self.time_discretization == 'semi-implicit':
      p, info = cg(
      self.Id + self.msh.dt * (-self.D + self.R + self.negative_divU_advection_term),
      self.Id + self.msh.dt * (-self.D + self.R),
      p_old
      - self.msh.dt
      * (
      self.A.rmatvec(p_old)
      + self.positive_divU_advection_term(p_old)
      + adjoint_derivative_obs_operator(self.msh.t, cost.gradient_objectif_wrt_d())
      ),
      * (self.A_adjoint.rmatvec(p_old) + adjoint_derivative_obs_operator(self.msh.t, cost.gradient_objectif_wrt_d())),
      x0=p_old,
      tol=self.tol_inversion,
      )
      # p = self.transpose_inv_matrix_semi_implicit.dot(
      # p_old
      # - self.msh.dt
      # * (-self.A_adjoint.matvec(p_old) + adjoint_derivative_obs_operator(self.msh.t, cost.gradient_objectif_wrt_d()))
      # )
      # inverse the linear system resulting the implicit time discretization using a gmres method for the current time step
      # inverse the linear system resulting the semi-implicit time discretization
      # using the pre-computed inverse matrix
      elif self.time_discretization == 'semi-implicit with stationnary matrix inversion':
      RHS = p_old - self.msh.dt * (
      self.A_adjoint.rmatvec(p_old) + adjoint_derivative_obs_operator(self.msh.t, cost.gradient_objectif_wrt_d())
      )
      if not np.linalg.norm(RHS, ord=np.inf) < 1e-16:
      LHS = (self.Id + self.msh.dt * (-self.D + self.R)).matvec(p_old)
      flag_residu = not np.linalg.norm(RHS - LHS) < self.tol_inversion * np.linalg.norm(RHS)
      if flag_residu:
      p = np.dot(self.transpose_inv_matrix_implicit_part, RHS)
      else:
      p = np.zeros_like(p)
      info = 0
      # inverse the linear system resulting the implicit time discretization
      # using a gmres method for the current time step
      elif self.time_discretization == 'implicit':
      p, info = gmres(
      self.Id + self.msh.dt * (-self.D + self.R + self.A_adjoint),
      ......@@ -329,18 +282,15 @@ class AdjointDiffusionConvectionReaction2DEquation:
      x0=p_old,
      tol=self.tol_inversion,
      )
      elif self.time_discretization == 'implicit with matrix inversion':
      p = self.transpose_inv_matrix_implicit_scheme[it, :, :].dot(
      p_old - self.msh.dt * adjoint_derivative_obs_operator(self.msh.t, cost.gradient_objectif_wrt_d())
      )
      info = 0
      # inverse the linear system resulting the implicit time discretization
      # using the pre-computed inverse matrix
      elif self.time_discretization == 'implicit with stationnary matrix inversion':
      RHS = p_old - self.msh.dt * adjoint_derivative_obs_operator(self.msh.t, cost.gradient_objectif_wrt_d())
      if not np.linalg.norm(RHS, ord=np.inf) < 1e-16:
      LHS = (self.Id + self.msh.dt * (-self.D + self.R + self.A_adjoint)).matvec(p_old)
      flag_residu = not np.linalg.norm(RHS - LHS) < self.tol_inversion * np.linalg.norm(RHS)
      if flag_residu:
      p = np.dot(self.transpose_inv_matrix_implicit_scheme, RHS)
      p = np.dot(self.transpose_inv_matrix_implicit_part, RHS)
      else:
      p = np.zeros_like(p)
      info = 0
      ......
      tests/test_2D_adjoint_model.py
      + 2
      − 14
      View file @ 50c7cb26
      ......@@ -8,7 +8,6 @@ from pytest_bdd import then
      from pytest_bdd import when
      from scipy.sparse.linalg import LinearOperator as LinOp
      from pheromone_dispersion.advection_operator import Advection
      from pheromone_dispersion.advection_operator import AdvectionAdjoint
      from pheromone_dispersion.diffusion_operator import Diffusion
      from pheromone_dispersion.diffusion_tensor import DiffusionTensor
      ......@@ -115,19 +114,8 @@ def the_diffusion_term_of_the_PDE_is_a_linear_operator(PDE_semi_imp, PDE_imp):
      @then("the convection term of the semi implicit PDE is a linear operator")
      def the_convection_term_of_the_semi_implicit_PDE_is_a_linear_operator(PDE_semi_imp):
      """the convection term of the semi implicit PDE is a linear operator"""
      assert isinstance(PDE_semi_imp.negative_divU_advection_term, LinOp) and isinstance(PDE_semi_imp.negative_divU_advection_term, Reaction)
      assert isinstance(PDE_semi_imp.positive_divU_advection_term, LinOp) and isinstance(PDE_semi_imp.positive_divU_advection_term, Reaction)
      assert isinstance(PDE_semi_imp.A, LinOp) and isinstance(PDE_semi_imp.A, Advection)
      @then("the div parts of the convection term of the semi implicit PDE have the expected values")
      def the_div_parts_of_the_convection_term_of_the_semi_implicit_PDE_have_the_expected_values(PDE_semi_imp):
      """the div parts of the convection term of the semi implicit PDE have the expected values"""
      x = np.random.random((PDE_semi_imp.msh.y.shape[0] * PDE_semi_imp.msh.x.shape[0],))
      res_exp = PDE_semi_imp.positive_divU_advection_term(x)
      res_imp = PDE_semi_imp.negative_divU_advection_term(x)
      assert (res_exp == 2 * x).all()
      assert (res_imp == 0.0).all()
      assert isinstance(PDE_semi_imp.A_adjoint, LinOp)
      assert isinstance(PDE_semi_imp.A_adjoint, AdvectionAdjoint)
      @then("the convection term of the implicit PDE is a linear operator")
      ......
      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