bc1.f.html | ![]() |
Source file: bc1.f | |
Directory: /Users/rjl/git/rjleveque/clawpack-4.6.3/book/chap9/acoustics/layered | |
Converted: Mon Jan 21 2013 at 20:15:48 using clawcode2html | |
This documentation file will not reflect any later changes in the source file. |
c c c ================================================================= subroutine bc1(maxmx,meqn,mbc,mx,xlower,dx,q,maux,aux,t,dt,mthbc) c ================================================================= c c # Standard boundary condition choices for claw2 c c # Modified for layered medium: c # For small t oscillating solid wall BCs are used to generate pulse c # For t>50, switch to periodic boundary conditions. The code c # can then be run to much larger t to observe how the pulse c # behaves as it loops around. c c # At each boundary k = 1 (left), 2 (right): c # mthbc(k) = 0 for user-supplied BC's (must be inserted!) c # = 1 for zero-order extrapolation c # = 2 for periodic boundary coniditions c # = 3 for solid walls, assuming this can be implemented c # by reflecting the data about the boundary and then c # negating the 2'nd component of q. c ------------------------------------------------ c c # Extend the data from the computational region c # i = 1, 2, ..., mx2 c # to the virtual cells outside the region, with c # i = 1-ibc and i = mx+ibc for ibc=1,...,mbc c implicit double precision (a-h,o-z) dimension q(1-mbc:maxmx+mbc, meqn) dimension aux(1-mbc:maxmx+mbc, *) dimension mthbc(2) common /comwall/ pi,t1,a1,tw1,t2,a2,tw2 c c # switch to periodic boundary conditions once pulse is generated: if (t .gt. 50.d0 .and. t .lt. 51.d0) then mthbc(1) = 2 mthbc(2) = 2 do m=1,meqn do ibc=1,mbc aux(1-ibc,m) = aux(mx+1-ibc,m) enddo enddo endif c c------------------------------------------------------- c # left boundary: c------------------------------------------------------- go to (100,110,120,130) mthbc(1)+1 c 100 continue c # user-specified boundary conditions c # oscillating wall do 105 m=1,meqn do 105 ibc=1,mbc q(1-ibc,m) = q(ibc,m) 105 continue c # wall velocity: vwall = a1*g0((t-t1)/tw1) c # adjust the normal velocity: do 106 ibc=1,mbc q(1-ibc,2) = 2.0d0*vwall - q(ibc,2) 106 continue go to 199 c 110 continue c # zero-order extrapolation: do 115 m=1,meqn do 115 ibc=1,mbc q(1-ibc,m) = q(1,m) 115 continue go to 199 120 continue c # periodic: do 125 m=1,meqn do 125 ibc=1,mbc q(1-ibc,m) = q(mx+1-ibc,m) 125 continue go to 199 130 continue c # solid wall (assumes 2'nd component is velocity or momentum in x): do 135 m=1,meqn do 135 ibc=1,mbc q(1-ibc,m) = q(ibc,m) 135 continue c # negate the normal velocity: do 136 ibc=1,mbc q(1-ibc,2) = -q(ibc,2) 136 continue go to 199 199 continue c c------------------------------------------------------- c # right boundary: c------------------------------------------------------- go to (200,210,220,230) mthbc(2)+1 c 200 continue c # user-specified boundary conditions c # oscillating wall do 205 m=1,meqn do 205 ibc=1,mbc q(mx+ibc,m) = q(mx+1-ibc,m) 205 continue c # wall velocity: vwall = a2*g0((t-t2)/tw2) c # adjust the normal velocity: do 206 ibc=1,mbc q(mx+ibc,2) = 2.0d0*vwall - q(mx+1-ibc,2) 206 continue go to 299 210 continue c # zero-order extrapolation: do 215 m=1,meqn do 215 ibc=1,mbc q(mx+ibc,m) = q(mx,m) 215 continue go to 299 220 continue c # periodic: do 225 m=1,meqn do 225 ibc=1,mbc q(mx+ibc,m) = q(ibc,m) 225 continue go to 299 230 continue c # solid wall (assumes 2'nd component is velocity or momentum in x): do 235 m=1,meqn do 235 ibc=1,mbc q(mx+ibc,m) = q(mx+1-ibc,m) 235 continue do 236 ibc=1,mbc q(mx+ibc,2) = -q(mx+1-ibc,2) 236 continue go to 299 299 continue c return end c =============================== double precision function g0(t) c =============================== implicit double precision (a-h,o-z) common /comwall/ pi,t1,a1,t2,a2,tw1,tw2 if (dabs(t) .lt. 1.d0) then g0 = 1.d0 + dcos(pi*t) else g0 = 0.d0 endif return end