File contents
c=======================================================================
c
c \\\\\\\\\\ B E G I N S U B R O U T I N E //////////
c ////////// G C O L L A P S E \\\\\\\\\\
c
c Developed by
c Laboratory of Computational Astrophysics
c University of Illinois at Urbana-Champaign
c
c=======================================================================
c
subroutine gcollapse
c
c gcollapse <----------- initialises spherical gas cloud for
c self gravitational collapse, november, 1999
c
c written by: PSLi
c
c PURPOSE: Sets up a spherical region at a specified point on
c the grid (x10, x20, x30) with the specified radius (r) whose density
c differ from the rest of the initial grid.
c
c LOCAL VARIABLES:
c r initial radius of overpressured region
c x10,x20,x30 coordinates of centre of overpressured region.
c drat ratio of density across blast front
c prat ratio of pressure across blast front
c d0 density in ambient medium (default = 1.0)
c p0 pressure in ambient medium (default = 0.6)
c e0 internal energy in ambient medium (default = 0.9)
c er0 radiation energy in ambient medium (default = 1.0)
c v10 1-velocity in ambient medium
c v20 2-velocity in ambient medium
c v30 3-velocity in ambient medium
c b10 1-magnetic field on entire grid
c b20 2-magnetic field on entire grid
c b30 3-magnetic field on entire grid
c d1 density in central region (default = 1.0)
c p1 pressure in central region (default = 0.6)
c e1 internal energy in central region (default = 0.9)
c er1 radiation energy in central region (default = 30.)
c v11 1-velocity in central region
c v21 2-velocity in central region
c v31 3-velocity in central region
c m,drs,drc parameters for specifying a sphere whose surface is
c sinusoidally perturbed (spherical coordinates only
c For an unperturbed sphere, set all values to zero
c (default).
c
c EXTERNALS:
c OVERLAP [XYZ only]
c BNDYALL
c BSETMAG
c
c-----------------------------------------------------------------------
c
use real_prec
use config
use param
use field
use bndry
use grid
use root
use scratch
use cons
#ifdef MPI_USED
use mpiyes
#else
use mpino
#endif
use mpipar
c
integer :: i, j, k, ip1, jp1, kp1, m
c
real(rl) :: r , x10, x20 , x30 , drat,
. prat , d0 , p0 , e0 , v10,
. v20 , v30, b10 , b20 , b30,
. d1 , p1 , e1 , v11 , v21,
. v31 , drs, cofrac , mass ,
. drc , rsq, rin , rout , frac,
. er0 , er1, ros_mfp, dx_min, flx_lim,
. dmc_max
c
integer :: iin (ijkn), iout(ijkn), jin (ijkn),
1 jout(ijkn), kin (ijkn), kout(ijkn)
c
real(rl) :: massk(ijkn), sasum, overlap
c
implicit NONE
c
namelist / pgen /
. r , x10, x20 , x30 , drat,
. prat, d0 , p0 , e0 ,er0,
. v10 , v20, v30,
. b10 , b20, b30,
. d1 , p1 , e1 , er1,
. v11 , v21, v31,
. drs , drc, m
c
c-----------------------------------------------------------------------
c
r = 1.0
x10 = 0.0
x20 = 0.0
x30 = 0.0
drat = 0.0
prat = 0.0
d0 = 1.0
p0 = 0.6
e0 = 0.0
v10 = 0.0
v20 = 0.0
v30 = 0.0
b10 = 0.0
b20 = 0.0
b30 = 0.0
d1 = 1.0
p1 = 0.6
e1 = 0.0
v11 = 0.0
v21 = 0.0
v31 = 0.0
drs = 0.0
drc = 0.0
m = 0
er0 = 0.0
er1 = 0.0
c
if (myid .eq. 0) then
read (1, pgen)
write (2, pgen)
#ifdef MPI_USED
buf_in( 1) = r
buf_in( 2) = x10
buf_in( 3) = x20
buf_in( 4) = x30
buf_in( 5) = drat
buf_in( 6) = prat
buf_in( 7) = d0
buf_in( 8) = p0
buf_in( 9) = e0
buf_in(10) = v10
buf_in(11) = v20
buf_in(12) = v30
buf_in(13) = b10
buf_in(14) = b20
buf_in(15) = b30
buf_in(16) = d1
buf_in(17) = p1
buf_in(18) = e1
buf_in(19) = v11
buf_in(20) = v21
buf_in(21) = v31
buf_in(22) = drs
buf_in(23) = drc
buf_in(24) = er0
buf_in(25) = er1
ibuf_in( 1) = m
#endif
endif ! myid
#ifdef MPI_USED
call MPI_BCAST( buf_in, 25, MPI_FLOAT
& , 0, comm3d, ierr )
call MPI_BCAST( ibuf_in, 1, MPI_INTEGER
& , 0, comm3d, ierr )
if (myid .ne. 0) then
r = buf_in( 1)
x10 = buf_in( 2)
x20 = buf_in( 3)
x30 = buf_in( 4)
drat = buf_in( 5)
prat = buf_in( 6)
d0 = buf_in( 7)
p0 = buf_in( 8)
e0 = buf_in( 9)
v10 = buf_in(10)
v20 = buf_in(11)
v30 = buf_in(12)
b10 = buf_in(13)
b20 = buf_in(14)
b30 = buf_in(15)
d1 = buf_in(16)
p1 = buf_in(17)
e1 = buf_in(18)
v11 = buf_in(19)
v21 = buf_in(20)
v31 = buf_in(21)
drs = buf_in(22)
drc = buf_in(23)
er0 = buf_in(24)
er1 = buf_in(25)
m = ibuf_in( 1)
endif ! myid
#endif
c
c Set up atmosphere.
c
if(e0 .ne. 0.0) p0 = e0 * gamm1
if(e0 .eq. 0.0) e0 = p0 / gamm1
if(lrad .ne. 0) er0 = rad_con * (gamm1*e0*mmw*mh/(d0*boltz))**4
C print *, myid_w,x1a(is),x1a(ie+1)
do 30 k=ks,ke
do 20 j=js,je
do 10 i=is,ie
d (i,j,k) = d0
v1(i,j,k) = v10
v2(i,j,k) = v20
v3(i,j,k) = v30
if(xiso .eqv. .false.) e (i,j,k) = e0
if(lrad .ne. 0) er(i,j,k) = er0
10 continue
20 continue
30 continue
if(xmhd) then
do 31 k=ks-2,ke+3
do 21 j=js-2,je+3
do 11 i=is-2,ie+3
b1(i,j,k) = b10
b2(i,j,k) = b20
b3(i,j,k) = b30
11 continue
21 continue
31 continue
endif ! xmhd
c
c Set up central region.
c
do 40 i=is,ie
ip1 = i + 1
if ( abs(x1a(i)-x10) .lt. abs(x1a(ip1)-x10) ) then
iin (i) = i
iout(i) = ip1
else
iin (i) = ip1
iout(i) = i
endif
40 continue
c
do 50 j=js,je
jp1 = j + 1
if ( abs(x2a(j)-x20) .lt. abs(x2a(jp1)-x20) ) then
jin (j) = j
jout(j) = jp1
else
jin (j) = jp1
jout(j) = j
endif
50 continue
c
do 60 k=ks,ke
kp1 = k + 1
if ( abs(x3a(k)-x30) .lt. abs(x3a(kp1)-x30) ) then
kin (k) = k
kout(k) = kp1
else
kin (k) = kp1
kout(k) = k
endif
massk(k) = 0.0
60 continue
c
if (drat .ne. 0.0) d1 = d0 * drat
if (prat .ne. 0.0) p1 = p0 * prat
if (e1 .ne. 0.0) then
p1 = e1 * gamm1
else
e1 = p1 / gamm1
endif
if(lrad .ne. 0) er1 = rad_con * (gamm1*e1*mmw*mh/(d1*boltz))**4
c
do 90 k=ks,ke
do 80 j=js,je
do 70 i=is,ie
if(lgeom .eq. 1) then
rsq = r**2
rin = ( x1a(iin (i)) - x10 )**2
1 + ( x2a(jin (j)) - x20 )**2
2 + ( x3a(kin (k)) - x30 )**2
rout = ( x1a(iout(i)) - x10 )**2
1 + ( x2a(jout(j)) - x20 )**2
2 + ( x3a(kout(k)) - x30 )**2
endif ! CARTESIAN
if(lgeom .eq. 2) then
rsq = r**2
rin = ( x1a(iin (i)) - x10 )**2
1 + ( x2a(jin (j)) - x20 )**2
rout = ( x1a(iout(i)) - x10 )**2
1 + ( x2a(jout(j)) - x20 )**2
endif ! CYLINDRICAL
if(lgeom .eq. 3) then
rsq = r**2 * ( 1.0 + drs * sin (m * x2a(j))
1 + drc * cos (m * x2a(j)) )**2
rin = ( x1a(iin (i)) - x10 )**2
rout = ( x1a(iout(i)) - x10 )**2
endif ! SPHERICAL
if ( (rin .lt. rsq) .and. (rout .le. rsq) ) then
d (i,j,k) = d1
v1(i,j,k) = v11
v2(i,j,k) = v21
v3(i,j,k) = v31
if(xiso .eqv. .false.) e (i,j,k) = e1
if(lrad .ne. 0) er(i,j,k) = er1
massk(k) = massk(k) + d1 * dvl1a(i) * dvl2a(j) * dvl3a(k)
endif
if ( (rin .lt. rsq) .and. (rout .gt. rsq) ) then
if(lgeom .eq. 1) then
frac = overlap ( 1, r, x10, x20, x30
1 , x1a(iin (i)), x2a(jin (j)), x3a(kin (k))
2 , x1a(iout(i)), x2a(jout(j)), x3a(kout(k)) )
else ! lgeom
frac = ( rsq - rin ) / ( rout - rin )
endif ! lgeom
cofrac = 1.0 - frac
d(i,j,k) = d1 * frac + d0 * cofrac
if(xiso .eqv. .false.) e(i,j,k) = e1 * frac + e0 * cofrac
if(lrad .ne. 0) er(i,j,k) = er1 * frac + er0 * cofrac
massk(k) = massk(k)
1 + d1 * frac * dvl1a(i) * dvl2a(j) * dvl3a(k)
endif
70 continue
80 continue
90 continue
mass = SASUM ( nx3z, massk(ks), 1 )
C write(*,"('total mass = ',1pd12.4)")mass
c
call gravity
c
return
end
c
c=======================================================================
c
c \\\\\\\\\\ E N D S U B R O U T I N E //////////
c ////////// B L A S T \\\\\\\\\\
c
c=======================================================================
c
c
c=======================================================================
c
c \\\\\\\\\\ B E G I N F U N C T I O N //////////
c ////////// O V E R L A P \\\\\\\\\\
c
c=======================================================================
c
real(rl) function overlap ( ishp, rad, x0, y0, z0, xin, yin
1 , zin, xout, yout, zout )
c
c dac:zeus3d.overlap <--------- overlap of region over Cartesian zone
c april, 1990
c
c written by: David Clarke
c modified 1:
c
c PURPOSE: Determines the fraction of a Cartesian zone that overlaps
c the specified geometrical region (sphere or right cylinder). This
c is done by dividing the zone into 20**3 "subzones", and finding the
c fraction of subzone centres lying inside the surface of the region.
c
c INPUT VARIABLES:
c
c ishp =1 => sphere
c =2 => right cylinder
c rad radius of region
c x0,y0,z0 coordinates of centre of curvature.
c xin,yin,zin coordinates of zone corner known to lie inside
c region.
c xout,yout,zout coordinates of zone corner diametrically opposed to
c zone corner known to lie inside region.
c
c OUTPUT VARIABLES:
c
c LOCAL VARIABLES:
c
c EXTERNALS: [NONE]
c
c-----------------------------------------------------------------------
c
use real_prec
use config
use param
c
integer :: i, j, k, nx, ny, nz, ishp
c
real(rl) :: rad , x0 , y0 , z0 , xin ,
. yin , zin , xout, yout, zout,
. delx , dely, delz, r , fact,
. scount
c
real(rl) :: xsq(20), ysq(20), zsq(20), count(20)
c
implicit NONE
c-----------------------------------------------------------------------
c
if(lgeom .eq. 1) then
c Number of subzones in the x-direction is "nx", etc. for "ny" and
c "nz". Increment between subzones in x-direction is "delx", etc. for
c "dely" and "delz".
c
nx = 20
ny = 20
nz = 20
delx = ( xout - xin ) / real( nx )
dely = ( yout - yin ) / real( ny )
delz = ( zout - zin ) / real( nz )
c
c Set up subgrid inside zone.
c
do 10 i=1,nx
xsq (i) = ( xin + ( 0.5 + real(i-1) ) * delx - x0 )**2
count(i) = 0.0
10 continue
do 20 j=1,ny
ysq (j) = ( yin + ( 0.5 + real(j-1) ) * dely - y0 )**2
20 continue
do 30 k=1,nz
zsq (k) = ( zin + ( 0.5 + real(k-1) ) * delz - z0 )**2
30 continue
c
c Count the number of subzones lying inside the surface of the
c region which passes through the zone.
c
fact = 1.0
if (ishp .eq. 2) fact = 0.0
do 60 k=1,nz
do 50 j=1,ny
do 40 i=1,nx
r = sqrt ( fact * xsq(i) + ysq(j) + zsq(k) )
if (r .le. rad) count(i) = count(i) + 1.0
40 continue
50 continue
60 continue
scount = 0.0
do 70 i=1,nx
scount = scount + count(i)
70 continue
scount = max ( one, scount )
c
c Set the fraction of the zone which overlaps the region.
c
overlap = scount / real( nx * ny * nz )
else ! lgeom
overlap = 1.0
endif ! lgeom
c
return
end
c
c=======================================================================
c
c \\\\\\\\\\ E N D F U N C T I O N //////////
c ////////// O V E R L A P \\\\\\\\\\
c
c=======================================================================
c
c
Forgot your password?
Click here to get the file
The Computational Astrophysics Data Analysis Center