Blog: 2.5D Pulsed Jet: problem.f90

File problem.f90, 16.1 KB (added by Baowei Liu, 13 years ago)

problem.f90

Line 
1!> @dir PulsedJet
2!! @brief Contains files necessary for module PulsedJet
3
4!> @file problem.f90
5!! @brief Main file for module PulsedJet
6
7!> @defgroup PulsedJet PulsedJet Module
8!! @brief Module for simulating a pulsed 2.5D Jet
9!! @ingroup Modules
10
11!> Pulsed Jet problem module
12!! @ingroup 1DJet
13MODULE Problem
14 USE GlobalDeclarations
15 USE PhysicsDeclarations
16 USE DataDeclarations
17 USE CoolingSrc
18 USE ZCooling
19 USE Chemistry
20 USE Refinements
21 USE Shapes
22 USE Fields
23 USE Emissions
24 USE Projections
25 USE ProcessingDeclarations
26 USE Ambients
27 IMPLICIT NONE
28 PUBLIC ProblemModuleInit, ProblemGridInit, &
29 ProblemBeforeStep, ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep
30
31 REAL(KIND=qPREC) :: namb, njet, tamb, vjet, Rjet, Bm, vrange, newvjet, period, pamb, alpha, &
32 Betam, Rm, freq, amp, yupper, theta, ycoeff, xzcoeff, precfreq, precperiod, CoolingRes
33 REAL(KIND=qPREC), ALLOCATABLE :: randoms(:), vtimes(:)
34 INTEGER :: Bgeom, nrandoms, seed, vtype, simtype, N
35 INTEGER, PARAMETER :: toroidal = 1, helical = 2
36 INTEGER, PARAMETER :: constant = 0, sinusoidal = 1, random = 2
37 INTEGER, PARAMETER :: single = 0, colliding_same = 1, colliding_opp = 2
38 LOGICAL :: ljet
39
40
41CONTAINS
42
43 !> Initializes module variables
44 SUBROUTINE ProblemModuleInit()
45 TYPE(AmbientDef), POINTER :: Ambient
46 TYPE(RefinementDef), POINTER :: InjectRefine, TopInjectRefine, CoolingLengthRefinement
47 TYPE(ProjectionDef), POINTER :: Projection
48 REAL(KIND=qPREC) :: rand
49 INTEGER :: i, jj
50
51 NAMELIST/problemdata/ namb, tamb, njet, vjet, vrange, Rjet, Rm, Betam, Bgeom, vtype, freq, &
52 theta, precfreq, seed, simtype, N, CoolingRes
53
54 OPEN(PROBLEM_DATA_HANDLE,FILE='problem.data',STATUS='old')
55 READ(PROBLEM_DATA_HANDLE,NML=problemdata)
56 CLOSE(PROBLEM_DATA_HANDLE)
57
58 ! CALL AddDiagnosticVar(OI_Field, 'OI') ! 6300
59 ! CALL AddDiagnosticVar(NII_Field, 'NII') ! 6583
60 ! CALL AddDiagnosticVar(SII_6716_Field, 'SII_6716')
61 ! CALL AddDiagnosticVar(SII_6731_Field, 'SII_6731')
62 ! CALL AddDiagnosticVar(Halpha_Field, 'Halpha')
63 ! CALL AddDiagnosticVar(CellSize_Field, 'dx')
64 ! CALL AddDiagnosticVar(CoolingLength_Field, 'CoolingLength')
65
66 ! Set up some parameters for the velocity pulses
67 SELECT CASE(vtype)
68 CASE(constant)
69 CASE(sinusoidal)
70 period = final_time/freq
71 amp = vrange/vjet
72 CASE(random)
73 ! set up an array of pseudo-random numbers to use for random jet velocities
74 nrandoms = NINT(freq)
75 period = final_time/REAL(nrandoms,qPrec)
76 ALLOCATE(randoms(nrandoms))
77 ALLOCATE(vtimes(nrandoms))
78 CALL RANDOM_SEED(seed)
79 DO jj=1, nrandoms
80 CALL RANDOM_NUMBER(rand)
81 randoms(jj) = rand
82 vtimes(jj) = jj*period
83 END DO
84 END SELECT
85
86 pamb = namb*Boltzmann*tamb
87 IF(lMHD) THEN
88 Bm = SQRT(8d0*Pi*pamb/Betam)/Bscale
89 alpha = 1d0 - (Rm/Rjet)**2d0/Betam
90 ELSE
91 Bm = 0d0
92 alpha = 1d0
93 Rm = -10d0
94 END IF
95
96 ! Set up parameters related to a precessing injection velocity
97 theta = theta*Pi/180d0
98 ycoeff = SQRT(1d0/(1d0 + 2d0*TAN(theta)**2d0))
99 xzcoeff = SQRT(1d0 - ycoeff**2d0)
100 precperiod = final_time/precfreq
101
102 ! Set up ambient
103 CALL CreateAmbient(Ambient)
104 Ambient%density= namb/nScale
105 Ambient%pressure= pamb/pScale
106
107 ! Set up refinement
108 CALL ClearAllRefinements()
109
110 CALL CreateRefinement(CoolingLengthRefinement)
111 CoolingLengthRefinement%field=CoolingTime_Field
112 CoolingLengthRefinement%limit = LESSTHAN
113 CoolingLengthRefinement%Threshold(0:MaxLevel)=CoolingRes*levels(0:MaxLevel)%dx/(0.1*vjet/VelScale)
114
115 CALL CreateRefinement(InjectRefine)
116 CALL CreateShape(InjectRefine%Shape)
117 InjectRefine%Shape%type = RECTANGULAR_PRISM
118 InjectRefine%Shape%position = (/ 0d0, 0d0, 0d0 /)
119 InjectRefine%Shape%size_param = (/Rjet, 1*levels(0)%dx, Rjet/)
120 InjectRefine%Maxlevel = 3
121 CALL SetShapeBounds(InjectRefine%Shape)
122
123 IF (simtype /= single) THEN
124 yupper = GxBounds(2,2)
125 CALL CreateRefinement(TopInjectRefine)
126 CALL CreateShape(TopInjectRefine%Shape)
127 TopInjectRefine%Shape%type = RECTANGULAR_PRISM
128 TopInjectRefine%Shape%position = (/ 0d0, yupper, 0d0 /)
129 TopInjectRefine%Shape%size_param = Rjet
130 InjectRefine%Maxlevel = 3
131 CALL SetShapeBounds(TopInjectRefine%Shape)
132 END IF
133
134 ! Set up projections to make emission maps
135 CALL CreateProjection(Projection)
136 IF (iCylindrical == 0) THEN
137 Projection%dim = 3
138 ELSE
139 Projection%dim = CYLINDRICAL_PROJECTION
140 END IF
141 Projection%PlotLevel = MIN(MaxLevel, 4)
142 Projection%field(1)%id = Halpha_Field
143 Projection%field(1)%name = "Halpha"
144 Projection%field(2)%id = SII_6716_Field
145 Projection%field(2)%name = "SII_6716"
146 Projection%field(3)%id = SII_6731_Field
147 Projection%field(3)%name = "SII_6731"
148 ! Projection%field(4)%id = NII_Field
149 ! Projection%field(4)%name = "NII"
150
151 END SUBROUTINE ProblemModuleInit
152
153 !> Applies initial conditions
154 !! @param Info Info object
155 SUBROUTINE ProblemGridInit(Info)
156 TYPE(InfoDef) :: Info
157 CALL PlaceJet(Info, levels(Info%level)%dx)
158 END SUBROUTINE ProblemGridInit
159
160 !> Applies Boundary conditions
161 !! @param Info Info object
162 SUBROUTINE ProblemBeforeStep(Info)
163 TYPE(InfoDef) :: Info
164 CALL PlaceJet(Info, 0d0)
165 END SUBROUTINE ProblemBeforeStep
166
167 SUBROUTINE PlaceJet(Info, yjet)
168 TYPE(InfoDef) :: Info
169 REAL(KIND=qPREC), POINTER, DIMENSION(:,:,:,:) :: q
170 INTEGER :: nn, i, j, k, rmbc, zrmbc, level, mx, my, mz, ii, jj, kk, Nz
171 REAL(KIND=qPREC) :: x, y , z, r, phi, dx, xlower, ylower, zlower, dt, time,&
172 Bx, By, Bz, randv, pjet, Bjet, nH, nHe, Ay, ptot, padd, ddx, xx, yy, zz, rr, yjet, hdx, pB, sumpB, ddz
173 REAL(KIND=qPREC),ALLOCATABLE :: A(:,:,:,:), BzAux(:,:)
174 REAL(KIND=qPREC), DIMENSION(1:NrHydroVars) :: newq, sumq
175
176
177 level=Info%level
178 q=>Info%q
179 rmbc = levels(level)%gmbc(levels(level)%step)
180 zrmbc = 0
181 mx = Info%mX(1)
182 my = Info%mX(2)
183 mz = Info%mX(3)
184 dx = levels(level)%dX
185 ddx = dx/N
186 ddz = 0
187 Nz = 1
188 IF (nDim == 3) THEN
189 zrmbc = rmbc
190 Nz=N
191 ddz=ddx
192 END IF
193
194 xlower = Info%xBounds(1,1)
195 ylower = Info%xBounds(2,1)
196 zlower = Info%xBounds(3,1)
197 hdx=half*dx
198
199 dt = levels(level)%dt
200 time = levels(level)%tnow
201 ! IF (time <= dt) RETURN ! don't want to overwrite initial conditions
202
203 phi = 0.0
204
205 ! sets new jet velocity depending on time and type of velocity changes
206 SELECT CASE(vtype)
207 CASE(constant)
208 newvjet = vjet
209 CASE(sinusoidal)
210 newvjet = vjet*(1d0 + amp*SIN(2d0*Pi*time/period))
211 phi = 2d0*Pi*time/precperiod
212 CASE(random)
213 DO nn=1, nrandoms
214 IF (time == vtimes(nn)) THEN
215 randv = ((randoms(nn)*2d0)-1d0)*vrange
216 newvjet = vjet + randv
217 ELSE IF (time < vtimes(1)) THEN
218 newvjet = vjet
219 END IF
220 END DO
221 END SELECT
222
223
224 DO j=1-rmbc, my+rmbc
225 y = ylower + (REAL(j,xPrec) - 1d0)*dx
226
227 IF (simtype /= single) THEN
228 IF (y >= 0d0+yjet .AND. y <= yupper-yjet-dx) cycle
229 ELSE
230 IF (y >= 0d0+yjet) THEN
231! write(*,*) j, y
232 cycle
233 END IF
234 END IF
235
236 DO i=1-rmbc, mx+rmbc
237 x = xlower + (REAL(i,xPrec) - 1d0)*dx
238 DO k=1-zrmbc, mz+zrmbc
239 z = zlower + (REAL(k,xPrec) - 1d0)*dx
240
241 sumq=0d0
242 sumpB=0d0
243 CALL Cons_To_Prim(q(i,j,k,:))
244
245 DO jj=1, N
246 yy = y + (REAL(jj,xPrec) - half)*ddx
247 ! ONLY DO FOLLOWING INSIDE THE -y BOUNDARY. ALSO,
248 ! ABOVE THE +y BOUNDARY IF THIS IS A COLLIDING JETS RUN
249 IF (simtype /= single) THEN
250 ljet = (yy <= 0d0+yjet .AND. yy >= yupper-yjet)
251 ELSE
252 ljet = (yy <= 0d0+yjet)
253! write(*,'(2I6,2E15.4,L2)') j,jj,y,yy,ljet
254 END IF
255 DO ii=1, N
256 xx = x + (REAL(ii,xPrec) - half)*ddx
257 DO kk=1, Nz
258
259 IF (ljet) THEN
260 IF (nDim == 2) THEN
261 rr = xx
262 ELSE IF (nDim == 3) THEN
263 zz = z + (REAL(kk,xPrec) - half)*ddz
264 rr = SQRT(xx**2d0+zz**2d0)
265 END IF
266
267 IF (rr < Rjet) THEN
268 newq(1) = njet/nScale
269 newq(ivx) = xzcoeff*COS(phi)*newvjet/velScale
270 newq(ivy) = ycoeff*newvjet/velScale
271 ! Second jet injection region, if this is a colliding jets run
272 IF (ivz /= 0) THEN
273 IF (simtype /= single .AND. y > yupper-yjet) newq(ivy) = -newq(ivy)
274 newq(ivz) = xzcoeff*SIN(phi)*newvjet/velScale
275 END IF
276 IF (rr < Rm) THEN
277 pjet = pamb*(alpha + 2d0/Betam*(1d0 - (rr/Rm)**2d0))/pScale
278 ELSE
279 pjet = alpha*pamb/pScale
280 END IF
281 newq(iE) = pjet !/gamma1 + half*SUM(newq(ivx:ivz)**2d0)/newq(1)
282 IF (rr < Rm) THEN
283 pB= .5*(Bm*rr/Rm)**2d0
284 ELSE
285 pB = .5*(Bm*Rm/rr)**2d0
286 END IF
287 IF (lTrackHydrogen) THEN
288 nH = newq(1)/(1d0+xH+HeFrac*(1d0+xHe+2d0*xHeII)+ZFrac*(1d0+xZ))
289 newq(iH) = (1d0-xH)*nH
290 newq(iHII) = xH*nH
291 IF (lTrackHelium) THEN
292 nHe = HeFrac*nH
293 newq(iHe) = (1d0-xHe-xHeII)*nHe
294 newq(iHeII) = xHe*nHe
295 newq(iHeIII) = xHeII*nHe
296 END IF
297 END IF
298
299 ELSE
300 newq=Info%q(i,j,k,1:NrHydroVars)
301 ! newq(1) = newq(1) + Info%q(i,j,k,1,1) !namb/nScale
302 ! newq(ivx:ivz) = newq(ivx:ivz) + 0d0
303 ! newq(iE) = newq(iE) + pamb/pScale/gamma1
304 pB=0d0
305 END IF
306 ELSE
307 newq=Info%q(i,j,k,1:NrHydroVars)
308 pB=0d0
309 END IF
310 sumq=sumq+newq
311 sumpB=sumpB+pB
312 END DO
313 END DO
314 END DO
315 q(i,j,k,1:NrHydroVars) = sumq/N**nDim
316 CALL Prim_To_Cons(q(i,j,k,:))
317 q(i,j,k,iE)=q(i,j,k,iE)+sumpB/N**nDim
318 END DO
319 END DO
320 END DO
321
322 IF (lMHD) THEN
323 ALLOCATE(A(1-rmbc:mx+rmbc+1,1-rmbc:my+rmbc+1,1-zrmbc:mz+zrmbc+1,3))
324 IF (nDim == 2) ALLOCATE(BzAux(1-rmbc:mx+rmbc+1,1-rmbc:my+rmbc+1))
325 A(:,:,:,1) = 0d0
326 A(:,:,:,3) = 0d0
327 DO j=1-rmbc, my+rmbc+1
328 y = ylower + (REAL(j,xPrec) - half)*dx
329 ljet = .false.
330 ! Check if inside jet injection region(s)
331 IF (simtype /= single) THEN
332 IF (y < 0d0+yjet .OR. y > yupper-yjet) ljet = .true.
333 ELSE
334 IF (y < 0d0+yjet) ljet = .true.
335 END IF
336 DO i=1-rmbc, mx+rmbc+1
337 x = xlower + (REAL(i,xPrec) - 1d0)*dx
338 DO k=1-zrmbc, mz+zrmbc+1
339 ! A M B I E N T
340 Ay = -Bm*Rm*(LOG(Rjet/Rm)+half)
341
342 IF (ljet) THEN
343 IF (nDim == 2) THEN
344 r = x
345 ELSE IF (nDim == 3) THEN
346 z = zlower + (REAL(k,xPrec) - 1d0)*dx
347 r = SQRT(x**2d0 + z**2d0)
348 END IF
349 IF (r < Rjet) THEN
350 IF (r < Rm) THEN
351 Ay = -half*Bm*r**2d0/Rm
352 ELSE
353 Ay = -Bm*Rm*(LOG(r/Rm)+half)
354 END IF
355 END IF
356 END IF
357
358 IF (simtype == colliding_opp .AND. y > yupper-dx) Ay = -Ay
359 A(i,j,k,2) = Ay
360 END DO
361 END DO
362 END DO
363
364 !Then take the curl of A to get B
365 IF (nDim == 2) THEN
366 Info%aux(:,:,:,:) = 0d0
367 ELSE IF (nDim == 3) THEN
368 DO j=1-rmbc, my+rmbc
369 y = ylower + (REAL(j,xPrec) - half)*dx
370 IF (simtype /= single) THEN
371 IF (y >= 0d0+yjet .AND. y <= yupper-yjet) CYCLE
372 ELSE
373 IF (y >= 0d0+yjet) CYCLE
374 END IF
375 DO i=1-rmbc, mx+rmbc+1
376 DO k=1-zrmbc, mz+zrmbc
377 Info%aux(i,j,k,1) = (A(i,j+1,k,3)-A(i,j,k,3))/dx - (A(i,j,k+1,2)-A(i,j,k,2))/dx
378 END DO
379 END DO
380 END DO
381
382 DO j=1-rmbc, my+rmbc+1
383 y = ylower + (REAL(j,xPrec) - 1d0)*dx
384 IF (simtype /= single) THEN
385 IF (y >= 0d0+yjet .AND. y <= yupper+yjet) CYCLE
386 ELSE
387 IF (y >= 0d0+yjet) CYCLE
388 END IF
389 DO i=1-rmbc, mx+rmbc
390 DO k=1-zrmbc, mz+zrmbc
391 Info%aux(i,j,k,2) = (A(i,j,k+1,1)-A(i,j,k,1))/dx - (A(i+1,j,k,3)-A(i,j,k,3))/dx
392 END DO
393 END DO
394 END DO
395 END IF
396 DO j=1-rmbc, my+rmbc
397 y = ylower + (REAL(j,xPrec) - half)*dx
398 IF (simtype /= single) THEN
399 IF (y >= 0d0+yjet .AND. y <= yupper-yjet) CYCLE
400 ELSE
401 IF (y >= 0d0+yjet) CYCLE
402 END IF
403 DO i=1-rmbc, mx+rmbc
404 DO k=1-zrmbc, mz+zrmbc+1
405 IF (nDim == 2) THEN
406 BzAux(i,j) = (A(i+1,j,1,2)-A(i,j,1,2))/dx - (A(i,j+1,1,1)-A(i,j,1,1))/dx
407 ELSE IF (nDim == 3) THEN
408 Info%aux(i,j,k,3) = (A(i+1,j,k,2)-A(i,j,k,2))/dx - (A(i,j+1,k,1)-A(i,j,k,1))/dx
409 END IF
410 END DO
411 END DO
412 END DO
413 DEALLOCATE(A)
414
415 ! Average aux fields to get B-fields for q, update energy
416 DO j=1-rmbc, my+rmbc
417 y = ylower + (REAL(j,xPrec) - half)*dx
418 IF (simtype /= single) THEN
419 IF (y >= 0d0+yjet .AND. y <= yupper-yjet) CYCLE
420 ELSE
421 IF (y >= 0d0+yjet) CYCLE
422 END IF
423 DO i=1-rmbc, mx+rmbc
424 DO k=1-zrmbc, mz+zrmbc
425 IF (nDim == 2) THEN
426 q(i,j,k,iBx:iBy) = 0d0
427 q(i,j,k,iBz) = BzAux(i,j)
428 ELSE IF (nDim == 3) THEN
429 q(i,j,k,iBx) = half*(Info%aux(i,j,k,1)+Info%aux(i+1,j,k,1))
430 q(i,j,k,iBy) = half*(Info%aux(i,j,k,2)+Info%aux(i,j+1,k,2))
431 q(i,j,k,iBz) = half*(Info%aux(i,j,k,3)+Info%aux(i,j,k+1,3))
432 END IF
433! q(i,j,k,iE) = q(i,j,k,iE) + half*SUM(q(i,j,k,iBx:iBz)**2d0)
434
435 ! Pressure Check
436 ! ptot = (gamma1*(q(i,j,k,iE) - half*SUM(q(i,j,k,2:m_high)**2d0)/q(i,j,k,1)) + &
437 ! gamma13*half*SUM(q(i,j,k,iBx:iBz)**2d0))*pScale
438 ! IF (ptot < pamb) THEN
439 ! padd = pamb - ptot
440 ! q(i,j,k,iE) = q(i,j,k,iE) + padd/gamma1/pScale
441 ! END IF
442 END DO
443 END DO
444 END DO
445 IF (nDim == 2) DEALLOCATE(BzAux)
446 END IF
447
448 END SUBROUTINE PlaceJet
449
450 !> Could be used to update grids pre-output
451 !! @param Info Info Object
452 SUBROUTINE ProblemAfterStep(Info)
453 TYPE(InfoDef) :: Info
454 END SUBROUTINE ProblemAfterStep
455
456 !> Could be used to set force refinement
457 !! @param Info Info object
458 SUBROUTINE ProblemSetErrFlag(Info)
459 TYPE(InfoDef) :: Info
460 END SUBROUTINE ProblemSetErrFlag
461
462 SUBROUTINE ProblemBeforeGlobalStep(n)
463 INTEGER :: n
464 END SUBROUTINE ProblemBeforeGlobalStep
465
466
467END MODULE Problem
468