-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathsetup_dustybox.f90
132 lines (124 loc) · 3.13 KB
/
setup_dustybox.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
!----------------------------------------------------------------
! Set up a box for a dustybox test
! should work in 1, 2 and 3 dimensions
!----------------------------------------------------------------
subroutine setup
!
!--include relevant global variables
!
use dimen_mhd
use debug
use loguns
use bound
use eos
use options
use part
use setup_params
use uniform_distributions, only:set_uniform_cartesian
use mem_allocation, only:alloc
!
!--define local variables
!
implicit none
integer :: i
integer :: ntypes
integer :: ngas,npdust,jtype
real, dimension(ndimV) :: Bzero
real :: massp,masspdust
real :: spsoundi
real :: denszero,uuzero,przero,denszerodust
real :: dust_to_gas_ratio
!
!--setup parameters (could read in from a file)
!
Bzero = 0.
denszero = 1.0
uuzero = 1.0
if (idust.eq.1) then
ntypes = 1
else
ntypes = 2
endif
!
!--set boundaries
!
ibound = 3 ! periodic boundaries
nbpts = 0 ! use ghosts not fixed
xmin(:) = 0. ! set position of boundaries
xmax(1) = 1.0
if (ndim.GE.2) then
xmax(2:ndim) = 11.*psep ! would need to adjust this depending on grid setup
endif
!
!--initially set up a uniform density grid (also determines npart)
! (the call to set_uniform_cartesian means this works in 1,2 and 3D)
!
ngas = 0
npdust = 0
do jtype=1,ntypes
call set_uniform_cartesian(1,psep,xmin,xmax,adjustbound=.true.)
if (jtype.eq.1) then
ngas = npart
itype(1:ngas) = itypegas
elseif (jtype.eq.2) then
npdust = npart - ngas
itype(ngas+1:ngas+npdust) = itypedust
endif
enddo
massp = 1.0/FLOAT(ngas) ! average particle mass
dust_to_gas_ratio = 1.
if (idust.eq.1) then ! one fluid dust
massp = massp*(1. + dust_to_gas_ratio)
endif
! two fluid dust (if npdust > 0)
masspdust = 0.
if (npdust.gt.0) masspdust = dust_to_gas_ratio*1.0/FLOAT(npdust) ! average particle mass
denszerodust = dust_to_gas_ratio*denszero
if (ntypes.gt.1) print*,' ngas = ',ngas,' npdust = ',npdust
!
!--allocate memory here
!
call alloc(npart)
!
!--setup uniform density grid of particles
!
do i=1,npart
vel(:,i) = 0.
if (itype(i).eq.itypedust) then
dens(i) = denszerodust
pmass(i) = masspdust
uu(i) = 0.
else
if (idust.eq.1) then
dustfrac(1,i) = (dust_to_gas_ratio/(1. + dust_to_gas_ratio))
deltav(1,i) = -1. ! deltav = vdust - vgas
vel(1,i) = 0.5 ! v = vg + rhod/rho*deltav
else
vel(1,i) = 1. !--gas velocity
endif
dens(i) = denszero
pmass(i) = massp
uu(i) = uuzero
endif
if (imhd.GT.0) then
Bfield(:,i) = Bzero
else
Bfield(:,i) = 0.
endif
enddo
ntotal = npart
!
!--get sound speed from equation of state (want average sound speed, so
! before the density is perturbed)
!
call equation_of_state1(przero,spsoundi,uuzero,denszero)
print*,' gamma = ',gamma
print*,' pr = ',przero,' cs = ',spsoundi,' u = ',uu(1),' dens = ',dens(1)
return
end subroutine setup
!
! use this routine to modify the dump upon code restart
!
subroutine modify_dump()
implicit none
end subroutine modify_dump