Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2012-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 diffusion_procs
21 :
22 : use const_def, only: dp, ln10, pi4
23 : use chem_def
24 : use diffusion_support, only: tiny_mass
25 : use star_private_def
26 :
27 : implicit none
28 :
29 : private
30 : public :: fixup
31 : public :: set_new_xa
32 : public :: update_rad_accel_face
33 : public :: get_limit_coeffs
34 : public :: setup_struct_info
35 :
36 : integer, parameter :: ngp = 2
37 : real(dp), public, save :: fk_gam_old(ngp,17)
38 : logical :: initialize_gamma_grid = .true.
39 :
40 : contains
41 :
42 0 : subroutine fixup( &
43 : s, nz, nc, m, nzlo, nzhi, total_num_iters, &
44 : min_X_hard_limit, X_total_atol, X_total_rtol, &
45 0 : cell_dm, mtotal, xtotal_init, X, &
46 0 : lnT, sum_mass, mass, dX_dm, &
47 : bad_j, bad_k, bad_X, bad_sum, bad_Xsum, ierr)
48 :
49 : type (star_info), pointer :: s
50 : integer, intent(in) :: &
51 : nz, nc, m, nzlo, nzhi, total_num_iters
52 : real(dp), intent(in) :: mtotal, min_X_hard_limit, &
53 : X_total_atol, X_total_rtol
54 : real(dp), intent(in), dimension(:) :: cell_dm, xtotal_init, lnT
55 : real(dp), intent(inout), dimension(:,:) :: X
56 : real(dp), intent(inout), dimension(:) :: sum_mass
57 : real(dp), intent(inout), dimension(:,:) :: mass, dX_dm
58 : real(dp), intent(out) :: bad_X, bad_sum, bad_Xsum
59 : integer, intent(out) :: bad_j, bad_k, ierr
60 :
61 : integer :: j, k, op_err, nsmooth_x_in_fixup
62 : real(dp) :: max_lnT_for_smooth
63 :
64 : logical :: dbg
65 :
66 : include 'formats'
67 :
68 0 : ierr = 0
69 0 : nsmooth_x_in_fixup = 0
70 0 : max_lnT_for_smooth = ln10*99 ! max_logT_for_smooth -- always on for now.
71 0 : bad_j = 0
72 0 : bad_k = 0
73 0 : bad_X = 0d0
74 0 : bad_sum = 0d0
75 0 : bad_Xsum = 0d0
76 :
77 0 : dbg = .false. ! total_num_iters == 21
78 :
79 : ! set mass(j,k) using cell_dm and X
80 : ! note that can have bad X's, so can have bad masses too.
81 0 : do k = nzlo,nzhi
82 0 : do j=1,nc
83 0 : if (X(j,k) < min_X_hard_limit) then
84 0 : bad_j = j
85 0 : bad_k = k
86 0 : bad_X = X(j,k)
87 : if (dbg) write(*,3) 'min_X_hard_limit', j, k, X(j,k), min_X_hard_limit
88 : if (dbg) call mesa_error(__FILE__,__LINE__,'fixup')
89 0 : ierr = -1
90 0 : return
91 : end if
92 0 : mass(j,k) = cell_dm(k)*X(j,k)
93 : end do
94 : end do
95 :
96 0 : call fix_negative_masses(nz, nc, nzlo, nzhi, cell_dm, mass, ierr)
97 0 : if (ierr /= 0) return
98 :
99 : call fix_species_conservation( &
100 : nz, nc, nzlo, nzhi, mtotal, xtotal_init, X_total_atol, X_total_rtol, &
101 0 : mass, sum_mass, bad_Xsum, bad_j, ierr)
102 0 : if (ierr /= 0) return
103 :
104 : call fix_single_point_extremes( &
105 0 : nz, nc, nzlo, nzhi, max_lnT_for_smooth, lnT, sum_mass, mass, ierr)
106 0 : if (ierr /= 0) return
107 :
108 : if (nsmooth_x_in_fixup > 0) then
109 : call smooth_x( &
110 : nz, nc, nzlo, nzhi, nsmooth_x_in_fixup, &
111 : max_lnT_for_smooth, lnT, sum_mass, mass, ierr)
112 : if (ierr /= 0) return
113 : end if
114 :
115 0 : !$OMP PARALLEL DO PRIVATE(k, op_err) SCHEDULE(dynamic,2)
116 : do k = nzlo, nzhi
117 : call get1_dX_dm( &
118 : k, nz, nc, nzlo, nzhi, &
119 : mass, sum_mass, &
120 : dX_dm, .false., op_err)
121 : if (op_err /= 0) then
122 : bad_k = k
123 : ierr = op_err
124 : end if
125 : end do
126 : !$OMP END PARALLEL DO
127 :
128 0 : if (ierr /= 0) then
129 : if (dbg) write(*,2) 'failed in get1_dX_dm', bad_k
130 : if (dbg) call mesa_error(__FILE__,__LINE__,'fixup')
131 : return
132 : end if
133 :
134 : if (dbg) write(*,*) 'call redistribute_mass'
135 : call redistribute_mass( &
136 : s, nc, nzlo, nzhi, nz, total_num_iters, &
137 : sum_mass, mass, dX_dm, &
138 0 : xtotal_init, cell_dm, X, .false., ierr)
139 0 : if (ierr /= 0) then
140 : if (dbg) write(*,2) 'failed in redistribute_mass'
141 : if (dbg) call mesa_error(__FILE__,__LINE__,'fixup')
142 : return
143 : end if
144 :
145 : end subroutine fixup
146 :
147 :
148 0 : subroutine fix_single_point_extremes( &
149 0 : nz, nc, nzlo, nzhi, max_lnT_for_smooth, lnT, sum_mass, mass, ierr)
150 :
151 : integer, intent(in) :: nz, nc, nzlo, nzhi
152 : real(dp), intent(in) :: max_lnT_for_smooth, lnT(:)
153 : real(dp), intent(inout) :: sum_mass(:), mass(:,:)
154 : integer, intent(out) :: ierr
155 :
156 : integer :: k, j, k1
157 0 : real(dp) :: xm1, x00, xp1, x1, m0, sum0, m1, sum1, dm
158 :
159 : include 'formats'
160 :
161 0 : ierr = 0
162 :
163 : ! move mass to remove single point extremes in X
164 0 : do k=nzlo+1,nzhi-1
165 0 : if (lnT(k) > max_lnT_for_smooth) exit
166 0 : do j=1,nc
167 0 : xm1 = mass(j,k-1)/sum_mass(k-1)
168 0 : x00 = mass(j,k)/sum_mass(k)
169 0 : xp1 = mass(j,k+1)/sum_mass(k+1)
170 0 : if ((xm1-x00)*(x00-xp1) < -1d-24) then
171 0 : if (abs(xm1-x00) < abs(x00-xp1)) then
172 0 : k1 = k-1; x1 = xm1
173 : else
174 0 : k1 = k+1; x1 = xp1
175 : end if
176 : ! make X(j,k)==X(j,k1) while conserving total j mass
177 0 : m0 = mass(j,k)
178 0 : sum0 = sum_mass(k)
179 0 : m1 = mass(j,k1)
180 0 : sum1 = sum_mass(k1)
181 : ! find dm s.t. xnew = (m0+dm)/(sum0+dm) = (m1-dm)/(sum1-dm)
182 0 : dm = (m0*sum1 - m1*sum0)/(m0 + m1 - sum0 - sum1)
183 0 : if (dm > 0) then ! moving mass from k1
184 0 : dm = min(dm,0.99999d0*sum1)
185 : else ! moving mass from k
186 0 : dm = -min(-dm,0.99999d0*sum0)
187 : end if
188 0 : mass(j,k) = m0 + dm
189 0 : sum_mass(k) = sum(mass(1:nc,k)) ! sum0 + dm
190 0 : if (sum_mass(k) < 0d0) then
191 0 : write(*,2) 'sum_mass(k)', k, sum_mass(k)
192 0 : write(*,1) 'sum0', sum0
193 0 : write(*,1) 'sum1', sum1
194 0 : write(*,1) 'm0', m0
195 0 : write(*,1) 'm1', m1
196 0 : write(*,1) 'dm', dm
197 0 : call mesa_error(__FILE__,__LINE__,'fixup extremes')
198 : end if
199 0 : mass(j,k1) = m1 - dm
200 0 : sum_mass(k1) = sum(mass(1:nc,k1)) ! sum1 - dm
201 0 : if (sum_mass(k1) < 0d0) then
202 0 : write(*,2) 'sum_mass(k1)', k1, sum_mass(k1)
203 0 : write(*,1) 'sum0', sum0
204 0 : write(*,1) 'sum1', sum1
205 0 : write(*,1) 'm0', m0
206 0 : write(*,1) 'm1', m1
207 0 : write(*,1) 'dm', dm
208 0 : call mesa_error(__FILE__,__LINE__,'fixup extremes')
209 : end if
210 : end if
211 : end do
212 : end do
213 :
214 0 : do k=nzlo,nzhi
215 0 : if (sum_mass(k) < 0d0) then
216 0 : write(*,2) 'sum_mass(k)', k, sum_mass(k)
217 0 : call mesa_error(__FILE__,__LINE__,'fixup 3a')
218 : end if
219 0 : do j=1,nc
220 0 : if (mass(j,k) > sum_mass(k)) then
221 0 : write(*,3) 'sum_mass(k)', j, k, mass(j,k), sum_mass(k)
222 0 : call mesa_error(__FILE__,__LINE__,'fixup 3b')
223 : end if
224 : end do
225 : end do
226 :
227 0 : end subroutine fix_single_point_extremes
228 :
229 :
230 0 : subroutine fix_negative_masses( &
231 0 : nz, nc, nzlo, nzhi, cell_dm, mass, ierr)
232 :
233 : integer, intent(in) :: nz, nc, nzlo, nzhi
234 : real(dp), intent(in), dimension(:) :: cell_dm
235 : real(dp), intent(inout), dimension(:,:) :: mass
236 : integer, intent(out) :: ierr
237 :
238 : integer :: k, j, cnt, maxcnt, k_hi, k_lo, kk, jj
239 0 : real(dp) :: dm, source_mass, frac, sum_m
240 :
241 : include 'formats'
242 :
243 0 : ierr = 0
244 :
245 0 : do k = nzlo,nzhi
246 :
247 0 : fix1: do j=1,nc
248 :
249 0 : if (mass(j,k) >= 0d0) cycle fix1
250 0 : if (mass(j,k) >= -1d-13*cell_dm(k)) then
251 0 : mass(j,k) = 0d0
252 0 : cycle fix1
253 : end if
254 :
255 0 : k_hi = min(k+1,nzhi)
256 0 : k_lo = max(k-1,nzlo)
257 0 : maxcnt = 2
258 0 : do cnt = 1, maxcnt
259 0 : sum_m = sum(mass(j,k_lo:k_hi))
260 0 : if (sum_m >= tiny_mass) exit
261 0 : if (cnt == maxcnt .or. mass(j,k_lo) < 0d0 .or. mass(j,k_hi) < 0d0) then
262 0 : mass(j,k) = 0d0
263 0 : cycle fix1
264 : end if
265 0 : k_hi = min(k_hi+1,nzhi)
266 0 : k_lo = max(k_lo-1,nzlo)
267 : end do
268 :
269 0 : dm = -mass(j,k) ! dm > 0
270 0 : mass(j,k) = 0d0
271 : ! remove dm from neighbors
272 0 : source_mass = sum_m + dm
273 0 : frac = sum_m/source_mass
274 0 : do kk = k_lo, k_hi
275 0 : mass(j,kk) = mass(j,kk)*frac
276 : end do
277 0 : if (abs(sum_m - sum(mass(j,k_lo:k_hi))) > 1d-12*sum_m) then
278 :
279 0 : write(*,5) 'bad (sum_m - sum mass(j,:))/sum_m', j, k, k_lo, k_hi, &
280 0 : (sum_m - sum(mass(j,k_lo:k_hi)))/sum_m
281 0 : write(*,1) 'sum_m - sum(mass(j,k_lo:k_hi))', sum_m - sum(mass(j,k_lo:k_hi))
282 0 : write(*,1) 'sum(mass(j,k_lo:k_hi))', sum(mass(j,k_lo:k_hi))
283 0 : write(*,1) 'sum_m', sum_m
284 0 : write(*,1) 'dm', dm
285 :
286 0 : write(*,1) 'sum_m + dm', sum_m + dm
287 0 : write(*,1) 'frac', frac
288 :
289 0 : write(*,1) 'dm/cell_dm(k)', dm/cell_dm(k)
290 0 : write(*,2) 'nzlo', nzlo
291 0 : write(*,2) 'k_lo', k_lo
292 0 : write(*,2) 'k', k
293 0 : write(*,2) 'k_hi', k_hi
294 0 : do jj=k_lo,k_hi
295 0 : write(*,3) 'mass', j, jj, mass(j,jj), mass(j,jj)/cell_dm(jj)
296 : end do
297 0 : call mesa_error(__FILE__,__LINE__,'fixup')
298 : end if
299 :
300 : end do fix1
301 :
302 : end do
303 :
304 0 : do k=nzlo,nzhi
305 0 : do j=1,nc
306 0 : if (mass(j,k) < 0d0) then
307 0 : write(*,3) 'mass(j,k)', j, k, mass(j,k)
308 0 : call mesa_error(__FILE__,__LINE__,'fix_negative_masses')
309 : end if
310 : end do
311 : end do
312 :
313 0 : end subroutine fix_negative_masses
314 :
315 :
316 0 : subroutine fix_species_conservation( &
317 0 : nz, nc, nzlo, nzhi, mtotal, xtotal_init, X_total_atol, X_total_rtol, &
318 0 : mass, sum_mass, bad_Xsum, bad_j, ierr)
319 :
320 : integer, intent(in) :: nz, nc, nzlo, nzhi
321 : real(dp), intent(in) :: &
322 : mtotal, xtotal_init(:), X_total_atol, X_total_rtol
323 : real(dp), intent(inout) :: mass(:,:), sum_mass(:)
324 : real(dp), intent(out) :: bad_Xsum
325 : integer, intent(out) :: bad_j, ierr
326 :
327 : integer :: k, j
328 0 : real(dp) :: xtotal_new, frac, err
329 :
330 : logical, parameter :: dbg = .false.
331 :
332 : include 'formats'
333 :
334 0 : ierr = 0
335 0 : bad_j = 0
336 0 : bad_Xsum = 0d0
337 :
338 0 : do j=1,nc
339 0 : if (xtotal_init(j) < 1d-50) cycle
340 0 : xtotal_new = sum(mass(j,nzlo:nzhi))/mtotal
341 0 : frac = xtotal_new/xtotal_init(j)
342 : err = abs(xtotal_new - xtotal_init(j)) / (X_total_atol + &
343 0 : X_total_rtol*max(xtotal_new, xtotal_init(j)))
344 0 : if (err > 1d0) then
345 : if (dbg) write(*,2) 'fixup err', j, err
346 : if (dbg) write(*,2) 'xtotal_new', j, xtotal_new
347 : if (dbg) write(*,2) 'xtotal_init(j)', j, xtotal_init(j)
348 : if (dbg) write(*,2) 'X_total_atol', j, X_total_atol
349 : if (dbg) write(*,2) 'X_total_rtol', j, X_total_rtol
350 : if (dbg) write(*,2) 'frac', j, frac
351 : if (dbg) call mesa_error(__FILE__,__LINE__,'fixup 1')
352 0 : bad_j = j
353 0 : bad_Xsum = err
354 0 : ierr = -1
355 0 : return
356 : end if
357 :
358 : ! don't try to apply rescaling corrections to anything that
359 : ! ended up completely set to zero by fix_negative_masses
360 0 : if (frac <= 0d0) cycle
361 :
362 0 : do k=nzlo,nzhi
363 0 : mass(j,k) = mass(j,k)/frac
364 : end do
365 : end do
366 :
367 0 : do k=nzlo,nzhi
368 0 : sum_mass(k) = sum(mass(1:nc,k))
369 0 : if (sum_mass(k) < 0d0) then
370 0 : write(*,2) 'sum_mass(k)', k, sum_mass(k)
371 0 : call mesa_error(__FILE__,__LINE__,'fixup 2')
372 : end if
373 : end do
374 :
375 : end subroutine fix_species_conservation
376 :
377 :
378 : subroutine smooth_x( &
379 : nz, nc, nzlo, nzhi, nsmooth_x_in_fixup, &
380 : max_lnT_for_smooth, lnT, sum_mass, mass, ierr)
381 :
382 : integer, intent(in) :: nz, nc, nzlo, nzhi, nsmooth_x_in_fixup
383 : real(dp), intent(in) :: max_lnT_for_smooth, lnT(:)
384 : real(dp), intent(inout) :: sum_mass(:), mass(:,:)
385 : integer, intent(out) :: ierr
386 :
387 : integer :: rep, i, k, j
388 : real(dp) :: sum_m1, sum_00, sum_p1, m_m1, m_00, m_p1, &
389 : xavg, dm, dm_m1, dm_p1, max_lnT
390 :
391 : logical, parameter :: dbg = .false.
392 :
393 : include 'formats'
394 :
395 : ierr = 0
396 :
397 : max_lnT = max_lnT_for_smooth
398 :
399 : do rep = 1, nsmooth_x_in_fixup
400 :
401 : do i = 0, 2
402 :
403 : do k = nzlo+1+i, nzhi-1, 3
404 :
405 : if (lnT(k) > max_lnT) exit
406 :
407 : sum_m1 = sum_mass(k-1)
408 : sum_00 = sum_mass(k)
409 : sum_p1 = sum_mass(k+1)
410 :
411 : do j=1,nc
412 :
413 : m_m1 = mass(j,k-1)
414 : m_00 = mass(j,k)
415 : m_p1 = mass(j,k+1)
416 : if (m_m1 + m_p1 < tiny_mass) cycle
417 : xavg = (m_m1 + m_00 + m_p1)/(sum_m1 + sum_00 + sum_p1)
418 : if (abs(xavg - 1d0) < 1d-10) cycle
419 :
420 : ! find dm s.t. xavg = (m_00 + dm)/(sum_00 + dm)
421 : dm = (sum_00*xavg - m_00)/(1d0 - xavg)
422 : if (dm >= 0d0) then
423 : dm = min(dm, 0.999d0*(m_m1 + m_p1))
424 : else ! dm < 0
425 : dm = -min(-dm, 0.999d0*m_00)
426 : end if
427 :
428 : mass(j,k) = m_00 + dm
429 : if (mass(j,k) < 0d0) then
430 : write(*,3) 'mass00', j, k, mass(j,k), m_00, dm
431 : call mesa_error(__FILE__,__LINE__,'fixup 5')
432 : end if
433 : sum_00 = sum(mass(1:nc,k)) ! sum_00 + dm
434 : if (sum_00 < 0d0 .or. is_bad_num(sum_00)) then
435 : write(*,3) 'sum_00', j, k, sum_00
436 : write(*,3) 'dm', j, k, dm
437 : write(*,3) 'sum_00 - dm', j, k, sum_00 - dm
438 : write(*,3) 'sum_mass(k)', j, k, sum_mass(k)
439 : write(*,'(A)')
440 : write(*,3) 'sum_m1', j, k, sum_m1
441 : write(*,3) 'sum_p1', j, k, sum_p1
442 : write(*,'(A)')
443 : write(*,3) 'm_m1', j, k, m_m1
444 : write(*,3) 'm_00', j, k, m_00
445 : write(*,3) 'm_p1', j, k, m_p1
446 : write(*,'(A)')
447 : write(*,3) 'xavg', j, k, xavg
448 : write(*,'(A)')
449 :
450 : call mesa_error(__FILE__,__LINE__,'fixup')
451 : end if
452 : sum_mass(k) = sum_00
453 :
454 : dm_m1 = dm*m_m1/(m_m1 + m_p1)
455 : mass(j,k-1) = m_m1 - dm_m1
456 : if (mass(j,k-1) < 0d0) then
457 : write(*,3) 'massm1', j, k, mass(j,k-1)
458 : write(*,3) 'm_m1', j, k-1, m_m1
459 : write(*,3) 'm_00', j, k, m_00
460 : write(*,3) 'm_p1', j, k+1, m_p1
461 : write(*,3) 'xavg', j, k, xavg
462 : write(*,3) 'dm', j, k, dm
463 : write(*,3) 'm_00 + dm', j, k, m_00 + dm
464 : write(*,3) 'dm_m1', j, k, dm_m1
465 : write(*,3) 'dm/(m_m1 + m_p1)', j, k, dm/(m_m1 + m_p1)
466 : write(*,3) 'm_m1/(m_m1 + m_p1)', j, k, m_m1/(m_m1 + m_p1)
467 : write(*,3) 'm_m1 - dm_m1', j, k, m_m1 - dm_m1
468 : call mesa_error(__FILE__,__LINE__,'fixup 6')
469 : end if
470 : sum_m1 = sum(mass(1:nc,k-1)) ! sum_m1 - dm_m1
471 : sum_mass(k-1) = sum_m1
472 : if (sum_m1 < 0d0 .or. is_bad_num(sum_m1)) then
473 : write(*,2) 'sum_m1', j, k, sum_m1, dm_m1
474 : call mesa_error(__FILE__,__LINE__,'fixup')
475 : end if
476 :
477 : dm_p1 = dm*m_p1/(m_m1 + m_p1)
478 : mass(j,k+1) = m_p1 - dm_p1
479 : if (mass(j,k+1) < 0d0) then
480 : write(*,3) 'massp1', j, k, mass(j,k+1)
481 : write(*,3) 'massm1', j, k, mass(j,k-1)
482 : write(*,3) 'm_m1', j, k-1, m_m1
483 : write(*,3) 'm_00', j, k, m_00
484 : write(*,3) 'm_p1', j, k+1, m_p1
485 : write(*,3) 'xavg', j, k, xavg
486 : write(*,3) 'dm', j, k, dm
487 : write(*,3) 'm_00 + dm', j, k, m_00 + dm
488 : write(*,3) 'dm_p1', j, k, dm_m1
489 : write(*,3) 'dm/(m_m1 + m_p1)', j, k, dm/(m_m1 + m_p1)
490 : write(*,3) 'm_p1/(m_m1 + m_p1)', j, k, m_p1/(m_m1 + m_p1)
491 : write(*,3) 'm_p1 - dm_p1', j, k, m_p1 - dm_p1
492 : call mesa_error(__FILE__,__LINE__,'fixup 7')
493 : end if
494 : sum_p1 = sum(mass(1:nc,k+1)) ! sum_p1 - dm_p1
495 : sum_mass(k+1) = sum_p1
496 : if (sum_p1 < 0d0 .or. is_bad_num(sum_p1)) then
497 : write(*,2) 'sum_p1', j, k, sum_p1, dm_p1
498 : call mesa_error(__FILE__,__LINE__,'fixup')
499 : end if
500 :
501 : if (abs(xavg - sum(mass(j,k-1:k+1))/sum(sum_mass(k-1:k+1))) > 1d-12*xavg) then
502 : write(*,3) 'bad new xavg', j, k, xavg, &
503 : sum(mass(j,k-1:k+1))/sum(sum_mass(k-1:k+1))
504 : call mesa_error(__FILE__,__LINE__,'fixup 8')
505 : end if
506 :
507 : end do
508 :
509 : end do
510 :
511 : end do
512 :
513 : max_lnT = max_lnT - 0.1d0
514 :
515 : end do
516 :
517 : do k=nzlo,nzhi
518 : if (sum_mass(k) < 0d0) then
519 : write(*,2) 'sum_mass(k)', k, sum_mass(k)
520 : call mesa_error(__FILE__,__LINE__,'fixup 3')
521 : end if
522 : end do
523 :
524 : end subroutine smooth_x
525 :
526 :
527 0 : subroutine get1_dX_dm( &
528 0 : k, nz, nc, nzlo, nzhi, mass, sum_mass, &
529 0 : dX_dm, dbg, ierr)
530 : integer, intent(in) :: k, nz, nc, nzlo, nzhi
531 : real(dp), intent(in) :: sum_mass(:), mass(:,:)
532 : real(dp), intent(inout) :: dX_dm(:,:)
533 : logical :: dbg
534 : integer, intent(out) :: ierr
535 :
536 0 : real(dp) :: slope, sm1, s00, xface_00, xface_p1, &
537 0 : dm_half, dm_00, dm_p1, dm_m1, dmbar_p1, dmbar_00, &
538 0 : x00, xm1, xp1
539 : integer :: j
540 : real(dp), parameter :: tiny_slope = 1d-10
541 :
542 : include 'formats'
543 :
544 0 : ierr = 0
545 :
546 0 : dm_00 = sum_mass(k)
547 0 : if (dm_00 < tiny_mass) then
548 0 : dX_dm(1:nc,k) = 0d0
549 0 : return
550 : end if
551 0 : dm_half = 0.5d0*dm_00
552 :
553 0 : if (nzlo < k .and. k < nzhi) then
554 :
555 0 : dm_m1 = sum_mass(k-1)
556 0 : dm_p1 = sum_mass(k+1)
557 0 : if (dm_m1 < tiny_mass .or. dm_p1 < tiny_mass) then
558 0 : dX_dm(1:nc,k) = 0d0
559 : return
560 : end if
561 :
562 0 : dmbar_00 = 0.5d0*(dm_00 + dm_m1)
563 0 : dmbar_p1 = 0.5d0*(dm_00 + dm_p1)
564 :
565 0 : do j=1,nc
566 0 : xm1 = mass(j,k-1)/dm_m1
567 0 : x00 = mass(j,k)/dm_00
568 0 : xp1 = mass(j,k+1)/dm_p1
569 0 : sm1 = (xm1 - x00)/dmbar_00
570 0 : s00 = (x00 - xp1)/dmbar_p1
571 0 : slope = 0.5d0*(sm1 + s00)
572 0 : if (sm1*s00 <= 0 .or. abs(slope) < tiny_slope) then
573 0 : dX_dm(j,k) = 0d0
574 : else
575 0 : dX_dm(j,k) = slope
576 0 : xface_00 = x00 + slope*dm_half ! value at face(k)
577 0 : xface_p1 = x00 - slope*dm_half ! value at face(k+1)
578 : if (xface_00 > 1d0 .or. xface_00 < 0d0 .or. &
579 : (xm1 - xface_00)*(xface_00 - x00) < 0 .or. &
580 0 : xface_p1 > 1d0 .or. xface_p1 < 0d0 .or. &
581 : (x00 - xface_p1)*(xface_p1 - xp1) < 0) then
582 0 : if (abs(sm1) <= abs(s00)) then
583 0 : dX_dm(j,k) = sm1
584 : else
585 0 : dX_dm(j,k) = s00
586 : end if
587 : end if
588 : end if
589 : end do
590 :
591 0 : else if (k == nzlo) then
592 :
593 0 : dm_p1 = sum_mass(k+1)
594 0 : if (dm_p1 < tiny_mass) then
595 0 : dX_dm(1:nc,k) = 0d0
596 : return
597 : end if
598 :
599 0 : dmbar_p1 = 0.5d0*(dm_00 + dm_p1)
600 :
601 0 : do j=1,nc
602 0 : x00 = mass(j,k)/dm_00
603 0 : xp1 = mass(j,k+1)/dm_p1
604 0 : slope = (x00 - xp1)/dmbar_p1
605 0 : if (abs(slope) < tiny_slope) then
606 0 : dX_dm(j,k) = 0d0
607 : else
608 0 : dX_dm(j,k) = slope
609 0 : xface_00 = x00 + slope*dm_half ! value at face(k)
610 0 : xface_p1 = x00 - slope*dm_half ! value at face(k+1)
611 0 : if (xface_00 > 1d0 .or. xface_00 < 0d0 .or. &
612 : (x00 - xface_p1)*(xface_p1 - xp1) < 0) then
613 0 : dX_dm(j,k) = 0d0
614 : end if
615 : end if
616 : end do
617 :
618 0 : else if (k == nzhi) then
619 :
620 0 : dm_m1 = sum_mass(k-1)
621 0 : if (dm_m1 < tiny_mass) then
622 0 : dX_dm(1:nc,k) = 0d0
623 : return
624 : end if
625 :
626 0 : dmbar_00 = 0.5d0*(dm_00 + dm_m1)
627 :
628 0 : do j=1,nc
629 0 : xm1 = mass(j,k-1)/dm_m1
630 0 : x00 = mass(j,k)/dm_00
631 0 : slope = (xm1 - x00)/dmbar_00
632 0 : if (abs(slope) < tiny_slope) then
633 0 : dX_dm(j,k) = 0d0
634 : else
635 0 : dX_dm(j,k) = slope
636 0 : xface_00 = x00 + slope*dm_half ! value at face(k)
637 0 : xface_p1 = x00 - slope*dm_half ! value at face(k+1)
638 0 : if (xface_p1 > 1d0 .or. xface_p1 < 0d0 .or. &
639 : (xm1 - xface_00)*(xface_00 - x00) < 0) then
640 0 : dX_dm(j,k) = 0d0
641 : end if
642 : end if
643 : end do
644 :
645 : else
646 :
647 0 : write(*,2) 'k bad', k
648 0 : call mesa_error(__FILE__,__LINE__,'get1_dX_dm')
649 :
650 : end if
651 :
652 : ! adjust so that sum(dX_dm) = 0
653 0 : if (sum(dX_dm(1:nc,k)) > 0d0) then
654 0 : j = maxloc(dX_dm(1:nc,k),dim=1)
655 : else
656 0 : j = minloc(dX_dm(1:nc,k),dim=1)
657 : end if
658 0 : dX_dm(j,k) = 0d0 ! remove from sum
659 0 : dX_dm(j,k) = -sum(dX_dm(1:nc,k))
660 :
661 : ! recheck for valid values at faces
662 0 : do j=1,nc
663 0 : x00 = mass(j,k)/dm_00
664 0 : slope = dX_dm(j,k)
665 0 : xface_00 = x00 + slope*dm_half ! value at face(k)
666 0 : xface_p1 = x00 - slope*dm_half ! value at face(k+1)
667 : if (xface_00 > 1d0 .or. xface_00 < 0d0 .or. &
668 0 : xface_p1 > 1d0 .or. xface_p1 < 0d0) then
669 0 : if (dbg) then ! .and. abs(slope) > 1d-10) then
670 0 : write(*,3) 'give up on dX_dm', j, k
671 0 : write(*,1) 'slope', slope
672 0 : write(*,1) 'dm_half', dm_half
673 0 : write(*,1) 'xface_00', xface_00
674 0 : write(*,1) 'xface_p1', xface_p1
675 0 : if (k > nzlo) then
676 0 : dm_m1 = sum_mass(k-1)
677 0 : write(*,1) 'xm1', mass(j,k-1)/dm_m1
678 : end if
679 0 : write(*,1) 'x00', x00
680 0 : if (k < nzhi) then
681 0 : dm_p1 = sum_mass(k+1)
682 0 : write(*,1) 'xp1', mass(j,k+1)/dm_p1
683 : end if
684 0 : write(*,'(A)')
685 : end if
686 0 : dX_dm(1:nc,k) = 0d0
687 : exit
688 : end if
689 : end do
690 :
691 : end subroutine get1_dX_dm
692 :
693 :
694 0 : subroutine redistribute_mass( &
695 0 : s, nc, nzlo, nzhi, nz, total_num_iters, sum_mass, mass, dX_dm, &
696 0 : xtotal_init, cell_dm, X_new, dbg_in, ierr)
697 : type (star_info), pointer :: s
698 : integer, intent(in) :: nc, nzlo, nzhi, nz, total_num_iters
699 : real(dp), intent(inout), dimension(:) :: sum_mass
700 : real(dp), intent(in), dimension(:) :: cell_dm, xtotal_init
701 : real(dp), intent(in), dimension(:,:) :: mass, dX_dm
702 : real(dp), intent(inout) :: X_new(:,:)
703 : logical, intent(in) :: dbg_in
704 : integer, intent(out) :: ierr
705 :
706 : integer :: k_source, max_iters, k, i, j
707 0 : real(dp) :: source_cell_mass, remaining_source_mass, cell_dm_k, &
708 0 : remaining_needed_mass, sumX, total_source, total_moved, &
709 0 : dm0, dm1, old_sum, new_sum, mtotal, err, dm, diff_dm
710 : logical :: okay, dbg
711 :
712 : include 'formats'
713 :
714 0 : ierr = 0
715 0 : total_moved = 0d0
716 0 : dbg = dbg_in
717 :
718 : ! redistribute mass to make sum(X_new) = 1 for all cells
719 : ! this is done serially from nzlo to nzhi
720 0 : k_source = nzlo
721 0 : source_cell_mass = sum_mass(k_source)
722 0 : remaining_source_mass = source_cell_mass
723 0 : max_iters = 100
724 :
725 0 : if (dbg) write(*,2) 'nzlo', nzlo
726 0 : if (dbg) write(*,2) 'nzhi', nzhi
727 :
728 0 : cell_loop: do k=nzlo,nzhi
729 :
730 : !dbg = total_num_iters == 26 .and. k == 535
731 :
732 0 : cell_dm_k = cell_dm(k)
733 0 : remaining_needed_mass = cell_dm_k
734 0 : X_new(1:nc,k) = 0d0
735 0 : if (dbg) write(*,*)
736 :
737 0 : fill_loop: do i=1,max_iters
738 :
739 0 : if (dbg) write(*,4) 'remaining_needed_mass/cell_dm_k', &
740 0 : i, k, k_source, remaining_needed_mass/cell_dm_k, &
741 0 : remaining_needed_mass, cell_dm_k, remaining_source_mass
742 0 : if (is_bad_num(remaining_needed_mass)) call mesa_error(__FILE__,__LINE__,'redistribute_mass')
743 0 : if (is_bad_num(remaining_source_mass)) then
744 0 : write(*,4) 'remaining_source_mass', &
745 0 : i, k, k_source, remaining_source_mass, &
746 0 : remaining_needed_mass, cell_dm_k
747 0 : call mesa_error(__FILE__,__LINE__,'redistribute_mass')
748 : end if
749 :
750 0 : if (remaining_needed_mass <= remaining_source_mass) then
751 0 : if (dbg) write(*,1) 'remaining_needed_mass <= remaining_source_mass', &
752 0 : remaining_needed_mass, remaining_source_mass
753 0 : diff_dm = remaining_needed_mass
754 0 : dm0 = remaining_source_mass - diff_dm
755 0 : dm1 = remaining_source_mass
756 0 : if (dm0 < 0 .or. dm1 < 0) then
757 0 : write(*,2) 'dm0', k, dm0
758 0 : write(*,2) 'dm1', k, dm1
759 0 : write(*,2) 'diff_dm', k, diff_dm
760 0 : write(*,2) 'remaining_source_mass', k, remaining_source_mass
761 0 : write(*,2) 'remaining_needed_mass', k, remaining_needed_mass
762 0 : call mesa_error(__FILE__,__LINE__,'redistribute_mass')
763 : end if
764 0 : total_moved = total_moved + remaining_needed_mass
765 0 : if (dbg) then
766 0 : write(*,2) 'cell_dm_k', k, cell_dm_k
767 0 : write(*,2) 'remaining_needed_mass', k, remaining_needed_mass
768 : write(*,2) &
769 0 : 'cell_dm_k - remaining_needed_mass', k, cell_dm_k - remaining_needed_mass
770 0 : write(*,2) 'sum(X_new)*cell_dm_k', k, &
771 0 : sum(X_new(1:nc,k))*cell_dm_k
772 : end if
773 0 : do j=1,nc
774 0 : if (dbg) write(*,3) 'init X_new(j,k)', j, k, X_new(j,k)
775 0 : dm = integrate_mass(j,dm0,dm1,diff_dm)
776 0 : if (dm < 0d0) then
777 0 : write(*,3) 'dm', j, k, dm
778 0 : write(*,'(A)')
779 0 : call mesa_error(__FILE__,__LINE__,'redistribute_mass')
780 : end if
781 0 : X_new(j,k) = X_new(j,k) + dm/cell_dm_k
782 0 : if (dbg .and. (X_new(j,k) > 1d0 .or. X_new(j,k) < 0d0)) then
783 :
784 0 : write(*,3) 'bad X_new(j,k)', j, k, X_new(j,k)
785 0 : write(*,3) 'dm/cell_dm_k', j, k, dm/cell_dm_k
786 0 : write(*,3) 'dm', j, k, dm
787 0 : write(*,3) 'cell_dm_k', j, k, cell_dm_k
788 0 : write(*,3) 'remaining_needed_mass', j, k, remaining_needed_mass
789 0 : write(*,3) 'diff_dm', j, k, diff_dm
790 0 : write(*,3) 'dm0', j, k, dm0
791 0 : write(*,3) 'dm1', j, k, dm1
792 0 : write(*,3) 'remaining_source_mass', j, k, remaining_source_mass
793 0 : write(*,'(A)')
794 0 : call mesa_error(__FILE__,__LINE__,'redistribute_mass')
795 : end if
796 0 : remaining_needed_mass = remaining_needed_mass - dm
797 : end do
798 0 : if (dbg) write(*,1) 'final remaining_needed_mass/cell_dm_k', &
799 0 : remaining_needed_mass/cell_dm_k, remaining_needed_mass, cell_dm_k
800 0 : remaining_source_mass = dm0
801 0 : remaining_needed_mass = 0d0
802 0 : exit fill_loop
803 : end if
804 :
805 0 : if (dbg) write(*,1) 'use all remaining source cell mass'
806 : ! use all remaining source cell mass
807 0 : diff_dm = remaining_source_mass
808 0 : dm0 = 0d0
809 0 : dm1 = diff_dm
810 0 : do j=1,nc
811 0 : dm = integrate_mass(j,dm0,dm1,diff_dm)
812 0 : X_new(j,k) = X_new(j,k) + dm/cell_dm_k
813 0 : if (dbg .and. X_new(j,k) > 1d0 .or. X_new(j,k) < 0d0) then
814 0 : write(*,3) 'bad X_new(j,k)', j, k, X_new(j,k)
815 0 : write(*,1) 'dm', dm
816 0 : write(*,1) 'dm0', dm0
817 0 : write(*,1) 'dm1', dm1
818 0 : write(*,1) 'remaining_source_mass', remaining_source_mass
819 0 : write(*,1) 'remaining_needed_mass', remaining_needed_mass
820 0 : write(*,1) 'cell_dm_k', cell_dm_k
821 0 : call mesa_error(__FILE__,__LINE__,'redistribute_mass')
822 : end if
823 : end do
824 0 : if (is_bad_num(remaining_needed_mass)) then
825 0 : write(*,4) 'remaining_needed_mass', i, k, k_source, remaining_needed_mass
826 0 : call mesa_error(__FILE__,__LINE__,'redistribute_mass')
827 : end if
828 0 : total_moved = total_moved + remaining_source_mass
829 0 : remaining_needed_mass = remaining_needed_mass - remaining_source_mass
830 0 : if (remaining_needed_mass > cell_dm_k) then
831 0 : write(*,2) 'remaining_source_mass', k, remaining_source_mass
832 0 : write(*,2) 'remaining_needed_mass', k, remaining_needed_mass
833 0 : write(*,2) 'cell_dm_k', k, cell_dm_k
834 0 : call mesa_error(__FILE__,__LINE__,'redistribute_mass')
835 : end if
836 :
837 : ! go to next source cell
838 0 : k_source = k_source + 1 ! okay to allow k_source > nzhi; see integrate_mass
839 0 : source_cell_mass = sum_mass(min(nzhi,k_source))
840 0 : remaining_source_mass = source_cell_mass
841 :
842 0 : if (source_cell_mass < 0d0 .or. is_bad_num(source_cell_mass)) then
843 0 : ierr = -1
844 0 : s% retry_message = 'bad source cell mass in element diffusion'
845 0 : if (s% report_ierr) then
846 0 : !$OMP critical (diffusion_source_cell_mass)
847 0 : write(*,4) 'source_cell_mass', &
848 0 : i, k, k_source, remaining_source_mass, source_cell_mass
849 : !$OMP end critical (diffusion_source_cell_mass)
850 : end if
851 0 : return
852 :
853 : write(*,4) 'source_cell_mass', &
854 : i, k, k_source, remaining_source_mass, source_cell_mass
855 : call mesa_error(__FILE__,__LINE__,'redistribute_mass')
856 : end if
857 :
858 : end do fill_loop
859 0 : if (dbg) write(*,1) 'finished fill_loop'
860 :
861 0 : if (remaining_needed_mass > 0d0) then
862 0 : ierr = -1
863 0 : s% retry_message = 'bad remaining mass in element diffusion'
864 0 : if (s% report_ierr) then
865 0 : !$OMP critical (diffusion_dist_cell_mass)
866 0 : write(*,1) 'redistribute_mass: remaining_needed_mass', &
867 0 : remaining_needed_mass
868 : !$OMP end critical (diffusion_dist_cell_mass)
869 : end if
870 0 : return
871 :
872 : write(*,5) 'remaining_needed_mass > 0d0', k, k_source, nzhi, nz, &
873 : remaining_needed_mass
874 : write(*,1) 'source_cell_mass', source_cell_mass
875 : write(*,1) 'remaining_source_mass', remaining_source_mass
876 : call mesa_error(__FILE__,__LINE__,'redistribute_mass')
877 : end if
878 :
879 0 : sumX = sum(X_new(1:nc,k))
880 0 : if (abs(sumX - 1d0) > 1d-10) then
881 0 : write(*,1) 'sum(X(k)) - 1d0', sumX - 1d0
882 0 : write(*,1) 'sum mass/source_cell_mass', sum(mass(1:nc,k_source))/source_cell_mass
883 0 : write(*,1) 'sum dX_dm k_source', sum(dX_dm(1:nc,k_source))
884 0 : write(*,2) 'k', k
885 0 : write(*,2) 'k_source', k_source
886 0 : write(*,2) 'nzlo', nzlo
887 0 : write(*,2) 'nzhi', nzhi
888 0 : write(*,2) 'total_num_iters', total_num_iters
889 0 : call mesa_error(__FILE__,__LINE__,'redistribute_mass')
890 : end if
891 :
892 0 : do j=1,nc
893 0 : X_new(j,k) = X_new(j,k)/sumX
894 : end do
895 0 : sum_mass(k) = sum_mass(k)/sumX
896 :
897 : end do cell_loop
898 :
899 0 : if (dbg) write(*,1) 'finished cell_loop'
900 :
901 0 : total_source = sum(sum_mass(nzlo:nzhi))
902 0 : if (abs(total_moved/total_source - 1d0) > 1d-12) then
903 0 : write(*,1) 'total_moved/total_source - 1', total_moved/total_source - 1d0
904 0 : write(*,1) 'total_source', total_source
905 0 : write(*,1) 'total_moved', total_moved
906 0 : write(*,1) 'sum cell_dm', sum(cell_dm(nzlo:nzhi))
907 0 : write(*,2) 'k_source', k_source
908 0 : write(*,2) 'nzlo', nzlo
909 0 : write(*,2) 'nzhi', nzhi
910 0 : write(*,2) 'nz', nz
911 0 : write(*,'(A)')
912 0 : call mesa_error(__FILE__,__LINE__,'redistribute_mass')
913 : end if
914 :
915 : ! check cell sums
916 0 : okay = .true.
917 0 : do k=nzlo,nzhi
918 0 : new_sum = sum(X_new(1:nc,k))
919 0 : if (abs(new_sum - 1d0) > 1d-14 .or. is_bad_num(new_sum)) then
920 0 : write(*,2) 'redistribute_mass: new_sum-1', k, new_sum-1d0
921 0 : okay = .false.
922 : end if
923 : end do
924 0 : if (.not. okay) call mesa_error(__FILE__,__LINE__,'redistribute_mass')
925 :
926 : ! recheck conservation
927 0 : mtotal = sum(cell_dm(nzlo:nzhi))
928 0 : okay = .true.
929 0 : do j=1,nc
930 0 : if (xtotal_init(j) < 1d-20) cycle
931 0 : old_sum = mtotal*xtotal_init(j)
932 0 : new_sum = dot_product(cell_dm(nzlo:nzhi),X_new(j,nzlo:nzhi))
933 0 : err = (new_sum - old_sum)/max(old_sum,new_sum)
934 0 : if (abs(err) > 1d-12 .or. is_bad_num(err)) then
935 0 : write(*,2) 'redistribute_mass err', j, err, old_sum, new_sum
936 0 : okay = .false.
937 : end if
938 : end do
939 :
940 :
941 : contains
942 :
943 :
944 0 : real(dp) function integrate_mass(j,dm0,dm1,diff)
945 : integer, intent(in) :: j
946 : real(dp), intent(in) :: dm0, dm1, diff
947 :
948 0 : real(dp) :: dm, x, slope, x0, x1, half_dm, xavg
949 : integer :: k
950 :
951 : include 'formats'
952 :
953 0 : if (k_source > nzhi) then ! reuse last source cell
954 : k = nzhi
955 : slope = 0d0
956 : else
957 0 : k = k_source
958 0 : slope = dX_dm(j,k)
959 : end if
960 0 : dm = sum_mass(k)
961 0 : if (dm < tiny_mass) then
962 0 : integrate_mass = 0d0
963 : return
964 : end if
965 0 : x = mass(j,k)/dm
966 0 : half_dm = 0.5d0*dm
967 0 : x0 = x + slope*(dm0 - half_dm)
968 0 : x1 = x + slope*(dm1 - half_dm)
969 0 : if (dm0 > dm1) then
970 0 : write(*,1) 'x0', x0
971 0 : write(*,1) 'x1', x1
972 0 : write(*,1) 'dm0', dm0
973 0 : write(*,1) 'dm1', dm1
974 0 : call mesa_error(__FILE__,__LINE__,'integrate_mass')
975 : end if
976 0 : xavg = min(1d0, max(0d0, 0.5d0*(x0 + x1)))
977 0 : integrate_mass = xavg*diff
978 :
979 0 : end function integrate_mass
980 :
981 :
982 : end subroutine redistribute_mass
983 :
984 :
985 0 : subroutine get_limit_coeffs( &
986 : s, nz, nzlo, nzhi, &
987 : gamma_full_on, gamma_full_off, &
988 : T_full_on, T_full_off, &
989 0 : gamma, T, limit_coeffs_face, k_max)
990 :
991 : type (star_info), pointer :: s
992 : integer, intent(in) :: nz, nzlo, nzhi
993 : real(dp), intent(in) :: gamma_full_on, gamma_full_off, &
994 : T_full_on, T_full_off
995 : real(dp), dimension(:), intent(in) :: gamma, T
996 : real(dp), intent(inout) :: limit_coeffs_face(:)
997 : integer, intent(out) :: k_max
998 :
999 0 : real(dp) :: lim, gamma_term, T_term, gamma_max, T_min
1000 :
1001 : integer :: k
1002 :
1003 : include 'formats'
1004 :
1005 0 : k_max = nzhi
1006 :
1007 0 : do k=nzhi,nzlo+1,-1
1008 :
1009 0 : gamma_max = max(gamma(k),gamma(k-1))
1010 0 : if (gamma_max >= gamma_full_off) then
1011 : gamma_term = 0
1012 0 : else if (gamma_max <= gamma_full_on) then
1013 : gamma_term = 1
1014 : else
1015 0 : gamma_term = (gamma_full_off - gamma_max) / (gamma_full_off - gamma_full_on)
1016 : end if
1017 :
1018 0 : T_min = min(T(k),T(k-1))
1019 0 : if (T_min >= T_full_on) then
1020 : T_term = 1
1021 0 : else if (T_min <= T_full_off) then
1022 : T_term = 0
1023 : else
1024 0 : T_term = (T_full_off - T_min) / (T_full_off - T_full_on)
1025 : end if
1026 :
1027 0 : lim = gamma_term*T_term*(1d0 - s% phase(k))
1028 0 : if (lim <= 0d0) then
1029 0 : limit_coeffs_face(k) = 0d0
1030 0 : k_max = k-1
1031 0 : else if (lim >= 1d0) then
1032 0 : limit_coeffs_face(k) = 1d0
1033 : else
1034 0 : limit_coeffs_face(k) = 0.5d0*(1d0 - cospi(lim))
1035 : end if
1036 :
1037 : end do
1038 :
1039 0 : end subroutine get_limit_coeffs
1040 :
1041 :
1042 0 : subroutine setup_struct_info( &
1043 0 : s, nz, nzlo, nzhi, species, nc, m, X, A, tiny_X, &
1044 0 : dlnP_dm_face, dlnT_dm_face, dlnRho_dm_face, cell_dm, dm_in, &
1045 0 : abar, free_e, T, lnT, rho, lnd, L_face, r_face, alfa_face, &
1046 0 : class, class_chem_id, calculate_ionization, nsmooth_typical_charge, &
1047 : min_T_for_radaccel, max_T_for_radaccel, &
1048 : min_Z_for_radaccel, max_Z_for_radaccel, &
1049 0 : screening, rho_face, T_face, four_pi_r2_rho_face, &
1050 0 : dlnP_dr_face, dlnT_dr_face, dlnRho_dr_face, &
1051 0 : Z, typical_charge, xm_face, &
1052 0 : rad_accel_face, log10_g_rad, g_rad, &
1053 : kmax_rad_accel, ierr)
1054 :
1055 : type (star_info), pointer :: s
1056 : integer, intent(in) :: nz, nzlo, nzhi, species, nc, m, &
1057 : min_Z_for_radaccel, max_Z_for_radaccel
1058 : real(dp), intent(in) :: X(:,:), A(:)
1059 : real(dp), intent(in) :: tiny_X, &
1060 : min_T_for_radaccel, max_T_for_radaccel
1061 : real(dp), dimension(:), intent(in) :: &
1062 : dlnP_dm_face, dlnT_dm_face, dlnRho_dm_face, cell_dm, dm_in, &
1063 : abar, free_e, T, lnT, rho, lnd, L_face, r_face, alfa_face
1064 : integer, dimension(:), intent(in) :: class, class_chem_id
1065 : logical, intent(in) :: calculate_ionization, screening
1066 : integer, intent(in) :: nsmooth_typical_charge
1067 : real(dp), dimension(:), intent(out) :: &
1068 : rho_face, T_face, four_pi_r2_rho_face, &
1069 : dlnP_dr_face, dlnT_dr_face, dlnRho_dr_face
1070 : real(dp), dimension(:,:), intent(out) :: Z
1071 : real(dp), dimension(:,:), intent(out) :: &
1072 : typical_charge, rad_accel_face, log10_g_rad, g_rad
1073 : real(dp), dimension(:), intent(out) :: xm_face
1074 : integer, intent(out) :: kmax_rad_accel, ierr
1075 :
1076 : integer :: i, k, j, kmax
1077 :
1078 : logical, parameter :: dbg = .false.
1079 :
1080 : include 'formats'
1081 :
1082 0 : ierr = 0
1083 :
1084 : if (dbg) write(*,*) 'call do1 for each zone'
1085 :
1086 0 : if (calculate_ionization) then
1087 0 : !$OMP PARALLEL DO PRIVATE(k) SCHEDULE(dynamic,2)
1088 : do k=nzlo,nzhi
1089 : call do1(k)
1090 : end do
1091 : !$OMP END PARALLEL DO
1092 : else
1093 0 : do k=nzlo,nzhi
1094 0 : call do1(k)
1095 : end do
1096 : end if
1097 :
1098 0 : if (nsmooth_typical_charge > 0) then
1099 0 : do j=1,nsmooth_typical_charge
1100 0 : do i=1,nc
1101 : typical_charge(i,nzlo) = &
1102 0 : (2*typical_charge(i,nzlo) + typical_charge(i,nzlo+1))/3
1103 0 : do k = nzlo+1, nzhi-2
1104 : typical_charge(i,k) = &
1105 0 : (typical_charge(i,k-1) + typical_charge(i,k) + typical_charge(i,k+1))/3
1106 : end do
1107 : typical_charge(i,nzhi-1) = &
1108 0 : (2*typical_charge(i,nzhi-1) + typical_charge(i,nzhi-2))/3
1109 0 : do k = nzhi-2, nzlo+1, -1
1110 : typical_charge(i,k) = &
1111 0 : (typical_charge(i,k-1) + typical_charge(i,k) + typical_charge(i,k+1))/3
1112 : end do
1113 0 : typical_charge(i,nzhi) = typical_charge(i,nzhi-1)
1114 : end do
1115 : end do
1116 : end if
1117 :
1118 0 : xm_face(nzlo) = 0d0
1119 0 : do k=nzlo,nzhi
1120 0 : do i=1, nc
1121 0 : Z(i,k) = typical_charge(i,k)
1122 : end do
1123 0 : Z(m,k) = -1
1124 0 : if (k < nzhi) xm_face(k+1) = xm_face(k) + cell_dm(k)
1125 : end do
1126 :
1127 0 : if (T_face(nzlo+1) > max_T_for_radaccel) then
1128 0 : kmax_rad_accel = 0
1129 0 : return
1130 : end if
1131 :
1132 0 : kmax = nzhi
1133 0 : do k=nzlo+1,nzhi
1134 0 : if (T_face(k) > max_T_for_radaccel) then
1135 0 : kmax = k-1
1136 0 : exit
1137 : end if
1138 : end do
1139 0 : kmax_rad_accel = kmax
1140 :
1141 : if (dbg) write(*,*) 'call calc_g_rad'
1142 0 : if (s% op_mono_method == 'mombarg') then
1143 : call calc_g_rad_mombarg( &
1144 : s, nz, nzlo, nzhi, nc, m, kmax_rad_accel, X, A, &
1145 : class_chem_id, s% net_iso, s% op_mono_factors, &
1146 : L_face, rho_face, r_face, T_face, alfa_face, &
1147 : min_T_for_radaccel, max_T_for_radaccel, &
1148 : min_Z_for_radaccel, max_Z_for_radaccel, &
1149 : screening, log10_g_rad, g_rad, &
1150 0 : rad_accel_face, ierr)
1151 0 : else if (s% op_mono_method == 'hu') then
1152 : call calc_g_rad_hu( &
1153 : nz, nzlo, nzhi, nc, m, kmax_rad_accel, X, A, &
1154 : class_chem_id, s% net_iso, s% op_mono_factors, &
1155 : L_face, rho_face, r_face, T_face, alfa_face, &
1156 : min_T_for_radaccel, max_T_for_radaccel, &
1157 : min_Z_for_radaccel, max_Z_for_radaccel, &
1158 : screening, log10_g_rad, g_rad, &
1159 0 : rad_accel_face, ierr)
1160 : else
1161 0 : write(*,*) 'Invalid argument for op_mono_method.'
1162 0 : stop
1163 : end if
1164 :
1165 : if (dbg) write(*,*) 'done calc_g_rad'
1166 :
1167 : return
1168 : do k=nzlo,nzlo+9
1169 : write(*,2) 'alfa_face(k)', k, alfa_face(k)
1170 : write(*,2) 'rho_face(k)', k, rho_face(k)
1171 : write(*,2) 'r_face(k)', k, r_face(k)
1172 : write(*,2) 'four_pi_r2_rho_face(k)', k, four_pi_r2_rho_face(k)
1173 : write(*,2) 'dlnRho_dr_face(k)', k, dlnRho_dr_face(k)
1174 : write(*,2) 'dlnT_dr_face(k)', k, dlnT_dr_face(k)
1175 : end do
1176 :
1177 :
1178 : contains
1179 :
1180 :
1181 0 : subroutine do1(k)
1182 : use mod_typical_charge, only: eval_typical_charge
1183 : integer, intent(in) :: k
1184 : integer :: i
1185 :
1186 0 : if (k > nzlo) then
1187 0 : T_face(k) = alfa_face(k)*T(k) + (1d0-alfa_face(k))*T(k-1)
1188 0 : rho_face(k) = alfa_face(k)*rho(k) + (1d0-alfa_face(k))*rho(k-1)
1189 0 : four_pi_r2_rho_face(k) = pi4*r_face(k)*r_face(k)*rho_face(k)
1190 0 : dlnP_dr_face(k) = four_pi_r2_rho_face(k)*dlnP_dm_face(k)
1191 0 : dlnT_dr_face(k) = four_pi_r2_rho_face(k)*dlnT_dm_face(k)
1192 0 : dlnRho_dr_face(k) = four_pi_r2_rho_face(k)*dlnRho_dm_face(k)
1193 : end if
1194 :
1195 0 : if (calculate_ionization) then
1196 0 : do i=1, nc
1197 : typical_charge(i,k) = eval_typical_charge( &
1198 : class_chem_id(i), abar(k), free_e(k), &
1199 0 : T(k), lnT(k)/ln10, rho(k), lnd(k)/ln10)
1200 : end do
1201 : end if
1202 :
1203 0 : do i=1,nc
1204 0 : rad_accel_face(i,k) = 0
1205 0 : g_rad(i,k) = 0
1206 : end do
1207 0 : rad_accel_face(m,k) = 0
1208 :
1209 0 : end subroutine do1
1210 :
1211 : end subroutine setup_struct_info
1212 :
1213 :
1214 0 : subroutine calc_g_rad_mombarg( &
1215 0 : s, nz, nzlo, nzhi, nc, m, kmax_rad_accel, X, A, &
1216 0 : class_chem_id, net_iso, op_mono_factors, &
1217 0 : L_face, rho_face, r_face, T_face, alfa_face, &
1218 : min_T_for_radaccel, max_T_for_radaccel, &
1219 : min_Z_for_radaccel, max_Z_for_radaccel, &
1220 0 : screening, log10_g_rad, g_rad, &
1221 0 : rad_accel_face, ierr)
1222 :
1223 : use kap_lib, only: get_op_mono_params, call_load_op_master, call_compute_gamma_grid_mombarg
1224 : use utils_lib, only: utils_OMP_GET_MAX_THREADS, utils_OMP_GET_THREAD_NUM
1225 : type (star_info), pointer :: s
1226 :
1227 : integer, intent(in) :: &
1228 : nz, nzlo, nzhi, nc, m, kmax_rad_accel, class_chem_id(:), net_iso(:), &
1229 : min_Z_for_radaccel, max_Z_for_radaccel
1230 : real(dp), intent(in) :: &
1231 : min_T_for_radaccel, max_T_for_radaccel
1232 : real(dp), dimension(:), intent(in) :: &
1233 : A, L_face, rho_face, r_face, T_face, alfa_face, op_mono_factors
1234 : real(dp), dimension(:,:), intent(in) :: X
1235 : logical, intent(in) :: screening
1236 : real(dp), dimension(:,:), intent(out) :: &
1237 : log10_g_rad, g_rad, rad_accel_face
1238 : integer, intent(out) :: ierr
1239 :
1240 0 : integer :: iZ(nc), iZ_rad(nc), i, j, k, kk, kmax, op_err, k1, k2, dk
1241 0 : real(dp) :: alfa, beta, X_face(nc), blend_fac(-15:15)
1242 :
1243 :
1244 0 : real(dp) :: fk(17), delta, blend
1245 : character(len=4) :: e_name
1246 :
1247 : logical, parameter :: dbg = .false.
1248 :
1249 : include 'formats'
1250 :
1251 0 : blend_fac = [(1._dp - FLOAT(i)/31._dp, i=1,31)]
1252 0 : ierr = 0
1253 :
1254 0 : kmax = kmax_rad_accel
1255 :
1256 0 : g_rad(1:nc,kmax+1:nzhi) = 0
1257 0 : log10_g_rad(1:nc,kmax+1:nzhi) = 1
1258 :
1259 0 : g_rad(1:nc,nzlo) = 0
1260 0 : log10_g_rad(1:nc,nzlo) = 1
1261 :
1262 0 : iZ(1:nc) = chem_isos% Z(class_chem_id(1:nc))
1263 :
1264 0 : kk = 0
1265 0 : do i = 1, nc
1266 0 : if (iZ(i) >= min_Z_for_radaccel .and. &
1267 0 : iZ(i) <= max_Z_for_radaccel) then
1268 0 : kk = kk+1
1269 0 : iZ_rad(kk) = iZ(i)
1270 : end if
1271 : end do
1272 :
1273 0 : do j = 1, ngp
1274 0 : k1 = INT((kmax - (nzlo+1))/ ngp * (j-1) + (nzlo+1))
1275 0 : if (j == ngp) THEN
1276 : k2 = kmax
1277 : else
1278 0 : k2 = INT((kmax - (nzlo+1))/ ngp * j + (nzlo))
1279 : end if
1280 0 : fk = 0
1281 0 : do i=1, s% species
1282 0 : e_name = chem_isos% name(s% chem_id(i))
1283 0 : if (e_name == 'h1') fk(1) = SUM(s% xa(i,k1:k2))/ (chem_isos% W(s% chem_id(i)) * (k2-k1)) !s% xa(i,kmax)/ chem_isos% W(s% chem_id(i))
1284 0 : if (e_name == 'he4') fk(2) = SUM(s% xa(i,k1:k2))/ (chem_isos% W(s% chem_id(i)) * (k2-k1))
1285 0 : if (e_name == 'c12') fk(3) = SUM(s% xa(i,k1:k2))/ (chem_isos% W(s% chem_id(i)) * (k2-k1))
1286 0 : if (e_name == 'n14') fk(4) = SUM(s% xa(i,k1:k2))/ (chem_isos% W(s% chem_id(i)) * (k2-k1))
1287 0 : if (e_name == 'o16') fk(5) = SUM(s% xa(i,k1:k2))/ (chem_isos% W(s% chem_id(i)) * (k2-k1))
1288 0 : if (e_name == 'ne20')fk(6) = SUM(s% xa(i,k1:k2))/ (chem_isos% W(s% chem_id(i)) * (k2-k1))
1289 0 : if (e_name == 'na23')fk(7) = SUM(s% xa(i,k1:k2))/ (chem_isos% W(s% chem_id(i)) * (k2-k1))
1290 0 : if (e_name == 'mg24')fk(8) = SUM(s% xa(i,k1:k2))/ (chem_isos% W(s% chem_id(i)) * (k2-k1))
1291 0 : if (e_name == 'al27')fk(9) = SUM(s% xa(i,k1:k2))/ (chem_isos% W(s% chem_id(i)) * (k2-k1))
1292 0 : if (e_name == 'si28')fk(10) = SUM(s% xa(i,k1:k2))/ (chem_isos% W(s% chem_id(i)) * (k2-k1))
1293 0 : if (e_name == 's32') fk(11) = SUM(s% xa(i,k1:k2))/ (chem_isos% W(s% chem_id(i)) * (k2-k1))
1294 0 : if (e_name == 'ar40')fk(12) = SUM(s% xa(i,k1:k2))/ (chem_isos% W(s% chem_id(i)) * (k2-k1))
1295 0 : if (e_name == 'ca40')fk(13) = SUM(s% xa(i,k1:k2))/ (chem_isos% W(s% chem_id(i)) * (k2-k1))
1296 0 : if (e_name == 'cr52')fk(14) = SUM(s% xa(i,k1:k2))/ (chem_isos% W(s% chem_id(i)) * (k2-k1))
1297 0 : if (e_name == 'mn55')fk(15) = SUM(s% xa(i,k1:k2))/ (chem_isos% W(s% chem_id(i)) * (k2-k1))
1298 0 : if (e_name == 'fe56')fk(16) = SUM(s% xa(i,k1:k2))/ (chem_isos% W(s% chem_id(i)) * (k2-k1))
1299 0 : if (e_name == 'ni58')fk(17) = SUM(s% xa(i,k1:k2))/ (chem_isos% W(s% chem_id(i)) * (k2-k1))
1300 : end do
1301 0 : fk = fk / sum(fk)
1302 0 : delta = MAXVAL(ABS(fk - fk_gam_old(j,:))/fk_gam_old(j,:))
1303 0 : if (initialize_gamma_grid) then
1304 0 : call call_load_op_master(s% emesh_data_for_op_mono_path, ierr)
1305 :
1306 0 : write(*,*) 'Precompute gamma grid initialize.'
1307 0 : call call_compute_gamma_grid_mombarg(j, fk, ierr)
1308 : !write(*,*) 'Done precomputing gamma grid initialize.'
1309 0 : fk_gam_old(j,:) = fk
1310 0 : if (j == ngp) initialize_gamma_grid = .false.
1311 0 : else if (delta > 1d-4) then
1312 0 : write(*,*) 'Precompute gamma grid.'
1313 0 : call call_compute_gamma_grid_mombarg(j, fk, ierr)
1314 : !write(*,*) 'Done precomputing gamma grid.'
1315 0 : fk_gam_old(j,:) = fk
1316 : end if
1317 : end do
1318 :
1319 : !! PARALLEL DO PRIVATE(i,k,thread_num,op_err,j,alfa,beta,X_face,umesh,ff,ta,rs,sz,offset) SCHEDULE(guided)
1320 0 : !$OMP PARALLEL DO PRIVATE(i,k,op_err,j,alfa,beta,X_face,blend,dk) SCHEDULE(guided)
1321 : do k = nzlo+1, kmax
1322 : !if (k > nzlo .and. k <= kmax) then
1323 : if (T_face(k) < min_T_for_radaccel) then
1324 : do j = 1, nc
1325 : rad_accel_face(j,k) = 0d0
1326 : end do
1327 : rad_accel_face(m,k) = 0d0
1328 : cycle
1329 : end if
1330 : i = k - nzlo
1331 : alfa = alfa_face(k)
1332 : beta = 1d0 - alfa
1333 : do j = 1, nc
1334 : X_face(j) = alfa*X(j,k) + beta*X(j,k-1)
1335 : end do
1336 : op_err = 0
1337 :
1338 : if (dbg) write(*,2) 'call set1_g_rad', k
1339 : if (dbg) stop
1340 :
1341 : dk = k - INT((kmax - (nzlo+1))/ ngp + (nzlo+1))
1342 : if (ABS(dk) <= 15) then
1343 : j = -1
1344 : blend = blend_fac(dk)
1345 : else if (k < INT((kmax - (nzlo+1))/ ngp + (nzlo+1))) then
1346 : j = 1
1347 : blend = 0d0
1348 : else
1349 : j = 2
1350 : blend = 0d0
1351 : end if
1352 :
1353 : call set1_g_rad_mombarg( &
1354 : s, k, nc, j, blend, iZ, kk, T_face(k), rho_face(k), &
1355 : L_face(k), r_face(k), A, X_face, X, &
1356 : log10_g_rad(:,k), g_rad(:,k), &
1357 : op_err)
1358 :
1359 :
1360 : if (op_err /= 0) ierr = op_err
1361 : do j = 1, nc
1362 : rad_accel_face(j,k) = g_rad(j,k)
1363 : end do
1364 : rad_accel_face(m,k) = 0d0
1365 :
1366 : end do
1367 : !$OMP END PARALLEL DO
1368 :
1369 0 : k = nzlo
1370 0 : do j = 1, m
1371 0 : rad_accel_face(j,k) = rad_accel_face(j,k+1) ! for plotting
1372 : end do
1373 : !write(*,*) 'grad done'
1374 0 : end subroutine calc_g_rad_mombarg
1375 :
1376 :
1377 0 : subroutine set1_g_rad_mombarg( &
1378 0 : s, k, nc, j, blend, iZ, kk, T, rho, L, r, A, X, X_c,&
1379 0 : log10_g_rad, g_rad, ierr)
1380 0 : use kap_lib, only : call_compute_grad_mombarg
1381 : use chem_def
1382 :
1383 : type (star_info), pointer :: s
1384 :
1385 : integer, intent(in) :: k, nc, kk, j
1386 : integer, dimension(:), intent(in) :: iZ(:) ! (nc)
1387 : real(dp), intent(in) :: T, rho, L, r, blend
1388 : real(dp), dimension(:), intent(in) :: X, A
1389 : real(dp), dimension(:,:), intent(in) :: X_c
1390 :
1391 : real(dp), dimension(:), intent(out) :: &
1392 : log10_g_rad, g_rad
1393 0 : real(dp) :: logKappa
1394 : integer, intent(out) :: ierr
1395 :
1396 : integer :: i, ii
1397 0 : real(dp) :: logT, logRho
1398 0 : real(dp), dimension(17) ::lgrad, fk
1399 : integer :: iZ_rad2(17)
1400 : character(len=4) :: e_name
1401 : logical, parameter :: dbg = .false.
1402 :
1403 : include 'formats'
1404 :
1405 0 : ierr = 0
1406 0 : iZ_rad2 = [1,2,6,7,8,10,11,12,13,14,16,18,20,24,25,26,28]
1407 :
1408 : if (dbg) write(*,*) 'call op_mono_get_radacc'
1409 :
1410 :
1411 0 : fk = 0
1412 0 : do i=1, s% species
1413 0 : e_name = chem_isos% name(s% chem_id(i))
1414 0 : if (e_name == 'h1') fk(1) = s% xa(i,k)/ chem_isos% W(s% chem_id(i))
1415 0 : if (e_name == 'he4') fk(2) = s% xa(i,k)/ chem_isos% W(s% chem_id(i))
1416 0 : if (e_name == 'c12') fk(3) = s% xa(i,k)/ chem_isos% W(s% chem_id(i))
1417 0 : if (e_name == 'n14') fk(4) = s% xa(i,k)/ chem_isos% W(s% chem_id(i))
1418 0 : if (e_name == 'o16') fk(5) = s% xa(i,k)/ chem_isos% W(s% chem_id(i))
1419 0 : if (e_name == 'ne20')fk(6) = s% xa(i,k)/ chem_isos% W(s% chem_id(i))
1420 0 : if (e_name == 'na23')fk(7) = s% xa(i,k)/ chem_isos% W(s% chem_id(i))
1421 0 : if (e_name == 'mg24')fk(8) = s% xa(i,k)/ chem_isos% W(s% chem_id(i))
1422 0 : if (e_name == 'al27')fk(9) = s% xa(i,k)/ chem_isos% W(s% chem_id(i))
1423 0 : if (e_name == 'si28')fk(10) = s% xa(i,k)/ chem_isos% W(s% chem_id(i))
1424 0 : if (e_name == 's32') fk(11) = s% xa(i,k)/ chem_isos% W(s% chem_id(i))
1425 0 : if (e_name == 'ar40')fk(12) = s% xa(i,k)/ chem_isos% W(s% chem_id(i))
1426 0 : if (e_name == 'ca40')fk(13) = s% xa(i,k)/ chem_isos% W(s% chem_id(i))
1427 0 : if (e_name == 'cr52')fk(14) = s% xa(i,k)/ chem_isos% W(s% chem_id(i))
1428 0 : if (e_name == 'mn55')fk(15) = s% xa(i,k)/ chem_isos% W(s% chem_id(i))
1429 0 : if (e_name == 'fe56')fk(16) = s% xa(i,k)/ chem_isos% W(s% chem_id(i))
1430 0 : if (e_name == 'ni58')fk(17) = s% xa(i,k)/ chem_isos% W(s% chem_id(i))
1431 : end do
1432 0 : fk = fk / sum(fk)
1433 : call call_compute_grad_mombarg(k, j, blend, fk, T, rho,&
1434 : L, r,&
1435 0 : logKappa, lgrad, ierr)
1436 :
1437 : if (dbg) write(*,*) 'done op_mono_get_radacc'
1438 :
1439 0 : do i=1,nc
1440 0 : g_rad(i) = 0d0
1441 0 : log10_g_rad(i) = 1d0
1442 : end do
1443 :
1444 0 : if (ierr == 101) then ! logRho out of range
1445 0 : ierr = 0
1446 0 : write(*,2) 'op_mono_get_radacc bad logT logRho', k, logT, logRho
1447 : return
1448 : end if
1449 :
1450 0 : if (ierr /= 0) stop 'set1_g_rad' !return
1451 :
1452 0 : do ii = 1, 17 !kk
1453 0 : do i = 1, nc
1454 0 : if (iZ(i) == iZ_rad2(ii)) then
1455 0 : log10_g_rad(i) = lgrad(ii)
1456 0 : g_rad(i) = exp10(lgrad(ii))
1457 0 : exit
1458 : end if
1459 : end do
1460 : end do
1461 :
1462 0 : end subroutine set1_g_rad_mombarg
1463 :
1464 0 : subroutine calc_g_rad_hu( &
1465 0 : nz, nzlo, nzhi, nc, m, kmax_rad_accel, X, A, &
1466 0 : class_chem_id, net_iso, op_mono_factors, &
1467 0 : L_face, rho_face, r_face, T_face, alfa_face, &
1468 : min_T_for_radaccel, max_T_for_radaccel, &
1469 : min_Z_for_radaccel, max_Z_for_radaccel, &
1470 0 : screening, log10_g_rad, g_rad, &
1471 0 : rad_accel_face, ierr)
1472 :
1473 0 : use kap_lib, only: get_op_mono_params
1474 : use utils_lib, only: utils_OMP_GET_MAX_THREADS, utils_OMP_GET_THREAD_NUM
1475 :
1476 : integer, intent(in) :: &
1477 : nz, nzlo, nzhi, nc, m, kmax_rad_accel, class_chem_id(:), net_iso(:), &
1478 : min_Z_for_radaccel, max_Z_for_radaccel
1479 : real(dp), intent(in) :: &
1480 : min_T_for_radaccel, max_T_for_radaccel
1481 : real(dp), dimension(:), intent(in) :: &
1482 : A, L_face, rho_face, r_face, T_face, alfa_face, op_mono_factors
1483 : real(dp), dimension(:,:), intent(in) :: X
1484 : logical, intent(in) :: screening
1485 : real(dp), dimension(:,:), intent(out) :: &
1486 : log10_g_rad, g_rad, rad_accel_face
1487 : integer, intent(out) :: ierr
1488 :
1489 0 : integer :: iZ(nc), iZ_rad(nc), n, i, j, k, kk, kmax, op_err, sz, offset, &
1490 : nptot, ipe, nrad, thread_num
1491 0 : real(dp) :: alfa, beta, X_face(nc)
1492 0 : real, pointer, dimension(:) :: umesh1, semesh1, ff1, ta1, rs1
1493 : real, pointer :: &
1494 0 : umesh(:), semesh(:), ff(:,:,:,:), ta(:,:,:,:), rs(:,:,:)
1495 :
1496 : logical, parameter :: dbg = .false.
1497 :
1498 : include 'formats'
1499 :
1500 0 : ierr = 0
1501 0 : kmax = kmax_rad_accel
1502 :
1503 0 : g_rad(1:nc,kmax+1:nzhi) = 0
1504 0 : log10_g_rad(1:nc,kmax+1:nzhi) = 1
1505 :
1506 0 : g_rad(1:nc,nzlo) = 0
1507 0 : log10_g_rad(1:nc,nzlo) = 1
1508 :
1509 0 : iZ(1:nc) = chem_isos% Z(class_chem_id(1:nc))
1510 :
1511 0 : kk = 0
1512 0 : do i = 1, nc
1513 0 : if (iZ(i) >= min_Z_for_radaccel .and. &
1514 0 : iZ(i) <= max_Z_for_radaccel) then
1515 0 : kk = kk+1
1516 0 : iZ_rad(kk) = iZ(i)
1517 : end if
1518 : end do
1519 :
1520 : if (dbg) write(*,*) 'call get_op_mono_params'
1521 0 : call get_op_mono_params(nptot, ipe, nrad)
1522 : if (dbg) write(*,*) 'done get_op_mono_params'
1523 0 : n = utils_OMP_GET_MAX_THREADS()
1524 : if (dbg) write(*,2) 'nptot', nptot
1525 : if (dbg) write(*,2) 'ipe', ipe
1526 : if (dbg) write(*,2) 'nrad', nrad
1527 : if (dbg) write(*,2) 'kmax', kmax
1528 : if (dbg) write(*,2) 'nzlo', nzlo
1529 : if (dbg) write(*,2) 'n', n
1530 : if (dbg) write(*,2) 'nptot*ipe*4*4*n', nptot*ipe*4*4*n
1531 0 : if (nptot*ipe*4*4*n < 0) call mesa_error(__FILE__,__LINE__,'integer overflow for array size in calc_g_rad')
1532 : allocate( &
1533 : umesh1(nptot*n), semesh1(nptot*n), ff1(nptot*ipe*4*4*n), &
1534 : ta1(nptot*nrad*4*4*n), rs1(nptot*4*4*n), &
1535 0 : stat=ierr)
1536 0 : if (ierr /= 0) return
1537 :
1538 0 : !$OMP PARALLEL DO PRIVATE(i,k,thread_num,op_err,j,alfa,beta,X_face,umesh,ff,ta,rs,sz,offset) SCHEDULE(dynamic,2)
1539 : do k=nzlo+1, kmax
1540 : if (T_face(k) < min_T_for_radaccel) then
1541 : do j = 1, nc
1542 : rad_accel_face(j,k) = 0d0
1543 : end do
1544 : rad_accel_face(m,k) = 0d0
1545 : cycle
1546 : end if
1547 : i = k - nzlo
1548 : alfa = alfa_face(k)
1549 : beta = 1d0 - alfa
1550 : do j = 1, nc
1551 : X_face(j) = alfa*X(j,k) + beta*X(j,k-1)
1552 : end do
1553 : op_err = 0
1554 : thread_num = utils_OMP_GET_THREAD_NUM()
1555 : sz = nptot; offset = thread_num*sz
1556 : umesh(1:nptot) => umesh1(offset+1:offset+sz)
1557 : semesh(1:nptot) => semesh1(offset+1:offset+sz)
1558 : sz = nptot*ipe*4*4; offset = thread_num*sz
1559 : ff(1:nptot,1:ipe,1:4,1:4) => ff1(offset+1:offset+sz)
1560 : sz = nptot*nrad*4*4; offset = thread_num*sz
1561 : ta(1:nptot,1:nrad,1:4,1:4) => ta1(offset+1:offset+sz)
1562 : sz = nptot*4*4; offset = thread_num*sz
1563 : rs(1:nptot,1:4,1:4) => rs1(offset+1:offset+sz)
1564 : sz = nptot*nrad*4*4; offset = thread_num*sz
1565 :
1566 : if (dbg) write(*,2) 'call set1_g_rad', k
1567 : if (dbg) stop
1568 : call set1_g_rad_hu( &
1569 : k, nc, iZ, kk, iZ_rad, T_face(k), rho_face(k), &
1570 : L_face(k), r_face(k), A, X_face, screening, &
1571 : min_Z_for_radaccel, max_Z_for_radaccel, &
1572 : class_chem_id, net_iso, op_mono_factors, &
1573 : umesh, semesh, ff, ta, rs, &
1574 : log10_g_rad(:,k), g_rad(:,k), &
1575 : op_err)
1576 : if (op_err /= 0) ierr = op_err
1577 : do j = 1, nc
1578 : rad_accel_face(j,k) = g_rad(j,k)
1579 : end do
1580 : rad_accel_face(m,k) = 0d0
1581 : end do
1582 : !$OMP END PARALLEL DO
1583 :
1584 0 : deallocate(umesh1, semesh1, ff1, ta1, rs1)
1585 :
1586 0 : k = nzlo
1587 0 : do j = 1, m
1588 0 : rad_accel_face(j,k) = rad_accel_face(j,k+1) ! for plotting
1589 : end do
1590 :
1591 0 : end subroutine calc_g_rad_hu
1592 :
1593 :
1594 0 : subroutine set1_g_rad_hu( &
1595 0 : k, nc, iZ, kk, iZ_rad, T, rho, L, r, A, X, screening, &
1596 : min_Z_for_radaccel, max_Z_for_radaccel, &
1597 0 : class_chem_id, net_iso, op_mono_factors, &
1598 : umesh, semesh, ff, ta, rs, &
1599 0 : log10_g_rad, g_rad, ierr)
1600 :
1601 0 : use kap_lib, only: op_mono_get_radacc
1602 :
1603 : integer, intent(in) :: k, nc, kk, class_chem_id(:), net_iso(:), &
1604 : min_Z_for_radaccel, max_Z_for_radaccel
1605 : integer, dimension(:), intent(in) :: iZ(:) ! (nc)
1606 : integer, dimension(:), intent(in) :: iZ_rad(:) ! (kk)
1607 : real(dp), intent(in) :: T, rho, L, r
1608 : real(dp), dimension(:), intent(in) :: A, X, op_mono_factors
1609 : logical, intent(in) :: screening
1610 : real, pointer :: &
1611 : umesh(:), semesh(:), ff(:,:,:,:), ta(:,:,:,:), rs(:,:,:)
1612 : real(dp), dimension(:), intent(out) :: &
1613 : log10_g_rad, g_rad
1614 : integer, intent(out) :: ierr
1615 :
1616 : integer :: i, j, ii
1617 0 : real(dp) :: tot, logT, logRho, flux, logKappa
1618 : real(dp), dimension(nc) :: &
1619 0 : fa, fac, lgrad
1620 :
1621 : logical, parameter :: dbg = .false.
1622 :
1623 : include 'formats'
1624 :
1625 0 : ierr = 0
1626 :
1627 0 : tot = 0
1628 0 : fa = 0
1629 0 : do i=1,nc
1630 0 : fa(i) = X(i)/A(i)
1631 0 : tot = tot + fa(i)
1632 0 : j = net_iso(class_chem_id(i))
1633 0 : if (j /= 0) then
1634 0 : fac(i) = op_mono_factors(j)
1635 : else
1636 0 : write(*,*) 'bad class_chem_id? in set1_g_rad'
1637 0 : call mesa_error(__FILE__,__LINE__,'set1_g_rad')
1638 0 : fac(i) = 1
1639 : end if
1640 : end do
1641 0 : do i=1,nc
1642 0 : fa(i) = fa(i)/tot ! number fractions
1643 : end do
1644 :
1645 0 : logT = log10(T)
1646 0 : logRho = log10(rho)
1647 0 : flux = L/(pi4*r*r)
1648 :
1649 : if (dbg) write(*,*) 'call op_mono_get_radacc'
1650 : call op_mono_get_radacc( &
1651 : ! input
1652 : kk, iZ_rad, nc, iZ, fa, fac, &
1653 : flux, logT, logRho, screening, &
1654 : ! output
1655 : logKappa, &
1656 : lgrad, &
1657 : ! work arrays
1658 : umesh, semesh, ff, ta, rs, &
1659 0 : ierr)
1660 : if (dbg) write(*,*) 'done op_mono_get_radacc'
1661 :
1662 0 : do i=1,nc
1663 0 : g_rad(i) = 0d0
1664 0 : log10_g_rad(i) = 1d0
1665 : end do
1666 :
1667 0 : if (ierr == 101) then ! logRho out of range
1668 0 : ierr = 0
1669 0 : write(*,2) 'op_mono_get_radacc bad logT logRho', k, logT, logRho
1670 : return
1671 : end if
1672 :
1673 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'set1_g_rad') !return
1674 :
1675 0 : do ii = 1, kk
1676 0 : do i = 1, nc
1677 0 : if (iZ(i) == iZ_rad(ii)) then
1678 0 : log10_g_rad(i) = lgrad(ii)
1679 0 : g_rad(i) = exp10(lgrad(ii))
1680 0 : exit
1681 : end if
1682 : end do
1683 : end do
1684 :
1685 0 : end subroutine set1_g_rad_hu
1686 :
1687 0 : subroutine update_rad_accel_face( &
1688 : nzlo, nzhi, nc, m, A, X_init, X, &
1689 0 : log10_g_rad, g_rad, &
1690 0 : rad_accel_face, kmax_rad_accel)
1691 : integer, intent(in) :: nzlo, nzhi, nc, m, kmax_rad_accel
1692 : real(dp), intent(in) :: A(:) ! (nc)
1693 : real(dp), dimension(:,:), intent(in) :: &
1694 : X_init, X, log10_g_rad, g_rad
1695 : real(dp), dimension(:,:), intent(inout) :: rad_accel_face
1696 :
1697 : integer :: j, k, i
1698 :
1699 : include 'formats'
1700 :
1701 0 : do k=nzlo+1,kmax_rad_accel
1702 0 : do i=1,nc
1703 0 : rad_accel_face(i,k) = pow(10d0, log10_g_rad(i,k))
1704 : end do
1705 0 : rad_accel_face(m,k) = 0d0
1706 : end do
1707 :
1708 0 : k = nzlo
1709 0 : do j = 1, m
1710 0 : rad_accel_face(j,k) = rad_accel_face(j,k+1) ! for plotting
1711 : end do
1712 :
1713 0 : end subroutine update_rad_accel_face
1714 :
1715 :
1716 0 : subroutine set_new_xa( &
1717 0 : nz, nzlo, nzhi, species, nc, m, class, X_init, X, cell_dm, xa)
1718 : integer, intent(in) :: nz, nzlo, nzhi, species, nc, m, class(:)
1719 : real(dp), intent(in) :: X_init(:,:) ! (nc,nz)
1720 : real(dp), intent(in) :: X(:,:) ! (m,nz)
1721 : real(dp), intent(in) :: cell_dm(:) ! (nz)
1722 : real(dp), intent(inout) :: xa(:,:) ! (species,nz)
1723 : integer :: j, k, i
1724 0 : real(dp) :: tmp
1725 : include 'formats'
1726 0 : do k=nzlo,nzhi
1727 0 : do j=1,species
1728 0 : i = class(j)
1729 0 : if (X_init(i,k) <= 0 .or. is_bad_num(X_init(i,k))) then
1730 0 : write(*,3) 'X_init(i,k)', i, k, X_init(i,k)
1731 0 : call mesa_error(__FILE__,__LINE__,'set_new_xa')
1732 : end if
1733 0 : xa(j,k) = xa(j,k)*X(i,k)/X_init(i,k)
1734 : end do
1735 0 : tmp = sum(xa(1:species,k))
1736 0 : if (tmp <= 0 .or. is_bad_num(tmp)) then
1737 0 : write(*,2) 'tmp', k, tmp
1738 0 : call mesa_error(__FILE__,__LINE__,'set_new_xa')
1739 : end if
1740 0 : do j=1,species
1741 0 : xa(j,k) = xa(j,k)/tmp
1742 : end do
1743 : end do
1744 0 : end subroutine set_new_xa
1745 :
1746 : end module diffusion_procs
|
