Blog: meeting notes: source_control.f90

File source_control.f90, 22.3 KB (added by Erica Kaminski, 6 years ago)
Line 
1!#########################################################################
2!
3! Copyright (C) 2003-2012 Department of Physics and Astronomy,
4! University of Rochester,
5! Rochester, NY
6!
7! source_control.f90 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 source
24!! @brief directory containing modules for handling source terms
25
26!> @defgroup Source Source Terms
27!! @brief Group containing modules for handling source terms
28
29!> @file source_control.f90
30!! @brief Main file for module SourceControl
31
32!> @defgroup SourceControl Source Control
33!! @brief Module for managing various source terms
34!! @ingroup Source
35
36!> Module for managing various source terms
37!! @ingroup SourceControl
38MODULE SourceControl
39
40 USE DataDeclarations
41 USE PhysicsDeclarations
42 USE EOS
43 USE CoolingSrc
44 USE CylindricalSrc
45 USE RadForceSrc
46 USE UniformGravitySrc
47 USE PointGravitySrc
48 USE SelfGravitySrc
49 USE RotatingSrc
50 USE Chemistry
51 USE Outflows
52 USE ParticleDeclarations
53
54
55 IMPLICIT NONE
56 PRIVATE
57 PUBLIC Src, SrcInit, SrcFinalizeInit, ReadSourceObjects, SourceCell !CreateSrc, CoolingDef, CreateCoolingObject
58 SAVE
59 TYPE(LevelDef) :: currlevel
60
61 ! sources verbosity shorthand
62 INTEGER :: vb=0
63
64 !> Generic interface for creating sources
65 INTERFACE CreateSrc
66 MODULE PROCEDURE CreatePointGravityObject!, CreateCylindricalObject, &
67 !, CreateUniformGravityObject
68 END INTERFACE
69
70
71CONTAINS
72
73 !> Read in and initialize parameters for source terms
74 SUBROUTINE SrcInit(isrcsolvetype,iverbosity)
75 INTEGER,OPTIONAL :: isrcsolvetype,iverbosity
76 ! allocate sources & zero out components
77 CALL PointGravityInit()
78 IF (lRestart) CALL ReadSourceObjects(restart_frame)
79 END SUBROUTINE SrcInit
80
81 !> Finalize source term(s) initialization
82 !! @details All source terms should be set up by user at this point, and
83 !! all hyperbolic, elliptic, etc., variables defined. Now,
84 !! set up any additional tables, etc.
85 !! [BDS][20110106]: Do not initialize cooling source objects on restarts--use the I/O routines for that.
86 SUBROUTINE SrcFinalizeInit
87 END SUBROUTINE SrcFinalizeInit
88
89 !> Calculate source terms and update q.
90 !! @details Take in an Info structure, range, and hydro time step to
91 !! integrate source term(s).
92 !! @param Info Info structure
93 !! @param mb Restricted range to calculate source on grid
94 !! @param hdt Timestep from the hydro calculation
95 SUBROUTINE Src(Info, mb, tnow, hdt,divv_opt)
96 TYPE(InfoDef) :: Info
97 INTEGER,INTENT(IN) :: mb(3,2)
98 REAL(KIND=qPrec),INTENT(IN) :: hdt, tnow
99 REAL(KIND=qPREC), DIMENSION(:,:,:), OPTIONAL :: divv_opt
100 CALL Protect(Info,mb, 'before source')
101 IF (PRESENT(divv_opt)) THEN
102 CALL ExplicitSrc(Info, mb, tnow, hdt, divv_opt)
103 ELSE
104 CALL ExplicitSrc(Info, mb, tnow, hdt)
105 END IF
106 CALL Protect(Info,mb, 'after source')
107 END SUBROUTINE Src
108
109
110 ! ==================================================================
111 ! = Explicit Scheme Section =
112 ! ==================================================================
113
114
115 !> Use subcycling algorithm to calculate source terms.
116 !! @details When there are no strong source terms present,
117 !! algorithm uses 2nd order RK to quickly update q.
118 !! When stiff source terms present, it switches to
119 !! 4th/5th order RK.
120 !! @param Info Info structure
121 !! @param mb Restricted range to calculate source on grid
122 !! @param hdt Timestep from the hydro calculation
123 SUBROUTINE ExplicitSrc(Info, mb, tnow, hdt, divv_opt)
124 ! Interface declarations
125 TYPE(InfoDef) :: Info
126 INTEGER,INTENT(IN) :: mb(3,2)
127 REAL(KIND=qPrec),INTENT(IN) :: hdt
128 REAL(KIND=qPREC), DIMENSION(:,:,:), OPTIONAL :: divv_opt
129 REAL(KIND=qPREC) :: divv
130 ! Internal declarations
131 REAL(KIND=qPrec) :: tnow, dv
132 INTEGER :: i,j,k,ip(3)
133 REAL(KIND=qPrec) :: volumefactor
134 LOGICAL :: ghost, ghostz, ghostyz
135 REAL(KIND=qPREC), DIMENSION(:), ALLOCATABLE :: q
136 dv = Levels(Info%level)%dx**ndim
137 ALLOCATE(q(NrHydroVars))
138 DO k=mb(3,1),mb(3,2)
139 Ghostz = .NOT. (k >= 1 .AND. k <= Info%mx(3))
140 DO j=mb(2,1),mb(2,2)
141 Ghostyz = Ghostz .OR. .NOT. (j >= 1 .AND. j <= Info%mx(2))
142 DO i=mb(1,1),mb(1,2)
143 Ghost = Ghostyz .OR. .NOT. (i >= 1 .AND. i <= Info%mx(1))
144 IF (Ghost) THEN
145 VolumeFactor = 0d0
146 ELSE
147 IF (Info%Level==MaxLevel) THEN
148 VolumeFactor = 1d0
149 ELSE
150 IF (Info%ChildMask(i,j,k) >= 1) THEN
151 VolumeFactor = 0d0
152 ELSE
153 VolumeFactor = 1d0
154 END IF
155 END IF
156 END IF
157
158 ip=(/i,j,k/)
159
160 ! Method:
161 ! 1) Try 2nd order Runge-Kutta to get estimate of error for timestep.
162 ! Should speed things up where source terms are weak.
163 ! 2) If error is above tolerance, reduce timestep and switch to
164 ! 5th/4th order RK to do the integration.
165 ! 3) Iterate at 5th/4th until full timestep reached.
166
167 q=Info%q(i,j,k,1:NrHydroVars)
168 IF (PRESENT(divv_opt)) THEN
169 divv=divv_opt(i+1-mb(1,1), j+1-mb(2,1), k+1-mb(3,1))
170 ELSE
171 divv=0
172 END IF
173 CALL Cons_To_Source(q)
174 CALL SourceCell(q, Info, ip, tnow, hdt, dv*VolumeFactor,divv)!where it gets complicated
175 CALL Source_To_Cons(q)
176
177 IF(ANY(q(:).ne.q(:) .OR. abs(q(:)).gt.huge(abs(q(:)))) .AND. .NOT. lRequestRestart) THEN
178 PRINT*, 'Src routine produced a NAN'
179 CALL PrintQ(Info, q, tnow, i, j, k)
180 lRequestRestart=.true.
181 ! STOP
182 DEALLOCATE(q)
183 RETURN
184 END IF
185 Info%q(i,j,k,1:NrHydroVars)=q
186
187 END DO
188 END DO
189 END DO
190 DEALLOCATE(q)
191 !IF(vb>0) THEN
192
193 !END IF
194 END SUBROUTINE ExplicitSrc
195
196 SUBROUTINE SourceCell(q, Info, ip, tnow, hdt, dv, divv)
197 ! Interface declarations
198 TYPE(InfoDef) :: Info
199 REAL(KIND=qPrec),INTENT(IN) :: hdt
200 ! Internal declarations
201 REAL(KIND=qPrec) :: q(:), neutrinocoolingrate
202 REAL(KIND=qPrec) :: pos(3)
203 REAL(KIND=qPrec) :: dq(NrHydroVars),tnow,tnext, dv, divv
204 INTEGER :: ip(3)
205 REAL(KIND=qPrec) :: error
206 LOGICAL :: success
207
208
209! NeutrinoCoolingRate = GetCoolingStrength(q)
210! IF (q(1)>0.0001) THEN
211! IF (MPI_ID==0) PRINT *, 'iE_init=', q(iE), 'hdt=', hdt
212! IF (MPI_ID==0) PRINT *, 'NeutrinoCoolingRate (de/dt)=', NeutrinoCoolingRate, '*hdt=', NeutrinoCoolingRate*hdt
213! IF (MPI_ID==0) PRINT *, 'Predicted change in iE over timestep=', q(iE)-(NeutrinoCoolingRate*hdt)
214! END IF
215
216 tnext=tnow+hdt
217 pos(1:nDim)=Info%xBounds(1:nDim,1)+(REAL(ip(1:nDim))-half)*levels(Info%level)%dx
218 IF (nDim == 2) pos(3)=Info%xBounds(3,1)
219 IF (iCylindrical == NOCYL) CALL PointGravity(q,hdt,pos,tnow,dv,info%Level)
220 ! CALL PointGravity(q,hdt,pos,tnow,dv,info%Level)
221 ! 2nd order RK, uses hdt to get full timestep error
222 CALL RKOrder2(q,dq,pos,tnow,hdt,error,Info,ip, divv)
223 IF(vb>1) PRINT*,'::ExplicitSrc:RKOrder2 -- dq,error',dq,error
224
225 IF(Error > SrcPrecision) THEN
226 IF(vb>1) PRINT*,'::ExplicitSrc:RKOrder45, qbefore, tnow, tnext',q,tnow,tnext
227 ! iterate with 5th/4th order RK to get error & adjust timestep
228 ! until full timestep is completed
229 ! The subroutine updates q
230 CALL RKOrder45(q,pos,hdt,tnow,tnext,Info,ip,divv,success)
231 IF (.NOT. success .AND. .NOT. lRequestRestart) THEN
232 write(*,*) 'RK Failed'
233 CALL PrintQ(Info, q, tnow, ip(1), ip(2), ip(3))
234 lRequestRestart=.true.
235 ! CALL OutputDoubleArray(Info%q(i-1:i+1,j-1:j+1,k-1,iPhi))
236 ! CALL OutputDoubleArray(Info%q(i-1:i+1,j-1:j+1,k,iPhi))
237 ! CALL OutputDoubleArray(Info%q(i-1:i+1,j-1:j+1,k+1,iPhi))
238
239 ! CALL OutputDoubleArray(Info%q(i-1:i+1,j-1:j+1,k-1,iPhiDot))
240 ! CALL OutputDoubleArray(Info%q(i-1:i+1,j-1:j+1,k,iPhiDot))
241 ! CALL OutputDoubleArray(Info%q(i-1:i+1,j-1:j+1,k+1,iPhiDot))
242 ! STOP
243 RETURN
244 END IF
245 IF(vb>1) PRINT*,'::ExplicitSrc:RKOrder45, qafter, tnow, tnext',q,tnow,tnext
246
247 ELSE
248 ! good enough at 2nd order, update q
249 q=q+dq
250 END IF
251
252! IF (q(1)>0.0001) THEN
253! IF (MPI_ID==0) PRINT *, 'Actual change in q(iE) over timestep=', q(iE)
254! END IF
255 END SUBROUTINE SourceCell
256
257 !> 2nd order Runge-Kutta scheme
258 !! @params q Fluid variables for this cell
259 !! @params dq Change in q from source terms
260 !! @params pos Location of cell
261 !! @params dt Timestep to integrate over
262 SUBROUTINE RKOrder2(q,dq,pos,t,dt,error,Info,ip,divv)
263 ! Interface declarations
264 TYPE(InfoDef) :: Info
265 INTEGER :: ip(3)
266 REAL(KIND=qPrec) :: q(:)
267 REAL(KIND=qPrec) :: dq(NrHydroVars), pos(3), dt, error,t
268 ! Internal declarations
269 REAL(KIND=qPrec), ALLOCATABLE :: qstar(:)
270 REAL(KIND=qPrec) :: k1(NrHydroVars),k2(NrHydroVars), qerr(NrHydroVars), minscales(NrHydroVars)
271 REAL(KIND=qPREC) :: divv
272 dq=0d0
273
274 ALLOCATE(qStar(NrHydroVars))
275 ! Runge-Kutta 2nd order: for a function y(t), going from t=t0:t1 with dt=t1-t0,
276 ! the second order integration is
277 ! k1 = d(y0)/dt
278 ! y* = y0 + dt k1
279 ! k2 = d(y*)/dt
280 ! y1 = y0 + dt/2 (k1+k2)
281 ! where in our case y(t) -> q(t) and dq/dt is calculated from the source terms.
282 ! The error is
283 ! error = |y*/y1|
284
285 k1 = SrcDerivs(q,pos,t,Info,ip,divv)
286 WHERE(k1/=k1)k1=0d0
287 qstar = q + dt*k1
288 k2 = SrcDerivs(qstar,pos,t+dt,Info,ip,divv)
289 WHERE(k2/=k2)k2=0d0
290
291 dq=5d-1*dt*(k1+k2)
292
293 IF(ALL(ABS(dq(:))<TINY(dq(:)))) THEN
294 dq=0d0
295 error=0d0
296 ELSE
297 minscales(1)=minDensity
298 IF (lMHD) minscales(iBx:iBz) = max(sqrt(max(minDensity, q(1)))*max(SoundSpeed(q), sqrt(MinTemp/TempScale)),maxval(abs(q(iBx:iBz))))
299 IF (nTracerLo /= 0) minscales(nTracerLo:nTracerHi) = minDensity
300 minscales(m_low:m_high)=PrimSoundSpeed(q)
301 IF (iE .ne. 0) minscales(iE)=max(minDensity, q(1))*max(SoundSpeed(q)**2,MinTemp/TempScale)
302 qerr=q+dq-qstar
303 error=MAXVAL(ABS(qErr(:))/max(abs(q(:)), minscales(:)))
304 ! error = MAXVAL(ABS(qstar(:))/ABS(q(:)+dq(:))) !5d-1*dt*(k1+k2)))
305 END IF
306
307
308
309
310
311 DEALLOCATE(qStar)
312 END SUBROUTINE RKOrder2
313
314
315 !> Combination 5th/4th order Runge-Kutta scheme
316 !! @details Algorithm adjusts timestep to subcycle until full timestep is integrated over.
317 !! Method is based on the Cash-Karp RK algorithm in Numerical Recipes (1992), 710ff.
318 !! It takes advantage of the fact that six function calls results in formulas
319 !! for both 4th and 5th order integration, allowing quick ability to measure error
320 !! and adjust timestep as appropriate.
321 !! Once the error is acceptable, it uses the 5th order result to update q and continue.
322 !! @params q Fluid variables for this cell
323 !! @params pos Location of cell
324 !! @params hdt Timestep to integrate over
325 !! @params htnow Time at beginning of timestep
326 !! @params htnext Time at end of timestep
327 SUBROUTINE RKOrder45(q,pos,hdt,htnow,htnext,Info,ip,divv,success)
328 ! Interface declarations
329 TYPE(InfoDef) :: Info
330 INTEGER :: ip(3)
331 REAL(KIND=qPrec) :: q(:)
332 REAL(KIND=qPrec) :: pos(3),hdt,htnow,htnext
333 REAL(KIND=qPREC) :: divv
334 ! Internal declarations
335 REAL(KIND=qPrec) :: dt,tnow,tnext
336 INTEGER,PARAMETER :: MaxIters=1d4
337 INTEGER :: niter
338 REAL(KIND=qPrec),ALLOCATABLE :: qStar(:)
339 ! CKRK stuff, see above reference for table
340 REAL(KIND=qPrec) :: k1(NrHydroVars),k2(NrHydroVars),k3(NrHydroVars),k4(NrHydroVars),k5(NrHydroVars),k6(NrHydroVars), qTemp(NrHydroVars), qErr(NrHydroVars), maxerror, minscales(NrHydroVars)
341 REAL(KIND=qPrec),PARAMETER :: A1=0.0, A2=0.2, A3=0.3, A4=0.6, A5=1.0, A6=0.875 &
342 , B21=0.2 &
343 , B31=3./40., B32=9./40. &
344 , B41=0.3, B42=-0.9, B43=1.2 &
345 , B51=-11./54., B52=2.5, B53=-70./27., B54=35./27. &
346 , B61=1631./55296., B62=175./512., B63=575./13824., B64=44275./110592., B65=253./4096. &
347 , C1=37./378., C3=250./621., C4=125./594., C6=512./1771. &
348 , DC1=C1-2825./27648., DC3=C3-18575./48384., DC4=C4-13525./55296., DC5=-277./14336., DC6=C6-0.25
349 REAL(KIND=qPrec),PARAMETER :: SAFETY=0.9, PSHRNK=-0.25, PGROW=-0.2
350 LOGICAL :: success
351 !
352 ALLOCATE(qStar(NrHydroVars))
353 dt=hdt
354 tnow=htnow
355 tnext=htnext
356
357 k1=0d0;k2=0d0;k3=0d0;k4=0d0;k5=0d0;k6=0d0
358
359 DO niter=1,MaxIters
360 !
361 ! Cash-Karp RK functions, see above reference for details
362 !
363 k1 = SrcDerivs(q,pos,tnow,Info,ip, divv)
364 WHERE(k1/=k1)k1=0
365 qStar = q + dt*( B21*k1 )
366 k2 = SrcDerivs(qStar,pos,tnow+dt*A2,Info,ip, divv)
367 WHERE(k2/=k2)k2=0
368 qStar = q + dt*( B31*k1 + B32*k2 )
369 k3 = SrcDerivs(qStar,pos,tnow+dt*(A3),Info,ip, divv)
370 WHERE(k3/=k3)k3=0
371 qStar = q + dt*( B41*k1 + B42*k2 + B43*k3 )
372 k4 = SrcDerivs(qStar,pos,tnow+dt*(A4),Info,ip, divv)
373 WHERE(k4/=k4)k4=0
374 qStar = q + dt*( B51*k1 + B52*k2 + B53*k3 + B54*k4 )
375 k5 = SrcDerivs(qStar,pos,tnow+dt*(A5),Info,ip, divv)
376 WHERE(k5/=k5)k5=0
377 qStar = q + dt*( B61*k1 + B62*k2 + B63*k3 + B64*k4 + B65*k5 )
378 k6 = SrcDerivs(qStar,pos,tnow+dt*(A6),Info,ip, divv)
379 WHERE(k6/=k6)k6=0
380 qTemp = q + dt*( C1*k1 + C3*k3 + C4*k4 + C6*k6 )
381 qErr = dt*( DC1*k1 + DC3*k3 + DC4*k4 + DC5*k5 + DC6*k6 )
382
383 !
384 ! Error evaluation & timestep adjustment
385 !
386! IF(vb>1) WRITE(*,'(A,10E20.12)')'qtemp ',qtemp(1)
387! IF(vb>1) WRITE(*,'(A,10E20.12)')'qerr ',qerr(1)
388 IF(vb>1) PRINT*,'qtemp ',qtemp
389 IF(vb>1) PRINT*,'qerr ',qerr
390
391 minscales(1)=minDensity
392 IF (lMHD) minscales(iBx:iBz) = max(sqrt(max(minDensity, q(1)))*max(SoundSpeed(q), sqrt(MinTemp/TempScale)),maxval(abs(q(iBx:iBz))))
393 IF (nTracerLo /= 0) minscales(nTracerLo:nTracerHi) = minDensity
394 minscales(m_low:m_high)=PrimSoundSpeed(q)
395 IF (iE .ne. 0) minscales(iE)=max(minDensity, q(1))*max(SoundSpeed(q)**2,MinTemp/TempScale)
396 maxerror=MAXVAL(ABS(qErr(:))/max(abs(q(:)), minscales(:)))
397 IF (vb>1) print*,maxerror
398 IF(maxerror > SrcPrecision) THEN
399 ! reduce timestep and retry
400 dt = MAX( SAFETY*dt*(maxerror/SrcPrecision)**PSHRNK, dt/MaxIters)
401 IF(dt <= TINY(dt) .AND. .NOT. lRequestRestart) THEN
402 PRINT*,'Src Error::RKOrder45 -- timestep insignificant. Halting.'
403 ! PRINT*,' Iterations/MaxIters: ',niter, MaxIters
404 ! PRINT*,' Last timestep, full timestep: ',dt, hdt
405 ! PRINT*,' Current time, hydro begin/end time: ',tnow,htnow,htnext
406 ! PRINT*,' Time left to go: ',htnow-tnow
407 CALL PrintQ(Info, q, htnow, ip(1),ip(2),ip(3))
408 lRequestRestart=.true.
409 success=.false.
410 EXIT
411 ! STOP
412 ! RETURN
413
414 END IF
415 ! tnow remains the same; update tnext
416 tnext=tnow+dt
417 ELSE
418 q = qTemp
419 tnow=tnow+dt
420 dt = MAX( MIN( SAFETY*dt*(maxerror/SrcPrecision)**PGROW, 5*dt, htnext-tnow), 0d0)
421 tnext=tnow+dt
422 IF(dt <= TINY(dt)) THEN
423 ! end of full timestep reached
424 success=.true.
425 EXIT
426 END IF
427 END IF
428 END DO
429
430 IF(niter>MaxIters .AND. .NOT. lRequestRestart) THEN
431 PRINT*,'Src Error::RKOrder45 -- maximum number of iterations reached.'
432 ! PRINT*,' Max iterations: ',MaxIters
433 ! PRINT*,' Last timestep, full timestep: ',dt, hdt
434 ! PRINT*,' Current time, hydro begin/end time: ',tnow,htnow,htnext
435 ! PRINT*,' Time left to go: ',htnow-tnow
436 CALL PrintQ(Info, q, htnow, ip(1), ip(2), ip(3))
437 ! STOP
438 lRequestRestart=.true.
439 success=.false.
440 END IF
441 DEALLOCATE(qStar)
442 END SUBROUTINE RKOrder45
443
444
445
446 ! ==================================================================
447 ! = ~~~ Wrapper Subroutines/Functions ~~~ =
448 ! ==================================================================
449
450 !> Changes to q (dqdt) for all source terms
451 !! @params q Fluid variables for this cell
452 !! @params pos Location of cell
453 FUNCTION SrcDerivs(q,pos,t,Info,ip, divv)
454 ! Interface declarations
455 TYPE(InfoDef) :: Info
456 INTEGER :: ip(3)
457 REAL(KIND=qPrec) :: q(:), ne, Temp, r
458 REAL(KIND=qPrec) :: pos(3), SrcDerivs(NrHydroVars), t, particlepos(3), offsetpos(3)
459 REAL(KIND=qPREC) :: divv
460 ! Internal declarations
461 REAL(KIND=qPrec) :: dqdt(NrHydroVars)
462 TYPE(particleListDef), Pointer :: particlelist
463 TYPE(particleDef), Pointer :: particle
464 dqdt=0d0
465 IF (iCooling /= 0 .OR. lChemistry) THEN
466 Temp=SourceTemperature(q)
467 IF (lTrackHydrogen .OR. lTrackHelium) ne=get_ne(q)
468 IF (lChemistry) CALL ChemReactions(q, dqdt, ne, Temp)
469 IF (Temp > FloorTemp) THEN
470 IF (iCooling /= 0) CALL Cooling(q, dqdt, ne, Temp, divv, pos)
471 end if
472 END IF
473
474 IF (iCylindrical.ge.NoCyl .and. iCylindrical.le.WithAngMom) THEN
475 IF (iCylindrical.ne.NoCyl) CALL Cylindrical(q,dqdt,pos)!,currlevel%dx)
476 ELSE
477 print*,'';print*,'iCylindrical= ',iCylindrical
478 print*,'but should be either 0, 1 or 2. Aborting.'
479 stop
480 END IF
481
482 IF (luniformgravity) CALL UniformGravity_src(q,dqdt,t)
483 IF (OmegaRot /= 0d0) CALL Rotating(q, dqdt, pos)
484 IF (iCylindrical /= NOCYL) CALL PointGravity_inst(q,dqdt,pos,t) !moved to source cell for improved mom. cons.
485 IF (CentralLuminosity /= 0) CALL RadForce(q,dqdt,pos,t)
486
487 SrcDerivs=dqdt
488 END FUNCTION SrcDerivs
489
490
491
492 !> Reads source-term objects from a Chombo file. Currently, only cooling is implemented.
493 !! @param nframe The number of the restart frame.
494 SUBROUTINE ReadSourceObjects(nframe)
495
496#if defined _HDF5
497 USE ChomboDeclarations, ONLY: ChomboHandle, CreateChomboHandle, CloseChomboHandle, CHOMBO_HANDLE_READ, &
498 Chombo_OpenSourceGroup, Chombo_CloseSourceGroup
499
500 USE PointGravitySrc, ONLY: PointGravityDef, CreatePointGravityObject, PointGravity_ReadObjectFromChombo
501 USE Outflows, ONLY: OutflowDef, CreateOutflow, Outflow_ReadObjectFromChombo
502#endif
503
504 INTEGER :: nframe
505
506#if defined _HDF5
507 ! Variable declarations
508 TYPE(ChomboHandle), POINTER :: chandle
509 TYPE(PointGravityDef), POINTER :: gravity_obj
510 TYPE(OutflowDef), POINTER :: outflow_obj
511 INTEGER :: nr_pointgravityObj, nr_outflowobj
512 CHARACTER(LEN=23) :: s_filename
513
514
515 ! Open a reading handle for the specified frame.
516 WRITE(s_filename, '(A10,I5.5,A4)') 'out/chombo', nframe, '.hdf'
517 CALL CreateChomboHandle(s_filename, chandle, CHOMBO_HANDLE_READ)
518
519
520!!! Begin Check For Point Gravity Terms !!!
521 ! Open the 'point_gravity_objects' group and save the handle.
522 nr_pointgravityObj = Chombo_OpenSourceGroup(chandle, "point_gravity_objects")
523 chandle%source_offset = 0
524
525! IF (MPI_ID == 0) write(*,*) "Reading in point graity objects", nr_pointgravityObj
526 ! Read source objects from the data file until the expected number of objects is read.
527 DO WHILE (chandle%source_offset < nr_pointgravityObj)
528
529 ! Create a new gravity object.
530 NULLIFY(gravity_obj)
531 CALL CreatePointGravityObject(gravity_obj)
532
533 ! Read in the point gravity data from the Chombo file. This subroutine also advances the
534 ! chandle%source_offset variable.
535 CALL PointGravity_ReadObjectFromChombo(chandle, gravity_obj)
536! IF (MPI_ID == 0) write(*,*) 'read in gravity_obj', gravity_obj%id
537 END DO
538 ! Close the source term group. This also clears the source term offset.
539 CALL Chombo_CloseSourceGroup(chandle)
540
541 nr_outflowObj = Chombo_OpenSourceGroup(chandle, "outflow_objects")
542 chandle%source_offset = 0
543
544! IF (MPI_ID == 0) write(*,*) "Reading in outflow objects", nr_outflowObj
545 ! Read source objects from the data file until the expected number of objects is read.
546 DO WHILE (chandle%source_offset < nr_outflowobj)
547
548 ! Create a new gravity object.
549 NULLIFY(outflow_obj)
550 CALL CreateOutflow(outflow_obj)
551
552 ! Read in the point gravity data from the Chombo file. This subroutine also advances the
553 ! chandle%source_offset variable.
554 CALL Outflow_ReadObjectFromChombo(chandle, outflow_obj)
555! IF (MPI_ID == 0) write(*,*) "read in outflow_obj with id", outflow_obj%objid
556 CALL UpdateOutflow(outflow_obj)
557
558 END DO
559 ! Close the source term group. This also clears the source term offset.
560 CALL Chombo_CloseSourceGroup(chandle)
561
562 CALL CloseChomboHandle(chandle)
563
564#else
565 IF (MPI_ID == 0) THEN
566 WRITE(*,*) "Unable to Read Source Objects without chombo support"
567 STOP
568 END IF
569#endif
570
571 END SUBROUTINE ReadSourceObjects
572END MODULE SourceControl
573