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 .
Type | Intent | Optional | 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
|
||||||||||||||||||||||||||||||||||||||||||||||||
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.
|
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