c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
c     This file is part of the Test Set for IVP solvers
c     http://www.dm.uniba.it/~testset/
c
c        Water tube system
c        index 2 DAE of dimension 49
c
c     DISCLAIMER: see
c     http://www.dm.uniba.it/~testset/disclaimer.php
c
c     The most recent version of this source file can be found at
c     http://www.dm.uniba.it/~testset/src/problems/caraxis.f
c
c     This is revision
c     $Id: water.F,v 1.2 2006/10/02 10:29:14 testset Exp $
c
c-----------------------------------------------------------------------
      integer function pidate()
      pidate = 20060828
      return
      end
c-----------------------------------------------------------------------
      subroutine prob(fullnm,problm,type,
     +                neqn,ndisc,t,
     +                numjac,mljac,mujac,
     +                nummas,mlmas,mumas,
     +                ind)
      character*(*) fullnm, problm, type
      integer neqn,ndisc,mljac,mujac,mlmas,mumas,ind(*)
      double precision t(0:*)
      logical numjac, nummas

      integer i

      fullnm = 'Water tube system'
      problm = 'water'
      type   = 'DAE'
      neqn   = 49
      ndisc  = 0
      t(0)   = 0d0
      t(1)   = 17d0*3600d0
      numjac = .true.
      mljac  = neqn
      mujac  = neqn
      mlmas  = 0
      mumas  = 0
      do 10 i=1,38
         ind(i) = 1
   10 continue
      do 20 i=39,neqn
         ind(i) = 2
   20 continue

      return
      end
c-----------------------------------------------------------------------
      subroutine init(neqn,t,y,yprime,consis)
      integer neqn
      double precision t,y(neqn),yprime(neqn)
      logical consis

      integer i

      do 10 i=1,neqn
         y(i) = 0d0
         yprime(i) = 0d0
   10 continue
      do 20 i=19,36
         y(i) = 0.47519404529185289807d-1
   20 continue
      do 30 i=37,49
         y(i) = 109800d0
   30 continue

      consis = .true.

      return
      end
c-----------------------------------------------------------------------
      subroutine settolerances(neqn,rtol,atol,tolvec)
      integer neqn 
      logical tolvec
      double precision rtol(neqn), atol(neqn)

      integer i

      do 30 i=2,49
            atol(i) = atol(1)
            rtol(i) = rtol(1)
   30    continue
c
c        pressures are much larger:
c
      do 40 i=37,49
            atol(i) = 1d6*atol(1)
   40    continue

      tolvec = .true.       
      
      

      return
      end

c-----------------------------------------------------------------------
      subroutine setoutput(neqn,solref,printsolout,
     +                    nindsol,indsol)

      logical solref, printsolout
      integer neqn, nindsol
      integer indsol(neqn)

c the reference solution is available
      solref = .true.  

c output file is required
      printsolout = .true.

c default values if printsolout is .true.
      nindsol = neqn
c only nindsol component of indsol are referenced
      do i=1,nindsol
          indsol(i) = i
      end do  

 

      return
      end

c-----------------------------------------------------------------------
      subroutine feval(neqn,t,y,yprime,f,ierr,rpar,ipar)
      integer neqn,ierr,ipar(*)
      double precision t,y(neqn),yprime(neqn),f(neqn),rpar(*)

      integer nnodes
      parameter(nnodes=13)
      double precision phi(nnodes,nnodes),lambda(nnodes,nnodes),
     +                 p(nnodes),
     +                 fdba(nnodes,nnodes),netflo(nnodes),
     +                 ein(nnodes),eout(nnodes),rghres(nnodes,nnodes)

      double precision nu,g,rho,rcrit,length,k,d,b,pi
      parameter(nu    = 1.31d-6,
     +          g     = 9.8d0,
     +          rho   = 1.0d3,
     +          rcrit = 2.3d3,
     +          length= 1.0d3,
     +          k     = 2.0d-4,
     +          d     = 1.0d0,
     +          b     = 2.0d2,
     +          pi    = 3.141592653589793238462643383d0 )

      double precision a,mu
      integer i,j,n

      double precision rtla,r,that,that2

      a = pi*d**2/4d0
      mu = nu*rho

      do 5 j=1,nnodes
         ein(j) = 0d0
         eout(j) = 0d0
    5 continue

      that=t/3600d0
      that2=that*that

      ein(1) = (1d0-cos(exp(-that)-1d0))/200d0
      ein(13) = (1d0-cos(exp(-that)-1d0))/80d0
      eout(10) = that2*(3d0*that2-92d0*that+720d0)/1d6

      do 8 j=1,nnodes
         do 7 i=1,nnodes
            phi(i,j) = 0d0
            lambda(i,j) = 1d0
    7    continue
    8 continue

      phi( 1, 2) = y( 1)
      phi( 2, 3) = y( 2)
      phi( 2, 6) = y( 3)
      phi( 3, 4) = y( 4)
      phi( 3, 5) = y( 5)
      phi( 4, 5) = y( 6)
      phi( 5,10) = y( 7)
      phi( 6, 5) = y( 8)
      phi( 7, 4) = y( 9)
      phi( 7, 8) = y(10)
      phi( 8, 5) = y(11)
      phi( 8,10) = y(12)
      phi( 9, 8) = y(13)
      phi(11, 9) = y(14)
      phi(11,12) = y(15)
      phi(12, 7) = y(16)
      phi(12, 8) = y(17)
      phi(13,11) = y(18)

      lambda( 1, 2) = y(19)
      lambda( 2, 3) = y(20)
      lambda( 2, 6) = y(21)
      lambda( 3, 4) = y(22)
      lambda( 3, 5) = y(23)
      lambda( 4, 5) = y(24)
      lambda( 5,10) = y(25)
      lambda( 6, 5) = y(26)
      lambda( 7, 4) = y(27)
      lambda( 7, 8) = y(28)
      lambda( 8, 5) = y(29)
      lambda( 8,10) = y(30)
      lambda( 9, 8) = y(31)
      lambda(11, 9) = y(32)
      lambda(11,12) = y(33)
      lambda(12, 7) = y(34)
      lambda(12, 8) = y(35)
      lambda(13,11) = y(36)

      p( 5) = y(37)
      p( 8) = y(38)
      p( 1) = y(39)
      p( 2) = y(40)
      p( 3) = y(41)
      p( 4) = y(42)
      p( 6) = y(43)
      p( 7) = y(44)
      p( 9) = y(45)
      p(10) = y(46)
      p(11) = y(47)
      p(12) = y(48)
      p(13) = y(49)

      do 20 j=1,nnodes
         do 10 i=1,nnodes
            if (lambda(i,j) .lt. 0d0) then
               ierr = -1
               return
            endif
            rtla=sqrt(lambda(i,j))
            r = abs(phi(i,j)*d/(nu*a))
            if (r.gt.rcrit) then
               rghres(i,j) = 1/rtla - 1.74d0 +
     +                       2d0*log10(2d0*k/d + 18.7d0/(r*rtla))
               fdba(i,j) = p(i) - p(j) -
     +                     lambda(i,j)*rho*length*phi(i,j)**2/(a**2*d)
            else
               rghres(i,j) = 1/rtla - 1.74d0 +
     +                       2d0*log10(2d0*k/d + 18.7d0/(rcrit*rtla))
c              rghres(i,j) = lambda(i,j) - 0.47519404529185289807d-1
               fdba(i,j) = p(i) - p(j) -
     +                     32d0*mu*length*phi(i,j)/(a*d**2)
            endif
   10    continue
   20 continue

      do 50 n=1,nnodes
         netflo(n) = ein(n)-eout(n)
         do 30 i=1,nnodes
            netflo(n) = netflo(n)+phi(i,n)
   30    continue
         do 40 j=1,nnodes
            netflo(n) = netflo(n)-phi(n,j)
   40    continue
   50 continue

      f( 1) = fdba( 1, 2)
      f( 2) = fdba( 2, 3)
      f( 3) = fdba( 2, 6)
      f( 4) = fdba( 3, 4)
      f( 5) = fdba( 3, 5)
      f( 6) = fdba( 4, 5)
      f( 7) = fdba( 5,10)
      f( 8) = fdba( 6, 5)
      f( 9) = fdba( 7, 4)
      f(10) = fdba( 7, 8)
      f(11) = fdba( 8, 5)
      f(12) = fdba( 8,10)
      f(13) = fdba( 9, 8)
      f(14) = fdba(11, 9)
      f(15) = fdba(11,12)
      f(16) = fdba(12, 7)
      f(17) = fdba(12, 8)
      f(18) = fdba(13,11)

      f(19) = rghres( 1, 2)
      f(20) = rghres( 2, 3)
      f(21) = rghres( 2, 6)
      f(22) = rghres( 3, 4)
      f(23) = rghres( 3, 5)
      f(24) = rghres( 4, 5)
      f(25) = rghres( 5,10)
      f(26) = rghres( 6, 5)
      f(27) = rghres( 7, 4)
      f(28) = rghres( 7, 8)
      f(29) = rghres( 8, 5)
      f(30) = rghres( 8,10)
      f(31) = rghres( 9, 8)
      f(32) = rghres(11, 9)
      f(33) = rghres(11,12)
      f(34) = rghres(12, 7)
      f(35) = rghres(12, 8)
      f(36) = rghres(13,11)

      f(37) = netflo( 5)
      f(38) = netflo( 8)
      f(39) = netflo( 1)
      f(40) = netflo( 2)
      f(41) = netflo( 3)
      f(42) = netflo( 4)
      f(43) = netflo( 6)
      f(44) = netflo( 7)
      f(45) = netflo( 9)
      f(46) = netflo(10)
      f(47) = netflo(11)
      f(48) = netflo(12)
      f(49) = netflo(13)

      return
      end
c-----------------------------------------------------------------------
      subroutine jeval(ldim,neqn,t,y,yprime,dfdy,ierr,rpar,ipar)
      integer ldim,neqn,ierr,ipar(*)
      double precision t,y(neqn),yprime(neqn),dfdy(ldim,neqn),rpar(*)
c
c     dummy subroutine
c
      return
      end
c-----------------------------------------------------------------------
      subroutine meval(ldim,neqn,t,y,yprime,dfddy,ierr,rpar,ipar)
      integer ldim,neqn,ierr,ipar(*)
      double precision t,y(neqn),yprime(neqn),dfddy(ldim,neqn),rpar(*)

      double precision nu,g,rho,rcrit,length,k,d,b,pi
      parameter(nu    = 1.31d-6,
     +          g     = 9.8d0,
     +          rho   = 1.0d3,
     +          rcrit = 2.3d3,
     +          length= 1.0d3,
     +          k     = 2.0d-4,
     +          d     = 1.0d0,
     +          b     = 2.0d2,
     +          pi    = 3.141592653589793238462643383d0 )

      double precision a,c,v
      integer i

      a = pi*d**2/4d0
      c = b/(rho*g)
      v = rho*length/a

      do 10 i=1,18
         dfddy(1,i) = v
   10 continue

      do 20 i=19,36
         dfddy(1,i) = 0d0
   20 continue

      do 30 i=37,38
         dfddy(1,i) = c
   30 continue

      do 40 i=39,neqn
         dfddy(1,i) = 0d0
   40 continue

      return
      end
c-----------------------------------------------------------------------
      subroutine solut(neqn,t,y)
      integer neqn
      double precision t,y(neqn)

c Test Set for IVP Solvers (release 2.0)
c
c Solving Water tube system using PSIDE
c
c User input:
c
c give relative error tolerance: 1d-14
c give absolute error tolerance: 1d-14
c
c Numerical solution:
c
c                                               scd
c        solution component               ----------------   ignore
c                                          abs        rel
c ---------------------------------       -----      -----   ------
      y(  1) =  0.2298488296477430d-002
      y(  2) =  0.1188984650746585d-002
      y(  3) =  0.1109503645730845d-002
      y(  4) =  0.1589620100314825d-003
      y(  5) =  0.1030022640715102d-002
      y(  6) =  0.8710606306836165d-003
      y(  7) =  0.3243571480903489d-002
      y(  8) =  0.1109503645730845d-002
      y(  9) =  0.7120986206521341d-003
      y( 10) =  0.6414613963833099d-003
      y( 11) =  0.9416978549524347d-003
      y( 12) =  0.3403428519096511d-002
      y( 13) =  0.2397639310739395d-002
      y( 14) =  0.2397639310739395d-002
      y( 15) =  0.3348581430454180d-002
      y( 16) =  0.1353560017035444d-002
      y( 17) =  0.1995021413418736d-002
      y( 18) =  0.5746220741193575d-002
      y( 19) =  0.4751940452918529d-001
      y( 20) =  0.4751940452918529d-001
      y( 21) =  0.4751940452918529d-001
      y( 22) =  0.4751940452918529d-001
      y( 23) =  0.4751940452918529d-001
      y( 24) =  0.4751940452918529d-001
      y( 25) =  0.4311196778792902d-001
      y( 26) =  0.4751940452918529d-001
      y( 27) =  0.4751940452918529d-001
      y( 28) =  0.4751940452918529d-001
      y( 29) =  0.4751940452918529d-001
      y( 30) =  0.4249217433601160d-001
      y( 31) =  0.4732336439609648d-001
      y( 32) =  0.4732336439609648d-001
      y( 33) =  0.4270002118868241d-001
      y( 34) =  0.4751940452918529d-001
      y( 35) =  0.4751940452918529d-001
      y( 36) =  0.3651427026675656d-001
      y( 37) =  0.1111268591478108d+006
      y( 38) =  0.1111270045592387d+006
      y( 39) =  0.1111271078730254d+006
      y( 40) =  0.1111269851929858d+006
      y( 41) =  0.1111269255355337d+006
      y( 42) =  0.1111269322658045d+006
      y( 43) =  0.1111269221703983d+006
      y( 44) =  0.1111270121140691d+006
      y( 45) =  0.1111274419515807d+006
      y( 46) =  0.1111255158881087d+006
      y( 47) =  0.1111278793439227d+006
      y( 48) =  0.1111270995171642d+006
      y( 49) =  0.1111298338971779d+006
c
c Integration characteristics:
c
c    number of integration steps        4059
c    number of accepted steps           3746
c    number of f evaluations           66593
c    number of Jacobian evaluations       85
c    number of LU decompositions        4692

      return
      end