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 6 : 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 6 : allocate(h,stat=ierr)
42 6 : if (ierr /= 0) return
43 :
44 6 : h% imax = imax
45 6 : h% jmax = jmax
46 6 : h% with_coulomb_corrections = .true.
47 :
48 6 : call alloc_1d_array(h% d, imax)
49 6 : call alloc_1d_array(h% t, jmax)
50 :
51 : !..for the helmholtz free energy tables
52 6 : call alloc_2d_array(h% f, imax, jmax)
53 6 : call alloc_2d_array(h% fd, imax, jmax)
54 6 : call alloc_2d_array(h% ft, imax, jmax)
55 6 : call alloc_2d_array(h% fdd, imax, jmax)
56 6 : call alloc_2d_array(h% ftt, imax, jmax)
57 6 : call alloc_2d_array(h% fdt, imax, jmax)
58 6 : call alloc_2d_array(h% fddt, imax, jmax)
59 6 : call alloc_2d_array(h% fdtt, imax, jmax)
60 6 : call alloc_2d_array(h% fddtt, imax, jmax)
61 :
62 : !..for the pressure derivative with density tables
63 6 : call alloc_2d_array(h% dpdf, imax, jmax)
64 6 : call alloc_2d_array(h% dpdfd, imax, jmax)
65 6 : call alloc_2d_array(h% dpdft, imax, jmax)
66 6 : call alloc_2d_array(h% dpdfdt, imax, jmax)
67 :
68 : !..for chemical potential tables
69 6 : call alloc_2d_array(h% ef, imax, jmax)
70 6 : call alloc_2d_array(h% efd, imax, jmax)
71 6 : call alloc_2d_array(h% eft, imax, jmax)
72 6 : call alloc_2d_array(h% efdt, imax, jmax)
73 :
74 : !..for the number density tables
75 6 : call alloc_2d_array(h% xf, imax, jmax)
76 6 : call alloc_2d_array(h% xfd, imax, jmax)
77 6 : call alloc_2d_array(h% xft, imax, jmax)
78 6 : call alloc_2d_array(h% xfdt, imax, jmax)
79 :
80 : !..for storing the differences
81 6 : call alloc_1d_array(h% dt_sav, jmax)
82 6 : call alloc_1d_array(h% dt2_sav, jmax)
83 6 : call alloc_1d_array(h% dti_sav, jmax)
84 6 : call alloc_1d_array(h% dt2i_sav, jmax)
85 6 : call alloc_1d_array(h% dt3i_sav, jmax)
86 :
87 6 : call alloc_1d_array(h% dd_sav, imax)
88 6 : call alloc_1d_array(h% dd2_sav, imax)
89 6 : call alloc_1d_array(h% ddi_sav, imax)
90 6 : call alloc_1d_array(h% dd2i_sav, imax)
91 6 : call alloc_1d_array(h% dd3i_sav, imax)
92 :
93 : contains
94 :
95 72 : subroutine alloc_1d_array(ptr,sz)
96 : real(dp), dimension(:), pointer :: ptr
97 : integer, intent(in) :: sz
98 72 : allocate(ptr(sz),stat=ierr)
99 6 : end subroutine alloc_1d_array
100 :
101 126 : subroutine alloc_2d_array(ptr,sz1,sz2)
102 : real(dp), dimension(:,:), pointer :: ptr
103 : integer, intent(in) :: sz1,sz2
104 252 : allocate(ptr(sz1,sz2),stat=ierr)
105 126 : end subroutine alloc_2d_array
106 :
107 :
108 : end subroutine alloc_helm_table
109 :
110 :
111 6 : 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 : real(dp) :: dth,dt2,dti,dt2i,dt3i,dd,dd2,ddi,dd2i,dd3i
117 : !..construct the temperature and density deltas and their inverses
118 6006 : do j=1,jmax-1
119 6000 : dth = h% t(j+1) - h% t(j)
120 6000 : dt2 = dth * dth
121 6000 : dti = 1.0d0/dth
122 6000 : dt2i = 1.0d0/dt2
123 6000 : dt3i = dt2i*dti
124 6000 : h% dt_sav(j) = dth
125 6000 : h% dt2_sav(j) = dt2
126 6000 : h% dti_sav(j) = dti
127 6000 : h% dt2i_sav(j) = dt2i
128 6006 : h% dt3i_sav(j) = dt3i
129 : end do
130 16206 : do i=1,imax-1
131 16200 : dd = h% d(i+1) - h% d(i)
132 16200 : dd2 = dd * dd
133 16200 : ddi = 1.0d0/dd
134 16200 : dd2i = 1.0d0/dd2
135 16200 : dd3i = dd2i*ddi
136 16200 : h% dd_sav(i) = dd
137 16200 : h% dd2_sav(i) = dd2
138 16200 : h% ddi_sav(i) = ddi
139 16200 : h% dd2i_sav(i) = dd2i
140 16206 : h% dd3i_sav(i) = dd3i
141 : end do
142 6 : end subroutine setup_td_deltas
143 :
144 :
145 6 : subroutine read_helm_table(h, data_dir, cache_dir, temp_cache_dir, use_cache, ierr)
146 6 : 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 : real(dp), target :: vec_ary(20)
163 : real(dp), pointer :: vec(:)
164 : integer :: i,j,k,ios,imax,jmax,n
165 : real(dp) :: tsav,dsav
166 : logical, parameter :: dmp = .false.
167 :
168 6 : ierr = 0
169 6 : vec => vec_ary
170 :
171 : !..read the normal helmholtz free energy table
172 6 : h% logtlo = 3.0d0
173 6 : h% logthi = 13.0d0
174 6 : h% logdlo = -12.0d0
175 6 : h% logdhi = 15.0d0
176 :
177 6 : h% templo = exp10(h% logtlo)
178 6 : h% temphi = exp10(h% logthi)
179 6 : h% denlo = exp10(h% logdlo)
180 6 : h% denhi = exp10(h% logdhi)
181 :
182 6 : imax = h% imax
183 6 : jmax = h% jmax
184 6 : h% logtstp = (h% logthi - h% logtlo)/real(jmax-1,kind=dp)
185 6 : h% logtstpi = 1.0d0/h% logtstp
186 6 : h% logdstp = (h% logdhi - h% logdlo)/real(imax-1,kind=dp)
187 6 : h% logdstpi = 1.0d0/h% logdstp
188 :
189 6 : ios = -1
190 6 : if (use_cache) then
191 5 : write(filename,'(2a)') trim(cache_dir), '/helm_table.bin'
192 : open(unit=19,file=trim(filename), &
193 5 : action='read',status='old',iostat=ios,form='unformatted')
194 : end if
195 :
196 6 : if (ios == 0) then
197 :
198 1 : read(19) imax
199 1 : read(19) jmax
200 :
201 1 : if (imax /= h% imax .or. jmax /= h% jmax) then
202 0 : ios = 1 ! wrong cached info
203 : else
204 1 : read(19) h% f(1:imax,1:jmax)
205 1 : read(19) h% fd(1:imax,1:jmax)
206 1 : read(19) h% ft(1:imax,1:jmax)
207 1 : read(19) h% fdd(1:imax,1:jmax)
208 1 : read(19) h% ftt(1:imax,1:jmax)
209 1 : read(19) h% fdt(1:imax,1:jmax)
210 1 : read(19) h% fddt(1:imax,1:jmax)
211 1 : read(19) h% fdtt(1:imax,1:jmax)
212 1 : read(19) h% fddtt(1:imax,1:jmax)
213 1 : read(19) h% dpdf(1:imax,1:jmax)
214 1 : read(19) h% dpdfd(1:imax,1:jmax)
215 1 : read(19) h% dpdft(1:imax,1:jmax)
216 1 : read(19) h% dpdfdt(1:imax,1:jmax)
217 1 : read(19) h% ef(1:imax,1:jmax)
218 1 : read(19) h% efd(1:imax,1:jmax)
219 1 : read(19) h% eft(1:imax,1:jmax)
220 1 : read(19) h% efdt(1:imax,1:jmax)
221 1 : read(19) h% xf(1:imax,1:jmax)
222 1 : read(19) h% xfd(1:imax,1:jmax)
223 1 : read(19) h% xft(1:imax,1:jmax)
224 1 : read(19) h% xfdt(1:imax,1:jmax)
225 :
226 1002 : do j=1,jmax
227 1001 : tsav = h% logtlo + (j-1)*h% logtstp
228 1002 : h% t(j) = exp10(tsav)
229 : end do
230 2702 : do i=1,imax
231 2701 : dsav = h% logdlo + (i-1)*h% logdstp
232 2702 : h% d(i) = exp10(dsav)
233 : end do
234 : end if
235 :
236 1 : close(unit=19)
237 : end if
238 :
239 6 : if (ios /= 0) then
240 :
241 5 : write(filename,'(2a)') trim(data_dir), '/helm_table.dat'
242 5 : write(*,*) 'read ', trim(filename)
243 :
244 5 : ios = 0
245 5 : open(unit=19,file=trim(filename),action='read',status='old',iostat=ios)
246 5 : 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 5010 : do j=1,jmax
253 5005 : tsav = h% logtlo + (j-1)*h% logtstp
254 5005 : h% t(j) = exp10(tsav)
255 13523515 : do i=1,imax
256 13518505 : dsav = h% logdlo + (i-1)*h% logdstp
257 13518505 : h% d(i) = exp10(dsav)
258 13518505 : read(19,'(a)',iostat=ierr) buf
259 13518505 : if (ierr == 0) call str_to_vector(buf, vec, n, ierr)
260 13518505 : 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 13518505 : h% f(i,j) = vec(1)
266 13518505 : h% fd(i,j) = vec(2)
267 13518505 : h% ft(i,j) = vec(3)
268 13518505 : h% fdd(i,j) = vec(4)
269 13518505 : h% ftt(i,j) = vec(5)
270 13518505 : h% fdt(i,j) = vec(6)
271 13518505 : h% fddt(i,j) = vec(7)
272 13518505 : h% fdtt(i,j) = vec(8)
273 13518505 : h% fddtt(i,j) = vec(9)
274 5005 : 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 5010 : do j=1,jmax
285 13523515 : do i=1,imax
286 13518505 : read(19,'(a)',iostat=ierr) buf
287 13518505 : if (ierr == 0) call str_to_vector(buf, vec, n, ierr)
288 13518505 : 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 13518505 : h% dpdf(i,j) = vec(1)
294 13518505 : h% dpdfd(i,j) = vec(2)
295 13518505 : h% dpdft(i,j) = vec(3)
296 13518505 : h% dpdfdt(i,j) = vec(4)
297 5005 : 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 5010 : do j=1,jmax
308 13523515 : do i=1,imax
309 13518505 : read(19,'(a)',iostat=ierr) buf
310 13518505 : if (ierr == 0) call str_to_vector(buf, vec, n, ierr)
311 13518505 : 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 13518505 : h% ef(i,j) = vec(1)
317 13518505 : h% efd(i,j) = vec(2)
318 13518505 : h% eft(i,j) = vec(3)
319 13518505 : h% efdt(i,j) = vec(4)
320 5005 : 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 5010 : do j=1,jmax
331 13523515 : do i=1,imax
332 13518505 : read(19,'(a)',iostat=ierr) buf
333 13518505 : if (ierr == 0) call str_to_vector(buf, vec, n, ierr)
334 13518505 : 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 13518505 : h% xf(i,j) = vec(1)
340 13518505 : h% xfd(i,j) = vec(2)
341 13518505 : h% xft(i,j) = vec(3)
342 13518505 : h% xfdt(i,j) = vec(4)
343 5005 : 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 5 : close(unit=19)
353 : !..write cachefile
354 :
355 : if (dmp) call mesa_error(__FILE__,__LINE__,'helm_alloc')
356 :
357 5 : ios = -1
358 5 : if (use_cache) then
359 4 : write(filename,'(2a)') trim(cache_dir), '/helm_table.bin'
360 4 : write(temp_filename,'(2a)') trim(temp_cache_dir), '/helm_table.bin'
361 4 : write(*,*) 'write ', trim(filename)
362 : open(unit=19,file=trim(switch_str(temp_filename, filename, use_mesa_temp_cache)), &
363 4 : status='replace', iostat=ios,action='write',form='unformatted')
364 : end if
365 :
366 5 : if (ios == 0) then
367 4 : write(19) imax
368 4 : write(19) jmax
369 4 : write(19) h% f(1:imax,1:jmax)
370 4 : write(19) h% fd(1:imax,1:jmax)
371 4 : write(19) h% ft(1:imax,1:jmax)
372 4 : write(19) h% fdd(1:imax,1:jmax)
373 4 : write(19) h% ftt(1:imax,1:jmax)
374 4 : write(19) h% fdt(1:imax,1:jmax)
375 4 : write(19) h% fddt(1:imax,1:jmax)
376 4 : write(19) h% fdtt(1:imax,1:jmax)
377 4 : write(19) h% fddtt(1:imax,1:jmax)
378 4 : write(19) h% dpdf(1:imax,1:jmax)
379 4 : write(19) h% dpdfd(1:imax,1:jmax)
380 4 : write(19) h% dpdft(1:imax,1:jmax)
381 4 : write(19) h% dpdfdt(1:imax,1:jmax)
382 4 : write(19) h% ef(1:imax,1:jmax)
383 4 : write(19) h% efd(1:imax,1:jmax)
384 4 : write(19) h% eft(1:imax,1:jmax)
385 4 : write(19) h% efdt(1:imax,1:jmax)
386 4 : write(19) h% xf(1:imax,1:jmax)
387 4 : write(19) h% xfd(1:imax,1:jmax)
388 4 : write(19) h% xft(1:imax,1:jmax)
389 4 : write(19) h% xfdt(1:imax,1:jmax)
390 4 : close(unit=19)
391 :
392 4 : if (use_mesa_temp_cache) call mv(temp_filename,filename,.true.)
393 :
394 : end if
395 :
396 : end if
397 :
398 6 : call setup_td_deltas(h, imax, jmax)
399 :
400 6 : end subroutine read_helm_table
401 :
402 :
403 2 : subroutine free_helm_table(h)
404 6 : use eos_def
405 :
406 : type (Helm_Table), pointer :: h
407 :
408 4 : call do_free(h% d)
409 4 : call do_free(h% t)
410 :
411 4 : call do_free2(h% f)
412 4 : call do_free2(h% fd)
413 4 : call do_free2(h% ft)
414 4 : call do_free2(h% fdd)
415 4 : call do_free2(h% ftt)
416 4 : call do_free2(h% fdt)
417 4 : call do_free2(h% fddt)
418 4 : call do_free2(h% fdtt)
419 4 : call do_free2(h% fddtt)
420 :
421 : !..for the pressure derivative with density tables
422 4 : call do_free2(h% dpdf)
423 4 : call do_free2(h% dpdfd)
424 4 : call do_free2(h% dpdft)
425 4 : call do_free2(h% dpdfdt)
426 :
427 : !..for chemical potential tables
428 4 : call do_free2(h% ef)
429 4 : call do_free2(h% efd)
430 4 : call do_free2(h% eft)
431 4 : call do_free2(h% efdt)
432 :
433 : !..for the number density tables
434 4 : call do_free2(h% xf)
435 4 : call do_free2(h% xfd)
436 4 : call do_free2(h% xft)
437 4 : call do_free2(h% xfdt)
438 :
439 : !..for storing the differences
440 4 : call do_free(h% dt_sav)
441 4 : call do_free(h% dt2_sav)
442 4 : call do_free(h% dti_sav)
443 4 : call do_free(h% dt2i_sav)
444 4 : call do_free(h% dt3i_sav)
445 4 : call do_free(h% dd_sav)
446 4 : call do_free(h% dd2_sav)
447 4 : call do_free(h% ddi_sav)
448 4 : call do_free(h% dd2i_sav)
449 4 : call do_free(h% dd3i_sav)
450 :
451 2 : deallocate(h)
452 2 : nullify(h)
453 :
454 : contains
455 :
456 24 : subroutine do_free(array_ptr)
457 : real(dp), pointer :: array_ptr(:)
458 4 : if (associated(array_ptr)) deallocate(array_ptr)
459 2 : end subroutine do_free
460 :
461 42 : subroutine do_free2(array_ptr)
462 : real(dp), pointer :: array_ptr(:,:)
463 2 : 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
|
