Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2017-2019 Rich Townsend & 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 predictive_mix
21 :
22 : use const_def, only: dp, ln10, pi4, boltzm, amu, convective_mixing
23 : use star_private_def
24 : use chem_def
25 : use chem_lib
26 : use num_lib
27 :
28 : implicit none
29 :
30 : private
31 : public :: add_predictive_mixing
32 :
33 : contains
34 :
35 33 : subroutine add_predictive_mixing (s, ierr)
36 :
37 : type(star_info), pointer :: s
38 : integer, intent(out) :: ierr
39 :
40 : logical, parameter :: dbg = .false.
41 :
42 : integer :: i
43 : integer :: j
44 : logical :: is_core_zone
45 : logical :: is_surf_zone
46 : logical :: match_zone_type
47 : logical :: match_zone_loc
48 : logical :: match_bdy_loc
49 : integer :: k
50 :
51 66 : logical :: mix_mask(s%nz)
52 :
53 : include 'formats'
54 :
55 : ! Initialize
56 :
57 33 : ierr = 0
58 :
59 : if (dbg) then
60 : write(*, *) 'add_predictive_mixing; model, n_conv_bdy=', &
61 : s%model_number, s%num_conv_boundaries
62 : end if
63 :
64 : ! Loop over convective boundaries, from center to surface
65 :
66 40766 : mix_mask = .false.
67 :
68 138 : conv_bdy_loop : do i = 1, s%num_conv_boundaries
69 :
70 : ! Skip this boundary if it's at the surface, since we don't
71 : ! predictively mix there
72 :
73 105 : if (s%conv_bdy_loc(i) == 1) then
74 : if (dbg) then
75 : write(*,*) 'skip since s%conv_bdy_loc(i) == 1', i
76 : end if
77 : cycle conv_bdy_loop
78 : end if
79 :
80 : ! Loop over predictive mixing criteria
81 :
82 1818 : criteria_loop : do j = 1, NUM_PREDICTIVE_PARAM_SETS
83 :
84 1680 : if (.NOT. s% predictive_mix(j)) cycle criteria_loop
85 :
86 : ! Check if the criteria match the current boundary
87 :
88 0 : select case (s%predictive_zone_type(j))
89 : case ('burn_H')
90 0 : match_zone_type = s%burn_h_conv_region(i)
91 : case ('burn_He')
92 0 : match_zone_type = s%burn_he_conv_region(i)
93 : case ('burn_Z')
94 0 : match_zone_type = s%burn_z_conv_region(i)
95 : case ('nonburn')
96 : match_zone_type = .NOT. ( &
97 : s%burn_h_conv_region(i) .OR. &
98 : s%burn_he_conv_region(i) .OR. &
99 0 : s%burn_z_conv_region(i) )
100 : case ('any')
101 0 : match_zone_type = .true.
102 : case default
103 0 : write(*,*) 'Invalid predictive_zone_type: j, s%predictive_zone_type(j)=', j, s%predictive_zone_type(j)
104 0 : ierr = -1
105 0 : return
106 : end select
107 :
108 0 : is_core_zone = (i == 1 .AND. s%R_center == 0d0 .AND. s%top_conv_bdy(i))
109 :
110 0 : if (s%top_conv_bdy(i)) then
111 : is_surf_zone = s%conv_bdy_loc(i) == 1
112 : else
113 0 : is_surf_zone = s%conv_bdy_loc(i+1) == 1
114 : end if
115 :
116 : select case (s%predictive_zone_loc(j))
117 : case ('core')
118 0 : match_zone_loc = is_core_zone
119 : case ('shell')
120 0 : match_zone_loc = .NOT. (is_core_zone .OR. is_surf_zone)
121 : case ('surf')
122 0 : match_zone_loc = is_surf_zone
123 : case ('any')
124 0 : match_zone_loc = .true.
125 : case default
126 0 : write(*,*) 'Invalid predictive_zone_loc: j, s%predictive_zone_loc(j)=', j, s%predictive_zone_loc(j)
127 0 : ierr = -1
128 0 : return
129 : end select
130 :
131 0 : select case (s%predictive_bdy_loc(j))
132 : case ('bottom')
133 0 : match_bdy_loc = .NOT. s%top_conv_bdy(i)
134 : case ('top')
135 0 : match_bdy_loc = s%top_conv_bdy(i)
136 : case ('any')
137 0 : match_bdy_loc = .true.
138 : case default
139 0 : write(*,*) 'Invalid predictive_bdy_loc: j, s%predictive_bdy_loc(j)=', j, s%predictive_bdy_loc(j)
140 0 : ierr = -1
141 0 : return
142 : end select
143 :
144 0 : if (.NOT. (match_zone_type .AND. match_zone_loc .AND. match_bdy_loc)) cycle criteria_loop
145 :
146 0 : if (s%conv_bdy_q(i) < s%predictive_bdy_q_min(j) .OR. &
147 : s%conv_bdy_q(i) > s%predictive_bdy_q_max(j)) cycle criteria_loop
148 :
149 : if (dbg) then
150 : write(*,*) 'Predictive mixing at convective boundary: i, j=', i, j
151 : write(*,*) ' s%predictive_zone_type=', TRIM(s%predictive_zone_type(j))
152 : write(*,*) ' s%predictive_zone_loc=', TRIM(s%predictive_zone_loc(j))
153 : write(*,*) ' s%predictive_bdy_loc=', TRIM(s%predictive_bdy_loc(j))
154 : end if
155 :
156 : ! Perform the predictive mixing for this boundary
157 :
158 0 : if (s% do_conv_premix) then
159 0 : call mesa_error(__FILE__,__LINE__,'Predictive mixing and convective premixing cannot be enabled at the same time')
160 0 : stop
161 : end if
162 :
163 : !if (s% MLT_option == 'TDC') then
164 : ! call mesa_error(__FILE__,__LINE__,'Predictive mixing and TDC cannot be enabled at the same time')
165 : ! stop
166 : !end if
167 :
168 0 : call do_predictive_mixing(s, i, j, ierr, mix_mask)
169 0 : if (ierr /= 0) return
170 :
171 : ! Finish (we apply at most a single predictive mix to each boundary)
172 :
173 1785 : exit criteria_loop
174 :
175 : end do criteria_loop
176 :
177 : end do conv_bdy_loop
178 :
179 : ! Perform a sanity check on D_mix
180 :
181 40766 : check_loop : do k = 1, s%nz
182 :
183 40766 : if (is_bad_num(s%D_mix(k))) then
184 :
185 0 : if (s%stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'add_predictive_mixing')
186 :
187 : end if
188 :
189 : end do check_loop
190 :
191 : return
192 :
193 : end subroutine add_predictive_mixing
194 :
195 :
196 0 : subroutine do_predictive_mixing (s, i, j, ierr, mix_mask)
197 :
198 : type(star_info), pointer :: s
199 : integer, intent(in) :: i
200 : integer, intent(in) :: j
201 : integer, intent(out) :: ierr
202 : logical, intent(inout) :: mix_mask(:)
203 :
204 : logical, parameter :: dbg = .false.
205 : logical, parameter :: DUMP_PREDICTIONS = .false.
206 :
207 0 : real(dp) :: superad_thresh
208 0 : real(dp) :: ingest_factor
209 : integer :: iso_id
210 : integer :: iso_r
211 : integer :: iso_i
212 : integer :: k_bot_cz
213 : integer :: k_top_cz
214 : integer :: k_bot_ez
215 : integer :: k_top_ez
216 : integer :: k_bot_mz
217 : integer :: k_top_mz
218 0 : real(dp) :: xa_cz(s%species)
219 0 : real(dp) :: xa_cz_burn(s%species)
220 0 : real(dp) :: xa_ez(s%species)
221 0 : real(dp) :: xa_ez_burn(s%species)
222 0 : real(dp) :: xa_mz(s%species)
223 0 : real(dp) :: xa_mz_burn(s%species)
224 : logical :: outward
225 : logical :: ledoux_extension
226 : integer :: k_a
227 : integer :: k_b
228 0 : real(dp) :: D(s%nz)
229 0 : real(dp) :: vc(s%nz)
230 0 : real(dp) :: gradr(s%nz)
231 0 : real(dp) :: grada(s%nz)
232 0 : real(dp) :: m_ingest
233 0 : real(dp) :: m_ingest_limit
234 0 : real(dp) :: superad_min
235 : integer :: dk
236 : integer :: k
237 0 : real(dp) :: rho
238 0 : real(dp) :: cdc
239 0 : real(dp) :: dg0
240 0 : real(dp) :: dg1
241 : character(256) :: filename
242 : integer :: unit
243 :
244 : ! Perform predictive mixing at the i'th convective boundary,
245 : ! using the j'th set of predictive mixing parameters. This
246 : ! follows the scheme described in the MESA IV instrument paper
247 :
248 0 : ierr = 0
249 :
250 : ! Extract parameters
251 :
252 0 : superad_thresh = MAX(s%predictive_superad_thresh(j), 0._dp)
253 :
254 0 : if (s%predictive_avoid_reversal(j) /= '') then
255 0 : iso_id = chem_get_iso_id(s%predictive_avoid_reversal(j))
256 0 : if (iso_id == nuclide_not_found) then
257 0 : write(*,*) 'Invalid isotope name in predictive_avoid_reversal'
258 0 : ierr = -1
259 0 : return
260 : end if
261 0 : iso_r = s%net_iso(iso_id)
262 : else
263 : iso_r = 0
264 : end if
265 :
266 0 : if (s%predictive_limit_ingestion(j) /= '') then
267 0 : iso_id = chem_get_iso_id(s%predictive_limit_ingestion(j))
268 0 : if (iso_id == nuclide_not_found) then
269 0 : write(*,*) 'Invalid isotope name in predictive_limit_ingestion'
270 0 : ierr = -1
271 0 : return
272 : end if
273 0 : iso_i = s%net_iso(iso_id)
274 0 : ingest_factor = s%predictive_ingestion_factor(j)
275 : else
276 : iso_i = 0
277 : end if
278 :
279 : ! Determine cell indices spanning the initial convection zone,
280 : ! including the cells which contain the actual convective
281 : ! boundaries
282 :
283 0 : if (s%top_conv_bdy(i)) then
284 :
285 0 : if (i == 1) then
286 0 : k_bot_cz = s%nz
287 : else
288 0 : if (s%top_conv_bdy(i-1)) then
289 0 : write(*,*) 'Double top boundary in predictive_mix; i=', i
290 0 : ierr = -1
291 0 : return
292 : end if
293 0 : k_bot_cz = s%conv_bdy_loc(i-1) - 1
294 : end if
295 :
296 0 : k_top_cz = s%conv_bdy_loc(i)
297 :
298 : else
299 :
300 0 : k_bot_cz = s%conv_bdy_loc(i) - 1
301 :
302 0 : if (i == s%num_conv_boundaries) then
303 0 : k_top_cz = 1
304 : else
305 0 : if (.NOT. s%top_conv_bdy(i+1)) then
306 0 : write(*,*) 'Double bottom boundary in predictive_mix; i=', i
307 0 : ierr = -1
308 0 : return
309 : end if
310 0 : k_top_cz = s%conv_bdy_loc(i+1)
311 : end if
312 :
313 : end if
314 :
315 : if (dbg) then
316 : if (k_bot_cz < s%nz) then
317 : write(*,*) 'Predictive mixing: i, j, q_top, q_bot:', i, j, s%q(k_top_cz), s%q(k_bot_cz+1)
318 : else
319 : write(*,*) 'Predictive mixing: i, j, q_top, q_bot:', i, j, s%q(k_top_cz), 0._dp
320 : end if
321 : end if
322 :
323 : ! Determine average abundances of the initial convection zone
324 :
325 0 : call eval_abundances(s, k_bot_cz, k_top_cz, xa_cz, xa_cz_burn)
326 :
327 : ! Decide whether we are starting in the "Ledoux extension" phase,
328 : ! where the boundary moves to where it would be if the
329 : ! Schwarzschild (rather than Ledoux) criterion had been used in
330 : ! mix_info. During Ledoux extension, some of the predictive mixing
331 : ! checks are disabled
332 :
333 0 : ledoux_extension = s%use_ledoux_criterion
334 :
335 : ! Initialize the extended-zone data
336 :
337 0 : k_bot_ez = k_bot_cz
338 0 : k_top_ez = k_top_cz
339 :
340 0 : call eval_abundances(s, k_bot_ez, k_top_ez, xa_ez, xa_ez_burn)
341 :
342 : ! Begin the predictive mixing search, expanding the extent of the
343 : ! mixed zone until one of a number of criteria are met
344 :
345 0 : outward = s%top_conv_bdy(i)
346 :
347 0 : k_bot_mz = k_bot_cz
348 0 : k_top_mz = k_top_cz
349 :
350 : search_loop : do
351 :
352 : ! Grow the mixed zone by one cell
353 :
354 0 : if (outward) then
355 0 : k_top_mz = k_top_mz - 1
356 : else
357 0 : k_bot_mz = k_bot_mz + 1
358 : end if
359 :
360 : ! Exit search if the mixed region has gone out of bounds
361 :
362 0 : if (( outward .AND. k_top_mz < 1) .OR. &
363 : (.NOT. outward .AND. k_bot_mz > s%nz)) then
364 : exit search_loop
365 : end if
366 :
367 : ! Evaluate average abundance in the mixed zone
368 :
369 0 : call eval_abundances(s, k_bot_mz, k_top_mz, xa_mz, xa_mz_burn)
370 :
371 : ! Evaluate mixing coefficients D and vc, together with grad's,
372 : ! at faces inside the mixed zone
373 :
374 : call eval_mixing_coeffs(s, k_bot_mz, k_top_mz, xa_mz_burn, &
375 0 : k_a, k_b, D, vc, grada, gradr, ierr)
376 0 : if (ierr /= 0) return
377 :
378 : ! Update the Ledoux extension flag
379 :
380 0 : if (ledoux_extension) then
381 :
382 : ! Check if we've reached the Schwarzschild boundary
383 :
384 0 : if (.NOT. ALL(s%gradr(k_a:k_b) > s%grada_face(k_a:k_b))) then
385 :
386 : ledoux_extension = .false.
387 :
388 : else
389 :
390 : ! Otherwise, update the extended-zone data
391 :
392 0 : k_bot_ez = k_bot_mz
393 0 : k_top_ez = k_top_mz
394 :
395 0 : call eval_abundances(s, k_bot_ez, k_top_ez, xa_ez, xa_ez_burn)
396 :
397 : end if
398 :
399 : end if
400 :
401 : ! Perform checks that only occur after the Ledoux extension
402 : ! phase has completed
403 :
404 : if (.NOT. ledoux_extension) then
405 :
406 : ! Check whether the predictive mixing will lead to a
407 : ! reversal in the abundance evolution of isotope iso_r due
408 : ! to nuclear burning; if so, finish the search.
409 :
410 0 : if (iso_r /= 0) then
411 :
412 0 : if (SIGN(1._dp, xa_mz_burn(iso_r)-xa_ez(iso_r)) /= SIGN(1._dp, xa_ez_burn(iso_r)-xa_ez(iso_r))) then
413 : if (dbg) then
414 : write(*,*) 'Exiting predictive search due to abundance reversal'
415 : end if
416 : exit search_loop
417 : end if
418 :
419 : end if
420 :
421 : ! Check whether the predictive mixing will cause the
422 : ! ingestion rate for isotope iso_i to exceed the limit
423 :
424 0 : if (iso_i /= 0) then
425 :
426 : ! Calculate the mass ingested
427 :
428 : m_ingest = xa_mz(iso_i)*SUM(s%dm(k_top_mz:k_bot_mz)) - &
429 0 : xa_ez(iso_i)*SUM(s%dm(k_top_ez:k_bot_ez))
430 :
431 : ! Calculate the limiting ingestion mass
432 :
433 0 : if (outward) then
434 0 : call eval_ingest_limit(s, grada, gradr, ingest_factor, k_a, m_ingest_limit)
435 : else
436 0 : call eval_ingest_limit(s, grada, gradr, ingest_factor, k_b, m_ingest_limit)
437 : end if
438 :
439 : ! If the mass ingested exceeds the limit, finish the search
440 :
441 0 : if (m_ingest > m_ingest_limit) then
442 : if (dbg) then
443 : write(*,*) 'Exiting predictive search due to ingestion limit exceeded'
444 : end if
445 : exit search_loop
446 : end if
447 :
448 : end if
449 :
450 : end if
451 :
452 : ! See if the growing boundary face of the mixed region is
453 : ! non-convective; this signals that we've gone too far, and so
454 : ! finish the search
455 :
456 0 : if (( outward .AND. gradr(k_a) < grada(k_a)) .OR. &
457 : (.NOT. outward .AND. gradr(k_b) < grada(k_b))) then
458 : if (dbg) then
459 : write(*,*) 'Exiting predictive search due to non-convective growing boundary'
460 : end if
461 : exit search_loop
462 : end if
463 :
464 : ! See if any faces (apart from the growing boundary face) have
465 : ! a ratio gradr/grada below the superadiabaticity threshold;
466 : ! this signals we're too close to splitting, and so finish the
467 : ! search
468 :
469 : if (outward) then
470 0 : superad_min = MINVAL(gradr(k_a+1:k_b)/grada(k_a+1:k_b)) - 1._dp
471 : else
472 0 : superad_min = MINVAL(gradr(k_a:k_b-1)/grada(k_a:k_b-1)) - 1._dp
473 : end if
474 :
475 0 : if (superad_min <= superad_thresh) then
476 : if (dbg) then
477 : write(*,*) 'Exiting predictive search due to convection-zone split'
478 : end if
479 : exit search_loop
480 : end if
481 :
482 : end do search_loop
483 :
484 : ! If necessary, dump out the mixing prediction for this boundary
485 :
486 : if (DUMP_PREDICTIONS) then
487 : write(filename, 100) i, j, s%model_number
488 : 100 format('pred-mix.i',I0,'.j',I0,'.n',I0)
489 : open(NEWUNIT=unit, FILE=TRIM(filename), STATUS='REPLACE')
490 : dump_loop : do k = k_a, k_b
491 : write(unit, *) k, s%q(k), D(k), gradr(k), grada(k)
492 : end do dump_loop
493 : close(unit)
494 : print *,'Writing prediction data to file:',TRIM(filename)
495 : end if
496 :
497 : ! Back off the mixing by one zone
498 :
499 0 : if (outward) then
500 0 : k_top_mz = k_top_mz + 1
501 : else
502 0 : k_bot_mz = k_bot_mz - 1
503 : end if
504 :
505 : ! Check the mix mask
506 :
507 0 : if (ANY(mix_mask(k_top_mz:k_bot_mz))) then
508 0 : do k = 1, s%nz
509 0 : print *,k,mix_mask(k),k <= k_bot_mz .AND. k >= k_top_mz
510 : end do
511 0 : call mesa_error(__FILE__,__LINE__,'Double predictive')
512 : else
513 0 : mix_mask(k_top_mz:k_bot_mz) = .false.
514 : end if
515 :
516 : ! Return now if no additional mixing should occur
517 :
518 0 : if (outward .AND. k_top_mz == k_top_cz) then
519 : if (dbg) then
520 : write(*,*) 'No predictive mixing at top of zone; boundary i=', i
521 : end if
522 : return
523 0 : elseif (.NOT. outward .AND. k_bot_mz == k_bot_cz) then
524 : if (dbg) then
525 : write(*,*) 'No predictive mixing at bottom of zone; boundary i=', i
526 : end if
527 : return
528 : end if
529 :
530 : ! Re-calculate mixed abundances and mixing coefficients, since we
531 : ! backed off by one zone
532 :
533 0 : call eval_abundances(s, k_bot_mz, k_top_mz, xa_mz, xa_mz_burn)
534 :
535 : call eval_mixing_coeffs(s, k_bot_mz, k_top_mz, xa_mz_burn, &
536 0 : k_a, k_b, D, vc, grada, gradr, ierr)
537 0 : if (ierr /= 0) then
538 : if (dbg) write(*,*) 'Non-zero return from eval_mixing_coeffs in do_predictive_mixing/predictive_mix'
539 : return
540 : end if
541 :
542 : if (outward) then
543 0 : superad_min = MINVAL(gradr(k_a+1:k_b)/grada(k_a+1:k_b)) - 1._dp
544 : else
545 : superad_min = MINVAL(gradr(k_a:k_b-1)/grada(k_a:k_b-1)) - 1._dp
546 : end if
547 :
548 : ! Update the model with the new D and vc
549 :
550 0 : if (outward) then
551 0 : k_b = k_a
552 0 : k_a = k_top_cz
553 0 : dk = -1
554 : else
555 0 : k_a = k_bot_cz + 1
556 0 : dk = 1
557 : end if
558 :
559 0 : face_loop : do k = k_a, k_b, dk
560 :
561 : ! See if we would overwrite an adjacent convection zone; if so,
562 : ! stop
563 :
564 0 : if (s%mixing_type(k) == convective_mixing) then
565 0 : k_b = k - dk
566 0 : exit face_loop
567 : end if
568 :
569 0 : if (k > 1) then
570 : rho = (s%dq(k-1)*s%rho(k) + s%dq(k)*s%rho(k-1))/ &
571 0 : (s%dq(k-1) + s%dq(k))
572 : else
573 0 : rho = s%rho(k)
574 : end if
575 :
576 0 : cdc = (pi4*s%r(k)*s%r(k)*rho)*(pi4*s%r(k)*s%r(k)*rho)*D(k) ! gm^2/sec
577 :
578 0 : s%cdc(k) = cdc
579 0 : s%D_mix(k) = D(k)
580 0 : s%conv_vel(k) = vc(k)
581 0 : s%mixing_type(k) = convective_mixing
582 :
583 : end do face_loop
584 :
585 : ! Store the mass-fractional location of the new convective
586 : ! boundary, and reset the locations for the old convective
587 : ! boundary (cf. set_cz_boundary_info)
588 :
589 0 : if (outward) then
590 :
591 0 : s%cz_bdy_dq(k_top_cz) = 0._dp
592 :
593 0 : dg0 = s%grada_face(k_b-1) - s%gradr(k_b-1)
594 0 : dg1 = grada(k_b) - gradr(k_b)
595 :
596 0 : if (dg0*dg1 < 0) then
597 0 : s%cz_bdy_dq(k_top_mz) = find0(0._dp, dg0, s%dq(k_top_mz), dg1)
598 0 : if (s%cz_bdy_dq(k_top_mz) < 0._dp .OR. s%cz_bdy_dq(k_top_mz) > s%dq(k_top_mz)) then
599 0 : write(2,*) 'bad cz_bdy_dq in do_predictive_mixing', k_top_mz, s%cz_bdy_dq(k_top_mz), s%dq(k_top_mz)
600 0 : ierr = -1
601 0 : return
602 : end if
603 : end if
604 :
605 : else
606 :
607 0 : s%cz_bdy_dq(k_bot_cz) = 0._dp
608 :
609 0 : dg0 = grada(k_b) - gradr(k_b)
610 0 : dg1 = s%grada_face(k_b+1) - s%gradr(k_b+1)
611 :
612 0 : if (dg0*dg1 < 0) then
613 0 : s%cz_bdy_dq(k_bot_mz) = find0(0._dp, dg0, s%dq(k_bot_mz), dg1)
614 0 : if (s%cz_bdy_dq(k_bot_mz) < 0._dp .or. s%cz_bdy_dq(k_bot_mz) > s%dq(k_bot_mz)) then
615 0 : write(2,*) 'bad cz_bdy_dq in do_predictive_mixing', k_bot_mz, s%cz_bdy_dq(k_bot_mz), s%dq(k_bot_mz)
616 0 : ierr = -1
617 0 : return
618 : end if
619 : end if
620 :
621 : end if
622 :
623 : if (dbg) then
624 : write(*,*) 'Predictive mixing: i, k_a, k_b, q_a, q_b, superad_min=', i, k_a, k_b, s%q(k_a), s%q(k_b), &
625 : superad_min
626 : end if
627 :
628 : return
629 :
630 : end subroutine do_predictive_mixing
631 :
632 :
633 0 : subroutine eval_abundances (s, k_bot, k_top, xa, xa_burn)
634 :
635 : type(star_info), pointer :: s
636 : integer, intent(in) :: k_bot
637 : integer, intent(in) :: k_top
638 : real(dp), intent(out) :: xa(:)
639 : real(dp), intent(out) :: xa_burn(:)
640 :
641 0 : real(dp) :: mc
642 : integer :: l
643 :
644 : ! Evaluate average abundances over cells k_top:k_bot, without and
645 : ! with a burning correction
646 :
647 0 : mc = SUM(s%dm(k_top:k_bot))
648 :
649 0 : do l = 1, s%species
650 :
651 0 : xa(l) = SUM(s%dm(k_top:k_bot)*s%xa(l,k_top:k_bot))/mc
652 :
653 0 : xa_burn(l) = xa(l) + SUM(s%dm(k_top:k_bot)*s%dxdt_nuc(l,k_top:k_bot)*s%dt)/mc
654 :
655 : end do
656 :
657 : ! Apply limiting
658 :
659 0 : xa = MAX(MIN(xa, 1._dp), 0._dp)
660 0 : xa = xa/SUM(xa)
661 :
662 0 : xa_burn = MAX(MIN(xa_burn, 1._dp), 0._dp)
663 0 : xa_burn = xa_burn/SUM(xa_burn)
664 :
665 0 : return
666 :
667 : end subroutine eval_abundances
668 :
669 :
670 0 : subroutine eval_mixing_coeffs (s, k_bot_mz, k_top_mz, xa_mx, k_a, k_b, D, vc, grada, gradr, ierr)
671 :
672 : use eos_def
673 : use micro
674 : use turb_info, only: do1_mlt_2
675 :
676 : type(star_info), pointer :: s
677 : integer, intent(in) :: k_bot_mz
678 : integer, intent(in) :: k_top_mz
679 : real(dp), intent(in) :: xa_mx(:)
680 : integer, intent(out) :: k_a
681 : integer, intent(out) :: k_b
682 : real(dp), intent(out) :: D(:)
683 : real(dp), intent(out) :: vc(:)
684 : real(dp), intent(out) :: grada(:)
685 : real(dp), intent(out) :: gradr(:)
686 : integer, intent(out) :: ierr
687 :
688 : logical, parameter :: dbg = .false.
689 :
690 0 : real(dp) :: xh
691 0 : real(dp) :: xhe
692 0 : real(dp) :: z
693 : real(dp) :: abar
694 0 : real(dp) :: zbar, z53bar
695 0 : real(dp) :: z2bar
696 0 : real(dp) :: ye
697 0 : real(dp) :: mass_correction
698 0 : real(dp) :: sumx
699 : integer :: h1
700 0 : real(dp) :: lnd_save(s%nz)
701 0 : real(dp) :: Cp_save(s%nz)
702 0 : real(dp) :: Cv_save(s%nz)
703 0 : real(dp) :: gamma1_save(s%nz)
704 0 : real(dp) :: grada_save(s%nz)
705 0 : real(dp) :: chiRho_save(s%nz)
706 0 : real(dp) :: chiT_save(s%nz)
707 0 : real(dp) :: lnfree_e_save(s%nz)
708 0 : real(dp) :: d_eos_dlnd_save(num_eos_basic_results,s%nz)
709 0 : real(dp) :: d_eos_dlnT_save(num_eos_basic_results,s%nz)
710 0 : real(dp) :: xa_save(s%species,s%nz)
711 0 : real(dp) :: zbar_save(s%nz)
712 : integer :: k
713 0 : real(dp) :: w
714 0 : real(dp) :: rho_face_save(s%nz)
715 : integer :: op_err
716 : logical :: make_gradr_sticky_in_solver_iters
717 :
718 : ! Evaluate abundance data
719 :
720 : call basic_composition_info(s%species, s%chem_id, xa_mx, &
721 0 : xh, xhe, z, abar, zbar, z2bar, z53bar, ye, mass_correction, sumx)
722 :
723 0 : h1 = s%net_iso(ih1)
724 :
725 : ! Update EOS data for cells spanning the mixed region
726 :
727 0 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
728 : update_cell_loop_eos : do k = k_top_mz, k_bot_mz
729 :
730 : op_err = 0
731 :
732 : lnd_save(k) = s%lnd(k)
733 :
734 : Cp_save(k) = s%Cp(k)
735 : Cv_save(k) = s%Cv(k)
736 :
737 : gamma1_save(k) = s%gamma1(k)
738 : grada_save(k) = s%grada(k)
739 :
740 : chiRho_save(k) = s%chiRho(k)
741 : chiT_save(k) = s%chiT(k)
742 :
743 : lnfree_e_save(k) = s%lnfree_e(k)
744 :
745 : d_eos_dlnd_save(:,k) = s%d_eos_dlnd(:,k)
746 : d_eos_dlnT_save(:,k) = s%d_eos_dlnT(:,k)
747 :
748 : call eval_eos(s, k, z, xh, abar, zbar, xa_mx, &
749 : s%lnd(k), s%Cp(k), s%Cv(k), s%gamma1(k), s%grada(k), &
750 : s%chiRho(k), s%chiT(k), s%lnfree_e(k), &
751 : s%d_eos_dlnd(:,k), s%d_eos_dlnT(:,k), op_err)
752 : if (op_err /= 0) ierr = op_err
753 :
754 : xa_save(:,k) = s%xa(:,k)
755 : zbar_save(k) = s%zbar(k)
756 :
757 : s%xa(:,k) = xa_mx
758 : s%zbar(k) = zbar
759 :
760 : end do update_cell_loop_eos
761 : !$OMP END PARALLEL DO
762 0 : if (ierr /= 0) then
763 0 : if (s% report_ierr) write(*,*) 'Non-zero return from eval_eos in eval_mixing_coeffs/predictive_mix'
764 0 : return
765 : end if
766 :
767 :
768 0 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
769 : update_cell_loop_kap : do k = k_top_mz, k_bot_mz
770 : op_err = 0
771 : call do_kap_for_cell(s, k, op_err)
772 : if (op_err /= 0) ierr = op_err
773 : end do update_cell_loop_kap
774 : !$OMP END PARALLEL DO
775 0 : if (ierr /= 0) then
776 0 : if (s% report_ierr) write(*,*) 'Non-zero return from do_kap_for_cells in eval_mixing_coeffs/predictive_mix'
777 0 : print *,'xa:',xa_mx
778 0 : return
779 : end if
780 :
781 :
782 :
783 : ! Evaluate mixing coefficients D and vc, together with grad's,
784 : ! at faces inside the mixed region
785 :
786 0 : k_a = k_top_mz + 1
787 0 : k_b = k_bot_mz
788 :
789 0 : eval_face_loop: do k = k_a, k_b
790 :
791 : ! Update the face density
792 :
793 0 : rho_face_save(k) = s%rho_face(k)
794 :
795 0 : w = s%dq(k-1)/(s%dq(k-1) + s%dq(k))
796 :
797 0 : s%rho_face(k) = w*exp(s%lnd(k)) + (1._dp-w)*exp(s%lnd(k-1))
798 :
799 : ! Evaluate mixing coefficients etc.
800 : ! Explicitly set gradL_composition_term to 0 in this call.
801 : call do1_mlt_2(s, k, make_gradr_sticky_in_solver_iters, op_err, &
802 0 : s% alpha_mlt(k), 0._dp)
803 0 : if (op_err /= 0) call mesa_error(__FILE__,__LINE__,'non-zero op_err')
804 :
805 0 : D(k) = s%mlt_D(k)
806 0 : vc(k) = s%mlt_vc(k)
807 :
808 0 : grada(k) = s%grada_face(k)
809 0 : gradr(k) = s%gradr(k)
810 :
811 : end do eval_face_loop
812 :
813 : ! Restore EOS and MLT data
814 :
815 0 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
816 : restore_cell_loop : do k = k_top_mz, k_bot_mz
817 :
818 : op_err = 0
819 :
820 : s%lnd(k) = lnd_save(k)
821 :
822 : s%Cp(k) = Cp_save(k)
823 : s%Cv(k) = Cv_save(k)
824 :
825 : s%gamma1(k) = gamma1_save(k)
826 : s%grada(k) = grada_save(k)
827 :
828 : s%chiRho(k) = chiRho_save(k)
829 : s%chiT(k) = chiT_save(k)
830 :
831 : s%lnfree_e(k) = lnfree_e_save(k)
832 :
833 : s%d_eos_dlnd(:,k) = d_eos_dlnd_save(:,k)
834 : s%d_eos_dlnT(:,k) = d_eos_dlnT_save(:,k)
835 :
836 : s%xa(:,k) = xa_save(:,k)
837 : s%zbar(k) = zbar_save(k)
838 :
839 : call do_kap_for_cell(s, k, op_err)
840 : if (op_err /= 0) ierr = op_err
841 :
842 : end do restore_cell_loop
843 : !$OMP END PARALLEL DO
844 0 : if (ierr /= 0) then
845 0 : if (s% report_ierr) write(*,*) 'Non-zero return from do_kap_for_cells in eval_mixing_coeffs/predictive_mix'
846 0 : return
847 : end if
848 :
849 0 : restore_face_loop: do k = k_a, k_b
850 :
851 0 : s%rho_face(k) = rho_face_save(k)
852 0 : call do1_mlt_2(s, k, make_gradr_sticky_in_solver_iters, op_err)
853 0 : if (op_err /= 0) call mesa_error(__FILE__,__LINE__,'non-zero op_err')
854 :
855 : end do restore_face_loop
856 :
857 : return
858 :
859 0 : end subroutine eval_mixing_coeffs
860 :
861 :
862 0 : subroutine eval_eos (s, k, z, x, abar, zbar, xa, &
863 : lnRho, Cp, Cv, gamma1, grada, chiRho, chiT, &
864 0 : lnfree_e, d_dlnd, d_dlnT, ierr)
865 :
866 0 : use eos_def
867 : use eos_support
868 :
869 : type(star_info), pointer :: s
870 : integer, intent(in) :: k
871 : real(dp), intent(in) :: z
872 : real(dp), intent(in) :: x
873 : real(dp), intent(in) :: abar
874 : real(dp), intent(in) :: zbar
875 : real(dp), intent(in) :: xa(:)
876 : real(dp), intent(out) :: lnRho
877 : real(dp), intent(out) :: Cp
878 : real(dp), intent(out) :: Cv
879 : real(dp), intent(out) :: gamma1
880 : real(dp), intent(out) :: grada
881 : real(dp), intent(out) :: chiRho
882 : real(dp), intent(out) :: chiT
883 : real(dp), intent(out) :: lnfree_e
884 : real(dp), intent(out) :: d_dlnd(:)
885 : real(dp), intent(out) :: d_dlnT(:)
886 : integer, intent(out) :: ierr
887 :
888 : logical, parameter :: dbg = .false.
889 : real(dp), parameter :: LOGRHO_TOL = 1E-8_dp
890 : real(dp), parameter :: LOGPGAS_TOL = 1E-8_dp
891 :
892 0 : real(dp) :: logRho
893 0 : real(dp) :: res(num_eos_basic_results)
894 0 : real(dp) :: d_dxa(num_eos_d_dxa_results,s%species)
895 :
896 : ! Evaluate EOS data in cell k, assuming the cell's temperature and
897 : ! pressure are as specified in the model, but with abundances
898 : ! given by xa and other input abundance parameters
899 :
900 : ! (NEEDS FIXING TO HANDLE CASE WHEN LNPGAS_FLAG = .true.)
901 :
902 : call solve_eos_given_PgasT( &
903 : s, k, xa, &
904 : s%lnT(k)/ln10, s%lnPgas(k)/ln10, s%lnd(k)/ln10, LOGRHO_TOL, LOGPGAS_TOL, &
905 : logRho, res, d_dlnd, d_dlnT, d_dxa, &
906 0 : ierr)
907 0 : if (ierr /= 0) then
908 : if (dbg) write(*,*) 'Non-zero return from solve_eos_given_PgasT in eval_eos/predictive_mix'
909 : return
910 : end if
911 :
912 0 : lnRho = logRho*ln10
913 :
914 0 : Cp = res(i_Cp)
915 0 : Cv = res(i_Cv)
916 :
917 0 : gamma1 = res(i_gamma1)
918 0 : grada = res(i_grad_ad)
919 :
920 0 : chiRho = res(i_chiRho)
921 0 : chiT = res(i_chiT)
922 :
923 0 : lnfree_e = res(i_lnfree_e)
924 :
925 0 : return
926 :
927 0 : end subroutine eval_eos
928 :
929 :
930 0 : subroutine eval_ingest_limit (s, grada, gradr, ingest_factor, k, m_ingest_limit)
931 :
932 : type(star_info), pointer :: s
933 : real(dp), intent(in) :: grada(:)
934 : real(dp), intent(in) :: gradr(:)
935 : real(dp), intent(in) :: ingest_factor
936 : integer, intent(in) :: k
937 : real(dp), intent(out) :: m_ingest_limit
938 :
939 0 : real(dp) :: alfa
940 0 : real(dp) :: beta
941 0 : real(dp) :: T_face
942 :
943 : ! Evaluate the mass ingestion limit at face k following Spruit
944 : ! (2015)
945 :
946 : ! Interpolate T to the face
947 :
948 0 : if (k == 1) then
949 : alfa = 1._dp
950 : else
951 0 : alfa = s%dq(k-1)/(s%dq(k-1) + s%dq(k))
952 : end if
953 0 : beta = 1._dp - alfa
954 :
955 0 : T_face = alfa*s%T(k) + beta*s%T(k-1)
956 :
957 : ! Evaluate the limit
958 :
959 0 : m_ingest_limit = ingest_factor*(amu*s%L(k)/(boltzm*T_face))*(1._dp - grada(k)/gradr(k))*s%dt
960 :
961 0 : return
962 :
963 0 : end subroutine eval_ingest_limit
964 :
965 : end module predictive_mix
|
