coriolis_block Derived Type

  type, public :: coriolis_block
    !* Author: Chris MacMackin
    !  Date: January 2017
    !
    ! A data type representing a matrix operator for the momentum
    ! components of the linear parts plume equations, with the
    ! Coriolis force. It can be useful when preconditioning a the
    ! plume solver. It is inherently 1-D in its implementation, but
    ! has a transverse velocity component.
    !
    ! This type solves the linear, coupled, differential equation 
    ! \[ \frac{d}{dx} \bm{Y} = \bm{A}\bm{Y} + \bm{F} \]
    ! where \( \bm{Y} = [U, V, U', V']^{T} \),
    ! \[ \bm{A} = \begin{bmatrix}
    ! 0        & 0         & 1 & 0 \\
    ! 0        & 0         & 0 & 1 \\
    ! 0        & -\Phi/\nu & 0 & 0 \\
    ! \Phi/\nu & 0         & 0 & 0
    ! \end{bmatrix},
    ! \]
    ! and \( \bm{F} \) is the input value to which this block is 
    ! applied in the preconditioner. \(\Phi\) is the dimensionless
    ! coriolis parameter, while \(\nu\) is the eddy viscosity.
    !
    ! The matrix \(\bm{A}\) can be diagonalised by factoring it such
    ! that \(\bm{A} = \bm{V}\bm{D}\bm{V}^{-1}\), where \(\bm{V}\) is a
    ! change of basis matrix with columns made up of the eigenvectors 
    ! of \(\bm{A}\) and \(\bm{D}\) is a diagonal matrix made up of
    ! the corresponding eigenvalues. These have the values
    ! \[ \bm{V} = \begin{bmatrix}
    ! \tfrac{-1-i}{\alpha} & \tfrac{-1+i}{\alpha} & \tfrac{1-i}{\alpha} & \tfrac{1+i}{\alpha} \\
    ! \tfrac{-1+i}{\alpha} & \tfrac{-1-i}{\alpha} & \tfrac{1+i}{\alpha} & \tfrac{1-i}{\alpha} \\
    ! i                    & -i                   & -i                  & i                   \\
    ! 1                    & 1                    & 1                   & 1
    ! \end{bmatrix},
    ! \]
    ! \[ \bm{D} = \beta\begin{bmatrix}
    ! -1-i & 0    & 0   & 0   \\
    ! 0    & -1+i &     & 0   \\
    ! 0    & 0    & 1-i & 0   \\
    ! 0    & 0    & 0   & 1+i
    ! \end{bmatrix},
    ! \]
    ! \[ \bm{V}^{-1} = \frac{1}{4}\begin{bmatrix}
    ! (-1+i)\beta & -(1+i)\beta & -i & 1 \\
    ! (-1-i)\beta & -(1-i)\beta & i  & 1 \\
    ! (1+i)\beta  & (1-i)\beta  & i  & 1 \\
    ! (1-i)\beta  & (1+i)\beta  & -i & 1
    ! \end{bmatrix},
    ! \]
    ! \[ \alpha = \sqrt{\frac{2\Phi}{\nu}}, \quad \beta = \sqrt{\frac{\Phi}{2\nu}}. \]
    ! It can be shown that the solution to the differential equation is
    ! \[ \bm{Y} = \bm{V}\left(e^{\bm{D}x}\bm{B} + 
    !    e^{\bm{D}x}\int_0^{x}e^{\bm{D}x'}\bm{V}^{-1}\bm{F}dx' \right),\]
    ! where \(\bm{B}\in\mathbb{R}^{4}\) is a vector chosen to satisfy
    ! the boundary conditions on the system. It can be found by solving a
    ! 4×4 linear system.
    !
    ! This type inherits 
    !
    private
    real(r8), dimension(4)   :: D_r
      !! Real component of the diagonal matrix, \(\bm{D}\), with only
      !! diagonal values stored
    real(r8), dimension(4)   :: D_i
      !! Imaginary component of the diagonal matrix, \(\bm{D}\), with
      !! only diagonal values stored
    type(cheb1d_scalar_field), dimension(4,4) :: emDxVinv_r
      !! Real component of \(e^{-\bm{D}x}\bm{V}^{-1}\)
    type(cheb1d_scalar_field), dimension(4,4) :: emDxVinv_i
      !! Imaginary component of \(e^{-\bm{D}x}\bm{V}^{-1}\)
    type(cheb1d_scalar_field), dimension(4) :: eDx_r
      !! Real component of \(e^{\bm{D}x}\), with only diagonal values
      !! stored
    type(cheb1d_scalar_field), dimension(4) :: eDx_i
      !! Imaginary component of \(e^{\bm{D}x}\), with only diagonal
      !! values stored
    real(r8), dimension(4,4) :: V_r
      !! Real component of the change of basis matrix, \(\bm{V}\)
    real(r8), dimension(4,4) :: V_i
      !! Imaginary component of the change of basis matrix, \(\bm{V}\)
    type(pseudospec_block)   :: integrator
      !! A pseudospectral differentiation block which can be used to
      !! perform integration
    integer                  :: vel_bound_loc
      !! Location code for the velocity's boundary condition
    integer                  :: dvel_bound_loc
      !! Location code for the velocity derivative's boundary
      !! condition
    integer                  :: integrate_bound
      !! Location from which to perform the integration
    real(r8), dimension(4)   :: xbounds
      !! Boundary location for each component of the solution vector
    complex(r8), dimension(4,4) :: bound_matrix
      !! Matrix for the system to solve in order to satisfy the
      !! boundary conditions
    complex(r8), dimension(4,4) :: bound_matrix_scaled
      !! Matrix for the system to solve in order to satisfy the
      !! boundary conditions, which has been scaled by LAPACK95 
      !! to improve conditioning.
    complex(r8), dimension(4,4) :: factored_matrix
      !! Factored matrix for the system to solve in order to satisfy
      !! the boundary conditions
    integer, dimension(4)    :: pivots
      !! The pivots used in the factorisation of the matrix used to
      !! satisfy boundary conditions
    real(r8), dimension(4)   :: r_scales
      !! Row scale factors from equilibrating the bound_matrix
    real(r8), dimension(4)   :: c_scales
      !! Column scale factors from equilibrating the bound_matrix
    character(len=1)         :: equed
      !! The method used to equilibrate bound_matrix
    integer :: int
  contains
    private
    procedure, public :: solve_for
    procedure         :: assign
    generic, public   :: assignment(=) => assign
  end type coriolis_block