ls_timmes_eos_module.f
Go to the documentation of this file.
1 #include "../macros.h"
2  module ls_timmes_eos_module
3 
4  use units_module
5 
6  implicit none
7 
8 !.....LS & Timmes switch parameters.....................................
9  real(r_kind),parameter :: tref=5.00e9 ![K] temperature LS-Timmes switch
10  ! 5.80e9 K approx. 0.5 MeV
11  real(r_kind),parameter :: dref=1.00e7 ![g/cm3] density LS-Timmes switch
12 
13 !.....EoS input variable parameters.....................................
14  integer,parameter :: ink=0
15  integer,parameter :: ine=1
16  integer,parameter :: ins=2
17 
18 !.....nr & bisection parameters.........................................
19  real(r_kind),parameter :: t1=1.00e2 !lower temperature bracket [K]
20  real(r_kind),parameter :: t2=1.00e12 !upper temperature bracket [K]
21  ! 1.17e12 K approx. 100 MeV
22 
23  type timmes_eos_state
24  integer :: ef !eos flag
25  real(r_kind) :: t !Temperature [K]
26  real(r_kind) :: f !Helmholtz free energy [erg/g]
27  real(r_kind) :: s !Entropy [kB/baryon]
28  real(r_kind) :: p !Pressure [erg/cm3]
29  real(r_kind) :: d2fdtdd !d^2F/dTdrho [erg*cm3/(g2*K)]
30  real(r_kind) :: d2fdt2 !d^2F/dT^2 [erg/(g*K2)]
31  real(r_kind) :: cs !Sound speed [cm/s]
32  real(r_kind) :: e !Internal energy [erg/g]
33  real(r_kind) :: abar !Average mass number [-]
34  real(r_kind) :: abarxx !Average mass number [-]
35  real(r_kind) :: zbar !Average charge number [-]
36  real(r_kind) :: xa !Mass fraction of alpha particles [-]
37  real(r_kind) :: cv !Specific heat capacity at constant volume
38  real(r_kind) :: xp !Mass fraction of exterior protons [-]
39  real(r_kind) :: xn !Mass fraction of exterior neutrons [-]
40  real(r_kind) :: muhat !Neutron chemical potential minus
41  ! proton chemical potential [MeV]
42  real(r_kind) :: mun !Neutron chemical potential [MeV]
43  real(r_kind) :: mue !Electron chemical potential [MeV]
44  end type
45 
46  contains
47 
48  subroutine timmes_eos(var,vin,d,Ye,state,status)
49 
50  implicit none
51 
52  integer,intent(in) :: var
53  real(r_kind),intent(in) :: vin,d,Ye
54  type(timmes_eos_state),intent(inout) :: state
55  integer,intent(out) :: status
56 
57 !-----------------------------------------------------------------------
58 !
59 ! LS & Timmes EoS routine for temperature, int. energy or entropy
60 !
61 ! Input: var ... selected input variable (see definitions
62 ! up in module)
63 ! vin=T,e,s ... Temperature [K], internal energy [erg/g],
64 ! entropy [kB/baryon]
65 ! d ... Density [g/cm3]
66 ! Ye ... Electron abundance [-]
67 !
68 ! Output: state ... see timmes_eos_state type declaration
69 ! status ... 0 for zero problems and 1 for trouble
70 !
71 !-----------------------------------------------------------------------
72 
73  real(r_kind) :: T
74 
75 !.....call relevant routine.............................................
76  if ( var .eq. ink ) then
77  call timmes_eos_interface(vin,d,ye,state,status)
78  elseif ( var .eq. ine ) then
79  !find corresponding temperature
80  call timmes_eos_e_bisec(vin,d,ye,state%abar,t,status)
81  if ( status .eq. 0 ) call timmes_eos_interface(t,d,ye,state,
82  1 status)
83  elseif ( var .eq. ins ) then
84  !find corresponding temperature
85  call timmes_eos_s_bisec(vin,d,ye,state%abar,t,status)
86  if ( status .eq. 0 ) call timmes_eos_interface(t,d,ye,state,
87  1 status)
88  else
89  if (verbose_level .gt. 1) then
90  write(6,*) 'timmes_eos_module: Check input variable '
91  1 ,'setting...'
92  endif
93  status = 1
94  endif
95 
96  return
97 
98  end subroutine timmes_eos
99 
100  subroutine timmes_eos_interface(T,d,Ye,state,status)
101 
102  use timmes_eos_module
103 
104  implicit none
105 
106  real(r_kind),intent(in) :: T,d,Ye
107  type(timmes_eos_state),intent(inout) :: state
108  integer,intent(out) :: status
109 
110 !-----------------------------------------------------------------------
111 !
112 ! Timmes EoS interface
113 !
114 ! NOTE: state%abar is used as an input!!!
115 !
116 ! Input: T ... Temperature [K]
117 ! d ... Density [g/cm3]
118 ! Ye ... Electron abundance [-]
119 ! state%abar ... average mass number [-]
120 !
121 ! Output: state ... see timmes_eos_state type declaration
122 ! status ... 0 for zero problems and 1 for trouble
123 !
124 !-----------------------------------------------------------------------
125 
126  include 'vector_eos.dek'
127 
128 !.....set input for timmes eos..........................................
129  temp_row(1) = t
130  den_row(1) = d
131  abar_row(1) = state%abar
132  zbar_row(1) = state%abar*ye
133 
134 !.....call timmes eos...................................................
135  !here the pipeline is only 1 element long
136  jlo_eos = 1
137  jhi_eos = 1
138 
139  call eosfxt
140 
141 !.....set output variables..............................................
142  state%T = t ![K]
143  state%e = etot_row(1) ![erg/g]
144  state%S = stot_row(1)/(units%kB/units%mb) ![kB/baryon]
145  state%F = etot_row(1) - t*stot_row(1) ![erg/g]
146  state%p = ptot_row(1) ![erg/cm3]
147  state%d2FdTdd = dpt_row(1)/(d*d) ![erg*cm3/(g2*K)]
148  state%d2FdT2 = -dst_row(1) ![erg/(g*K2)]
149  state%cs = cs_row(1) ![cm/s]
150  state%zbar = state%abar*ye ![-]
151  state%xa = 0. ![-]
152  state%xp = 0. ![-]
153  state%xn = 0. ![-]
154  state%cv = cv_row(1)/units%Na/units%MeV ![Mev/(baryon*K)]
155  state%muhat = 0. ![MeV]
156  state%mun = 0. ![MeV]
157  state%mue = etaele_row(1) ![MeV]
158 ! state%mue = 0. ![MeV]
159 
160 !debug:
161 ! write(6,*) 'pion = ',pion_row(1)
162 ! write(6,*) 'prad = ',prad_row(1)
163 ! write(6,*) 'pele = ',pele_row(1) + ppos_row(1)
164 ! write(6,*) 'pcou = ',pcou_row(1)
165 ! write(6,*) 'ptot = ',ptot_row(1)
166 
167 !.....set status according EoS..........................................
168  if ( .not.eosfail ) then
169  status = 0
170  else
171  status = 1
172  endif
173 
174  end subroutine timmes_eos_interface
175 
176  subroutine timmes_eos_e_bisec(e,d,Ye,abar,T,status)
177 
178  implicit none
179 
180  real(r_kind),intent(in) :: e,d,Ye,abar
181  real(r_kind),intent(out) :: T
182  integer,intent(out) :: status
183 
184 !-----------------------------------------------------------------------
185 !
186 ! Find Temperature for given internal energy of LS-Timmes EoS
187 ! with bisection algorithm ... see Numerical Recipes
188 !
189 ! Input: e ... Internal energy [erg/g]
190 ! d ... Density [g/cm3]
191 ! Ye ... Electron abundance [-]
192 !
193 ! Output: T ... Temperature [K]
194 ! status ... 0 for zero problems and 1 for trouble
195 !
196 !-----------------------------------------------------------------------
197 
198  integer,parameter :: jmax=100 !max number of bisec iterations
199  real(r_kind),parameter :: precision=5.e-4 !eps*(T1+T2)/2
200  integer :: j
201  type(timmes_eos_state) :: state
202  real(r_kind) :: Tmid,dT,f,fmid
203 
204 !.....start bisection...................................................
205  state%abar = abar
206  call timmes_eos_interface(t1,d,ye,state,status)
207  if ( status .eq. 0 ) then
208  f = e - state%e
209  else
210  if (verbose_level .gt. 1) then
211  write(6,*) 'ls_timmes_eos_module: Check bracket T1 =',t1
212  1 ,'in eos_find_T.'
213  end if
214  status = 1
215  return
216  endif
217  call timmes_eos_interface(t2,d,ye,state,status)
218  if ( status .eq. 0 ) then
219  fmid = e - state%e
220  else
221  if (verbose_level .gt. 1) then
222  write(6,*) 'ls_timmes_eos_module: Check bracket T2 =',t2
223  1 ,'in eos_find_T.'
224  end if
225  status = 1
226  return
227  endif
228  if( f*fmid .ge. 0. ) then
229  if (verbose_level .gt. 1) then
230  write(6,*) 'ls_timmes_eos_module: Root not bracketed... '
231  1 ,'try other values of T1,T2 ...'
232  end if
233  status = 1
234  return
235  endif
236  if ( f .lt. 0. ) then
237  t = t1
238  dt = t2 - t1
239  else
240  t = t2
241  dt = t1 - t2
242  endif
243  do j=1,jmax
244  dt = 0.5*dt
245  tmid = t + dt
246  call timmes_eos_interface(tmid,d,ye,state,status)
247  if ( status .eq. 0 ) then
248  fmid = e - state%e
249  else
250  if (verbose_level .gt. 1) then
251  write(6,*) 'ls_timmes_eos_module: Error while bisecting...'
252  end if
253  status = 1
254  return
255  endif
256  if ( fmid .le. 0. ) t = tmid
257  if ( (abs(dt) .lt. precision) .or. (fmid .eq. 0.) ) then
258  status = 0
259  return
260  endif
261  enddo
262 
263 !.....too many iterations...............................................
264  status = 1
265  if (verbose_level .gt. 1) then
266  write(6,*) 'ls_timmes_eos_module: Too many iterations...'
267  end if
268 
269  end subroutine timmes_eos_e_bisec
270 
271  subroutine timmes_eos_e_nr(e,d,Ye,abar,T,status)
272 
273  implicit none
274 
275  real(r_kind),intent(in) :: e,d,Ye,abar
276  real(r_kind),intent(inout) :: T
277  integer,intent(out) :: status
278 
279 !-----------------------------------------------------------------------
280 !
281 ! Find Temperature for given internal energy of LS-Timmes EoS
282 ! with Newton-Raphson algorithm ... see Numerical Recipes
283 !
284 ! Input: e ... Internal energy [erg/g]
285 ! d ... Density [g/cm3]
286 ! Ye ... Electron abundance [-]
287 !
288 ! Output: T ... Temperature [K]
289 ! status ... 0 for zero problems and 1 for trouble
290 !
291 !-----------------------------------------------------------------------
292 
293  integer,parameter :: jmax=100 !max number of nr iterations
294  real(r_kind),parameter :: precision=5.e-4 !eps*(T1+T2)/2
295  integer :: j
296  type(timmes_eos_state) :: state
297  real(r_kind) :: f,df
298 
299 
300 !.....check the eventually given initial guess for T....................
301  state%abar = abar
302  call timmes_eos_interface(t,d,ye,state,status)
303  if ( status .eq. 1 ) then
304  t = t2
305  call timmes_eos_interface(t,d,ye,state,status)
306  endif
307  f = state%e - e
308  df = -t*state%d2FdT2
309 
310 !.....nr iterations.....................................................
311  do j=1,jmax
312  t = t - f/df
313  if ( t .le. 0. ) t = t2
314  call timmes_eos_interface(t,d,ye,state,status)
315  if ( status .ne. 0 ) then
316  if (verbose_level .gt. 1) then
317  write(6,*) 'ls_timmes_eos_module: NR iteration failed '
318  1 ,'with T =',t
319  end if
320  status = 1
321  return
322  endif
323  f = state%e - e
324  df = -t*state%d2FdT2
325  if ( (f .eq. 0.) .or. (abs(f/df) .le. precision) ) then
326  status = 0
327  return
328  endif
329  enddo
330 
331 !.....too many iterations...............................................
332  status = 1
333  if (verbose_level .gt. 1) then
334  write(6,*) 'ls_timmes_eos_module: Too many iterations...'
335  end if
336  end subroutine timmes_eos_e_nr
337 
338  subroutine timmes_eos_s_bisec(s,d,Ye,abar,T,status)
339 
340  implicit none
341 
342  real(r_kind),intent(in) :: s,d,Ye,abar
343  real(r_kind),intent(out) :: T
344  integer,intent(out) :: status
345 
346 !-----------------------------------------------------------------------
347 !
348 ! Find Temperature for given entropy of LS-Timmes EoS
349 ! with bisection algorithm ... see Numerical Recipes
350 !
351 ! Input: s ... Entropy [kB/baryon]
352 ! d ... Density [g/cm3]
353 ! Ye ... Electron abundance [-]
354 !
355 ! Output: T ... Temperature [K]
356 ! status ... 0 for zero problems and 1 for trouble
357 !
358 !-----------------------------------------------------------------------
359 
360  integer,parameter :: jmax=100 !max number of bisec iterations
361  real(r_kind),parameter :: precision=5.e-4 !eps*(T1+T2)/2
362  integer :: j
363  type(timmes_eos_state) :: state
364  real(r_kind) :: Tmid,dT,f,fmid
365 
366 !.....start bisection...................................................
367  state%abar = abar
368  call timmes_eos_interface(t1,d,ye,state,status)
369  if ( status .eq. 0 ) then
370  f = s - state%s
371  else
372  if (verbose_level .gt. 1) then
373  write(6,*) 'ls_timmes_eos_module: Check bracket T1 =',t1
374  1 ,'in eos_find_T.'
375  end if
376  status = 1
377  return
378  endif
379  call timmes_eos_interface(t2,d,ye,state,status)
380  if ( status .eq. 0 ) then
381  fmid = s - state%s
382  else
383  if (verbose_level .gt. 1) then
384  write(6,*) 'ls_timmes_eos_module: Check bracket T2 =',t2
385  1 ,'in eos_find_T.'
386  end if
387  status = 1
388  return
389  endif
390  if( f*fmid .ge. 0. ) then
391  if (verbose_level .gt. 1) then
392  write(6,*) 'ls_timmes_eos_module: Root not bracketed... '
393  1 ,'try other values of T1,T2 ...'
394  end if
395  status = 1
396  return
397  endif
398  if ( f .lt. 0. ) then
399  t = t1
400  dt = t2 - t1
401  else
402  t = t2
403  dt = t1 - t2
404  endif
405  do j=1,jmax
406  dt = 0.5*dt
407  tmid = t + dt
408  call timmes_eos_interface(tmid,d,ye,state,status)
409  if ( status .eq. 0 ) then
410  fmid = s - state%s
411  else
412  if (verbose_level .gt. 1) then
413  write(6,*) 'ls_timmes_eos_module: Error while bisecting...'
414  end if
415  status = 1
416  return
417  endif
418  if ( fmid .le. 0. ) t = tmid
419  if ( (abs(dt) .lt. precision) .or. (fmid .eq. 0.) ) then
420  status = 0
421  return
422  endif
423  enddo
424 
425 !.....too many iterations...............................................
426  status = 1
427  if (verbose_level .gt. 1) then
428  write(6,*) 'ls_timmes_eos_module: Too many iterations...'
429  end if
430  end subroutine timmes_eos_s_bisec
431 
432  subroutine timmes_eos_s_nr(s,d,Ye,abar,T,status)
433 
434  implicit none
435 
436  real(r_kind),intent(in) :: s,d,Ye,abar
437  real(r_kind),intent(inout) :: T
438  integer,intent(out) :: status
439 
440 !-----------------------------------------------------------------------
441 !
442 ! Find Temperature for given entropy of LS-Timmes EoS
443 ! with Newton-Raphson algorithm ... see Numerical Recipes
444 !
445 ! Input: s ... Entropy [kB/baryon]
446 ! d ... Density [g/cm3]
447 ! Ye ... Electron abundance [-]
448 !
449 ! Output: T ... Temperature [K]
450 ! status ... 0 for zero problems and 1 for trouble
451 !
452 !-----------------------------------------------------------------------
453 
454  integer,parameter :: jmax=100 !max number of nr iterations
455  real(r_kind),parameter :: precision=5.e-4 !eps*(T1+T2)/2
456  integer :: j
457  type(timmes_eos_state) :: state
458  real(r_kind) :: f,df
459 
460 
461 !.....check the eventually given initial guess for T....................
462  state%abar = abar
463  call timmes_eos_interface(t,d,ye,state,status)
464  if ( status .eq. 1 ) then
465  t = t2
466  call timmes_eos_interface(t,d,ye,state,status)
467  endif
468  f = (state%s - s)*(units%kB/units%mb)
469  df = -state%d2FdT2
470 
471 !.....nr iterations.....................................................
472  do j=1,jmax
473  t = t - f/df
474  if ( t .le. 0. ) t = t2
475  call timmes_eos_interface(t,d,ye,state,status)
476  if ( status .ne. 0 ) then
477  if (verbose_level .gt. 1) then
478  write(6,*) 'ls_timmes_eos_module: NR iteration failed '
479  1 ,'with T =',t
480  end if
481  status = 1
482  return
483  endif
484  f = (state%s - s)*(units%kB/units%mb)
485  df = -state%d2FdT2
486  if ( (f .eq. 0.) .or. (abs(f/df) .le. precision) ) then
487  status = 0
488  return
489  endif
490  enddo
491 
492 !.....too many iterations...............................................
493  status = 1
494  if (verbose_level .gt. 1) then
495  write(6,*) 'ls_timmes_eos_module: Too many iterations...'
496  end if
497  end subroutine timmes_eos_s_nr
498 
499  end module ls_timmes_eos_module
ls_timmes_eos_module::ink
integer, parameter ink
Definition: ls_timmes_eos_module.f:14
ls_timmes_eos_module
Definition: ls_timmes_eos_module.f:2
ls_timmes_eos_module::tref
real(r_kind), parameter tref
Definition: ls_timmes_eos_module.f:9
ls_timmes_eos_module::timmes_eos_s_bisec
subroutine timmes_eos_s_bisec(s, d, Ye, abar, T, status)
Definition: ls_timmes_eos_module.f:339
units_module
Definition: units_module.f:7
ls_timmes_eos_module::t2
real(r_kind), parameter t2
Definition: ls_timmes_eos_module.f:20
ls_timmes_eos_module::timmes_eos
subroutine timmes_eos(var, vin, d, Ye, state, status)
Definition: ls_timmes_eos_module.f:49
ls_timmes_eos_module::t1
real(r_kind), parameter t1
Definition: ls_timmes_eos_module.f:19
units_module::units
type(units_in_cgs), parameter units
Definition: units_module.f:16
ls_timmes_eos_module::ine
integer, parameter ine
Definition: ls_timmes_eos_module.f:15
ls_timmes_eos_module::timmes_eos_interface
subroutine timmes_eos_interface(T, d, Ye, state, status)
Definition: ls_timmes_eos_module.f:101
timmes_eos_module
Definition: timmes_eos_module.f:1
eosfxt
subroutine eosfxt
Definition: eosfxt.f:13
r_kind
#define r_kind
Definition: macros.h:46
ls_timmes_eos_module::timmes_eos_e_nr
subroutine timmes_eos_e_nr(e, d, Ye, abar, T, status)
Definition: ls_timmes_eos_module.f:272
ls_timmes_eos_module::ins
integer, parameter ins
Definition: ls_timmes_eos_module.f:16
ls_timmes_eos_module::timmes_eos_e_bisec
subroutine timmes_eos_e_bisec(e, d, Ye, abar, T, status)
Definition: ls_timmes_eos_module.f:177
ls_timmes_eos_module::timmes_eos_state
Definition: ls_timmes_eos_module.f:23
ls_timmes_eos_module::dref
real(r_kind), parameter dref
Definition: ls_timmes_eos_module.f:11
ls_timmes_eos_module::timmes_eos_s_nr
subroutine timmes_eos_s_nr(s, d, Ye, abar, T, status)
Definition: ls_timmes_eos_module.f:433