Changeset 1491

Timestamp:

Jun 1, 2005, 4:04:38 PM (20 years ago)

Author:

steve

Message:

 

Location:

inundation/ga/storm_surge

Files:

• Hobart/Landslide_1D.for (modified) (1 diff)
• Hobart/Landslide_Tadmor.for (modified) (16 diffs)
• Hobart/run_hobart_buildings.py (modified) (1 diff)
• zeus/pyvolution.zbi (modified) (previous)

inundation/ga/storm_surge/Hobart/Landslide_1D.for

@@ -348 +348 @@
 c
 c ---- Try for a flat subinterval 
-c
-      do i=md,nx+md+2
-         hmax = dmax1(h_l(i),h_r(i))
-         dz  = (z2(i) - z2(i-1))
-
-         if (h_l(i) .lt. 0.01d0*h_r(i) ) then
-           h_r(i)  = sqrt(abs(2.0d0*dz*h(i)))
-         endif
-
-         if (h_r(i) .lt. 0.01d0*h_l(i) ) then         
-           h_l(i)  = sqrt(abs(2.0d0*dz*h(i)))       
-         endif         
-      enddo
+cc
+c      do i=md,nx+md+2
+c         hmax = dmax1(h_l(i),h_r(i))
+c         dz  = (z2(i) - z2(i-1))
+c
+c         if (h_l(i) .lt. 0.01d0*h_r(i) ) then
+c           h_r(i)  = sqrt(abs(2.0d0*dz*h(i)))
+c         endif
+c
+c         if (h_r(i) .lt. 0.01d0*h_l(i) ) then         
+c           h_l(i)  = sqrt(abs(2.0d0*dz*h(i)))       
+c         endif         
+c      enddo
             
       return

inundation/ga/storm_surge/Hobart/Landslide_Tadmor.for

@@ -3 +3 @@
 c
 c     Stiff source term solution using the second order 
-c     central scheme of Tadmor and the modifications
+c     central scheme of Landslide and the modifications
 c     suggested by Liotta, Romano and Russo, SIAM
 c     J. Numer. Anal. 38(4), 1337-1356, 2000.
 c     This is the Randau scheme.
 c
-c     Solving Sampson's parabolic canal problem
 c============================================================
       implicit none
@@ -16 +15 @@
 
       integer nx, nxx, io, i, isteps
-      real*8 toll,cfl,theta,g,h0,dx,dtmax,Tout,Tfinal,tc,a,tau,B,t0,
-     .       p,s,zeta_0,zeta_1,zeta,u0
-
-      open(1,file='Stepflow5.out')
-      open(2,file='Stepflow5.t')
-      open(3,file='Stepflow5.h')
-      open(4,file='Stepflow5.uh')
-      open(5,file='Stepflow5.n')
-      open(8,file='Stepflow5.c')
+      real*8 toll,cfl,theta,g,dx,dtmax,Tout,Tfinal,tc,a,tau,B,t0,
+     .       p,s,zeta_0,zeta_1,zeta,u0,version,hl,hr
+
+      open(11,file='data.w')
+      open(12,file='data.t')
+      open(13,file='data.h')
+      open(14,file='data.uh')
+      open(15,file='data.n')
+      open(16,file='data.x')
+
+      open(1,file='MacDonald.bed')
 c
 c Parameters
-c      
-      nx = 501
+c   
+      version = 1.0   
+      nx = 101
       nxx = nx+2*md
       if (nx.gt.nxmax) write(*,*)'nx too large'
       toll = 1.d-6
-      cfl = 0.5d0
+      cfl = 0.2d0
       theta = 1.0d0
       g = 9.81d0
-      h0=10.0d0
+      hl=10.0d0
+      hr = 0.d0
       dx = 10000.d0/dble(nx-1)
       dtmax  = 0.5d0
-      Tout   = 20.0d0
-      Tfinal = 10000.0d0
-      a = 3000.d0
-      tau = 0.001d0
-      B = 5.d0
+      Tout   = 1.0d0
+      Tfinal = 1500.0d0
+
       t0 = 0.d0
 
@@ -48 +49 @@
 c Setup Terrain
 c
-      do io = 1,2
-         do i = 1,nxx
-            x(i,io) = dble(i-md-1)*dx+(io-1)*dx/2.0
-         enddo
-         call bed(x(1,io),z(1,io),dz(1,io),nxx)
+c ---- Non-uniform terrain from MacDonald
+      io = 1
+      do i = 1,nxx
+        read(1,*)x(i,io),z(i,io)
+        z(i,io) = z(i,io)*5.d0
+      end do
+      do i = 2,nxx-1
+        dz(i,io) = (z(i+1,io) - z(i-1,io))/2.d0/dx
+      end do
+      dz(1,io) = dz(2,io)
+      dz(nxx,io) = dz(nxx-1,io)
+
+      io = 2
+      do i = 1,nxx
+         x(i,io) = x(i,io-1)
+         z(i,io) = z(i,io-1)
+         dz(i,io) = dz(i,io-1)
       end do
       
@@ -58 +71 @@
 c Initial Conditions
 c 
-      p = dsqrt(8.d0*g*h0/a/a)
-      s = dsqrt(p*p-tau*tau)/2.d0
-      u0 = B*exp(-tau*t0/2.d0)*sin(s*t0)
-      zeta_0 = a*a*B*B*exp(-tau*t0)/(8.d0*g*g*h0)
-     .         *(-s*tau*cos(s*t0)+(tau*tau/4.d0-s*s)*cos(2.d0*s*t0))
-     .         -B*B*exp(-tau*t0)/4.d0/g
-      zeta_1 = -exp(-tau*t0/2.d0)/g*(B*s*cos(s*t0)
-     .         +tau*B/2.d0*sin(s*t0))
-      do i = 1,nxx
-        u(i,1)  = 20.d0 - h0*(1.d0 - (x(i,1) - 5000.d0)**2/a/a)
-        zeta = zeta_0 + (x(i,1) - 5000.d0)*zeta_1 + 20.d0
-        if(zeta.lt.u(i,1))then
-          u(i,1) = 0.d0
+c ---- Rectangular pulse initial condition
+      do i = md+1, nx+md
+        if (x(i,io).lt.100.d0.or.x(i,io).gt.200.d0) then
+          u(i,1) = hr
         else
-          u(i,1) = zeta - u(i,1)
-        end if
-        u(i,2) = 0.d0
-      end do
+          u(i,1) = hl
+        endif
+        u(i,2) = 0.0d0
+      end do  
 
 c
 c Setup for Time looping
 c          
-      isteps = 0
+
       tc = 0.0d0
       call dump_solution(u,z,nx,tc,isteps)
@@ -87 +91 @@
 c
       call momentum_mass(nx,dx,cfl,g,theta,Tfinal,
-     . u,dtmax,h0,toll,z,dz,Tout)
+     .           u,dtmax,toll,z,dz,Tout,isteps)
+
+      write(15,*)isteps,nxx,version
+      do i=1,nxx
+         write(16,*) x(i,1)
+      enddo
+
+
 c
 c Close down
@@ -104 +115 @@
 c The evolver
 c===============================================================
-      subroutine momentum_mass(nx,dx,cfl,g,theta,tf,u,dtmax,h0,toll,
-     .                         z,dz,Tout)
+      subroutine momentum_mass(nx,dx,cfl,g,theta,tf,u,dtmax,toll,
+     .                         z,dz,Tout,isteps)
 c
 c INPUT   nx # of cells in x-direction
@@ -129 +140 @@
      
       integer nx, nxx, io, i, isteps
-      real*8 toll,cfl,theta,g,h0,dx,dt,Tout,Tfinal,tc
+      real*8 toll,cfl,theta,g,dx,dt,Tout,Tfinal,tc
       
       integer istop,iout,nt,ii,oi,i2,m
@@ -173 +184 @@
 c see (2.1) and (4.3)
 c          
-          call swflux(f,u,1,nxx,em_x,g,toll)
-          
-          call print_twoarray(f,nxx,'flux')
-          call print_twoarray(u,nxx,'h and uh')          
+          call swflux(f,u,1,nxx,em_x,g,toll) 
 c
 c Compute numerical derivatives 'ux', 'fx'.
@@ -187 +195 @@
           call derivative(f(1,2),fx(1,2),nxx,theta)    
 
-          call print_twoarray(df,nxx,'df')
-          call print_twoarray(ux,nxx,'ux')    
 c
 c Compute time step size according to the input CFL #
@@ -194 +200 @@
           if(ii.eq.1) then
             dt = dmin1(cfl*dx/em_x,dtmax,Tout)
-c            write(*,*)em_x, dt
             if( ( tc + 2.d0*dt ) .ge. Tnext ) then
               dt = 0.5d0*(Tnext - tc )
@@ -204 +209 @@
             endif
           end if
-c          write(*,*)'dt = ',dt
-c          write(*,*)'Tnext = ',Tnext
+
           dtcdx2 = dt/dx/2.d0
           dtcdx3 = dt/dx/3.d0
@@ -226 +230 @@
           call swflux(f,u2,3,nxx-2,em_x,g,toll)
 
-          call print_twoarray(f,nxx,'flux half')
 c
 c Compute the values of 'u' at the next time level. see (2.16).
@@ -246 +249 @@
           end do
 
-          call print_twoarray(u,nxx,'u after h update')
 c
 c Momentum equation
@@ -267 +269 @@
 c
 c Friction term
-            u(i,m) = u(i,m)*exp(-0.001d0*dt)
-          end do
-          call print_twoarray(u,nxx,'u after uh update')          
+c            u(i,m) = u(i,m)*exp(-0.001d0*dt)
+          end do
 c
           tc = tc + dt
@@ -282 +283 @@
 c
  1001 write(*,*)isteps
-      write(5,*)isteps,nx
+
       return
  
       end
 
-c===========================================================
-c Setup Terrain with parabolic canal
-c===========================================================
-      subroutine bed(x,z,dz,n)
-
-      implicit none
-      integer md, nxmax,nxd,mn
-      parameter (md=3, nxmax=1601, nxd=nxmax+2*md, mn=2)
-      integer n
-      real*8 z(n),dz(n),x(n)
-      
-      integer i
-      real*8 h0, a
-
-      h0 = 10.d0
-      a = 3000.d0
-
-      do i=1,n
-         z(i)  = 20.d0 - h0*(1.d0 - (x(i) - 5000.d0)**2/a/a)
-         dz(i) = 2.d0*h0*(x(i) - 5000.d0)/a/a
-      enddo
-      return
-      end
 c===========================================================
 c Shallow Water Flux
@@ -393 +371 @@
       real*8 u(nxd,mn),z(nxd,2)
 
-      write(2,*)tc
+      write(12,*)tc
       io = 1
+      i = 2
       isteps = isteps + 1
       write(*,*)isteps, tc
-      do ii = 1,2
-        do i = md + 1, nx + md
-           write(3,200)u(i,1)
-           write(4,200)u(i,2)
-           write(1,200)z(i,io) + u(i,1)
+
+      do i = 1, nx + 2*md
+         write(13,200)u(i,1)
+         write(14,200)u(i,2)
+         write(11,200)z(i,io) + u(i,1)
   200 format(e20.8)
-        end do
-      end do
-
-      iflag = 0
-      do i = md + 1,nx + md
-        if(u(i,1).gt.1.d-04.and.iflag.eq.0)then
-c          write(8,200)u(i,1)
-          write(8,200)i*20.d0
-          iflag = 1
-        end if
-      end do
+      end do
+
       return
       end
-c===========================================================
-c Output Routine
-c===========================================================
-      subroutine print_twoarray(u,nxx,msg)
-      implicit real*8 (a-h)
-      implicit real*8 (o-z)
-      parameter (md=3, nxmax=1601, nxd=nxmax+2*md, mn=2)
-      real*8 u(nxd,mn)
-      character*(*) msg
-      return
-      write(*,*)msg
-
-      do i = 1,nxx
-         write(*,*)u(i,1),u(i,2)
-      end do     
-      return
-      end

inundation/ga/storm_surge/Hobart/run_hobart_buildings.py

@@ -47 +47 @@
 
 # Colours for Hobart
-domain.visualiser.stage_color = (0.0,0.38,0.1)
+domain.visualiser.stage_color    = (0.0,0.38,0.1)
 domain.visualiser.friction_color = (0.1,0.99,0.1)
-domain.visualiser.other_color = (0.99,0.4,0.1)
+domain.visualiser.other_color    = (0.99,0.4,0.1)