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
|
---|
38 | MODULE 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 |
|
---|
71 | CONTAINS
|
---|
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
|
---|
572 | END MODULE SourceControl
|
---|
573 |
|
---|