jacobian_block Derived Type

  type, public :: jacobian_block
    !* Author: Chris MacMackin
    !  Date: December 2016
    !
    ! A data type representing a submatrix of a
    ! Jacobian. Specifically, it represents the commonly occurring
    ! operation $$ \frac{\partial F}{\partial x_i} + F\Delta_i, $$
    ! where \( \Delta_i \) is the differentiation operator in the
    ! \(i\)-direction. Optionally, there can be an additional
    ! differentiation operator on the right-hand-side of this.
    !
    ! This data type is useful when constructing a preconditioner. A
    ! preconditioner, \(A\), should approximate the inverse Jacobian of
    ! a problem being solved. Rather than directly computing the
    ! inverse Jacobian, it is more efficient to approximate it. If
    ! \(d\) is the vector being preconditioned, and \(z\) is the
    ! result of applying the preconditioner, then $$ z = J^{-1}d
    ! \Rightarrow Jz = d.$$ Thus, the preconditioner can be applied
    ! by approximately solving this system for \(z\). Linearising
    ! \(J\), this system can be solved efficiently using Picard
    ! iteration.
    !
    ! Say that the Jacobian is an \(n\times n\) system of these blocks
    ! (each labeled as \(J_{j,k}\)) and that the vector being
    ! preconditioned constists of \(n\) scalar fields (\(d_j\)). Then
    ! the \(m^{th}\) estimate of the solution for the \(j^{th}\) field
    ! in the preconditioned vector (\(z^m_j\)) is the solution to $$
    ! J_{j,j}z^m_j = d_j - \sum_{\substack{k=1\ k\ne j}}^n
    ! J_{j,k}z^{m-1}_k. $$ Depending on the type of fields being used
    ! and the direction in which derivatives are being taken,
    ! \(J_{j,j}\) may be tridiaganol, meaning it can be solved
    ! efficiently. <!--Otherwise, it can be approximated that $$
    ! J_{j,j}z^m_j = \left(\frac{\partial F}{\partial x_i} +
    ! F\Delta_i\right)z^{m} \simeq \frac{\partial F}{\partial
    ! x_i}z^m_j + F\frac{\partial z^{m-1}_j}{\partial x_i}. $$ The
    ! first term in this approximation corresponds to a diagonal
    ! matrix (which can be solved trivially), while the second term is
    ! known and can be subtracted from the right-hand-side of the
    ! linear system.-->
    !
    ! For this purpose, a type-bound multiplication operator is
    ! provided for the Jacobian block type, which can be used when
    ! evaluating the right hand side of the linear system. There is
    ! also a `solve_for` method, which takes as an argument the field
    ! representing the right hand side of the system. <!--Before
    ! applying the latter operator, the `update_estimate` method must
    ! be called. This allows for the term involving derivative of the
    ! current estimate of the solution to be subtracted from the
    ! right-hand-side of the linear system, in cases where the system
    ! can not be expressed as a tridiagonal matrix.-->
    !
    ! Boundary conditions are a complicated issue. When constructing
    ! an instance of this type, locations in the raw field
    ! representation can be specified for which boundary conditions
    ! are given. By default these are all Dirichlet conditions, but
    ! this behaviour may be overridden to be Neumann conditions. Note
    ! that Neumann conditions result in approximately an order of
    ! magnitude lower accuracy. When using the Jacobian block for
    ! multiplation, the type of the boundary does not matter; the
    ! value at the boundary in the result will be set to 0 or to
    ! values produced by an optional user-provided function. When
    ! using the `solve_for` method, the tridiagonal matrix will be
    ! altered in the specified locations to reflect the correct
    ! boundary conditions. This effectively mimics the sorts of
    ! Jacobian rows typically produced by boundary conditions.
    !
    private
    integer                             :: direction
      !! The direction in which any derivatives are taken.
    integer                             :: extra_derivative = no_extra_derivative
      !! The direction in which to apply a differentiation on the
      !! right-hand-side of the Jacobian block operator. Defaults
      !! none.
    class(scalar_field), allocatable    :: contents
      !! The value, \(A\), to which the Jacobian block operation is
      !! beiing applied.
    real(r8)                            :: coef = 1._r8
      !! Optional coefficient by which the the \(\partial F/\partial x
      !! \) term in the operator will be multiplied.
    procedure(jacobian_block_bounds), pointer, nopass  :: &
                                           get_boundaries
      !! A subroutine which determines the expected boundary conditions
      !! (and their location in the raw array) for the solution of the
      !! Jacobian block.
    class(scalar_field), allocatable    :: derivative
      !! The cached derivative of `contents`
    real(r8), dimension(:), allocatable :: diagonal
      !! The diagonal of the tridiagonal matrix representation of this
      !! block.
    real(r8), dimension(:), allocatable :: super_diagonal
      !! The super-diagonal of the tridiagonal matrix representation
      !! of this block.
    real(r8), dimension(:), allocatable :: sub_diagonal
      !! The sub-diagonal of the tridiagonal matrix representation of
      !! this block.
    real(r8), dimension(:), allocatable :: l_multipliers
      !! Multipliers defining the L matrix in the LU factorisation of
      !! the tridiagonal matrix representation of this block.
    real(r8), dimension(:), allocatable :: u_diagonal
      !! The diagonal of the U matrix in the LU factorisation of
      !! the tridiagonal matrix representation of this block.
    real(r8), dimension(:), allocatable :: u_superdiagonal1
      !! The first superdiagonal of the U matrix in the LU
      !! factorisation of the tridiagonal matrix representation of
      !! this block.
    real(r8), dimension(:), allocatable :: u_superdiagonal2
      !! The second superdiagonal of the U matrix in the LU
      !! factorisation of the tridiagonal matrix representation of
      !! this block.
    integer, dimension(:), allocatable  :: pivots
      !! Pivot indicies from the LU factorisation of the tridiagonal
      !! matrix representation of this block.
    real(r8), dimension(:), allocatable :: boundary_vals
      !! Expected boundary values for the solution to the Jacobian
      !! block.
    integer, dimension(:), allocatable  :: boundary_locs
      !! Locations in the raw arrays which are used to specify
      !! boundary conditions.
    integer, dimension(:), allocatable  :: boundary_types
      !! The types of boundary conditions, specified using the
      !! parameters found in [[boundary_types_mod]].
    real(r8)                            :: real_increment
      !! A scalar value which is to be added to this Jacobian block
      !! (i.e. to the diagonal).
    class(scalar_field), allocatable    :: field_increment
      !! A scalar field which is to be added to this Jacobian block
      !! (i.e. to the diagonal).
    type(jacobian_block), pointer       :: block_increment => null()
      !! Another Jacobian block which is to be added to this one
    integer                             :: has_increment = 0
      !! Indicates whether or not there has been an increment added to
      !! this block. If not, then 0. If a scalar real value has been
      !! added, then 1. If a scalar value has been added, then 2.
  contains
    private
    procedure :: jacobian_block_multiply
    procedure :: jacobian_block_add_real
    procedure :: jacobian_block_add_field
    procedure :: jacobian_block_add_block
    procedure :: jacobian_block_assign
    procedure :: get_tridiag => jacobian_block_get_tridiag
    generic, public :: operator(*) => jacobian_block_multiply
    generic, public :: operator(+) => jacobian_block_add_real,  &
                                      jacobian_block_add_field, &
                                      jacobian_block_add_block
    generic, public :: assignment(=) => jacobian_block_assign
    procedure, public :: solve_for => jacobian_block_solve
  end type jacobian_block