example_web Program

program example_web
!! Example program for [[daskr]]:
!! DAE system derived from \(s\)-species interaction PDE in 2 dimensions.
!!
!! This program solves a DAE system that arises from a system of partial differential equations.
!! The PDE system is a food web population model, with predator-prey interaction and diffusion
!! on the unit square in two dimensions. The dependent variable vector is
!!
!! $$ c = [c_1, c_2 , ..., c_s] $$
!!
!! and the PDEs are as follows:
!!
!! $$\begin{aligned}
!! \frac{\partial c_i}{\partial t} &= d_i \left( \frac{\partial^2 c_i}{\partial x^2}
!!           + \frac{\partial^2 c_i}{\partial y^2} \right) + v_i \quad i=1,...,s/2 \\
!!                               0 &= d_i \left( \frac{\partial^2 c_i}{\partial x^2}
!!           + \frac{\partial^2 c_i}{\partial y^2} \right) + v_i \quad i=s/2+1,...,s
!! \end{aligned}$$
!!
!! where the rate of formation of species \(i\) is given by:
!!
!! $$ v_i(x,y,c) = c_i \left( b_i + \sum_{j=1}^s a_{ij} c_j \right) $$
!!
!! The number of species is \(s = 2 p\), with the first \(p\) being prey and the last \(p\)
!! being predators. The coefficients \(a_{ij}\), \(b_i\), \(d_i\) are:
!!
!! $$\begin{aligned}
!!   a_{ii} &= -a  \quad (\mathrm{all}\; i) \\
!!   a_{ij} &= -g  \quad (i \le p,\; j > p) \\
!!   a_{ij} &=  e  \quad (i > p,\; j \le p)
!! \end{aligned}$$
!!
!! $$\begin{aligned}
!!   b_i    &=  b (1 + \alpha x y + \beta \sin(4 \pi x) \sin(4 \pi y)) \quad (i \le p) \\
!!   b_i    &= -b (1 + \alpha x y + \beta \sin(4 \pi x) \sin(4 \pi y)) \quad (i > p)
!! \end{aligned}$$
!!
!! $$\begin{aligned}
!!   d_i    &= d_{prey} \quad (i \le p) \\
!!   d_i    &= d_{pred} \quad (i > p)
!! \end{aligned}$$
!!
!! The various scalar parameters are set in subroutine [[setpar]].
!!
!! The boundary conditions are of Neumann type (zero normal derivative) everywhere. A polynomial
!! in \(x\) and \(y\) is used to set the initial conditions. The PDEs are discretized by central
!! differencing on a \( M_x \times M_y\) mesh.
!!
!! The root function is:
!!
!! $$ r(t,c,c') = \iint c_1 dx dy - 20 $$
!!
!! The DAE system is solved by [[daskr]] with three different method options:
!!
!! 1. Direct band method for the linear systems (internal Jacobian),
!! 2. Preconditioned Krylov method for the linear systems, without block-grouping in the
!!    reaction-based factor, and
!! 3. Preconditioned Krylov method for the linear systems, with block-grouping in the
!!    reaction-based factor.
!!
!! In the Krylov cases, the preconditioner is the product of two factors:
!!
!! * The spatial factor uses a fixed number of Gauss-Seidel iterations based on the diffusion
!!   terms only.
!! * The reaction-based factor is a block-diagonal matrix based on the partial derivatives of
!!   the interaction terms `R` only.
!!
!! With block-grouping, only a subset of the \((s \times s)\) blocks are computed. An integer
!! flag, `jpre`, is set in the main program to specify whether the preconditioner is to use only
!! one of the two factors or both, and in which order.
!!
!! The reaction-based preconditioner factor is set up and solved in seven subroutines:
!!
!! * [[DMSET2]], [[DRBDJA]], [[DRBDPS]] in the case of no block-grouping, and
!! * [[DGSET2]], [[GSET1]], [[DRBGJA]], [[DRBGPS]]  in the case of block-grouping.
!!
!! These routines are provided separately for general use on problems arising from
!! reaction-transport systems.
!!
!! Two output files are written: one with the problem description and performance statistics on
!! and one with solution profiles at selected output times. The solution file is written only
!! in the case of the direct method.
!!
!! References:
!!
!! * Peter N. Brown and Alan C. Hindmarsh, "Reduced Storage Matrix Methods in Stiff ODE
!!   Systems", J. Appl. Math. & Comp., 31 (1989), pp. 40-91.
!! * Peter N. Brown, Alan C. Hindmarsh, and Linda R. Petzold, "Using Krylov Methods in the
!!   Solution of Large-Scale Differential-Algebraic Systems", SIAM J. Sci. Comput., 15 (1994),
!!   pp. 1467-1488.
!! * Peter N. Brown, Alan C. Hindmarsh, and Linda R. Petzold, "Consistent Initial Condition
!!   Calculation for Differential-Algebraic Systems", LLNL Report UCRL-JC-122175, August 1995;
!!   SIAM J. Sci. Comput., 19 (1998), pp. 1495 - 1512.

   use daskr_kinds, only: rk, zero, one
   use web_par, only: aa, alpha, bb, beta, dpred, dprey, ee, gg, mx, my, mxns, np, ns, setpar
   use web_m
   implicit none

   integer, parameter :: neq = ns*mx*my, maxm = max(mx, my)
   integer, parameter :: lrw = 63 + (3*ns*maxm + 11)*neq + maxm, liw = 40 + 2*neq
   integer :: idid, info(20), iout, ipar(2), iwork(liw), jbg, jpre, jroot, leniw, &
              lenrw, ldout, lcout, meth, ncfl, ncfn, ng, nli, nlidif, nni, nnidif, nout, &
              npe, nps, nqu, nre, nrt, nrte, nst, nxg, nyg
   real(rk) :: atol, avlin, c1ave, cc(neq), ccprime(neq), hu, pred_ic, rpar(neq), rtol, &
               rwork(lrw), t, tout

   ! Dimension solution arrays and work arrays.
   !
   ! When INFO(12) = 0, with INFO(5) = 0, INFO(6) = 1:
   !       The length required for RWORK is
   !       60 + 3*NRT + (2*ML+MU+11)*NEQ + 2*(NEQ/(ML+MU+1) + 1) .
   !       For MX = MY = (even number) and ML = MU = NS*MX, this length is
   !       60 + 3*NRT + (3*NS*MX + 11)*NEQ + MY .
   !       The length required for IWORK is  40 + 2*NEQ .
   !
   ! When INFO(12) = 1:
   !       The length required for RWORK is
   !       101 + 3*NRT + 19*NEQ + LENWP = 104 + 19*NEQ + NS*NS*NGRP .
   !       The length required for IWORK is
   !       40 + NEQ + LENIWP = 40 + NEQ + NS*NGRP .
   !
   ! The dimensions for the various arrays are set using parameters
   !       MAXN    which must be >= NEQ = NS*MX*MY,
   !       MAXM    which must be >= MAX(MX,MY).

   ! Open output files.
   open (newunit=ldout, file='./example/example_web_d.out', action="write", position="rewind")
   open (newunit=lcout, file='./example/example_web_c.out', action="write", position="rewind")

   ! Call SETPAR to set basic problem parameters.
   call setpar

   ! Set NRT = number of root functions.
   nrt = 1

   write (ldout, '(a, /)') 'DWEB: Example program for DASKR package'
   write (ldout, '(a, i4)') 'Food web problem with NS species, NS =', ns
   write (ldout, '(a, /)') 'Predator-prey interaction and diffusion on a 2-D square'
   write (ldout, '(a, e12.4, a, e12.4, a, e12.4)') &
      'Matrix parameters..  a =', aa, '  e =', ee, '  g =', gg
   write (ldout, '(21x, a, e12.4)') &
      'b =', bb
   write (ldout, '(a, e12.4, a, e12.4)') &
      'Diffusion coefficients: dprey =', dprey, '  dpred =', dpred
   write (ldout, '(a, e12.4, a, e12.4)') &
      'Rate parameters alphaa =', alpha, ' and beta =', beta
   write (ldout, '(a, 2i4, 5x, a, i7)') &
      'Mesh dimensions (MX,MY) =', mx, my, ' Total system size is NEQ =', neq
   write (ldout, '(a)') 'Root function is R(Y) = average(c1) - 20'

   ! Set the flat initial guess for the predators.
   pred_ic = 1e5_rk

   ! Set remaining method parameters for DDASKR.
   ! These include the INFO array and tolerances.
   info = 0

   ! Set INFO(11) = 1, indicating I.C. calculation requested.
   info(11) = 1

   ! Set INFO(14) = 1 to get the computed initial values.
   info(14) = 1

   ! Set INFO(15) = 1 to signal that a preconditioner setup routine is to be called in the
   ! Krylov case.
   info(15) = 1

   ! Set INFO(16) = 1 to get alternative error test (on the differential variables only).
   info(16) = 1

   ! Set the tolerances.
   rtol = 1e-5_rk
   atol = rtol

   write (ldout, '(/, a, e10.2, a, e10.2)') &
      'Tolerance parameters: RTOL =', rtol, '  ATOL =', atol
   write (ldout, '(a, i2, a)') &
      'Internal I.C. calculation flag INFO(11) =', info(11), '  (0 = off, 1 = on)'
   write (ldout, '(a, e10.2)') &
      'Predator I.C. guess =', pred_ic
   write (ldout, '(a, i2, a)') &
      'Alternate error test flag INFO(16) =', info(16), '  (0 = off, 1 = on)'

   ! Set NOUT = number of output times.
   nout = 18

   ! Loop over method options:
   ! METH = 0 means use INFO(12) = 0 (direct)
   ! METH = 1 means use INFO(12) = 1 (Krylov) without block-grouping in the reaction-based
   !          factor in the preconditioner.
   ! METH = 2 means use INFO(12) = 1 (Krylov) with block-grouping in the reaction-based factor
   !          in the preconditioner.
   ! A block-grouping flag JBG, communicated through IPAR, is set to 0 (no block-grouping) or
   ! 1 (use block-grouping) with METH = 1 or 2.
   ! Reset INFO(1) = 0 and INFO(11) = 1.

   do meth = 0, 2

      info(12) = min(meth, 1)
      info(1) = 0
      info(11) = 1
      jbg = meth - 1
      ipar(2) = jbg

      write (ldout, '(/, 80("."), //, a, i2, a)') &
         'Linear solver method flag INFO(12) =', info(12), '  (0 = direct, 1 = Krylov)'

      ! In the case of the direct method, set INFO(6) = 1 to signal a banded Jacobian, set
      ! IWORK(1) = IWORK(2) = MX*NS, the half-bandwidth, and call SETID to set the IWORK
      ! segment containing the block ID indicating the differential and algebraic components.
      if (info(12) == 0) then
         info(6) = 1
         iwork(1) = mxns
         iwork(2) = mxns
         call setid(mx, my, ns, np, 40, iwork)
         write (ldout, '(a, i4)') 'Difference-quotient banded Jacobian, half-bandwidths =', mxns
      end if

      ! In the case of the Krylov method, set and print various preconditioner parameters.
      if (info(12) == 1) then

         ! First set the preconditioner choice JPRE.
         !  JPRE = 1 means reaction-only (block-diagonal) factor A_R
         !  JPRE = 2 means spatial factor (Gauss-Seidel) A_S
         !  JPRE = 3 means A_S * A_R
         !  JPRE = 4 means A_R * A_S
         ! Use IPAR to communicate JPRE to the preconditioner solve routine.
         jpre = 3
         ipar(1) = jpre
         write (ldout, '(a, i3)') 'Preconditioner flag is JPRE =', jpre
         write (ldout, '(a)') &
            '(1 = reaction factor A_R, 2 = spatial factor A_S, 3 = A_S*A_R, 4 = A_R*A_S)'

         ! Call DMSET2 if JBG = 0, or DGSET2 if JBG = 1, to set the 2D mesh parameters and 
         ! block-grouping data, and the IWORK segment ID indicating the differential and 
         ! algebraic components.
         if (jbg == 0) then
            call dmset2(mx, my, ns, np, 40, iwork)
            write (ldout, '(a)') 'No block-grouping in reaction factor'
         end if
         if (jbg == 1) then
            nxg = 5
            nyg = 5
            ng = nxg*nyg
            call dgset2(mx, my, ns, np, nxg, nyg, 40, iwork)
            write (ldout, '(a)') 'Block-grouping in reaction factor'
            write (ldout, '(a, i5, a, i3, a, i3, a)') &
               'Number of groups =', ng, ' (NGX by NGY, NGX =', nxg, ',  NGY =', nyg, ')'
         end if
      end if

      ! Set the initial T and TOUT, and call CINIT to set initial values.
      t = zero
      tout = 1.0e-8_rk
      call cinit(cc, ccprime, pred_ic, rpar)

      nli = 0
      nni = 0

      ! Header for "ldout"
      write (ldout, '(/10x, a)') &
         't      <c1>  NSTEP   NRE   NNI   NLI   NPE    NQ          H    AVLIN'

      ! Loop over output times and call DASKR. At each output time, print average c1 value and
      ! performance data.
      ! The first call, with IOUT = 0, is to calculate initial values only.
      ! After the first call, reset INFO(11) = 0 and the initial TOUT.
      ! If a root was found, we flag this, and return to the DASKR call.
      do iout = 0, nout

         do
            call daskr(res, neq, t, cc, ccprime, tout, info, rtol, atol, &
                       idid, rwork, lrw, iwork, liw, rpar, ipar, jacrs, psolrs, &
                       rt, nrt, jroot)

            nst = iwork(11)
            nre = iwork(12)
            npe = iwork(13)
            nnidif = iwork(19) - nni
            nni = iwork(19)
            nlidif = iwork(20) - nli
            nli = iwork(20)
            nqu = iwork(8)
            hu = rwork(7)
            avlin = zero
            if (nnidif > 0) avlin = one*nlidif/nnidif

            if (meth == 0) then
               call out(t, cc, ns, mx, my, lcout)
            end if

            call c1_average(cc, c1ave)
            write (ldout, '(e11.5, f10.5, i7, i6, i6, i6, i6, i6, e11.2, f9.4)') &
               t, c1ave, nst, nre, nni, nli, npe, nqu, hu, avlin

            if (idid == 5) then
               write (ldout, '(15x, a, i3)') '*****   Root found, JROOT =', jroot
            else
               exit
            end if

         end do

         if (idid < 0) then
            write (ldout, '(//a, e12.4//)') 'Final time reached =', t
            exit
         end if

         if (tout > 0.9_rk) then
            tout = tout + one
         else
            tout = 10*tout
         end if

         if (iout == 0) then
            info(11) = 0
            tout = 1.0e-8_rk
            nli = 0
            nni = 0
         end if

      end do

      lenrw = iwork(18)
      leniw = iwork(17)
      nst = iwork(11)
      nre = iwork(12)
      npe = iwork(13)
      nni = iwork(19)
      nli = iwork(20)
      nps = iwork(21)
      if (nni > 0) avlin = real(nli)/real(nni)
      ncfn = iwork(15)
      ncfl = iwork(16)
      nrte = iwork(36)

      write (ldout, '(//a)') 'Final statistics for this run..'
      write (ldout, '(a, i8, a, i6)') 'RWORK size =', lenrw, '  IWORK size =', leniw
      write (ldout, '(a, i5)') 'Number of time steps              =', nst
      write (ldout, '(a, i5)') 'Number of residual evaluations    =', nre
      write (ldout, '(a, i5)') 'Number of root fn. evaluations    =', nrte
      write (ldout, '(a, i5)') 'Number of Jac. or prec. evals.    =', npe
      write (ldout, '(a, i5)') 'Number of preconditioner solves   =', nps
      write (ldout, '(a, i5)') 'Number of nonlinear iterations    =', nni
      write (ldout, '(a, i5)') 'Number of linear iterations       =', nli
      write (ldout, '(a, f8.4)') 'Average Krylov subspace dimension =', avlin
      write (ldout, '(i3, a, i5, a)') ncfn, ' nonlinear conv. failures,', ncfl, ' linear conv. failures'

   end do

   close (unit=ldout)
   close (unit=lcout)

end program example_web