This is a _temporary_ repository to be able to start to work on the

code right now. I have put in the public version of Whisky to start from.

Everybody with commit rights should get commit messages (and the other
way around). It should not be a problem to add people to that list, just
ask. I don't want to get into political problems because someone feels
excluded, but I also don't want to give everyone access per se.

Frank



git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/GRHydro_InitData/trunk@3 ac85fae7-cede-4708-beff-ae01c7fa1c26

1 files changed, 308 insertions, 0 deletions

diff --git a/src/tov_lsoda.F b/src/tov_lsoda.F
new file mode 100644
index 0000000..b7a1e9f
--- /dev/null
+++ b/src/tov_lsoda.F
@@ -0,0 +1,308 @@
+#ifdef HAVE_ODEPACK
+c===========================================================
+c     History: sode.f
+c
+c     Driver routine which integrates ODEs defining 
+c     TOV star
+c 
+c    Based on bstar.f by Matthew Choptuik
+c===========================================================
+      subroutine     tov_lsoda ( rho, p, m, alf,
+     &                          r,  nr, tol , n_succ, ssrc) 
+
+      implicit       none
+
+
+      integer       vsrc,  vsxynt
+      real*8         rho2p,         p2rho
+
+      integer        nr
+      real*8         rho(nr),p(nr), m(nr), alf(nr)
+      real*8         r(nr)
+      real*8         tol
+      integer        n_succ, ssrc
+
+      character*12    codenm
+      parameter    ( codenm = 'tov_lsoda' )
+
+      integer        iargc,         indlnb,         i4arg
+      real*8         r8arg
+
+      real*8         r8_never
+      parameter    ( r8_never = -1.0d-60 )
+
+c-----------------------------------------------------------
+c     Order of system.
+c-----------------------------------------------------------
+      integer        neq
+      parameter    ( neq = 3 )
+
+c-----------------------------------------------------------
+c     Storage for solution at requested output radii.
+c-----------------------------------------------------------
+      integer        maxout
+      parameter    ( maxout = 10 000 )
+
+      real*8         y0(neq)
+      real*8         vxout(maxout), vyout(maxout,neq),
+     &               work(maxout)
+      integer        nxout,         ixout,       nxout_succ
+
+      integer        ieq
+
+      logical        ltrace
+      parameter    ( ltrace = .false. )
+
+      integer        maxdump
+      parameter    ( maxdump = 50 )
+
+      logical        lsodatrace
+      parameter    ( lsodatrace = .false. )
+
+
+c-----------------------------------------------------------
+c     LSODA Variables.
+c-----------------------------------------------------------
+      external       tov_fcn,        tov_jac
+      real*8         y(neq)
+      real*8         tbgn,       tend
+      integer        itol
+      real*8         rtol,       atol
+      integer        itask,      istate,     iopt
+      integer        lrw
+
+c      parameter    ( lrw = 22 + neq * max(16 ,neq + 9) )
+      parameter    ( lrw = 22 + neq * (16) )
+      real*8         rwork(lrw)
+
+      integer        liw
+      parameter    ( liw = 20 + neq )
+      integer        iwork(liw)
+      integer        jt
+
+      integer        i
+
+c-----------------------------------------------------------
+c     Common communication with routine 'tov_fcn' in 'tov_fcn.f' ...
+c-----------------------------------------------------------
+c   rho0:  Value of rho at r=0
+c-----------------------------------------------------------
+      real*8  rho0
+      common / com_tov / rho0
+
+c-----------------------------------------------------------
+c     Parse arguments.  
+c-----------------------------------------------------------
+      rho0    = rho(1)
+
+      nxout   = nr
+      do ixout = 1, nxout
+         vxout(ixout) = r(ixout)
+      enddo
+
+c-----------------------------------------------------------
+c     Initialize those "boundary" conditions which are 
+c     fixed by regularity, or which are arbitrary.
+c-----------------------------------------------------------
+      y0(1) = rho2p(rho0)
+      y0(2) = 0.0d0
+      y0(3) = 1.0
+
+
+c-----------------------------------------------------------
+c     Echo command line arguments if "local tracing" 
+c     enabled ...
+c-----------------------------------------------------------
+      if( ltrace ) then
+         write(0,*) 'rho0: ', rho0
+         write(0,*) 'p(1): ', y0(1)
+         write(0,*) 'm(1): ', y0(2)
+         write(0,*) 'alf(1): ', y0(3)
+      end if
+
+c-----------------------------------------------------------
+c     Get output radii 
+c-----------------------------------------------------------
+      if( ltrace ) then
+         if( nxout .le. maxdump ) then
+            call dvdump(vxout,nxout,codenm//': output radii',0)
+         else 
+            call dvdmp1(vxout,work,
+     &                  int(1.0d0 * nxout / maxdump + 0.5d0),
+     &                  nxout,codenm//': selected output radii',0)
+         end if
+         write(0,*)
+         write(0,*) codenm//': Initial time: ', vxout(1)
+         write(0,*)
+      end if
+
+c-----------------------------------------------------------
+c     Set LSODA parameters ...
+c-----------------------------------------------------------
+      itol   = 1
+      rtol   = tol
+      atol   = tol
+      itask  = 1 
+      iopt   = 0
+      jt     = 2
+
+c-----------------------------------------------------------
+c     Initialize the solution ...
+c-----------------------------------------------------------
+      do ieq = 1 , neq
+         y(ieq)       = y0(ieq)
+         vyout(1,ieq) = y0(ieq)
+      end do
+
+c-----------------------------------------------------------
+c     Do the integration ...
+c-----------------------------------------------------------
+      ssrc = 0
+      do ixout = 2 , nxout
+         istate = 1
+c-----------------------------------------------------------
+c        Need these temporaries since lsoda overwrites 
+c        tend ... 
+c-----------------------------------------------------------
+100      continue
+         tbgn = vxout(ixout-1)
+         tend = vxout(ixout)
+         call lsoda(tov_fcn,neq,y,tbgn,tend,
+     &           itol,rtol,atol,itask,    
+     &           istate,iopt,rwork,lrw,iwork,liw,tov_jac,jt)
+         if( lsodatrace ) then
+            write(0,1000) codenm, ixout, vxout(ixout), 
+     &      vxout(ixout+1)
+1000        format(' ',a,': Step ',i4,'  t = ',1pe10.3,
+     &      ' .. ',1pe10.3)
+            write(0,*) codenm//': lsoda reurns ', istate
+         end if
+
+c         if (istate .eq. -2) then
+c            write(0,*) 'istate = -2, trying again...'
+c            itol   = 1
+c            rtol   = tol
+c            atol   = tol
+c            itask  = 1
+c            iopt   = 0
+c            jt     = 2
+c            istate = 3
+c            goto 100
+c         endif
+
+         if( istate .lt. -1 ) then
+            write(0,1500) codenm, istate, ixout, nxout, 
+     &                    vxout(ixout-1), vxout(ixout)
+1500        format(/' ',a,': Error return ',i2,
+     &              ' from integrator LSODA.'/
+     &              '       At output time ',i5,' of ',i5/
+     &              '       Interval ',1pe11.3,' .. ',
+     &             1pe11.3/)
+             nxout_succ = ixout - 1
+             ssrc = 1
+            go to 500
+
+         end if
+
+         do ieq = 1 , neq
+            vyout(ixout,ieq) = y(ieq)
+         end do
+
+      end do
+      nxout_succ = nxout
+
+500   continue
+
+c-----------------------------------------------------------
+c     OUTPUT
+c-----------------------------------------------------------
+      do ixout = 1 , nxout_succ
+         p(ixout)    = vyout(ixout, 1)
+         m(ixout)    = vyout(ixout, 2)
+         alf(ixout)    = vyout(ixout, 3)
+         rho(ixout)    = p2rho(p(ixout))
+      end do
+
+      n_succ = nxout_succ
+
+      return
+
+      end
+
+
+
+
+
+c===========================================================
+c     Some double precision vector utility routines. 
+c
+c     Originally from ~matt/vutil/dveclib.f
+c===========================================================
+ 
+c-----------------------------------------------------------
+c     Dumps vector labelled with LABEL on UNIT.
+c-----------------------------------------------------------
+      subroutine dvdump(v,n,label,unit)
+         implicit       none
+
+         real*8         v(*)
+         character*(*)  label
+         integer        i, n, st, unit
+ 
+         if( n .lt. 1 ) go to 130
+            write(unit,100) label
+ 100        format(/' <<< ',A,' >>>'/)
+            st = 1
+ 110        continue
+               write(unit,120) ( v(i) , i = st , min(st+3,n))
+ 120           FORMAT(' ',4(1PE19.10))
+               st = st + 4
+            if( st .le. n ) go to 110
+ 
+ 130     continue
+ 
+         return
+      end
+ 
+c-----------------------------------------------------------
+c     Extension of dvdump which dumps every 'inc'th element.
+c-----------------------------------------------------------
+      subroutine dvdmp1(v,w,inc,n,label,unit)
+         implicit       none
+
+         real*8         v(*),    w(*)
+         character*(*)  label
+         integer        inc,     n,      unit
+ 
+         call dvinj(v,w,inc,n)
+         call dvdump(w,1+(n-1)/inc,label,unit)
+
+         return
+
+      end
+ 
+c---------------------------------------------------------
+c     Injects every incth element of v1 into v2
+c---------------------------------------------------------
+      subroutine dvinj(v1,v2,inc,n)
+         implicit    none
+   
+         real*8      v1(*),     v2(*)
+         integer     i,         inc,     j,      n
+ 
+         j = 1
+         do i = 1 , n , inc
+            v2(j) = v1(i)
+            j = j + 1
+         end do
+
+ 
+         return
+      end
+
+#else
+      subroutine tov_lsoda()
+      return 
+      end
+ 
+#endif