DISPATCH
eos_mod.f90
1 !===============================================================================
2 !> EOS module, reading tabular EOS compatible with STAGGER code
3 !===============================================================================
4 MODULE eos_mod
5  USE io_mod
6  USE io_unit_mod
7  USE trace_mod
8  USE mpi_mod
9  USE units_mod
10  USE scaling_mod
11  USE kinds_mod
12  implicit none
13  private
14  !-----------------------------------------------------------------------------
15  ! Table related parameters
16  !-----------------------------------------------------------------------------
17  integer iupdte,nvar,md,mtable,njon
18  real(kind=4):: eosxmin,ur,ul,ut,eps,tff,grv,abnd,dbox
19  real(kind=4), allocatable, dimension(:):: tmean,tamp,rhm,xcorr,thmin,thmax, &
20  dth,eemin,eemax,deetab,tab
21  integer, allocatable, dimension(:):: itab,mtab
22  integer, save:: ncall=0, mbox
23  character(len=64):: top, tablefile, file
24  logical, save:: do_eos
25  logical, save:: first_time=.true.
26  !-----------------------------------------------------------------------------
27  ! Public part
28  !-----------------------------------------------------------------------------
29  integer, parameter:: n_lambda=4
30  type, public:: eos_t
31  integer :: n_lambda = n_lambda
32  real :: w_lambda(n_lambda) = 1.0
33  real :: gamma
34  real :: kappa
35  logical :: on = .false.
36  contains
37  procedure:: init
38  procedure:: lookup
39  procedure:: lookup_table
40  procedure:: pressure
41  procedure:: temperature
42  end type
43  integer, save:: verbose=0
44  type(eos_t), public:: eos
45 CONTAINS
46 
47 !===============================================================================
48 !> Initialize EOS module (only once)
49 !===============================================================================
50 SUBROUTINE init (self)
51  class(eos_t):: self
52  real, dimension(:,:,:), pointer:: d, ee, pg
53  real, save :: kappa = 1.0
54  integer:: ir
55  namelist /eos_params/ do_eos, top, tablefile, kappa, verbose
56  !-----------------------------------------------------------------------------
57  call trace%begin ('eos_t%init')
58  !$omp critical (eos_cr)
59  if (first_time) then
60  first_time = .false.
61  do_eos = .true.
62  top = '../../'
63  tablefile = 'data/eos/stagger/table.dat'
64  rewind(io%input); read (io%input,eos_params)
65  if (io%master) write (*,eos_params)
66  !-----------------------------------------------------------------------------
67  ! Read the table file
68  !-----------------------------------------------------------------------------
69  file = trim(top)//'/'//trim(tablefile)
70  open (io_unit%data, file=trim(file), form='unformatted', status='old', action='read')
71  read (io_unit%data) md,iupdte,nvar,mbox,eosxmin,dbox,ul,ut,ur,eps,tff,grv,abnd
72  njon = nvar-mbox
73  allocate (tmean(md),tamp(md),rhm(md),xcorr(md),thmin(md),thmax(md))
74  allocate (dth(md),eemin(md),eemax(md),deetab(md),itab(md),mtab(md))
75  if (io%do_trace.and.io%master) print *,mpi%rank, 'reading 12 times', md, ' values'
76  read (io_unit%data) tmean,tamp,rhm,xcorr,thmin,thmax,dth,eemin,eemax,deetab,itab,mtab
77  read (io_unit%data) mtable
78  allocate (tab(mtable))
79  read (io_unit%data) tab
80  close (io_unit%data)
81  if (mpi%master) then
82  print *, trim(tablefile)//': md, mtable =', md, mtable
83  allocate (d(1,1,1), ee(1,1,1), pg(1,1,1))
84  d = 1.0
85  ee = 5.0
86  !call self%lookup (shape(ee), ee=ee, d=d, pg=pg)
87  print *,'d, ee, pg =', d, ee, pg
88  deallocate (d, ee, pg)
89  if (verbose > 0) then
90  print *, ' density min(E) max(E) min(T) max(T)'
91  do ir=1,md
92  print '(i4,1p,e12.4,2(2x,2e12.4))', ir, rhm(ir), eemin(ir), eemax(ir), &
93  tab(itab(ir)+2*3*mtab(ir)), tab(itab(ir)+(2*3+1)*mtab(ir)-1)
94  end do
95  end if
96  end if
97  end if
98  self%on = do_eos
99  self%kappa = kappa
100  !$omp end critical (eos_cr)
101  !$omp barrier
102  call trace%end()
103 END SUBROUTINE init
104 
105 !-------------------------------------------------------------------------------
106 !-------------------------------------------------------------------------------
107 FUNCTION pressure (self, d, e) RESULT (pg)
108  class(eos_t):: self
109  real(kind=KindScalarVar), dimension(:,:,:), pointer:: d, e
110  real(kind=KindScalarVar), dimension(size(d,1),size(d,2),size(d,3)):: pg
111  !-----------------------------------------------------------------------------
112  pg = (self%gamma-1.0)*d*e
113 END FUNCTION
114 
115 !-------------------------------------------------------------------------------
116 !-------------------------------------------------------------------------------
117 FUNCTION temperature (self, d, e) RESULT (tt)
118  class(eos_t):: self
119  real, dimension(:,:,:), pointer:: d, e
120  real, dimension(size(d,1),size(d,2),size(d,3)):: tt
121  !-----------------------------------------------------------------------------
122  tt = (self%gamma-1.0)*e
123 END FUNCTION
124 
125 !===============================================================================
126 !> General EOS lookup routine, returning values in optional arguments
127 !===============================================================================
128 SUBROUTINE lookup (self, dim, lnd, ee, lnx, x, lny, y, pg, tt, ss, rk, src, gamma)
129  class(eos_t):: self
130  integer:: dim(3)
131  real, optional:: gamma
132  real, dimension(:,:,:), intent(in), pointer, optional :: lnx, x, lny, y, lnd, ee
133  real, dimension(:,:,:), optional :: pg, tt, ss
134  real, dimension(:,:,:,:), optional :: src, rk
135  real, dimension(:,:,:), pointer :: x_loc, y_loc, lnd_l, ee_l
136  integer:: i
137  integer, save:: itimer=0
138  real:: stefan
139  !-----------------------------------------------------------------------------
140  call trace%begin ('eos_t%lookup', itimer=itimer)
141  if (present(gamma)) then
142  self%gamma = gamma
143  end if
144  !-----------------------------------------------------------------------------
145  if (present(y)) then
146  y_loc => y
147  else if (present(lny)) then
148  allocate (y_loc(dim(1),dim(2),dim(3)))
149  y_loc = exp(lny)
150  else
151  call mpi%abort ('eos_t: missing 2nd argument')
152  end if
153  !-----------------------------------------------------------------------------
154  if (present(x)) then
155  x_loc => x
156  else if (present(lnx)) then
157  allocate (x_loc(dim(1),dim(2),dim(3)))
158  x_loc = exp(lnx)
159  else
160  call mpi%abort ('eos_t: missing 1st argument')
161  end if
162  !-----------------------------------------------------------------------------
163  if (present(lnd).and.present(ee)) then
164  if (present(src)) then
165  call lookup_table (self, dim, ee=ee, lnd=lnd, rk=rk, src=src)
166  else
167  call lookup_table (self, dim, ee=ee, lnd=lnd, rk=rk)
168  end if
169  else if (present(src).and.present(rk)) then
170  allocate (lnd_l(dim(1),dim(2),dim(3)))
171  allocate (ee_l(dim(1),dim(2),dim(3)))
172  lnd_l = log(x_loc)
173  ee_l = y_loc/x_loc
174  call lookup_table (self, dim, ee=ee_l, lnd=lnd_l, src=src, rk=rk)
175  deallocate (lnd_l, ee_l)
176  else
177  if (present(tt)) then
178  tt = (self%gamma-1.0)*x_loc
179  end if
180  if (present(src)) then
181  stefan = cgs%stefan/(scaling%p*scaling%u)*scaling%temp**4
182  do i=1,size(src,4)
183  src(:,:,:,i) = stefan*real(((self%gamma-1.0)*x_loc),kind=8)**4
184  end do
185  end if
186  if (present(pg)) then
187  pg = (self%gamma-1.0)*y_loc*x_loc
188  end if
189  if (present(rk)) then
190  do i=1,size(rk,4)
191  rk(:,:,:,i) = y_loc*self%kappa
192  end do
193  end if
194  end if
195  if (present(lnx)) deallocate(x_loc)
196  if (present(lny)) deallocate(y_loc)
197  call trace%end(itimer)
198 END SUBROUTINE lookup
199 
200 !===============================================================================
201 !> General EOS lookup routine, returning values in optional arguments
202 !===============================================================================
203 SUBROUTINE lookup_table (self, dim, e, ee, d, lnd, lne, pg, tt, ne, rk, src, &
204  ss, gamma)
205  class(eos_t):: self
206  integer:: dim(3)
207  real, dimension(:,:,:) , pointer, optional :: lnd,d,lne,e,ee,tt,ne,ss
208  real, dimension(:,:,:,:) , optional :: src,rk
209  real, dimension(:,:,:) , optional :: pg
210  real, optional:: gamma
211  !.............................................................................
212  real, dimension(:,:,:), pointer :: lnd_loc, ee_loc
213  integer, dimension(dim(1)) :: np,np1,ik,ntab,ik1,ntab1
214  real(kind=4), dimension(dim(1)) :: px,py,f00,f01,f10,f11, &
215  fx00,fx01,fx10,fx11,fy00,fy01,fy10,fy11
216  real :: rhm1,rhm2,drhm,algrk,eek,eek1,py1
217  real :: qx,qy,pxqx,pxpx,qxqx,pypy,pyqy,qyqy,pxpy,qxqy,pxqy,qxpy
218  logical :: newtable
219  integer :: ix,iy,iz,kee,kee1,j
220  integer :: mx,my,mz,mbox,nbelow_d, nbelow_e
221  integer, save:: itimer=0
222  !-----------------------------------------------------------------------------
223  ! Allocate temporary arrays
224  !-----------------------------------------------------------------------------
225  call trace%begin ('eos_t%lookup_table', itimer=itimer)
226  if (io%verbose>2) then
227  print *, 'present(lnd) =', present(lnd)
228  print *, 'present(ee) =', present(ee)
229  print *, 'present(pg) =', present(pg)
230  print *, 'present(src) =', present(src)
231  print *, 'present(rk) =', present(rk)
232  print *, 'dim:', dim
233  end if
234  mx=dim(1)
235  my=dim(2)
236  mz=dim(3)
237  newtable = .false.
238  if (present(lnd)) then
239  lnd_loc => lnd
240  else if (present(d)) then
241  allocate (lnd_loc(dim(1),dim(2),dim(3)))
242  lnd_loc = log(d)
243  else
244  call mpi%abort ('eos_t: missing density')
245  end if
246  if (present(ee)) then
247  ee_loc => ee
248  else if (present(e)) then
249  allocate (ee_loc(dim(1),dim(2),dim(3)))
250  if (present(d)) then
251  ee_loc = e/d
252  else
253  ee_loc = e/exp(lnd_loc)
254  end if
255  else if (present(lne)) then
256  allocate (ee_loc(dim(1),dim(2),dim(3)))
257  ee_loc = exp(lne)
258  else
259  call mpi%abort ('eos_t: missing internal energy')
260  end if
261  !
262  rhm1=log(rhm(1))
263  rhm2=log(rhm(md))
264  drhm=(rhm2-rhm1)/(md-1)
265  !-----------------------------------------------------------------------------
266  ! Loop over points
267  !-----------------------------------------------------------------------------
268  nbelow_d = 0
269  nbelow_e = 0
270  do iz=1,mz
271  do iy=1,my
272  !---------------------------------------------------------------------------
273  ! Density index (6 lops)
274  !---------------------------------------------------------------------------
275  do ix=1,mx
276  algrk=1.+(lnd_loc(ix,iy,iz)-rhm1)/drhm
277  nbelow_d = nbelow_d + merge(1,0,lnd_loc(ix,iy,iz)<rhm1)
278  np(ix)=max0(1,min0(md-1,int(algrk)))
279  px(ix)=algrk-np(ix)
280  end do
281  do ix=1,mx
282  !-------------------------------------------------------------------------
283  ! Internal energy index -- this loop has indirect referencing (9 flops)
284  !-------------------------------------------------------------------------
285  ntab(ix) = mtab(np(ix))
286  eek = 1.+(ee_loc(ix,iy,iz)-eemin(np(ix)))/deetab(np(ix))
287  nbelow_e = nbelow_e + merge(1,0,eek<1.0)
288  eek = max(eek,1.0)
289  kee = min(ntab(ix)-1,max(1,int(eek)))
290  py(ix) = eek-kee
291  ik(ix) = itab(np(ix))+kee-1
292  ntab1(ix) = mtab(np(ix)+1)
293  ! -- the integer offset below adjusts for a possible larger eemin at the
294  ! -- next table density. Zero for square tables (e.g. default table.dat)
295  kee1 = kee - nint((eemin(np(ix)+1)-eemin(np(ix))) &
296  /deetab(np(ix)))
297  kee1 = min0(ntab1(ix)-1,max0(1,kee1))
298  ik1(ix) = itab(np(ix)+1)+kee1-1
299  end do
300  !---------------------------------------------------------------------------
301  ! Separate loop for weights, vectorizes well
302  !---------------------------------------------------------------------------
303  do ix=1,mx
304 ! +5 = 20 flops
305  qx = 1. - px(ix)
306  pxqx = px(ix) * qx
307  pxpx = px(ix) * px(ix)
308  qxqx = qx * qx
309 ! +4 = 24 flops
310  qy = 1. - py(ix)
311  pyqy = py(ix) * qy
312  pypy = py(ix) * py(ix)
313  qyqy = qy * qy
314 ! +4 = 28 flops
315  pxqy = px(ix) * qy
316  pxpy = px(ix) * py(ix)
317  qxqy = qx * qy
318  qxpy = qx * py(ix)
319 ! +5 = 33 flops
320  f00(ix) = qxqy * (1. + pxqx - pxpx + pyqy - pypy)
321  f01(ix) = qxpy * (1. + pxqx - pxpx + pyqy - qyqy)
322  f10(ix) = pxqy * (1. - qxqx + pxqx + pyqy - pypy)
323  f11(ix) = pxpy * (1. - qxqx + pxqx + pyqy - qyqy)
324 ! +4 = 37 flops
325  fx00(ix) = qxqy * pxqx
326  fx01(ix) = qxpy * pxqx
327  fx10(ix) = - pxqy * pxqx
328  fx11(ix) = - pxpy * pxqx
329 ! +4 = 41 flops
330  fy00(ix) = qxqy * pyqy
331  fy01(ix) = - qxpy * pyqy
332  fy10(ix) = pxqy * pyqy
333  fy11(ix) = - pxpy * pyqy
334  end do
335  !---------------------------------------------------------------------------
336  ! Loop over table entries, nvar = 4. Would this be more efficent
337  ! if split into one gather loop and one calculation loop?
338  !---------------------------------------------------------------------------
339  if (present(pg)) then
340  do ix=1,mx
341  pg(ix,iy,iz) = exp( &
342  f00(ix) * tab(ik(ix))+ f01(ix) * tab(ik(ix) + 1) &
343  + f10(ix) * tab(ik1(ix))+ f11(ix) * tab(ik1(ix) + 1) &
344  + fx00(ix) * tab(ik(ix) + 2 * ntab(ix) ) &
345  + fx01(ix) * tab(ik(ix) + 2 * ntab(ix) + 1) &
346  + fx10(ix) * tab(ik1(ix) + 2 * ntab1(ix) ) &
347  + fx11(ix) * tab(ik1(ix) + 2 * ntab1(ix) + 1) &
348  + fy00(ix) * tab(ik(ix) + ntab(ix) ) &
349  + fy01(ix) * tab(ik(ix) + ntab(ix) + 1) &
350  + fy10(ix) * tab(ik1(ix) + ntab1(ix) ) &
351  + fy11(ix) * tab(ik1(ix) + ntab1(ix) + 1) &
352  )
353  end do
354  end if
355  ik = ik + 3*ntab
356  ik1 = ik1 + 3*ntab1
357  !---------------------------------------------------------------------------
358  if (present(rk)) then
359  do ix=1,mx
360  rk(ix,iy,iz,1) = exp( &
361  f00(ix) * tab(ik(ix))+ f01(ix) * tab(ik(ix) + 1) &
362  + f10(ix) * tab(ik1(ix))+ f11(ix) * tab(ik1(ix) + 1) &
363  + fx00(ix) * tab(ik(ix) + 2 * ntab(ix) ) &
364  + fx01(ix) * tab(ik(ix) + 2 * ntab(ix) + 1) &
365  + fx10(ix) * tab(ik1(ix) + 2 * ntab1(ix) ) &
366  + fx11(ix) * tab(ik1(ix) + 2 * ntab1(ix) + 1) &
367  + fy00(ix) * tab(ik(ix) + ntab(ix) ) &
368  + fy01(ix) * tab(ik(ix) + ntab(ix) + 1) &
369  + fy10(ix) * tab(ik1(ix) + ntab1(ix) ) &
370  + fy11(ix) * tab(ik1(ix) + ntab1(ix) + 1) &
371  )
372  end do
373  if (size(rk,4)>1) then
374  do j=2,size(rk,4)
375  rk(:,iy,iz,j) = rk(:,iy,iz,j-1)*10.0
376  end do
377  end if
378  end if
379  ik = ik + 3*ntab
380  ik1 = ik1 + 3*ntab1
381  !---------------------------------------------------------------------------
382  if (present(tt)) then
383  do ix=1,mx
384  tt(ix,iy,iz) = ( &
385  f00(ix) * tab(ik(ix))+ f01(ix) * tab(ik(ix) + 1) &
386  + f10(ix) * tab(ik1(ix))+ f11(ix) * tab(ik1(ix) + 1) &
387  + fx00(ix) * tab(ik(ix) + 2 * ntab(ix) ) &
388  + fx01(ix) * tab(ik(ix) + 2 * ntab(ix) + 1) &
389  + fx10(ix) * tab(ik1(ix) + 2 * ntab1(ix) ) &
390  + fx11(ix) * tab(ik1(ix) + 2 * ntab1(ix) + 1) &
391  + fy00(ix) * tab(ik(ix) + ntab(ix) ) &
392  + fy01(ix) * tab(ik(ix) + ntab(ix) + 1) &
393  + fy10(ix) * tab(ik1(ix) + ntab1(ix) ) &
394  + fy11(ix) * tab(ik1(ix) + ntab1(ix) + 1) &
395  )
396  end do
397  end if
398  ik = ik + 3*ntab
399  ik1 = ik1 + 3*ntab1
400  !---------------------------------------------------------------------------
401  if (present(ne)) then
402  do ix=1,mx
403  ne(ix,iy,iz) = exp( &
404  f00(ix) * tab(ik(ix))+ f01(ix) * tab(ik(ix) + 1) &
405  + f10(ix) * tab(ik1(ix))+ f11(ix) * tab(ik1(ix) + 1) &
406  + fx00(ix) * tab(ik(ix) + 2 * ntab(ix) ) &
407  + fx01(ix) * tab(ik(ix) + 2 * ntab(ix) + 1) &
408  + fx10(ix) * tab(ik1(ix) + 2 * ntab1(ix) ) &
409  + fx11(ix) * tab(ik1(ix) + 2 * ntab1(ix) + 1) &
410  + fy00(ix) * tab(ik(ix) + ntab(ix) ) &
411  + fy01(ix) * tab(ik(ix) + ntab(ix) + 1) &
412  + fy10(ix) * tab(ik1(ix) + ntab1(ix) ) &
413  + fy11(ix) * tab(ik1(ix) + ntab1(ix) + 1) &
414  )
415  end do
416  end if
417  !---------------------------------------------------------------------------
418  if (present(src)) then
419  do j=1,size(src,4)
420  do ix=1,mx
421  ik(ix) = ik(ix) + 3*ntab(ix)
422  ik1(ix) = ik1(ix) + 3*ntab1(ix)
423  src(ix,iy,iz,j) = exp( &
424  f00(ix) * tab(ik(ix))+ f01(ix) * tab(ik(ix) + 1) &
425  + f10(ix) * tab(ik1(ix))+ f11(ix) * tab(ik1(ix) + 1) &
426  + fx00(ix) * tab(ik(ix) + 2 * ntab(ix) ) &
427  + fx01(ix) * tab(ik(ix) + 2 * ntab(ix) + 1) &
428  + fx10(ix) * tab(ik1(ix) + 2 * ntab1(ix) ) &
429  + fx11(ix) * tab(ik1(ix) + 2 * ntab1(ix) + 1) &
430  + fy00(ix) * tab(ik(ix) + ntab(ix) ) &
431  + fy01(ix) * tab(ik(ix) + ntab(ix) + 1) &
432  + fy10(ix) * tab(ik1(ix) + ntab1(ix) ) &
433  + fy11(ix) * tab(ik1(ix) + ntab1(ix) + 1) &
434  )
435  end do
436  end do
437  end if
438  end do
439  end do
440  if (nbelow_d+nbelow_e > 0) then
441  write(io%output,*) &
442  'WARNING: number of points below table d & e:', nbelow_d, nbelow_e
443  end if
444  if (present(d)) then
445  deallocate (lnd_loc)
446  nullify (lnd_loc)
447  end if
448  if (present(e).or.present(lne)) then
449  deallocate (ee_loc)
450  nullify (ee_loc)
451  end if
452  call trace%end (itimer)
453 END SUBROUTINE lookup_table
454 
455 END MODULE eos_mod
Define code units, in terms of (here) CGS units.
Definition: scaling_mod.f90:4
Equation of state module for any sort of tables, provided by a reader.
Definition: eos_mod.f90:4
Fundamental constants in CGS and SI units.
Definition: units_mod.f90:4
Definition: io_mod.f90:4