Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2006, 2007-2019 Frank Timmes & The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module helm_alloc
21 :
22 : use const_def, only: dp, use_mesa_temp_cache
23 : use math_lib
24 :
25 : implicit none
26 :
27 : contains
28 :
29 13 : subroutine alloc_helm_table(h, imax, jmax, ierr)
30 : ! This routine allocates a Helm_Table and places pointer to it in h.
31 : ! It also allocates the arrays in the Helm_Table record.
32 :
33 : use eos_def
34 :
35 : type (Helm_Table), pointer :: h
36 : integer, intent(in) :: imax, jmax
37 : integer, intent(out) :: ierr ! 0 means AOK.
38 :
39 : ierr = 0
40 :
41 13 : allocate(h,stat=ierr)
42 13 : if (ierr /= 0) return
43 :
44 13 : h% imax = imax
45 13 : h% jmax = jmax
46 13 : h% with_coulomb_corrections = .true.
47 :
48 13 : call alloc_1d_array(h% d, imax)
49 13 : call alloc_1d_array(h% t, jmax)
50 :
51 : !..for the helmholtz free energy tables
52 13 : call alloc_2d_array(h% f, imax, jmax)
53 13 : call alloc_2d_array(h% fd, imax, jmax)
54 13 : call alloc_2d_array(h% ft, imax, jmax)
55 13 : call alloc_2d_array(h% fdd, imax, jmax)
56 13 : call alloc_2d_array(h% ftt, imax, jmax)
57 13 : call alloc_2d_array(h% fdt, imax, jmax)
58 13 : call alloc_2d_array(h% fddt, imax, jmax)
59 13 : call alloc_2d_array(h% fdtt, imax, jmax)
60 13 : call alloc_2d_array(h% fddtt, imax, jmax)
61 :
62 : !..for the pressure derivative with density tables
63 13 : call alloc_2d_array(h% dpdf, imax, jmax)
64 13 : call alloc_2d_array(h% dpdfd, imax, jmax)
65 13 : call alloc_2d_array(h% dpdft, imax, jmax)
66 13 : call alloc_2d_array(h% dpdfdt, imax, jmax)
67 :
68 : !..for chemical potential tables
69 13 : call alloc_2d_array(h% ef, imax, jmax)
70 13 : call alloc_2d_array(h% efd, imax, jmax)
71 13 : call alloc_2d_array(h% eft, imax, jmax)
72 13 : call alloc_2d_array(h% efdt, imax, jmax)
73 :
74 : !..for the number density tables
75 13 : call alloc_2d_array(h% xf, imax, jmax)
76 13 : call alloc_2d_array(h% xfd, imax, jmax)
77 13 : call alloc_2d_array(h% xft, imax, jmax)
78 13 : call alloc_2d_array(h% xfdt, imax, jmax)
79 :
80 : !..for storing the differences
81 13 : call alloc_1d_array(h% dt_sav, jmax)
82 13 : call alloc_1d_array(h% dt2_sav, jmax)
83 13 : call alloc_1d_array(h% dti_sav, jmax)
84 13 : call alloc_1d_array(h% dt2i_sav, jmax)
85 13 : call alloc_1d_array(h% dt3i_sav, jmax)
86 :
87 13 : call alloc_1d_array(h% dd_sav, imax)
88 13 : call alloc_1d_array(h% dd2_sav, imax)
89 13 : call alloc_1d_array(h% ddi_sav, imax)
90 13 : call alloc_1d_array(h% dd2i_sav, imax)
91 13 : call alloc_1d_array(h% dd3i_sav, imax)
92 :
93 : contains
94 :
95 156 : subroutine alloc_1d_array(ptr,sz)
96 : real(dp), dimension(:), pointer :: ptr
97 : integer, intent(in) :: sz
98 156 : allocate(ptr(sz),stat=ierr)
99 13 : end subroutine alloc_1d_array
100 :
101 273 : subroutine alloc_2d_array(ptr,sz1,sz2)
102 : real(dp), dimension(:,:), pointer :: ptr
103 : integer, intent(in) :: sz1,sz2
104 546 : allocate(ptr(sz1,sz2),stat=ierr)
105 273 : end subroutine alloc_2d_array
106 :
107 :
108 : end subroutine alloc_helm_table
109 :
110 :
111 13 : subroutine setup_td_deltas(h, imax, jmax)
112 : use eos_def
113 : type (Helm_Table), pointer :: h
114 : integer, intent(in) :: imax, jmax
115 : integer :: i, j
116 13 : real(dp) :: dth,dt2,dti,dt2i,dt3i,dd,dd2,ddi,dd2i,dd3i
117 : !..construct the temperature and density deltas and their inverses
118 13013 : do j=1,jmax-1
119 13000 : dth = h% t(j+1) - h% t(j)
120 13000 : dt2 = dth * dth
121 13000 : dti = 1.0d0/dth
122 13000 : dt2i = 1.0d0/dt2
123 13000 : dt3i = dt2i*dti
124 13000 : h% dt_sav(j) = dth
125 13000 : h% dt2_sav(j) = dt2
126 13000 : h% dti_sav(j) = dti
127 13000 : h% dt2i_sav(j) = dt2i
128 13013 : h% dt3i_sav(j) = dt3i
129 : end do
130 35113 : do i=1,imax-1
131 35100 : dd = h% d(i+1) - h% d(i)
132 35100 : dd2 = dd * dd
133 35100 : ddi = 1.0d0/dd
134 35100 : dd2i = 1.0d0/dd2
135 35100 : dd3i = dd2i*ddi
136 35100 : h% dd_sav(i) = dd
137 35100 : h% dd2_sav(i) = dd2
138 35100 : h% ddi_sav(i) = ddi
139 35100 : h% dd2i_sav(i) = dd2i
140 35113 : h% dd3i_sav(i) = dd3i
141 : end do
142 13 : end subroutine setup_td_deltas
143 :
144 :
145 13 : subroutine read_helm_table(h, data_dir, cache_dir, temp_cache_dir, use_cache, ierr)
146 13 : use eos_def
147 : use utils_lib, only: mv, switch_str
148 :
149 :
150 : type (Helm_Table), pointer :: h
151 : character(*), intent(IN) :: data_dir, cache_dir, temp_cache_dir
152 : logical, intent(IN) :: use_cache
153 : integer, intent(out) :: ierr
154 :
155 : !..this routine reads the helmholtz eos file, and
156 : !..must be called once before the helmeos routine is invoked.
157 :
158 : !..declare local variables
159 : character (len=256) :: filename, message, temp_filename
160 : character (len=500) :: buf
161 : character (len=26) :: s26
162 273 : real(dp), target :: vec_ary(20)
163 : real(dp), pointer :: vec(:)
164 : integer :: i,j,k,ios,imax,jmax,n
165 13 : real(dp) :: tsav,dsav
166 : logical, parameter :: dmp = .false.
167 :
168 13 : ierr = 0
169 13 : vec => vec_ary
170 :
171 : !..read the normal helmholtz free energy table
172 13 : h% logtlo = 3.0d0
173 13 : h% logthi = 13.0d0
174 13 : h% logdlo = -12.0d0
175 13 : h% logdhi = 15.0d0
176 :
177 13 : h% templo = exp10(h% logtlo)
178 13 : h% temphi = exp10(h% logthi)
179 13 : h% denlo = exp10(h% logdlo)
180 13 : h% denhi = exp10(h% logdhi)
181 :
182 13 : imax = h% imax
183 13 : jmax = h% jmax
184 13 : h% logtstp = (h% logthi - h% logtlo)/real(jmax-1,kind=dp)
185 13 : h% logtstpi = 1.0d0/h% logtstp
186 13 : h% logdstp = (h% logdhi - h% logdlo)/real(imax-1,kind=dp)
187 13 : h% logdstpi = 1.0d0/h% logdstp
188 :
189 13 : ios = -1
190 13 : if (use_cache) then
191 12 : write(filename,'(2a)') trim(cache_dir), '/helm_table.bin'
192 : open(unit=19,file=trim(filename), &
193 12 : action='read',status='old',iostat=ios,form='unformatted')
194 : end if
195 :
196 13 : if (ios == 0) then
197 :
198 11 : read(19) imax
199 11 : read(19) jmax
200 :
201 11 : if (imax /= h% imax .or. jmax /= h% jmax) then
202 0 : ios = 1 ! wrong cached info
203 : else
204 11 : read(19) h% f(1:imax,1:jmax)
205 11 : read(19) h% fd(1:imax,1:jmax)
206 11 : read(19) h% ft(1:imax,1:jmax)
207 11 : read(19) h% fdd(1:imax,1:jmax)
208 11 : read(19) h% ftt(1:imax,1:jmax)
209 11 : read(19) h% fdt(1:imax,1:jmax)
210 11 : read(19) h% fddt(1:imax,1:jmax)
211 11 : read(19) h% fdtt(1:imax,1:jmax)
212 11 : read(19) h% fddtt(1:imax,1:jmax)
213 11 : read(19) h% dpdf(1:imax,1:jmax)
214 11 : read(19) h% dpdfd(1:imax,1:jmax)
215 11 : read(19) h% dpdft(1:imax,1:jmax)
216 11 : read(19) h% dpdfdt(1:imax,1:jmax)
217 11 : read(19) h% ef(1:imax,1:jmax)
218 11 : read(19) h% efd(1:imax,1:jmax)
219 11 : read(19) h% eft(1:imax,1:jmax)
220 11 : read(19) h% efdt(1:imax,1:jmax)
221 11 : read(19) h% xf(1:imax,1:jmax)
222 11 : read(19) h% xfd(1:imax,1:jmax)
223 11 : read(19) h% xft(1:imax,1:jmax)
224 11 : read(19) h% xfdt(1:imax,1:jmax)
225 :
226 11022 : do j=1,jmax
227 11011 : tsav = h% logtlo + (j-1)*h% logtstp
228 11022 : h% t(j) = exp10(tsav)
229 : end do
230 29722 : do i=1,imax
231 29711 : dsav = h% logdlo + (i-1)*h% logdstp
232 29722 : h% d(i) = exp10(dsav)
233 : end do
234 : end if
235 :
236 11 : close(unit=19)
237 : end if
238 :
239 13 : if (ios /= 0) then
240 :
241 2 : write(filename,'(2a)') trim(data_dir), '/helm_table.dat'
242 2 : write(*,*) 'read ', trim(filename)
243 :
244 2 : ios = 0
245 2 : open(unit=19,file=trim(filename),action='read',status='old',iostat=ios)
246 2 : if (ios /= 0) then
247 0 : write(*,'(3a,i6)') 'failed to open ', trim(filename), ' : ios ', ios
248 0 : ierr = -1
249 0 : return
250 : end if
251 :
252 2004 : do j=1,jmax
253 2002 : tsav = h% logtlo + (j-1)*h% logtstp
254 2002 : h% t(j) = exp10(tsav)
255 5409406 : do i=1,imax
256 5407402 : dsav = h% logdlo + (i-1)*h% logdstp
257 5407402 : h% d(i) = exp10(dsav)
258 5407402 : read(19,'(a)',iostat=ierr) buf
259 5407402 : if (ierr == 0) call str_to_vector(buf, vec, n, ierr)
260 5407402 : if (ierr /= 0 .or. n /= 9) then
261 0 : write(*,'(a)') 'failed while reading ' // trim(filename)
262 0 : close(19)
263 0 : return
264 : end if
265 5407402 : h% f(i,j) = vec(1)
266 5407402 : h% fd(i,j) = vec(2)
267 5407402 : h% ft(i,j) = vec(3)
268 5407402 : h% fdd(i,j) = vec(4)
269 5407402 : h% ftt(i,j) = vec(5)
270 5407402 : h% fdt(i,j) = vec(6)
271 5407402 : h% fddt(i,j) = vec(7)
272 5407402 : h% fdtt(i,j) = vec(8)
273 5407402 : h% fddtt(i,j) = vec(9)
274 2002 : if (dmp) then
275 : do k=1,9
276 : write(*,'(1pd24.16)',advance='no') vec(k)
277 : end do
278 : write(*,'(A)')
279 : end if
280 : end do
281 : end do
282 :
283 : !..read the pressure derivative with density table
284 2004 : do j=1,jmax
285 5409406 : do i=1,imax
286 5407402 : read(19,'(a)',iostat=ierr) buf
287 5407402 : if (ierr == 0) call str_to_vector(buf, vec, n, ierr)
288 5407402 : if (ierr /= 0 .or. n /= 4) then
289 0 : write(*,'(a)') 'failed while reading ' // trim(filename)
290 0 : close(19)
291 0 : return
292 : end if
293 5407402 : h% dpdf(i,j) = vec(1)
294 5407402 : h% dpdfd(i,j) = vec(2)
295 5407402 : h% dpdft(i,j) = vec(3)
296 5407402 : h% dpdfdt(i,j) = vec(4)
297 2002 : if (dmp) then
298 : do k=1,4
299 : write(*,'(1pd24.16)',advance='no') vec(k)
300 : end do
301 : write(*,'(A)')
302 : end if
303 : end do
304 : end do
305 :
306 : !..read the electron chemical potential table
307 2004 : do j=1,jmax
308 5409406 : do i=1,imax
309 5407402 : read(19,'(a)',iostat=ierr) buf
310 5407402 : if (ierr == 0) call str_to_vector(buf, vec, n, ierr)
311 5407402 : if (ierr /= 0 .or. n /= 4) then
312 0 : write(*,'(a)') 'failed while reading ' // trim(filename)
313 0 : close(19)
314 0 : return
315 : end if
316 5407402 : h% ef(i,j) = vec(1)
317 5407402 : h% efd(i,j) = vec(2)
318 5407402 : h% eft(i,j) = vec(3)
319 5407402 : h% efdt(i,j) = vec(4)
320 2002 : if (dmp) then
321 : do k=1,4
322 : write(*,'(1pd24.16)',advance='no') vec(k)
323 : end do
324 : write(*,'(A)')
325 : end if
326 : end do
327 : end do
328 :
329 : !..read the number density table
330 2004 : do j=1,jmax
331 5409406 : do i=1,imax
332 5407402 : read(19,'(a)',iostat=ierr) buf
333 5407402 : if (ierr == 0) call str_to_vector(buf, vec, n, ierr)
334 5407402 : if (ierr /= 0 .or. n /= 4) then
335 0 : write(*,'(a)') 'failed while reading ' // trim(filename)
336 0 : close(19)
337 0 : return
338 : end if
339 5407402 : h% xf(i,j) = vec(1)
340 5407402 : h% xfd(i,j) = vec(2)
341 5407402 : h% xft(i,j) = vec(3)
342 5407402 : h% xfdt(i,j) = vec(4)
343 2002 : if (dmp) then
344 : do k=1,4
345 : write(*,'(1pd24.16)',advance='no') vec(k)
346 : end do
347 : write(*,'(A)')
348 : end if
349 : end do
350 : end do
351 :
352 2 : close(unit=19)
353 : !..write cachefile
354 :
355 : if (dmp) call mesa_error(__FILE__,__LINE__,'helm_alloc')
356 :
357 2 : ios = -1
358 2 : if (use_cache) then
359 1 : write(filename,'(2a)') trim(cache_dir), '/helm_table.bin'
360 1 : write(temp_filename,'(2a)') trim(temp_cache_dir), '/helm_table.bin'
361 1 : write(*,*) 'write ', trim(filename)
362 : open(unit=19,file=trim(switch_str(temp_filename, filename, use_mesa_temp_cache)), &
363 1 : status='replace', iostat=ios,action='write',form='unformatted')
364 : end if
365 :
366 2 : if (ios == 0) then
367 1 : write(19) imax
368 1 : write(19) jmax
369 1 : write(19) h% f(1:imax,1:jmax)
370 1 : write(19) h% fd(1:imax,1:jmax)
371 1 : write(19) h% ft(1:imax,1:jmax)
372 1 : write(19) h% fdd(1:imax,1:jmax)
373 1 : write(19) h% ftt(1:imax,1:jmax)
374 1 : write(19) h% fdt(1:imax,1:jmax)
375 1 : write(19) h% fddt(1:imax,1:jmax)
376 1 : write(19) h% fdtt(1:imax,1:jmax)
377 1 : write(19) h% fddtt(1:imax,1:jmax)
378 1 : write(19) h% dpdf(1:imax,1:jmax)
379 1 : write(19) h% dpdfd(1:imax,1:jmax)
380 1 : write(19) h% dpdft(1:imax,1:jmax)
381 1 : write(19) h% dpdfdt(1:imax,1:jmax)
382 1 : write(19) h% ef(1:imax,1:jmax)
383 1 : write(19) h% efd(1:imax,1:jmax)
384 1 : write(19) h% eft(1:imax,1:jmax)
385 1 : write(19) h% efdt(1:imax,1:jmax)
386 1 : write(19) h% xf(1:imax,1:jmax)
387 1 : write(19) h% xfd(1:imax,1:jmax)
388 1 : write(19) h% xft(1:imax,1:jmax)
389 1 : write(19) h% xfdt(1:imax,1:jmax)
390 1 : close(unit=19)
391 :
392 1 : if (use_mesa_temp_cache) call mv(temp_filename,filename,.true.)
393 :
394 : end if
395 :
396 : end if
397 :
398 13 : call setup_td_deltas(h, imax, jmax)
399 :
400 13 : end subroutine read_helm_table
401 :
402 :
403 5 : subroutine free_helm_table(h)
404 13 : use eos_def
405 :
406 : type (Helm_Table), pointer :: h
407 :
408 10 : call do_free(h% d)
409 10 : call do_free(h% t)
410 :
411 10 : call do_free2(h% f)
412 10 : call do_free2(h% fd)
413 10 : call do_free2(h% ft)
414 10 : call do_free2(h% fdd)
415 10 : call do_free2(h% ftt)
416 10 : call do_free2(h% fdt)
417 10 : call do_free2(h% fddt)
418 10 : call do_free2(h% fdtt)
419 10 : call do_free2(h% fddtt)
420 :
421 : !..for the pressure derivative with density tables
422 10 : call do_free2(h% dpdf)
423 10 : call do_free2(h% dpdfd)
424 10 : call do_free2(h% dpdft)
425 10 : call do_free2(h% dpdfdt)
426 :
427 : !..for chemical potential tables
428 10 : call do_free2(h% ef)
429 10 : call do_free2(h% efd)
430 10 : call do_free2(h% eft)
431 10 : call do_free2(h% efdt)
432 :
433 : !..for the number density tables
434 10 : call do_free2(h% xf)
435 10 : call do_free2(h% xfd)
436 10 : call do_free2(h% xft)
437 10 : call do_free2(h% xfdt)
438 :
439 : !..for storing the differences
440 10 : call do_free(h% dt_sav)
441 10 : call do_free(h% dt2_sav)
442 10 : call do_free(h% dti_sav)
443 10 : call do_free(h% dt2i_sav)
444 10 : call do_free(h% dt3i_sav)
445 10 : call do_free(h% dd_sav)
446 10 : call do_free(h% dd2_sav)
447 10 : call do_free(h% ddi_sav)
448 10 : call do_free(h% dd2i_sav)
449 10 : call do_free(h% dd3i_sav)
450 :
451 5 : deallocate(h)
452 5 : nullify(h)
453 :
454 : contains
455 :
456 60 : subroutine do_free(array_ptr)
457 : real(dp), pointer :: array_ptr(:)
458 10 : if (associated(array_ptr)) deallocate(array_ptr)
459 5 : end subroutine do_free
460 :
461 105 : subroutine do_free2(array_ptr)
462 : real(dp), pointer :: array_ptr(:,:)
463 5 : if (associated(array_ptr)) deallocate(array_ptr)
464 : end subroutine do_free2
465 :
466 : end subroutine free_helm_table
467 :
468 : end module helm_alloc
|
