|
|
b4step1.f.html
| |
|
Source file: b4step1.f
|
|
Directory: /var/www/html/clawpack/links/shockvacuum10/shockvacuum
|
|
Converted: Mon Aug 16 2010 at 14:13:31
using clawcode2html
|
|
This documentation file will
not reflect any later changes in the source file.
|
c ============================================
subroutine b4step1(maxmx,mbc,mx,meqn,q,
& xlower,dx,t,dt,maux,aux)
c ============================================
c
c # called from claw1 before each call to step1.
c # used to print out shock location, h, u, dh each time step.
c
c
implicit double precision (a-h,o-z)
dimension q(1-mbc:maxmx+mbc, meqn)
dimension aux(1-mbc:maxmx+mbc, *)
integer midpoint, local_mthbc(2), step
double precision x0, x1, t1, midx, frac
logical regrid
integer icount, iold, inewleft, inewright
double precision qiold, qioldm1, qioldp1
double precision leftslope, rightslope, myslope
double precision npoly, nstar
logical treset
common /ceos/ grav, gamma, gamma1, sK, npoly, nstar, Crho, ieos
common /comxt/ dtcom,dxcom,tcom
common /comxt/ treset
c if (xlower.gt.-7.d0) then
c write(6,*) 'Stop in b4step1'
c stop
c endif
c # check for NANs in solution:
call check4nans(maxmx,meqn,mbc,mx,q,t,1)
c # Check for negative depth:
do i=1-mbc,mx+mbc
if (q(i,1) .lt. 0.d0) then
write(6,*) 'Resetting q(i,1) to zero, i = ',i
write(6,*) ' q(i,1) = ',q(i,1)
q(i,1) = 0.d0
stop
endif
enddo
regrid = .true.
frac = 0.1d0
x0 = -10.d0
x1 = -10.d0
t1 = 4.d0
step = 30
startx = x0 + t/t1*(x1-x0)
startpoint = int((startx-xlower)/dx)
if (startpoint .gt. mx) startpoint = mx
if (startpoint .lt. 1+3+step) startpoint = 1+3+step
write(30,*) t,startx,startpoint
dhmax = 0.d0
do i=startpoint,mx
dh = dabs(q(i+1,1)- q(i,1))
if (dh.gt.dhmax) then
dhmax = dh
imax = i
endif
enddo
c # Modification by RJL to look at where dv/dx is most negative:
dvmax = 0.d0
imax = mx
do i=mx,1,-1
dv = dabs(q(i,2)/q(i,1)- q(i+1,2)/q(i+1,1))
if (dv.gt.dvmax) then
dvmax = dv
imax = i
endif
enddo
xs = xlower + imax*dx
c write(6,*) '+++ imax = ',imax, ' dvmax = ',dvmax
hs = q(imax-3,1)
h1 = (hs - q(imax-3-step,1)) / (step*dx)
us = q(imax-3,2)/hs
u1 = (us - q(imax-3-step,2)/q(imax-3-step,1)) / (step*dx)
dh = q(imax-3,1) - q(imax+3,1)
c # Changed to write t and xs at end only after t is reset.
c write(29,291) t,xs,hs,us,dh,h1,u1
c write(39,*) t,imax,xs,dvmax
291 format(7e25.15)
c Regridding
if (regrid .and. xs .gt. frac*xlower) then
if (xs.gt.1.d-20) then
write(6,*) '***Aborting at t = ',t,
& ' since xs is too large, xs =', xs
stop
endif
write(6,*) 'regridding at t, xs =', t, xs
write(6,*) 'xlower from', xlower, 'to', 0.5d0*xlower
xlower = 0.5d0*xlower
dx = 0.5d0*dx
dt = 0.5d0*dt
dxcom = dx !# used in pres_psi
dtcom = dt
C$$$ccccccccccccccccccccccccccccccccccccc
c do i=1,mx
c write(501,*) i, q(i,1)
c write(502,*) i, q(i,2)/q(i,1)
c end do
c write(501,*) imax
c write(502,*) imax
C$$$cccccccccccccccccccccccccccccccccccccc
c Regrid the conserved quantities.
c Rewritten by RJL (12/09) to assume mx even, which simplifies things
c from CG's original version of minmod slopes assuming mx odd.
c q(i,:) is the cell average in cell i.
c Say we have cells 1, 2, 3, 4, 5, 6 (mx is always even).
c New cell 1 is the left half of old cell 4.
c New cell 6 is the right half of old cell 6.
mx2 = mx/2
do m=1,3
do iold=mx2+1, mx
qiold = q(iold,m)
c # one-sided slopes:
slopeleft = qiold - q(iold-1,m)
sloperight = q(iold+1,m) - qiold
c # minmod limiter:
if (dabs(slopeleft).lt.dabs(sloperight)) then
slope = slopeleft
else
slope = sloperight
endif
if (slopeleft*sloperight .le. 0.d0) then
slope = 0.d0
endif
c # Or, use zero slopes:
c slope = 0.d0
ileft = 2*(iold - mx2 - 1) + 1
iright = ileft + 1
q(ileft,m) = qiold - 0.25d0*slope
q(iright,m) = qiold + 0.25d0*slope
enddo
enddo
c In front of the shock, nothing should be moving, but when the
c shock gets very close to the beach, numerical error accumulates.
c To avoid this, overwrite with the rest solution.
c Is this still needed with wb method?
c Yes, since otherwise averaging densities, etc onto finer grid gives
c states that are no longer exact equilibrium states.
do i=imax+3,mx
xc = xlower + (i-0.5d0)*dx
q(i,1) = Crho*(-xc)**nstar
q(i,2) = 0.d0
pressure = grav/(nstar+1.d0)*Crho * (-xc)**(nstar+1.d0)
q(i,3) = pressure / gamma1
c +++ old:
c q(i,1) = (-grav*xc / (npoly*gamma*sK))**npoly
c q(i,2) = 0.d0
c q(i,3) = sK * q(i,1)**gamma / gamma1 +
c & 0.5d0*q(i,2)**2/q(i,1)
c +++ debug:
c write(51,501) i,xc,q(i,1),q(i,3)
c 501 format( i5,3e20.10)
end do
c The parameter mthbc has not been passed to us, but we know what it is:
local_mthbc(1)=1
local_mthbc(2)=3
call bc1(maxmx,meqn,mbc,mx,xlower,dx,q,maux,aux,t,dt,
$ local_mthbc)
c +++ debug:
c do i=mx-2,mx+2
c write(51,511) i,xc,q(i,1),q(i,2),q(i,3)
c 511 format( i5,4e16.6)
c end do
C$$$ccccccccccccccccccccccccccccccccccccc
c do i=1,mx
c write(501,*) i, q(i,1)
c write(502,*) i, q(i,2)/q(i,1)
c end do
c write(501,*) imax
c write(502,*) imax
c stop 'stop for debugging after first regridding'
C$$$cccccccccccccccccccccccccccccccccccccc
endif
c # always reset to equilibrium soln near right boundary
c # see if this works for n=1.5, nstar=3. case?
do i=mx-20,mx
xc = xlower + (i-0.5d0)*dx
write(51,501) i,xc,q(i,1),q(i,3)
q(i,1) = Crho*(-xc)**nstar
q(i,2) = 0.d0
pressure = grav/(nstar+1.d0)*Crho * (-xc)**(nstar+1.d0)
q(i,3) = pressure / gamma1
c +++ debug:
write(51,501) i,xc,q(i,1),q(i,3)
501 format( i5,3e20.10)
end do
c xreset = -1.d-5 ! used for Euler
xreset = -1.d-4
if (xlower.gt.xreset .and. (.not. treset)) then
write(6,*) '==== Resetting t to 0 at time t = ',t
write(26,261) t
261 format(e26.16,' time tc reset to 0')
t = 0.d0
treset = .true.
endif
if (treset) then
write(29,291) t,xs
else
write(28,291) t,xs
endif
return
end