jac_rbdpre Subroutine

public subroutine jac_rbdpre(t, u, r0, rblock, r1, rewt, cj, bd, ipbd, ierr)

Uses

  • proc~~jac_rbdpre~~UsesGraph proc~jac_rbdpre jac_rbdpre module~linpack linpack proc~jac_rbdpre->module~linpack module~daskr_kinds daskr_kinds module~linpack->module~daskr_kinds iso_fortran_env iso_fortran_env module~daskr_kinds->iso_fortran_env

This routine generates and preprocesses a block-diagonal preconditioner matrix, based on the part of the Jacobian corresponding to the reaction terms of the problem. It generates a matrix of the form . It calls dgefa from LINPACK to do the LU decomposition of each diagonal block. The computation of the diagonal blocks uses the mesh information in the module variables. One block per spatial point is computed. The Jacobian elements are generated by difference quotients. This routine calls a user routine rblock to compute the block (jx,jy) of .

Arguments

Type IntentOptional Attributes Name
real(kind=rk), intent(in) :: t

Current independent variable.

real(kind=rk), intent(inout) :: u(*)

Current dependent variables.

real(kind=rk), intent(in) :: r0(*)

Current value of the vector .

private subroutine rblock(t, jx, jy, uxy, rxy)

User routine that computes a single block of .

Arguments
Type IntentOptional Attributes Name
real(kind=rk), intent(in) :: t

Current independent variable.

integer, intent(in) :: jx

Spatial index in x-direction.

integer, intent(in) :: jy

Spatial index in y-direction.

real(kind=rk), intent(in) :: uxy(*)

Block of ns dependent variables at spatial point (jx,jy).

real(kind=rk), intent(out) :: rxy(*)

Block (jx,jy) of .

real(kind=rk), intent(inout) :: r1(*)

Real work space available to this subroutine.

real(kind=rk), intent(in) :: rewt(*)

Reciprocal error weights.

real(kind=rk), intent(in) :: cj

Scalar used in forming the system Jacobian.

real(kind=rk), intent(out) :: bd(*)

LU factors of the diagonal blocks.

integer, intent(out) :: ipbd(*)

Pivots for the LU factorizations.

integer, intent(inout) :: ierr

Error flag. 0: no error occurred; k > 0: a zero pivot was found at the k-th stage in one of the LU factorizations.


Calls

proc~~jac_rbdpre~~CallsGraph proc~jac_rbdpre jac_rbdpre proc~dgefa dgefa proc~jac_rbdpre->proc~dgefa interface~idamax idamax proc~dgefa->interface~idamax

Called by

proc~~jac_rbdpre~~CalledByGraph proc~jac_rbdpre jac_rbdpre proc~jac jac proc~jac->proc~jac_rbdpre

Source Code

   subroutine jac_rbdpre(t, u, r0, rblock, r1, rewt, cj, bd, ipbd, ierr)
   !! This routine generates and preprocesses a block-diagonal preconditioner matrix, based on
   !! the part of the Jacobian corresponding to the reaction terms \(R\) of the problem. It 
   !! generates a matrix of the form \(c_J I_d - \partial R/ \partial u\). It calls [[dgefa]]
   !! from LINPACK to do the LU decomposition of each diagonal block. The computation of the
   !! diagonal blocks uses the mesh information in the module variables. One block per spatial
   !! point is computed. The Jacobian elements are generated by difference quotients. 
   !! This routine calls a user routine `rblock` to compute the block `(jx,jy)` of \(R\).

      use linpack, only: dgefa

      real(rk), intent(in) :: t
         !! Current independent variable.
      real(rk), intent(inout) :: u(*)
         !! Current dependent variables.
      real(rk), intent(in) :: r0(*)
         !! Current value of the vector \(R(t,u)\).
      interface
         subroutine rblock(t, jx, jy, uxy, rxy)
            !! User routine that computes a single block of \(R\).
            import :: rk
            real(rk), intent(in) :: t
               !! Current independent variable.
            integer, intent(in) :: jx
               !! Spatial index in x-direction.
            integer, intent(in) :: jy
               !! Spatial index in y-direction.
            real(rk), intent(in) :: uxy(*)
               !! Block of `ns` dependent variables at spatial point `(jx,jy)`.
            real(rk), intent(out) :: rxy(*)
               !! Block `(jx,jy)` of \(R\).
         end subroutine rblock
      end interface
      real(rk), intent(inout) :: r1(*)
         !! Real work space available to this subroutine.
      real(rk), intent(in) :: rewt(*)
         !! Reciprocal error weights.
      real(rk), intent(in) :: cj
         !! Scalar used in forming the system Jacobian.
      real(rk), intent(out) :: bd(*)
         !! LU factors of the diagonal blocks.
      integer, intent(out) :: ipbd(*)
         !! Pivots for the LU factorizations.
      integer, intent(inout) :: ierr
         !! Error flag.
         !! `0`: no error occurred;
         !! `k > 0`: a zero pivot was found at the `k`-th stage in one of the LU factorizations.
      
      real(rk) :: del, dfac, fac, uj
      integer :: i, ibd, idiag, iip, j, j0, js, jx, jy

      ! Make MP calls to RBLOCK to approximate each diagonal block of dR/du.
      dfac = 1e-2_rk
      ibd = 0
      j0 = 0
      do jy = 1, meshy
         do jx = 1, meshx
            ! If R0 has not been set previously as an array of length NEQ, it can
            ! be set here, as an array of length MP, with the call:
            !    CALL RBLOCK (T, JX, JY, U(J0+1), R0)
            ! In this case, change R0(J0+I) below to R0(I).
            do js = 1, mp
               j = j0 + js
               uj = u(j)
               del = max(srur*abs(uj), dfac/rewt(j))
               u(j) = u(j) + del
               fac = -one/del
               call rblock(t, jx, jy, u(j0 + 1), r1)
               do i = 1, mp
                  bd(ibd + i) = (r1(i) - r0(j0 + i))*fac
               end do
               u(j) = uj
               ibd = ibd + mp
            end do
            j0 = j0 + mp
         end do
      end do

      ! Add matrix CJ*Id, and do LU decomposition on blocks.
      ibd = 1
      iip = 1
      do j = 1, meshx*meshy
         idiag = ibd
         do i = 1, mp
            if (i <= mpd) bd(idiag) = bd(idiag) + cj
            idiag = idiag + (mp + 1)
         end do
         call dgefa(bd(ibd), mp, mp, ipbd(iip), ierr)
         if (ierr /= 0) exit
         ibd = ibd + mpsq
         iip = iip + mp
      end do

   end subroutine jac_rbdpre