| 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 |
|
|---|
