Blog: Binary simulation: problem.f90

File problem.f90, 10.0 KB (added by Zhuo Chen, 11 years ago)
Line 
1!#########################################################################
2!
3! Copyright (C) 2003-2012 Department of Physics and Astronomy,
4! University of Rochester,
5! Rochester, NY
6!
7! problem.f90 of module Binary is part of AstroBEAR.
8!
9! AstroBEAR is free software: you can redistribute it and/or modify
10! it under the terms of the GNU General Public License as published by
11! the Free Software Foundation, either version 3 of the License, or
12! (at your option) any later version.
13!
14! AstroBEAR is distributed in the hope that it will be useful,
15! but WITHOUT ANY WARRANTY; without even the implied warranty of
16! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17! GNU General Public License for more details.
18!
19! You should have received a copy of the GNU General Public License
20! along with AstroBEAR. If not, see <http://www.gnu.org/licenses/>.
21!
22!#########################################################################
23!> @dir OrbitingParticles
24!! @brief Contains files necessary for the Orbiting Particles test problem
25
26!> @file problem.f90
27!! @brief Main file for module Problem
28
29!> @defgroup OrbitingParticles Orbiting Particle Module
30!! @brief Module for setting up orbiting particles
31!! @ingroup Modules
32
33!> Orbiting Particle Module
34!! @ingroup OrbitingParticles
35MODULE Problem
36 USE DataDeclarations
37 USE ParticleDeclarations
38 USE Ambients
39 USE Disks
40 USE Winds
41 IMPLICIT NONE
42 SAVE
43 PUBLIC ProblemModuleInit, ProblemGridInit, &
44 ProblemBeforeStep, ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep
45 REAL(KIND=qPREC) :: alpha, densw, velw, t1, t2
46 INTEGER :: radiusw, nWinds
47 LOGICAL :: windsPresent=.false., disk=.false.
48 !TYPE(pWindDef), DIMENSION(:), ALLOCATABLE :: MyWinds
49
50CONTAINS
51
52 !> Initializes module variables
53 SUBROUTINE ProblemModuleInit()
54 TYPE(InfoDef) :: Info
55 TYPE(WindDef), POINTER :: Wind
56 INTEGER :: nParticles, edge
57 REAL(KIND=qPREC) :: mass=0
58 REAL(KIND=qPREC) :: xloc(3)
59 REAL(KIND=qPREC) :: vel(3)
60 TYPE(ParticleDef), POINTER :: Particle
61 TYPE(ParticleListDef), POINTER :: particlelist
62 INTEGER :: i, grav_soft_rad
63 INTEGER :: ids(2)
64 TYPE(AmbientDef), POINTER :: Ambient
65 TYPE(DiskDef), POINTER :: mydisk
66 REAL(KIND=qPREC) :: time, rhoOut, pOut, vxOut, vyOut, vzOut, BxOut, ByOut, BzOut, &
67 buff(17)
68 NAMELIST /AmbientData/ rhoOut, pOut, vxOut, vyOut, vzOut, BxOut, ByOut, BzOut
69 NAMELIST /ProblemData/ nParticles, densw, velw, windsPresent, disk, t1, t2
70 NAMELIST /ParticleData/ mass,xloc,vel,alpha,radiusw,grav_soft_rad, buff
71 NAMELIST /RestartData/ ids
72
73! time=levels(Info%level)%tnow
74
75 OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', STATUS="OLD")
76 READ(PROBLEM_DATA_HANDLE,NML=ProblemData)
77
78if (.not. lrestart) then
79 CALL CreateAmbient(Ambient)
80 READ(PROBLEM_DATA_HANDLE,NML=AmbientData)
81 if (disk) then
82 Ambient%density=1d0
83 else ; Ambient%density=rhoOut
84 end if
85 Ambient%pressure=Ambient%density
86 Ambient%B(:)=(/BxOut, ByOut, BzOut/)
87 Ambient%velocity(:)=(/vxOut, vyOut, vzOut/)
88end if
89
90if (windsPresent) then
91 DO i=1,nDim
92 DO edge=1,2
93 IF (Gmthbc(i,edge) == 1) THEN
94 NULLIFY(Wind)
95 CALL CreateWind(Wind)
96 Wind%dir=i
97 Wind%edge=edge
98 Wind%type=OUTFLOW_ONLY
99 END IF
100 END DO
101 END DO
102! ALLOCATE(MyWinds(6))
103! nWinds=0
104! DO i=1,nDim
105! DO edge=1,2
106! IF (Gmthbc(i,edge) == 1) THEN
107! nWinds=nWinds+1
108! CALL CreateWind(MyWinds(nWinds)%p)
109! MyWinds(nWinds)%p%dir=i
110! MyWinds(nWinds)%p%edge=edge
111! MyWinds(nWinds)%p%type=OUTFLOW_ONLY
112! END IF
113! END DO
114! END DO
115end if
116
117 IF (lRestart) THEN
118 particleList=>SinkParticles
119 i=1
120 DO WHILE (ASSOCIATED(particlelist))
121 particle=>particlelist%self
122 READ(PROBLEM_DATA_HANDLE,NML=ParticleData)
123
124 if (i==1) then
125 Particle%PointGravityObj%soft_length=REAL(grav_soft_rad,qPREC)*sink_dx
126 Particle%PointGravityObj%alpha = alpha
127 if (alpha==0d0) then
128 Particle%iAccrete=NOACCRETION
129 else
130 Particle%iAccrete=KRUMHOLZ_ACCRETION
131 ! see particle_declarations
132 end if
133 end if
134
135 !if (time.ge.t1 .and. time.le.t2) then
136 ! Particle%PointGravityObj%alpha = alpha*(time-t1)/(t2-t1) ! (0,1)
137 !else ; Particle%PointGravityObj%alpha = 0d0
138 !end if
139
140 IF (i==2) THEN
141 CALL CreateOutflowObject(Particle%OutflowObj)
142 Particle%OutflowObj%duration = 1.e30
143 Particle%OutflowObj%radius =REAL(radiusw,qPREC)*sink_dx
144 Particle%OutflowObj%thickness =REAL(radiusw,qPREC)*sink_dx
145 Particle%OutflowObj%open_angle = Pi
146 Particle%OutflowObj%density = densw
147 Particle%OutflowObj%temperature = 1d0
148 Particle%OutflowObj%velocity = velw
149 Particle%OutflowObj%source_vel(1:nDim) = Particle%Q(imom(1:nDim))
150 Particle%OutflowObj%position= Particle%xloc
151 CALL UpdateOutflow(Particle%OutflowObj)
152 END IF
153 particlelist=>particlelist%next
154 i=i+1
155 END DO
156
157 ELSE!lrestart
158
159 DO i=1,nParticles
160 READ(PROBLEM_DATA_HANDLE,NML=ParticleData)
161 NULLIFY(Particle)
162 CALL CreateParticle(Particle)
163 Particle%Q(1)=mass
164 Particle%xloc=xloc
165 Particle%Q(imom(1:nDim))=vel(1:nDim)
166 Particle%Buffer=buff
167 CALL AddSinkParticle(Particle)
168
169 if (i==1) then
170 Particle%iAccrete=KRUMHOLZ_ACCRETION
171 CALL CreatePointGravityObject(Particle%PointGravityObj)
172 Particle%PointGravityObj%alpha = alpha
173 Particle%PointGravityObj%soft_length = REAL(grav_soft_rad,qPREC)*sink_dx
174 Particle%PointGravityObj%soft_function =SPLINESOFT
175 Particle%PointGravityObj%Mass =Particle%Q(1)
176 Particle%PointGravityObj%v0(1:nDim) =Particle%Q(imom(1:nDim))
177 Particle%PointGravityObj%x0 =Particle%xloc
178 end if
179!
180!!!
181if (disk) then
182ALLOCATE(MyDisk)
183! mydisk%HeightProfile=0
184 !mydisk%HeightProfile=1 !flared
185 mydisk%soft_length=REAL(grav_soft_rad,qPREC)*sink_dx
186 mydisk%soft_function=SPLINESOFT
187 mydisk%density=densw
188 mydisk%pressure= &!mydisk%density
189 Ambient%pressure
190 mydisk%velocity=Particle%Q(imom(1:nDim))
191 !
192 mydisk%theta=0d0 !5feb12
193 !mydisk%theta=10d0*(Pi/180d0) !4feb12
194 !mydisk%theta=Pi*.5d0 !3feb12
195 !
196 mydisk%phi=0d0
197 mydisk%position=Particle%xloc
198 mydisk%thickness=0d0
199 mydisk%radius= 2d0 !6dec
200 !1d0
201 !6d0*REAL(radiusw,qPREC)*sink_dx
202 !24d0*sink_dx
203 ! if (.not.mydisk%HeightProfile)
204 mydisk%height=2d0*REAL(grav_soft_rad,qPREC)*sink_dx
205 !mydisk%height=mydisk%radius*.1d0
206 mydisk%central_mass=mass
207 CALL UpdateDisk(mydisk)
208end if!lock
209!!!
210!
211 IF (i.eq.2) THEN
212 CALL CreateOutflowObject(Particle%OutflowObj)
213 Particle%OutflowObj%duration = 1.e30
214 Particle%OutflowObj%radius =REAL(radiusw,qPREC)*sink_dx
215 Particle%OutflowObj%thickness =REAL(radiusw,qPREC)*sink_dx
216 Particle%OutflowObj%open_angle = Pi
217 Particle%OutflowObj%density = densw
218 Particle%OutflowObj%temperature = 1d0
219 Particle%OutflowObj%velocity = velw
220 Particle%OutflowObj%source_vel(1:nDim) = Particle%Q(imom(1:nDim))
221 Particle%OutflowObj%position= Particle%xloc
222 CALL UpdateOutflow(Particle%OutflowObj)
223 END IF
224 END DO
225 END IF!lrestart=0/1
226CLOSE(PROBLEM_DATA_HANDLE)
227
228 END SUBROUTINE ProblemModuleInit
229
230 !> Applies initial conditions
231 !! @param Info Info object
232 SUBROUTINE ProblemGridInit(Info)
233 TYPE(InfoDef) :: Info
234 END SUBROUTINE ProblemGridInit
235
236 !> Applies Boundary conditions
237 !! @param Info Info object
238 SUBROUTINE ProblemBeforeStep(Info)
239 TYPE(InfoDef) :: Info
240 INTEGER :: i, nParticles, grav_soft_rad
241 REAL(KIND=xprec) :: beta,time, mass, xloc(3), vel(3),buff(17)
242 TYPE(ParticleListDef), POINTER :: particlelist
243 TYPE(ParticleDef), POINTER :: Particle
244 LOGICAL,SAVE :: lock=.false.
245 NAMELIST /ProblemData/ nParticles, densw, velw, windsPresent, disk, t1, t2
246 NAMELIST /ParticleData/ mass,xloc,vel,alpha,radiusw,grav_soft_rad, buff
247
248 return
249
250 time=levels(Info%level)%tnow
251 if (time.lt.t1) then ; return
252 else if (time.ge.t1 .and. time.le.t2) then
253 beta = (time-t1)/(t2-t1) ! 0<= beta <=1
254 else ; beta=1d0
255 end if
256
257 OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', STATUS="OLD")
258 READ(PROBLEM_DATA_HANDLE,NML=ProblemData)
259 particleList=>SinkParticles
260 i=1
261 DO WHILE (ASSOCIATED(particlelist))
262 particle=>particlelist%self
263 READ(PROBLEM_DATA_HANDLE,NML=ParticleData)
264 if (i==1) then
265 Particle%PointGravityObj%alpha =alpha*beta
266 !print*,'t1,t2i,a*b=',t1,t2,i,alpha*beta
267 end if
268 particlelist=>particlelist%next
269 i=i+1
270 END DO
271 CLOSE(PROBLEM_DATA_HANDLE)
272
273 if (time.gt.t2 .and. Particle%PointGravityObj%alpha.eq.1d0) lock=.true.
274
275 END SUBROUTINE ProblemBeforeStep
276
277 !> Could be used to update grids pre-output
278 !! @param Info Info Object
279 SUBROUTINE ProblemAfterStep(Info)
280 TYPE(InfoDef) :: Info
281 END SUBROUTINE ProblemAfterStep
282
283 !> Could be used to set force refinement
284 !! @param Info Info object
285 SUBROUTINE ProblemSetErrFlag(Info)
286 TYPE(InfoDef) :: Info
287 END SUBROUTINE ProblemSetErrFlag
288
289 SUBROUTINE ProblemBeforeGlobalStep(n)
290 INTEGER :: n
291 END SUBROUTINE ProblemBeforeGlobalStep
292
293END MODULE Problem
294