IP_SSAHardavForwardProblem.hh

Go to the documentation of this file.

// Copyright (C) 2013, 2014, 2015, 2016, 2017, 2020, 2021, 2022, 2023  David Maxwell
//
// This file is part of PISM.
//
// PISM is free software; you can redistribute it and/or modify it under the
// terms of the GNU General Public License as published by the Free Software
// Foundation; either version 2 of the License, or (at your option) any later
// version.
//
// PISM is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
// details.
//
// You should have received a copy of the GNU General Public License
// along with PISM; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 
#ifndef IP_SSAHARDAVFORWARDPROBLEM_HH_HAD68BK7
#define IP_SSAHARDAVFORWARDPROBLEM_HH_HAD68BK7
 
#include "pism/stressbalance/ssa/SSAFEM.hh"
#include "pism/inverse/IPDesignVariableParameterization.hh"
#include "pism/util/petscwrappers/KSP.hh"
#include "pism/util/petscwrappers/Mat.hh"
 
namespace pism {
namespace inverse {
 
//! Implements the forward problem of the map taking \f$\tau_c\f$ to the corresponding solution of the %SSA.
/*! The class SSAFEM solves the %SSA, and the solution depends on a large number of parameters.  Considering
  all of these to be fixed except for \f$\tau_c\f$, we obtain a map \f$F_{\rm SSA}\f$ taking
  \f$\tau_c\f$ to the corresponding solution \f$u_{\rm SSA}\f$ of the %SSA.  This is a forward problem which
  we would like to invert; given \f$u_{\rm SSA}\f$ determine \f$\tau_c\f$  such that
  \f$F_{\rm SSA}(\tau_c) = u_{\rm SSA}\f$.  The forward problem actually implemented by
  IP_SSAHardavForwardProblem is a mild variation \f$F_{\rm SSA}\f$.
 
  First, given the constraint \f$\tau_c\ge 0\f$ it can be helpful to parameterize \f$\tau_c\f$ by some other
  parameter \f$\zeta\f$,
  \f[
  \tau_c = g(\zeta),
  \f]
  where the function \f$g\f$ is non-negative.  The function \f$g\f$ is specified by an instance
  of IPDesignVariableParameterization.
 
  Second, there may be locations where the value of \f$\tau_c\f$ (and hence \f$\zeta\f$)
  is known a-priori, and should not be adjusted.  Let \f$\pi\f$ be the map that replaces
  the values of \f$\zeta\f$ with known values at these locations.
 
  IP_SSAHardavForwardProblem implements the forward problem
  \f[
  F(\zeta) = F_{\rm SSA}(g(\pi(\zeta))).
  \f]
 
  Algorithms to solve the inverse problem make use of variations on the linearization
  of \f$F\f$, which are explained below.
 
  The solution of the %SSA in SSAFEM is implemented using SNES to solve
  a nonlinear residual function of the form
  \f[
  \mathcal{R}_{SSA}(u;\tau_c, \text{other parameters})= \vec 0,
  \f]
  usually using Newton's method. Let
  \f[
  \mathcal{R}(u, \zeta) = \mathcal{R}_{SSA}(u; g(\pi(\zeta)), \text{other parameters}).
  \f]
 
  The derivative of \f$\mathcal{R}\f$ with respect to \f$ u\f$ is called the state Jacobian,
  \f$J_{\rm State}\f$. Specifically, if \f$u=[U_j]\f$ then
  \f[
  (J_{\rm State})_{ij} = \frac{d \mathcal{R}_i}{dU_j}.
  \f]
  This is exactly the same Jacobian that is computed by SSAFEM for solving the %SSA via SNES. The
  matrix for the design Jacobian can be obtained with \ref assemble_jacobian_state.
 
  The derivative of \f$\mathcal{R}\f$ with respect to \f$ \zeta\f$ is called the design Jacobian,
  \f$J_{\rm Design}\f$. Specifically, if \f$\zeta=[Z_k]\f$ then
  \f[
  (J_{\rm Design})_{ik} = \frac{d \mathcal{R}_i}{dZ_k}.
  \f]
  The map \f$J_{\rm Design}\f$ can be applied to a vector \f$d\zeta\f$ with
  apply_jacobian_design.  For inverse methods using adjoints, one also
  needs to be able to apply the transpose of \f$J_{\rm Design}\f$,
  which is done using \ref apply_jacobian_design_transpose.
 
  The forward problem map \f$F\f$ solves the implicit equation
  \f[
  \mathcal{R}(F(\zeta), \zeta) = 0.
  \f]
  Its linearisation \f$DF\f$ then satisfies
  \f[
  J_{\rm State}\; DF\; d\zeta + J_{\rm Design}\; d\zeta = 0
  \f]
  for any perturbation \f$d\zeta\f$.  Hence
  \f[
  DF = -J_{\rm State}^{-1}\circ J_{\rm Design}.
  \f]
  This derivative is sometimes called the reduced gradient in the literature.  To
  apply \f$DF\f$ to a perturbation \f$d\zeta\f$, use \ref apply_linearization.  Adjoint
  methods require the transpose of this map; to apply \f$DF^t\f$ to \f$du\f$ use
  \ref apply_linearization_transpose.
*/
class IP_SSAHardavForwardProblem : public stressbalance::SSAFEM
{
public:
 
  /// The function space for the design variable, i.e. \f$\tau_c\f$.
  typedef array::Scalar DesignVec;
  typedef array::Scalar1 DesignVecGhosted;
 
  /// The function space for the state variable, \f$u_{\rm SSA}\f$.
  typedef array::Vector StateVec;
  typedef array::Vector1 StateVec1;
 
  //! Constructs from the same objects as SSAFEM, plus a specification of how \f$\tau_c\f$ is parameterized.
  IP_SSAHardavForwardProblem(std::shared_ptr<const Grid> g,
                             IPDesignVariableParameterization &tp);
 
  virtual ~IP_SSAHardavForwardProblem() = default;
 
  //! Selects nodes where \f$\tau_c\f$ (more specifically \f$\zeta\f$) should not be adjusted.
  /*! The paramter \a locations should be set to 1 at each node where \f$\tau_c\f$
    is fixed. The forward map then effectively treats the design space as the subspace
    of nodes where \a locations is 0. Tangent vectors to this subspace, as would be
    generated by, e.g., \f$J_{\rm Design}^t\f$ are represented as vectors in the full
    space with entries set to zero in the fixed locations.  These can safely be added
    to preexisting values of \f$\zeta\f$ without changing the entries of \f$\zeta\f$ at the
    fixed locations.  Inversion can be done by setting an initial value of \f$\zeta\f$
    having the desired values in the fixed locations, and using set_tauc_fixed_locations()
    to indicated the nodes that should not be changed.
  */
  virtual void set_design_fixed_locations(array::Scalar &locations)
  {
    m_fixed_design_locations = &locations;
  }
 
  //! Returns the last solution of the %SSA as computed by \ref linearize_at.
  virtual std::shared_ptr<array::Vector> solution() {
    m_velocity_shared->copy_from(m_velocity);
    return m_velocity_shared;
  }
 
  //! Exposes the design variable parameterization being used.
  virtual IPDesignVariableParameterization & design_param() {
    return m_design_param;
  }
 
  void init();
 
  virtual void set_design(array::Scalar &zeta);
 
  virtual std::shared_ptr<TerminationReason> linearize_at(array::Scalar &zeta);
 
  virtual void assemble_residual(array::Vector &u, array::Vector &R);
  virtual void assemble_residual(array::Vector &u, Vec R);
 
  virtual void assemble_jacobian_state(array::Vector &u, Mat J);
 
  virtual void apply_jacobian_design(array::Vector &u, array::Scalar &dzeta, array::Vector &du);
  virtual void apply_jacobian_design(array::Vector &u, array::Scalar &dzeta, Vec du);
  virtual void apply_jacobian_design(array::Vector &u, array::Scalar &dzeta, Vector2d **du_a);
 
  virtual void apply_jacobian_design_transpose(array::Vector &u, array::Vector &du, array::Scalar &dzeta);
  virtual void apply_jacobian_design_transpose(array::Vector &u, array::Vector &du, Vec dzeta);
  virtual void apply_jacobian_design_transpose(array::Vector &u, array::Vector &du, double **dzeta);
 
  virtual void apply_linearization(array::Scalar &dzeta, array::Vector &du);
  virtual void apply_linearization_transpose(array::Vector &du, array::Scalar &dzeta);
 
  //! Exposes the DMDA of the underlying grid for the benefit of TAO.
  petsc::DM& get_da() const {
    return *m_da;
  }
 
protected:
  const int m_stencil_width;
 
  /// Current value of zeta, provided from caller.
  array::Scalar   *m_zeta;
  /// Storage for d_zeta with ghosts, if needed when an argument d_zeta is ghost-less.
  array::Scalar1   m_dzeta_local;
 
  /// Locations where \f$\tau_c\f$ should not be adjusted.
  array::Scalar *m_fixed_design_locations;
 
  /// The function taking \f$\zeta\f$ to \f$\tau_c\f$.
  IPDesignVariableParameterization &m_design_param;
 
  std::shared_ptr<array::Vector> m_velocity_shared;
 
  /// Temporary storage when state vectors need to be used without ghosts.
  array::Vector  m_du_global;
  /// Temporary storage when state vectors need to be used with ghosts.
  array::Vector1  m_du_local;
  /// Vertically-averaged ice hardness.
  array::Scalar1  m_hardav;
 
  fem::ElementIterator m_element_index;
  fem::Q1Element2       m_element;
 
  /// KSP used in \ref apply_linearization and \ref apply_linearization_transpose
  petsc::KSP  m_ksp;
  /// Mat used in \ref apply_linearization and \ref apply_linearization_transpose
  petsc::Mat  m_J_state;
 
  SNESConvergedReason m_reason;
 
  /// Flag indicating that the state jacobian matrix needs rebuilding.
  bool m_rebuild_J_state;
};
 
} // end of namespace inverse
} // end of namespace pism
 
#endif /* end of include guard: IP_SSAHARDAVFORWARDPROBLEM_HH_HAD68BK7 */