Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 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 pulse_utils
21 :
22 : use star_private_def
23 : use const_def, only: dp, pi, pi4
24 : use num_lib
25 : use star_utils
26 :
27 : implicit none
28 :
29 : private
30 : public :: set_segment_indices
31 : public :: eval_face
32 : public :: eval_face_X
33 : public :: eval_face_A_ast
34 : public :: eval_face_rho
35 : public :: eval_center
36 : public :: eval_center_X
37 : public :: eval_center_rho
38 : public :: eval_center_d2
39 :
40 : contains
41 :
42 0 : subroutine set_segment_indices (s, k_a, k_b, include_last_face)
43 :
44 : type(star_info), intent(in) :: s
45 : integer, allocatable, intent(out) :: k_a(:)
46 : integer, allocatable, intent(out) :: k_b(:)
47 : logical, intent(in) :: include_last_face
48 :
49 : logical, parameter :: dbg = .false.
50 :
51 0 : real(dp) :: grad_mu(s%nz)
52 0 : logical :: mask(s%nz)
53 : integer :: k
54 0 : integer :: i(s%nz)
55 : integer :: n_mk
56 : integer :: n_sg
57 : integer :: sg
58 :
59 : ! Set up index ranges for segments which are delineated by
60 : ! composition jumps
61 :
62 0 : if (s%add_double_points_to_pulse_data) then
63 :
64 : ! Calculate grad_mu
65 :
66 0 : grad_mu(1) = 0d0
67 :
68 0 : do k = 2, s%nz-1
69 0 : grad_mu(k) = log(s%mu(k)/s%mu(k-1))/log(s%Peos(k)/s%Peos(k-1))
70 : end do
71 :
72 0 : if (include_last_face) then
73 0 : grad_mu(k) = log(s%mu(s%nz)/s%mu(s%nz-1))/log(s%Peos(s%nz)/s%Peos(s%nz-1))
74 : else
75 0 : grad_mu(k) = 0d0
76 : end if
77 :
78 : ! Set up the mask marking faces which will have a double point
79 :
80 0 : mask = ABS(grad_mu) > s%threshold_grad_mu_for_double_point
81 :
82 0 : if (s%max_number_of_double_points > 0) then
83 :
84 : ! Limit the number of marked faces
85 :
86 0 : n_mk = MIN(s%max_number_of_double_points, s%nz)
87 :
88 0 : call qsort(i, s%nz, -ABS(grad_mu))
89 :
90 0 : mask(i(n_mk+1:)) = .false.
91 :
92 : end if
93 :
94 : else
95 :
96 0 : mask = .false.
97 :
98 : end if
99 :
100 : ! Use the mask to set up the index ranges
101 :
102 0 : n_sg = COUNT(mask) + 1
103 :
104 0 : allocate(k_a(n_sg))
105 0 : allocate(k_b(n_sg))
106 :
107 0 : sg = 1
108 :
109 0 : k_a(sg) = 1
110 :
111 0 : do k = 1, s%nz
112 :
113 0 : if (mask(k)) then
114 :
115 0 : k_b(sg) = k - 1
116 0 : sg = sg + 1
117 0 : k_a(sg) = k
118 :
119 : if (dbg) then
120 : write(*, 100) k, grad_mu(k)
121 : 100 format('placing double point at k=', I6, 1X, '(grad_mu=', 1PE10.3, ')')
122 : end if
123 :
124 : end if
125 :
126 : end do
127 :
128 0 : k_b(sg) = s%nz
129 :
130 :
131 0 : return
132 :
133 : end subroutine set_segment_indices
134 :
135 :
136 0 : real(dp) function eval_face (dq, v, k, k_a, k_b, v_lo, v_hi) result (v_face)
137 :
138 : real(dp), intent(in) :: dq(:)
139 : real(dp), intent(in) :: v(:)
140 : integer, intent(in) :: k
141 : integer, intent(in) :: k_a
142 : integer, intent(in) :: k_b
143 : real(dp), intent(in), optional :: v_lo
144 : real(dp), intent(in), optional :: v_hi
145 :
146 : ! Evaluate v at face k, by interpolating (or extrapolating, if
147 : ! k==k_a or k==k_b+1) from cells k_a:k_b
148 :
149 0 : if (k_b < k_a) call mesa_error(__FILE__,__LINE__,'eval_face: invalid segment indices')
150 0 : if (k < k_a .OR. k > k_b+1) call mesa_error(__FILE__,__LINE__,'eval_face: out-of-bounds interpolation')
151 :
152 0 : if (k_b == k_a) then
153 :
154 : ! Using a single cell
155 :
156 0 : v_face = v(k_a)
157 :
158 : else
159 :
160 : ! Using multiple cells
161 :
162 0 : if (k == k_a) then
163 0 : v_face = v(k_a) - dq(k_a)*(v(k_a+1) - v(k_a))/(dq(k_a+1) + dq(k_a))
164 0 : elseif (k == k_a+1) then
165 0 : v_face = v(k_a) + dq(k_a)*(v(k_a+1) - v(k_a))/(dq(k_a+1) + dq(k_a))
166 0 : elseif (k == k_b) then
167 0 : v_face = v(k_b) - dq(k_b)*(v(k_b) - v(k_b-1))/(dq(k_b) + dq(k_b-1))
168 0 : elseif (k == k_b+1) then
169 0 : v_face = v(k_b) + dq(k_b)*(v(k_b) - v(k_b-1))/(dq(k_b) + dq(k_b-1))
170 : else
171 0 : v_face = interp_val_to_pt(v(k_a:k_b), k-k_a+1, k_b-k_a+1, dq(k_a:k_b), 'pulse_utils : eval_face')
172 : end if
173 :
174 : end if
175 :
176 : ! Apply limits
177 :
178 0 : if (PRESENT(v_lo)) then
179 0 : v_face = MAX(v_face, v_lo)
180 : end if
181 :
182 0 : if (PRESENT(v_hi)) then
183 0 : v_face = MIN(v_face, v_hi)
184 : end if
185 :
186 :
187 : return
188 :
189 : end function eval_face
190 :
191 :
192 0 : real(dp) function eval_face_X (s, i, k, k_a, k_b) result (X_face)
193 :
194 : type(star_info), intent(in) :: s
195 : integer, intent(in) :: i
196 : integer, intent(in) :: k
197 : integer, intent(in) :: k_a
198 : integer, intent(in) :: k_b
199 :
200 : ! Evaluate the abundance for species i at face k, by interpolating
201 : ! (or extrapolating, if k==k_a or k==k_b+1) from cells k_a:k_b
202 :
203 0 : if (k_b < k_a) call mesa_error(__FILE__,__LINE__,'eval_face_X: invalid segment indices')
204 0 : if (k < k_a .OR. k > k_b+1) call mesa_error(__FILE__,__LINE__,'eval_face_X: out-of-bounds interpolation')
205 :
206 0 : if (i >= 1) then
207 :
208 0 : if (k_b == k_a) then
209 :
210 : ! Using a single cell
211 :
212 0 : X_face = s%xa(i,k_a)
213 :
214 : else
215 :
216 : ! Using multiple cells
217 :
218 0 : if (k == k_a) then
219 0 : X_face = s%xa(i,k_a) - s%dq(k_a)*(s%xa(i,k_a+1) - s%xa(i,k_a))/(s%dq(k_a+1) + s%dq(k_a))
220 0 : elseif (k == k_a+1) then
221 0 : X_face = s%xa(i,k_a) + s%dq(k_a)*(s%xa(i,k_a+1) - s%xa(i,k_a))/(s%dq(k_a+1) + s%dq(k_a))
222 0 : elseif (k == k_b) then
223 0 : X_face = s%xa(i,k_b) - s%dq(k_b)*(s%xa(i,k_b) - s%xa(i,k_b-1))/(s%dq(k_b) + s%dq(k_b-1))
224 0 : elseif (k == k_b+1) then
225 0 : X_face = s%xa(i,k_b) + s%dq(k_b)*(s%xa(i,k_b) - s%xa(i,k_b-1))/(s%dq(k_b) + s%dq(k_b-1))
226 : else
227 0 : X_face = interp_val_to_pt(s%xa(i,k_a:k_b), k-k_a+1, k_b-k_a+1, s%dq(k_a:k_b), 'pulse_utils : eval_face_X')
228 : end if
229 :
230 : end if
231 :
232 : ! Apply limits
233 :
234 0 : X_face = MIN(1d0, MAX(0d0, X_face))
235 :
236 : else
237 :
238 : X_face = 0d0
239 :
240 : end if
241 :
242 :
243 : return
244 :
245 : end function eval_face_X
246 :
247 :
248 0 : real(dp) function eval_face_A_ast (s, k, k_a, k_b) result (A_ast_face)
249 :
250 : type(star_info), intent(in) :: s
251 : integer, intent(in) :: k
252 : integer, intent(in) :: k_a
253 : integer, intent(in) :: k_b
254 :
255 0 : real(dp) :: A_ast_1
256 0 : real(dp) :: A_ast_2
257 :
258 : ! Evaluate A*=N2*r/g (A_ast) at face k, using data from faces
259 : ! k_a:k_b+1. If k==k_a or k==k_b+1, then use extrapolation from
260 : ! neighboring faces
261 :
262 0 : if (k_b < k_a) call mesa_error(__FILE__,__LINE__,'eval_face_A_ast: invalid segment indices')
263 0 : if (k < k_a .OR. k > k_b+1) call mesa_error(__FILE__,__LINE__,'eval_face_A_ast: out-of-bounds interpolation')
264 :
265 0 : if (.not. s% calculate_Brunt_N2) call mesa_error(__FILE__,__LINE__,'eval_face_A_ast: must have calculate_Brunt_N2 = .true.')
266 :
267 0 : if (k_b == k_a) then
268 :
269 0 : A_ast_face = s%brunt_N2(k)*s%r(k)/s%grav(k)
270 :
271 : else
272 :
273 0 : if (k == k_a) then
274 :
275 0 : A_ast_1 = s%brunt_N2(k_a+1)*s%r(k_a+1)/s%grav(k_a+1)
276 0 : A_ast_2 = s%brunt_N2(k_a+2)*s%r(k_a+2)/s%grav(k_a+2)
277 :
278 0 : A_ast_face = A_ast_1 - s%dq(k_a)*(A_ast_2 - A_ast_1)/s%dq(k_a+1)
279 :
280 0 : elseif (k == k_b+1) then
281 :
282 0 : A_ast_1 = s%brunt_N2(k_b-1)*s%r(k_b-1)/s%grav(k_b)
283 0 : A_ast_2 = s%brunt_N2(k_b)*s%r(k_b)/s%grav(k_b)
284 :
285 0 : A_ast_face = A_ast_2 + s%dq(k_b)*(A_ast_2 - A_ast_1)/s%dq(k_b-1)
286 :
287 : else
288 :
289 0 : A_ast_face = s%brunt_N2(k)*s%r(k)/s%grav(k)
290 :
291 : end if
292 :
293 : end if
294 :
295 :
296 : return
297 :
298 : end function eval_face_A_ast
299 :
300 :
301 0 : real(dp) function eval_face_rho (s, k, k_a, k_b) result (rho_face)
302 :
303 :
304 : type(star_info), intent(in) :: s
305 : integer, intent(in) :: k
306 : integer, intent(in) :: k_a
307 : integer, intent(in) :: k_b
308 :
309 0 : real(dp) :: r
310 0 : real(dp) :: dm
311 0 : real(dp) :: dlnr
312 :
313 : ! Evaluate rho at face k, using data from cells k_a:k_b
314 :
315 0 : if (k_b < k_a) call mesa_error(__FILE__,__LINE__,'eval_face_rho: invalid segment indices')
316 0 : if (k < k_a .OR. k > k_b+1) call mesa_error(__FILE__,__LINE__,'eval_face_rho: out-of-bounds interpolation')
317 :
318 0 : if (k_b == k_a) then
319 :
320 0 : rho_face = s%rho(k)
321 :
322 : else
323 :
324 0 : if (k == k_a) then
325 :
326 0 : r = s%r(k_a)
327 :
328 0 : dm = 0.5d0*s%dm(k_a)
329 0 : dlnr = 1d0 - s%rmid(k_a)/r
330 :
331 0 : elseif (k == k_b+1) then
332 :
333 0 : r = s%r(k_b+1)
334 :
335 0 : dm = 0.5d0*s%dm(k_b)
336 0 : dlnr = s%rmid(k_b)/r - 1d0
337 :
338 : else
339 :
340 0 : r = s%r(k)
341 :
342 0 : dm = 0.5d0*(s%dm(k) + s%dm(k-1))
343 0 : dlnr = (s%rmid(k-1) - s%rmid(k))/r
344 :
345 : end if
346 :
347 0 : rho_face = dm/(pi4*r*r*r*dlnr)
348 :
349 : end if
350 :
351 0 : end function eval_face_rho
352 :
353 :
354 0 : real(dp) function eval_center (r, v, k_a, k_b, v_lo, v_hi) result (v_center)
355 :
356 : real(dp), intent(in) :: r(:)
357 : real(dp), intent(in) :: v(:)
358 : integer, intent(in) :: k_a
359 : integer, intent(in) :: k_b
360 : real(dp), intent(in), optional :: v_lo
361 : real(dp), intent(in), optional :: v_hi
362 :
363 0 : real(dp) :: r_1
364 0 : real(dp) :: r_2
365 0 : real(dp) :: v_1
366 0 : real(dp) :: v_2
367 :
368 : ! Evaluate v at the center, by extrapolating from cells/faces
369 : ! k_a:k_b
370 :
371 0 : if (k_b < k_a) call mesa_error(__FILE__,__LINE__,'eval_center: invalid segment indices')
372 :
373 0 : if (k_a == k_b) then
374 :
375 : ! Using a single cell/face
376 :
377 0 : v_center = v(k_a)
378 :
379 : else
380 :
381 : ! Using the innermost two cells/faces in k_a:k_b; fit a parabola,
382 : ! with dv/dr = 0 at the center
383 :
384 0 : r_1 = r(k_b)
385 0 : r_2 = r(k_b-1)
386 :
387 0 : v_1 = v(k_b)
388 0 : v_2 = v(k_b-1)
389 :
390 0 : v_center = (v_1*r_2*r_2 - v_2*r_1*r_1)/(r_2*r_2 - r_1*r_1)
391 :
392 : end if
393 :
394 : ! Apply limits
395 :
396 0 : if (PRESENT(v_lo)) then
397 0 : v_center = MAX(v_center, v_lo)
398 : end if
399 :
400 0 : if (PRESENT(v_hi)) then
401 0 : v_center = MIN(v_center, v_hi)
402 : end if
403 :
404 0 : end function eval_center
405 :
406 :
407 0 : real(dp) function eval_center_X (s, i, k_a, k_b) result (X_center)
408 :
409 : type(star_info), intent(in) :: s
410 : integer, intent(in) :: i
411 : integer, intent(in) :: k_a
412 : integer, intent(in) :: k_b
413 :
414 0 : real(dp) :: r_1
415 0 : real(dp) :: r_2
416 0 : real(dp) :: X_1
417 0 : real(dp) :: X_2
418 :
419 : ! Evaluate the abundance for species i at the center, by
420 : ! extrapolating from cells k_a:k_b
421 :
422 0 : if (i >= 1) then
423 :
424 0 : if (k_b < k_a) call mesa_error(__FILE__,__LINE__,'eval_center: invalid segment indices')
425 :
426 0 : if (k_a == k_b) then
427 :
428 : ! Using a single cell
429 :
430 0 : X_center = s%xa(i,k_a)
431 :
432 : else
433 :
434 : ! Using the innermost two cells/faces in k_a:k_b; fit a parabola,
435 : ! with dv/dr = 0 at the center
436 :
437 0 : r_1 = s%rmid(k_b)
438 0 : r_2 = s%rmid(k_b-1)
439 :
440 0 : X_1 = s%xa(i,k_b)
441 0 : X_2 = s%xa(i,k_b-1)
442 :
443 0 : X_center = (X_1*r_2*r_2 - X_2*r_1*r_1)/(r_2*r_2 - r_1*r_1)
444 :
445 : end if
446 :
447 : ! Apply limits
448 :
449 0 : X_center = MIN(1d0, MAX(0d0, X_center))
450 :
451 : else
452 :
453 : X_center = 0d0
454 :
455 : end if
456 :
457 :
458 : return
459 :
460 : end function eval_center_X
461 :
462 :
463 0 : real(dp) function eval_center_rho (s, k_b) result (rho_center)
464 :
465 : type(star_info), intent(in) :: s
466 : integer, intent(in) :: k_b
467 :
468 0 : real(dp) :: r_1
469 0 : real(dp) :: rho_1
470 0 : real(dp) :: M_1
471 :
472 : ! Evaluate rho at the center, by extrapolating from cell k_b
473 :
474 : ! Fit a parabola with drho/dr = 0 at the center, which conserves mass
475 :
476 0 : r_1 = s%rmid(k_b)
477 0 : rho_1 = s%rho(k_b)
478 0 : M_1 = s%m(k_b) - 0.5d0*s%dm(k_b)
479 :
480 0 : rho_center = 3d0*(5d0*M_1/(pi*r_1*r_1*r_1) - 4d0*rho_1)/8d0
481 :
482 0 : end function eval_center_rho
483 :
484 :
485 0 : real(dp) function eval_center_d2 (r, v, k_a, k_b) result (d2v_center)
486 :
487 : real(dp), intent(in) :: r(:)
488 : real(dp), intent(in) :: v(:)
489 : integer, intent(in) :: k_a
490 : integer, intent(in) :: k_b
491 :
492 0 : real(dp) :: r_1
493 0 : real(dp) :: r_2
494 0 : real(dp) :: v_1
495 0 : real(dp) :: v_2
496 :
497 : ! Evaluate d2v/dr2 at the center, by extrapolating from
498 : ! cells/faces k_a:k_b
499 :
500 0 : if (k_b < k_a) call mesa_error(__FILE__,__LINE__,'eval_center_d2: invalid segment indices')
501 :
502 0 : if (k_a == k_b) then
503 :
504 : ! Using a single cell/face
505 :
506 : d2v_center = 0d0
507 :
508 : else
509 :
510 : ! Using the innermost two cells/faces; fit a parabola, with
511 : ! dv/dq = 0 at the center
512 :
513 0 : r_1 = r(k_b)
514 0 : r_2 = r(k_b-1)
515 :
516 0 : v_1 = v(k_b)
517 0 : v_2 = v(k_b-1)
518 :
519 0 : d2v_center = 2d0*(v_2 - v_1)/(r_2*r_2 - r_1*r_1)
520 :
521 : end if
522 :
523 0 : end function eval_center_d2
524 :
525 : end module pulse_utils
|
