DISPATCH
eos_mod.f90
1 !===============================================================================
2 !> Equation of state module for any sort of tables, provided by a reader
3 !===============================================================================
4 MODULE eos_mod
5  USE io_mod
6  USE mpi_mod
7  USE trace_mod
8  USE kinds_mod
9  USE scaling_mod
10  USE units_mod
11  implicit none
12  private
13  type, public:: eos_t
14  real:: gamma = 1.4 ! in case eos%init is never called
15  integer :: eos_type = 0 ! 0 - ideal, x - any other
16  integer :: n_lambda = 1
17  real :: w_lambda = 1
18  logical :: do_rt = .false.
19  real :: kappa
20  contains
21  procedure:: init
22  procedure:: pressure
23  procedure:: temperature
24  procedure:: lookup
25  procedure:: lookup_table
26  end type
27  type(eos_t), public:: eos
28 CONTAINS
29 
30 !===============================================================================
31 !> The eos reader is supposed to deliver table date in the eos_reader data
32 !> type, with pressure and internal energy in an (nT,nd,2) table, with NT
33 !> temperature value T(nT), and nd density values d(nd).
34 !===============================================================================
35 SUBROUTINE init (self)
36  class(eos_t):: self
37  integer:: iostat
38  logical, save:: first_time=.true.
39  real, save:: gamma=1.4
40  real, save:: kappa=0.1
41  logical, save:: do_rt = .false.
42  integer, save:: itimer=0
43  namelist /eos_params/ gamma, do_rt, kappa
44  !-----------------------------------------------------------------------------
45  call trace%begin ('eos_t%init', itimer=itimer)
46  if (io%master) print *, &
47  '------------------------ ideal equation-of-state -----------------------------'
48  !$omp critical (init_cr)
49  if (first_time) then
50  gamma = self%gamma
51  do_rt = self%do_rt
52  kappa = self%kappa
53  rewind(io%input); read (io%input, eos_params, iostat=iostat)
54  if (io%master) write (io%output, eos_params)
55  end if
56  !$omp end critical (init_cr)
57  self%gamma = gamma
58  self%do_rt = do_rt
59  self%kappa = kappa * scaling%m/scaling%l**2
60  call trace%end(itimer)
61 END SUBROUTINE init
62 
63 !-------------------------------------------------------------------------------
64 !-------------------------------------------------------------------------------
65 FUNCTION pressure (self, d, e) RESULT (pg)
66  class(eos_t):: self
67  real(kind=KindScalarVar), dimension(:,:,:), pointer:: d, e
68  real(kind=KindScalarVar), dimension(size(d,1),size(d,2),size(d,3)):: pg
69  !-----------------------------------------------------------------------------
70  pg = (self%gamma-1.0)*d*e
71 END FUNCTION
72 
73 !-------------------------------------------------------------------------------
74 !-------------------------------------------------------------------------------
75 FUNCTION temperature (self, d, e) RESULT (tt)
76  class(eos_t):: self
77  real, dimension(:,:,:), pointer:: d, e
78  real, dimension(size(d,1),size(d,2),size(d,3)):: tt
79  !-----------------------------------------------------------------------------
80  tt = (self%gamma-1.0)*e
81 END FUNCTION
82 
83 !===============================================================================
84 !> General EOS lookup routine, returning values in optional arguments
85 !===============================================================================
86 SUBROUTINE lookup (self, dim, lnd, ee, lnx, x, lny, y, pg, tt, ss, rk, src, gamma)
87  class(eos_t):: self
88  integer:: dim(3)
89  real, optional:: gamma
90  real, dimension(:,:,:), intent(in), pointer, optional :: lnx, x, lny, y, lnd, ee
91  real, dimension(:,:,:), optional :: pg, tt, ss
92  real, dimension(:,:,:,:), optional :: src, rk
93  real, dimension(:,:,:), pointer :: x_loc, y_loc
94  integer:: i
95  integer, save:: itimer=0
96  real:: stefan
97  !-----------------------------------------------------------------------------
98  call trace%begin ('eos_t%lookup', itimer=itimer)
99  if (present(gamma)) then
100  self%gamma = gamma
101  end if
102  if (present(lnd)) then
103  call mpi%abort('trying to use lnd in eos/ideal/eos_mod.f90')
104  else if (present(ee)) then
105  call mpi%abort('trying to use ee in eos/ideal/eos_mod.f90')
106  end if
107  !-----------------------------------------------------------------------------
108  if (present(y)) then
109  y_loc => y
110  else if (present(lny)) then
111  allocate (y_loc(dim(1),dim(2),dim(3)))
112  y_loc = exp(lny)
113  else
114  call mpi%abort ('eos_t: missing 2nd argument')
115  end if
116  !-----------------------------------------------------------------------------
117  if (present(x)) then
118  x_loc => x
119  else if (present(lnx)) then
120  allocate (x_loc(dim(1),dim(2),dim(3)))
121  x_loc = exp(lnx)
122  else
123  call mpi%abort ('eos_t: missing 1st argument')
124  end if
125  !-----------------------------------------------------------------------------
126  if (present(tt)) then
127  tt = (self%gamma-1.0)*x_loc
128  end if
129  if (present(src)) then
130  stefan = cgs%stefan/(scaling%p*scaling%u)*scaling%temp**4
131  do i=1,size(src,4)
132  src(:,:,:,i) = stefan*real(((self%gamma-1.0)*x_loc),kind=8)**4
133  end do
134  end if
135  if (present(pg)) then
136  pg = (self%gamma-1.0)*y_loc*x_loc
137  end if
138  if (present(rk)) then
139  do i=1,size(rk,4)
140  rk(:,:,:,i) = y_loc*self%kappa
141  end do
142  end if
143  if (present(lnx)) deallocate(x_loc)
144  if (present(lny)) deallocate(y_loc)
145  call trace%end(itimer)
146 END SUBROUTINE lookup
147 
148 !===============================================================================
149 !> General EOS lookup routine, returning values in optional arguments
150 !===============================================================================
151 !===============================================================================
152 !> General EOS lookup routine, returning values in optional arguments
153 !===============================================================================
154 SUBROUTINE lookup_table (self, dim, e, ee, d, lnd, lne,lnee, pg, tt, ne, rk, src, gamma)
155 
156  class(eos_t):: self
157  integer:: dim(3)
158  real, optional:: gamma
159  real, dimension(:,:,:), intent(in), pointer, optional :: e, lne, d, lnee, lnd, ee
160  real, dimension(:,:,:), optional :: pg, tt, ne
161  real, dimension(:,:,:,:), optional :: src, rk
162  real, dimension(:,:,:), pointer :: d_loc, ee_loc
163  integer:: i
164  integer, save:: itimer=0
165  real:: stefan
166  !-----------------------------------------------------------------------------
167  call trace%begin ('eos_t%lookup', itimer=itimer)
168  if (present(gamma)) then
169  self%gamma = gamma
170  end if
171  if (present(lnd)) then
172  call mpi%abort('trying to use lnd in eos/ideal/eos_mod.f90')
173  else if (present(ee)) then
174  call mpi%abort('trying to use ee in eos/ideal/eos_mod.f90')
175  end if
176  !-----------------------------------------------------------------------------
177  if (present(d)) then
178  d_loc => d
179  else if (present(lnd)) then
180  allocate (d_loc(dim(1),dim(2),dim(3)))
181  d_loc = exp(lnd)
182  else
183  call mpi%abort ('eos_t: missing 1st argument')
184  end if
185  !-----------------------------------------------------------------------------
186  if (present(ee)) then
187  ee_loc => ee
188  else if (present(lne)) then
189  allocate (ee_loc(dim(1),dim(2),dim(3)))
190  ee_loc = exp(lne)/d_loc
191  else if (present(e)) then
192  allocate (ee_loc(dim(1),dim(2),dim(3)))
193  ee_loc = e/d_loc
194  else
195  call mpi%abort ('eos_t: missing 2nd argument')
196  end if
197  !-----------------------------------------------------------------------------
198  if (present(tt)) then
199  tt = (self%gamma-1.0)*ee_loc
200  end if
201  if (present(src)) then
202  stefan = cgs%stefan/(scaling%p*scaling%u)*scaling%temp**4
203  do i=1,size(src,4)
204  src(:,:,:,i) = stefan*real(((self%gamma-1.0)*ee_loc),kind=8)**4
205  end do
206  end if
207  if (present(pg)) then
208  pg = (self%gamma-1.0)*ee_loc*d_loc
209  end if
210  if (present(rk)) then
211  do i=1,size(rk,4)
212  rk(:,:,:,i) = d_loc*self%kappa
213  end do
214  end if
215  if (present(lne).or.present(e).or.present(lnee)) deallocate(ee_loc)
216  if (present(lnd)) deallocate(d_loc)
217  call trace%end(itimer)
218 END SUBROUTINE lookup_table
219 
220 END MODULE eos_mod
221 
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