sedov.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# sedov - blast wave solutions #
4!# module: sedov.f90 #
5!# #
6!# Copyright (C) 2019 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# #
9!# This program 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 2 of the License, or (at #
12!# your option) any later version. #
13!# #
14!# This program is distributed in the hope that it will be useful, but #
15!# WITHOUT ANY WARRANTY; without even the implied warranty of #
16!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
17!# NON INFRINGEMENT. See the GNU General Public License for more #
18!# details. #
19!# #
20!# You should have received a copy of the GNU General Public License #
21!# along with this program; if not, write to the Free Software #
22!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
23!# #
24!#############################################################################
25
26!----------------------------------------------------------------------------!
31!----------------------------------------------------------------------------!
32MODULE sedov_mod
33 IMPLICIT NONE
34 !--------------------------------------------------------------------------!
35! #define DEBUG_OUTPUT 1
36#undef DEBUG_OUTPUT
37#ifdef f2003
38 use, INTRINSIC :: iso_fortran_env, only : stdin=>input_unit, &
39 stdout=>output_unit, &
40 stderr=>error_unit
41#else
42#define STDIN 5
43#define STDOUT 6
44#define STDERR 0
45#endif
46 !--------------------------------------------------------------------------!
47 PRIVATE
48 REAL, PARAMETER :: pi = 3.1415926535897932384626433832795028842
49 !--------------------------------------------------------------------------!
50 TYPE sedov_typ
51 LOGICAL :: initialized = .false.
52 INTEGER :: ndim = 3
53 REAL :: gamma = 1.4
54 REAL :: omega = 0.0
55 REAL :: vmin = 2./7. ! = 2./(ndim+2-omega)/gamma
56 REAL :: vmax = 1./3. ! = 4./(ndim+2-omega)/(gamma+1.)
57 REAL :: rho0 = 1.0
58 REAL :: e0 = 1.0
59 REAL :: p0 = 1.0e-5
60 REAL :: alpha = 0.8511
61 REAL :: eps = 1.0e-12 ! tolerance for integration and root finding
62 ! smaller values yield better results, but
63 ! may run into convergence problems
64 REAL, ALLOCATABLE, DIMENSION(:) :: plist
65 CONTAINS
66 PROCEDURE :: initparams
67 PROCEDURE :: computeenergyintegral
68 PROCEDURE :: printconfiguration
69 PROCEDURE :: computesolution
70 PROCEDURE :: shockposition
71 PROCEDURE :: abort
72 PROCEDURE :: destroy
73 END TYPE sedov_typ
74 !--------------------------------------------------------------------------!
75 INTERFACE sedov_typ
76 MODULE PROCEDURE createsedov
77 END INTERFACE
78 !--------------------------------------------------------------------------!
79 PUBLIC :: sedov_typ
80
81CONTAINS
82
84 FUNCTION createsedov(ndim,gamma,omega,rho0,E0,p0) RESULT(this)
85 IMPLICIT NONE
86 !-------------------------------------------------------------------!
87 TYPE(sedov_typ) :: this
88 INTEGER, OPTIONAL, INTENT(IN) :: ndim
89 REAL, OPTIONAL, INTENT(IN) :: gamma,omega,rho0,e0,p0
90 !-------------------------------------------------------------------!
91 ! perform sanity checks and set basic parameters
92 IF (PRESENT(ndim)) THEN
93 IF (ndim.LT.1.AND.ndim.GT.3) &
94 CALL this%Abort("CreateSedov: ndim must be one of {1,2,3}")
95 this%ndim = ndim
96 END IF
97 IF (PRESENT(gamma)) THEN
98 IF (gamma.LE.1.0.OR.gamma.GT.5./3.) &
99 CALL this%Abort("CreateSedov: gamma not in ( 1 , 5/3 )")
100 this%gamma = gamma
101 END IF
102 IF (PRESENT(omega)) THEN
103 IF (omega.LT.0.0.OR.omega.GE.this%ndim) &
104 CALL this%Abort("CreateSedov: omega not in [ 0 , ndim )")
105 this%omega = omega
106 ! default precision does not work for omega > 0
107 IF (this%omega.GT.0.0) this%eps = 1.0e-10
108 END IF
109 IF (PRESENT(rho0)) THEN
110 IF (rho0.LE.tiny(rho0).OR.rho0.GT.huge(rho0)) &
111 CALL this%Abort("CreateSedov: rho0 not in ( 0 , huge )")
112 this%rho0 = rho0
113 END IF
114 IF (PRESENT(e0)) THEN
115 IF (e0.LE.tiny(e0).OR.e0.GT.huge(e0)) &
116 CALL this%Abort("CreateSedov: E0 not in ( 0 , huge )")
117 this%E0 = e0
118 END IF
119 IF (PRESENT(p0)) THEN
120 IF (p0.LE.tiny(p0).OR.p0.GT.huge(p0)) &
121 CALL this%Abort("CreateSedov: p0 not in ( 0 , huge )")
122 this%p0 = p0
123 END IF
124 CALL this%InitParams()
125 this%initialized = .true.
126 END FUNCTION createsedov
127
131 SUBROUTINE initparams(this)
132 IMPLICIT NONE
133 !-------------------------------------------------------------------!
134 CLASS(sedov_typ), INTENT(INOUT) :: this
135 !-------------------------------------------------------------------!
136 INTEGER :: err
137 REAL :: gam,gam_p1,gam_m1,j,w,w2,w3,J1,J2
138#ifdef DEBUG_OUTPUT
139 INTEGER :: i,inum
140 REAL :: V,dV,x1,x2,x3,x4,f,g,h,lambda,dlambda
141#endif
142 !-------------------------------------------------------------------!
143 ! some frequently used constants
144 gam = this%gamma
145 gam_p1 = gam + 1.0d+00
146 gam_m1 = gam - 1.0d+00
147 j = real(this%ndim)
148 w = this%omega
149 w2 = (2*gam_m1 + j)/gam
150 w3 = j*(2-gam)
151 this%Vmin = 2./(j+2-w)/gam
152 this%Vmax = 4./(j+2-w)/gam_p1
153
154 IF (ALLOCATED(this%plist)) DEALLOCATE(this%plist)
155 ALLOCATE(this%plist(19),stat=err)
156 IF (err.NE.0) CALL this%Abort("InitParams: memory allocation failed")
157
158 this%plist(1) = gam ! = gamma
159 this%plist(2) = 2.0/(j+2.0-w) ! = alpha_0
160 this%plist(4) = gam_m1/(gam*(w-w2)) ! = alpha_2
161 this%plist(3) = (j+2.0-w)*gam/(2.0+j*gam_m1) & ! = alpha_1
162 *(2*(j*(2.0-gam)-w)/(gam*(j+2.0-w)**2)-this%plist(4))
163 this%plist(5) = (j-w)/(gam*(w2-w)) ! = alpha_3
164 this%plist(6) = (j+2-w)*(j-w)/(w3-w)*this%plist(3) ! = alpha_4
165 this%plist(7) = (w*gam_p1-2*j)/(w3-w) ! = alpha_5
166 this%plist(8) = 0.25*(j+2-w)*gam_p1 ! = a = dx1/dV
167 this%plist(9) = gam_p1/gam_m1 ! = b
168 this%plist(10) = 0.5*gam*(j+2.0-w) ! = c = 1/Vmin
169 this%plist(12) = 0.5*(2+j*gam_m1) ! = e
170 this%plist(11) = 1.0 & ! = d = a/(a-e) = 1/(1-e/a)
171 /(1.0-this%plist(12)/this%plist(8))
172 this%plist(13) = this%plist(9)*this%plist(10) ! = dx2/dV
173 this%plist(14) = -this%plist(11)*this%plist(12) ! = dx3/dV
174 this%plist(15) = -this%plist(9)*this%plist(10) & ! = dx4/dV
175 /this%plist(1)
176 this%plist(16) = j ! = j = ndim : {1,2,3}
177 this%plist(17) = w ! = omega
178 ! these are set later, see ComputeSolution
179 this%plist(18) = 1.0
180 this%plist(19) = 1.0
181
182 this%Vmin = this%Vmin + 2*epsilon(this%Vmin) ! avoid lambda(Vmin) -> infty
183
184 CALL this%ComputeEnergyIntegral(j1,j2)
185#ifdef DEBUG_OUTPUT
186 print '(A)',"# " // repeat("=",77)
187 print '(A)',"# Parameter:"
188 print '(A)',"# gamma omega dim alpha J1 J2"
189 print '(A,6(F11.4))',"# ",this%gamma,this%omega,real(this%ndim),this%alpha,j1,j2
190 print '(A)',"# " // repeat("-",77)
191 print '(A,10(F10.4))',"# p1-p7 ",this%plist(1:7)
192 print '(A,10(F10.4))',"# p8-p13 ",this%plist(8:13)
193 print '(A,10(F10.4))',"# p14-p19",this%plist(14:19)
194 print '(A)',"# " // repeat("=",77)
195 inum = 20
196 dv = (this%Vmax-this%Vmin)/inum
197 print '(A)',"# lambda V f g h dJ1/dlam dJ2/dlam"
198 DO i=0,inum
199 v = this%Vmin+i*dv
200 CALL funcd_lambda(v,this%plist,x1,x2,x3,lambda,dlambda)
201 x4 = this%plist(9)*(1.0-this%plist(10)*v/this%plist(1))
202 f = lambda*x1
203 g = func_g(x1,x2,x3,x4,this%plist)
204 h = func_h(x1,x2,x3,x4,this%plist)
205 print '(A,7(F11.4))','# ',lambda,v,f,g,h,func_j1(v,this%plist),func_j2(v,this%plist)
206 END DO
207#endif
208 END SUBROUTINE initparams
209
211 SUBROUTINE computeenergyintegral(this,J1,J2)
212 USE integration
213 IMPLICIT NONE
214 !-------------------------------------------------------------------!
215 CLASS(sedov_typ), INTENT(INOUT) :: this
216 !-------------------------------------------------------------------!
217 REAL, OPTIONAL, INTENT(OUT) :: J1,J2
218 !-------------------------------------------------------------------!
219 ! compute energy integral
220 j1 = integrate(func_j1,this%Vmin,this%Vmax,this%eps,this%plist,1)
221 j2 = integrate(func_j2,this%Vmin,this%Vmax,this%eps,this%plist,1)
222 IF (this%ndim.EQ.1) THEN
223 ! 1D
224 this%alpha = 0.5*(j1+1./(this%gamma-1.0)*j2)
225 ELSE ! 2D and 3D
226 this%alpha = (this%ndim-1)*pi*(j1+2./(this%gamma-1.0)*j2)
227 END IF
228 END SUBROUTINE computeenergyintegral
229
231 SUBROUTINE computesolution(this,t,r,rho,v,p)
232 USE roots
233 IMPLICIT NONE
234 !-------------------------------------------------------------------!
235 CLASS(sedov_typ), INTENT(INOUT) :: this
236 REAL, INTENT(IN) :: t,r
237 REAL, INTENT(OUT) :: rho,v,p
238 !-------------------------------------------------------------------!
239 INTEGER :: error
240 REAL :: Rshock,Vshock,Vstar,rho2,v2,p2,x1,x2,x3,x4,lambda,dlambda
241 REAL :: lnx,lnlambda,lnx1,lnx2,lnx3,lnx4
242 !-------------------------------------------------------------------!
243 ! compute current shock position
244 rshock = this%ShockPosition(t)
245 IF (r.GT.0.0.AND.r.LT.rshock) THEN
246 ! estimate for Vstar valid for 0< Vstar/Vmin - 1 << 1
247 lnx = -(log(this%plist(9)) + (log(r/rshock) + this%plist(2) &
248 *log(this%plist(8)/this%plist(10)) + this%plist(3) &
249 *log(this%plist(11)*(1.0-this%plist(12)/this%plist(10))))/this%plist(4))
250 vstar = (1.0+exp(lnx)) / this%plist(10)
251 IF (lnx.LT.log(8*epsilon(lnx))) THEN
252 ! Vstar = Vmin*(1.0+EXP(lnx)) with EXP(lnx) << epsilon
253 ! => computation of Vstar and other quantities depending on Vstar fails
254 ! => use logarithmic quantities and approximate for Vstar ~ Vmin
255 vshock = this%plist(2) * rshock / t
256 rho2 = this%plist(9) * this%rho0
257 v2 = 2./(this%gamma+1.0) * vshock
258 p2 = this%rho0 * v2 * vshock
259 lnlambda = log(r/rshock)
260 lnx1 = log(this%plist(8)*this%Vmin) ! = LOG(a*Vmin)
261 lnx2 = log(this%plist(9)) + lnx ! = LOG(b*x)
262 lnx3 = log(this%plist(11)*(1.0-this%plist(12)/this%plist(10))) ! = LOG(d*(1-e/c))
263 lnx4 = log(this%plist(9)*(1.-this%plist(10)*this%Vmin/this%gamma)) ! = LOG(b*(1-c*Vmin/gamma))
264 rho = rho2 * exp(func_lng(lnx1,lnx2,lnx3,lnx4,this%plist))
265 v = v2 * exp(lnx1+lnlambda)
266 p = p2 * exp(func_lnh(lnx1,lnx2,lnx3,lnx4,this%plist))
267 ELSE
268 ! Newton-Raphson to solve the implicit equation
269 this%plist(18) = rshock
270 this%plist(19) = r
271 ! Newton-Raphson to solve the implicit equation for Vstar
272 ! see Kamm & Timmes, 2007 eq. (68)
273 CALL getroot_newton(funcd,this%Vmin,this%Vmax,vstar,error,this%plist)
274! CALL GetRoot(func,this%Vmin,this%Vmax,Vstar,error,this%plist)
275 CALL funcd_lambda(vstar,this%plist,x1,x2,x3,lambda,dlambda)
276 x4 = this%plist(9)*(1.0-this%plist(10)*vstar/this%gamma)
277 vshock = this%plist(2) * rshock / t
278 rho2 = this%plist(9) * this%rho0
279 v2 = 2./(this%gamma+1.0) * vshock
280 p2 = this%rho0 * v2 * vshock
281 rho = rho2 * func_g(x1,x2,x3,x4,this%plist)
282 v = v2 * lambda * x1
283 p = p2 * func_h(x1,x2,x3,x4,this%plist)
284 END IF
285 ELSE
286 ! return unperturbed state
287 rho = this%rho0/r**this%omega
288 v = 0.0
289 p = this%p0
290 END IF
291 END SUBROUTINE computesolution
292
293 FUNCTION shockposition(this,time) RESULT(Rshock)
294 IMPLICIT NONE
295 !-------------------------------------------------------------------!
296 CLASS(sedov_typ), INTENT(INOUT) :: this
297 REAL, INTENT(IN) :: time
298 REAL :: rshock
299 !-------------------------------------------------------------------!
300 IF (.NOT.this%initialized) CALL this%Abort("sedov uninitialzed, aborting ...")
301 rshock = (this%E0/this%rho0*time*time/this%alpha)**(0.5*this%plist(2))
302 END FUNCTION shockposition
303
305 SUBROUTINE printconfiguration(this)
306 IMPLICIT NONE
307 !-------------------------------------------------------------------!
308 CLASS(sedov_typ), INTENT(INOUT) :: this
309 !-------------------------------------------------------------------!
310 WRITE (stdout,"(A)") "# ==========================================="
311 WRITE (stdout,"(A)") "# << Solution of Sedov explosion problem >> "
312 WRITE (stdout,"(A)") "# ==========================================="
313 WRITE (stdout,"(A,I3)") "# Geometry / Dimensions ", this%ndim
314 WRITE (stdout,"(A,F10.6)") "# Ratio of specific heat cap. ", this%gamma
315 WRITE (stdout,"(A,ES10.2)") "# Initial energy input ", this%E0
316 WRITE (stdout,"(A,ES10.2)") "# Ambient density ", this%rho0
317 WRITE (stdout,"(A,ES10.2)") "# Peak density ", this%rho0*this%plist(9)
318 WRITE (stdout,"(A,F10.6)") "# Density power law exponent ", this%omega
319 WRITE (stdout,"(A,F10.6)") "# Shock position parameter ", this%alpha
320 END SUBROUTINE printconfiguration
321
322 SUBROUTINE abort(this,msg)
323 IMPLICIT NONE
324 !-------------------------------------------------------------------!
325 CLASS(sedov_typ), INTENT(INOUT) :: this
326 !-------------------------------------------------------------------!
327 CHARACTER(LEN=*) :: msg
328 !-------------------------------------------------------------------!
329 WRITE (stderr,'(A)') "ERROR in " // trim(msg)
330 stop
331 END SUBROUTINE abort
332
333 SUBROUTINE destroy(this)
334 IMPLICIT NONE
335 !------------------------------------------------------------------------!
336 CLASS(sedov_typ), INTENT(INOUT) :: this
337 !-------------------------------------------------------------------!
338 DEALLOCATE(this%plist)
339 this%initialized = .false.
340 END SUBROUTINE destroy
341
342 PURE SUBROUTINE funcd(y,fy,dfy,plist)
343 IMPLICIT NONE
344 !------------------------------------------------------------------------!
345 REAL, INTENT(IN) :: y
346 REAL, INTENT(OUT) :: fy,dfy
347 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
348 !------------------------------------------------------------------------!
349 REAL :: x1,x2,x3,lambda,dlambda
350 !------------------------------------------------------------------------!
351 CALL funcd_lambda(y,plist,x1,x2,x3,lambda,dlambda)
352 fy = plist(18)*lambda - plist(19)
353 dfy = plist(18)*dlambda
354 END SUBROUTINE funcd
355
356 PURE SUBROUTINE func(y,fy,plist)
357 IMPLICIT NONE
358 !------------------------------------------------------------------------!
359 REAL, INTENT(IN) :: y
360 REAL, INTENT(OUT) :: fy
361 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
362 !------------------------------------------------------------------------!
363 REAL :: x1,x2,x3,lambda,dlambda
364 !------------------------------------------------------------------------!
365 CALL funcd_lambda(y,plist,x1,x2,x3,lambda,dlambda)
366 fy = plist(18)*lambda - plist(19)
367 END SUBROUTINE func
368
369 FUNCTION func_j1(x,plist) RESULT(fx)
370 IMPLICIT NONE
371 !------------------------------------------------------------------------!
372 REAL, INTENT(IN) :: x
373 REAL, INTENT(INOUT), DIMENSION(:), OPTIONAL :: plist
374 REAL :: fx
375 !------------------------------------------------------------------------!
376 REAL :: x1,x2,x3,x4,lambda,dlambda,g
377 !------------------------------------------------------------------------!
378 CALL funcd_lambda(x,plist,x1,x2,x3,lambda,dlambda)
379 x4 = plist(9)*(1.0-plist(10)*x/plist(1))
380 g = func_g(x1,x2,x3,x4,plist)
381 fx = plist(9) * lambda**(plist(16)+1.0) * g * x**2 * dlambda
382 END FUNCTION func_j1
383
384 FUNCTION func_j2(x,plist) RESULT(fx)
385 IMPLICIT NONE
386 !------------------------------------------------------------------------!
387 REAL, INTENT(IN) :: x
388 REAL, INTENT(INOUT), DIMENSION(:), OPTIONAL :: plist
389 REAL :: fx
390 !------------------------------------------------------------------------!
391 REAL :: x1,x2,x3,x4,lambda,dlambda,h
392 !------------------------------------------------------------------------!
393 CALL funcd_lambda(x,plist,x1,x2,x3,lambda,dlambda)
394 x4 = plist(9)*(1.0-plist(10)*x/plist(1))
395 h = func_h(x1,x2,x3,x4,plist)
396 !!! ATTENTION: There is an error in eq. 56 of Kamm & Timmes (2007);
397 !!! it must be lambda^(gamma-1) and NOT ..^(gamma+1)
398 fx = 0.5*(plist(1)+1.0)/plist(8)**2 * lambda**(plist(16)-1.0) * h * dlambda
399 END FUNCTION func_j2
400
401 PURE SUBROUTINE funcd_lambda(V,plist,x1,x2,x3,lambda,dlambda)
402 IMPLICIT NONE
403 !------------------------------------------------------------------------!
404 REAL, INTENT(IN) :: v
405 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
406 REAL, INTENT(OUT) :: x1,x2,x3,lambda,dlambda
407 !------------------------------------------------------------------------!
408 x1 = plist(8)*v
409 x2 = plist(9)*max(0.0,plist(10)*v-1.0) ! should be >=0, since V*p10 >= 1
410 x3 = plist(11)*(1.0-plist(12)*v)
411 lambda = x1**(-plist(2)) * x2**(-plist(4)) * x3**(-plist(3))
412! dlambda = -lambda*(plist(2)/x1*plist(8) + plist(4)/x2*plist(13) + plist(3)/x3*plist(14))
413 dlambda = -plist(2)*plist(8) * x1**(-plist(2)-1.0) * x2**(-plist(4)) * x3**(-plist(3)) &
414 -plist(4)*plist(13) * x1**(-plist(2)) * x2**(-plist(4)-1.0) * x3**(-plist(3)) &
415 -plist(3)*plist(14) * x1**(-plist(2)) * x2**(-plist(4)) * x3**(-plist(3)-1.0)
416 END SUBROUTINE funcd_lambda
417
418 PURE FUNCTION func_g(x1,x2,x3,x4,plist) RESULT(g)
419 IMPLICIT NONE
420 !------------------------------------------------------------------------!
421 REAL, INTENT(IN) :: x1,x2,x3,x4
422 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
423 REAL :: g
424 !------------------------------------------------------------------------!
425 g = x1**(plist(2)*plist(17)) * x2**(plist(5)+plist(4)*plist(17)) &
426 *x3**(plist(6)+plist(3)*plist(17)) * x4**plist(7)
427 END FUNCTION func_g
428
429 PURE FUNCTION func_lng(lnx1,lnx2,lnx3,lnx4,plist) RESULT(lng)
430 IMPLICIT NONE
431 !------------------------------------------------------------------------!
432 REAL, INTENT(IN) :: lnx1,lnx2,lnx3,lnx4
433 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
434 REAL :: lng
435 !------------------------------------------------------------------------!
436 lng = (plist(2)*plist(17))*lnx1 + (plist(5)+plist(4)*plist(17))*lnx2 &
437 + (plist(6)+plist(3)*plist(17))*lnx3 + plist(7)*lnx4
438 END FUNCTION func_lng
439
440 PURE FUNCTION func_h(x1,x2,x3,x4,plist) RESULT(h)
441 IMPLICIT NONE
442 !------------------------------------------------------------------------!
443 REAL, INTENT(IN) :: x1,x2,x3,x4
444 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
445 REAL :: h
446 !------------------------------------------------------------------------!
447 h = x1**(plist(2)*plist(16)) * x3**(plist(6)+plist(3)*(plist(17)-2.0)) &
448 *x4**(1.0+plist(7))
449 END FUNCTION func_h
450
451 PURE FUNCTION func_lnh(lnx1,lnx2,lnx3,lnx4,plist) RESULT(lnh)
452 IMPLICIT NONE
453 !------------------------------------------------------------------------!
454 REAL, INTENT(IN) :: lnx1,lnx2,lnx3,lnx4
455 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
456 REAL :: lnh
457 !------------------------------------------------------------------------!
458 lnh = (plist(2)*plist(16))*lnx1 + (plist(6)+plist(3)*(plist(17)-2.0))*lnx3 &
459 + (1.0+plist(7))*lnx4
460 END FUNCTION func_lnh
461
462END MODULE sedov_mod
463
464
465
funcd
pure subroutine funcd(y, fy, dfy, plist)
Definition: bondi2d.f90:416
integration
Numerical integration.
Definition: integration.f90:33
integration::integrate
real function, public integrate(fkt, xl, xr, eps, plist, method)
Numerical integration function.
Definition: integration.f90:125
roots
root finding subroutines
Definition: roots.f90:45
roots::getroot_newton
subroutine, public getroot_newton(funcd, x1, x2, root, error, plist, xm, iterations)
Definition: roots.f90:247
sedov_mod
module providing sedov class with basic init and solve subroutines
Definition: sedov.f90:32
sedov_mod::pi
real, parameter pi
Definition: sedov.f90:48
sedov_mod::initparams
subroutine initparams(this)
sets constants according to Kamm & Timmes "On Efficient Generation of Numerically Robust Sedov Soluti...
Definition: sedov.f90:132
sedov_mod::func_h
pure real function func_h(x1, x2, x3, x4, plist)
Definition: sedov.f90:441
sedov_mod::func_lng
pure real function func_lng(lnx1, lnx2, lnx3, lnx4, plist)
Definition: sedov.f90:430
sedov_mod::abort
subroutine abort(this, msg)
Definition: sedov.f90:323
sedov_mod::destroy
subroutine destroy(this)
Definition: sedov.f90:334
sedov_mod::shockposition
real function shockposition(this, time)
Definition: sedov.f90:294
sedov_mod::func_j1
real function func_j1(x, plist)
Definition: sedov.f90:370
sedov_mod::func_lnh
pure real function func_lnh(lnx1, lnx2, lnx3, lnx4, plist)
Definition: sedov.f90:452
sedov_mod::printconfiguration
subroutine printconfiguration(this)
prints configuration on the screen (stdout)
Definition: sedov.f90:306
sedov_mod::func_g
pure real function func_g(x1, x2, x3, x4, plist)
Definition: sedov.f90:419
sedov_mod::computesolution
subroutine computesolution(this, t, r, rho, v, p)
compute the energy integral and set parameter alpha
Definition: sedov.f90:232
sedov_mod::createsedov
type(sedov_typ) function createsedov(ndim, gamma, omega, rho0, E0, p0)
constructor for sedov class
Definition: sedov.f90:85
sedov_mod::computeenergyintegral
subroutine computeenergyintegral(this, J1, J2)
compute the energy integral and set parameter alpha
Definition: sedov.f90:212
sedov_mod::funcd_lambda
pure subroutine funcd_lambda(V, plist, x1, x2, x3, lambda, dlambda)
Definition: sedov.f90:402
sedov_mod::func_j2
real function func_j2(x, plist)
Definition: sedov.f90:385
func
pure subroutine func(x, fx, plist)
Definition: rootstest.f90:192
omega
pure real function, dimension(size(ephir%data2d, dim=1)) omega(ephir, velocity)
Definition: sources_viscosity.f90:432
sedov_mod::sedov_typ
Definition: sedov.f90:50
Generated by   doxygen 1.9.4