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 mesh_adjust
21 :
22 : use const_def, only: dp, ln10, one_third, four_thirds_pi
23 : use star_private_def
24 : use chem_def
25 : use interp_1d_def, only: pm_work_size
26 : use utils_lib
27 : use star_utils, only: get_r_from_xh, store_r_in_xh, store_rho_in_xh, &
28 : get_lnT_from_xh, store_lnT_in_xh, store_lnd_in_xh
29 :
30 : implicit none
31 :
32 : private
33 : public :: do_mesh_adjust, do_prune_mesh_surface, &
34 : set_lnT_for_energy_with_tol, set_lnT_for_energy
35 :
36 : integer, parameter :: nwork = pm_work_size
37 :
38 : real(dp), parameter :: eta_limit = -1d-6
39 :
40 :
41 : logical, parameter :: dbg = .false.
42 :
43 : contains
44 :
45 10 : subroutine do_mesh_adjust( &
46 : s, nz, nz_old, xh_old, xa_old, &
47 : energy_old, eta_old, lnd_old, lnPgas_old, &
48 : j_rot_old, i_rot_old, omega_old, D_omega_old, &
49 : mlt_vc_old, lnT_old, w_old, specific_PE_old, specific_KE_old, &
50 : old_m, old_r, old_rho, dPdr_dRhodr_info_old, D_mix_old, &
51 10 : cell_type, comes_from, dq_old, xq_old, xh, xa, dq, xq, ierr)
52 : use chem_lib, only: basic_composition_info
53 : use interp_1d_def
54 : use interp_1d_lib
55 : use star_utils, only: set_m_grav_and_grav
56 : use auto_diff_support
57 : type (star_info), pointer :: s
58 : integer, intent(in) :: nz, nz_old
59 : integer, dimension(:) :: cell_type, comes_from
60 : real(dp), dimension(:), pointer :: &
61 : dq_old, xq_old, dq, xq, energy_old, eta_old, &
62 : lnd_old, lnPgas_old, mlt_vc_old, lnT_old, w_old, &
63 : specific_PE_old, specific_KE_old, &
64 : old_m, old_r, old_rho, dPdr_dRhodr_info_old, &
65 : j_rot_old, omega_old, D_omega_old, D_mix_old
66 : type(auto_diff_real_star_order1), dimension(:), pointer :: i_rot_old
67 : real(dp), dimension(:,:), pointer :: xh_old, xa_old
68 : real(dp), dimension(:,:), pointer :: xh, xa
69 : integer, intent(out) :: ierr
70 :
71 10 : real(dp) :: dxa, xmstar, mstar, sumx, &
72 10 : total_internal_energy1, total_internal_energy2, err
73 : character (len=strlen) :: message
74 : integer :: k, j, op_err, nzlo, nzhi, nzlo_old, nzhi_old, species
75 10 : real(dp), pointer :: work(:)
76 : real(dp), dimension(:), allocatable :: &
77 10 : dqbar, dqbar_old, new_r, Vol_new, xq_old_plus1, &
78 10 : xout_old, xout_new, xq_new, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, &
79 10 : energy_new, density_new
80 10 : real(dp), dimension(:,:), allocatable :: xa_c0, xa_c1, xa_c2
81 :
82 : include 'formats'
83 :
84 10 : ierr = 0
85 10 : species = s% species
86 :
87 10 : xmstar = s% xmstar
88 10 : mstar = xmstar + s% M_center
89 :
90 : ! check xq's
91 10633 : do k=1,nz
92 10633 : if (xq(k) < 0 .or. xq(k) > 1) then
93 0 : ierr = -1
94 0 : return
95 :
96 : write(*,*) 'k', k
97 : write(*,*) 'xq(k)', xq(k)
98 : call mesa_error(__FILE__,__LINE__,'debug: do_mesh_adjust')
99 : end if
100 : end do
101 :
102 : if (dbg) write(*,*) 'enter do_mesh_adjust'
103 :
104 10 : nzlo = 0
105 8446 : do k = 1, nz
106 8446 : if (cell_type(k) /= unchanged_type) then
107 : if (dbg) write(*,2) 'nzlo changed', k
108 4 : nzlo = k; exit
109 : end if
110 : end do
111 10 : if (nzlo == 0) then
112 : if (dbg) write(*,2) 'no cells changed'
113 6 : nzlo = nz
114 : end if
115 :
116 10 : nzhi = nzlo
117 252 : do k = nz, nzlo, -1
118 252 : if (cell_type(k) /= unchanged_type) then
119 : if (dbg) write(*,2) 'nzhi changed', k
120 4 : nzhi = k; exit
121 : end if
122 : end do
123 :
124 : ! extend range for purposes of interpolation
125 10 : if (nzhi < nz) nzhi = nzhi+1
126 10 : if (nzlo > 1) nzlo = nzlo-1
127 :
128 10 : nzlo_old = comes_from(nzlo)
129 10 : if (nzhi == nz) then
130 6 : nzhi_old = nz_old
131 : else
132 4 : nzhi_old = comes_from(nzhi+1)
133 : end if
134 :
135 10 : call do_alloc(ierr)
136 10 : if (ierr /= 0) return
137 :
138 10613 : do k=2,nz-1
139 10613 : dqbar(k) = 0.5d0*(dq(k-1) + dq(k))
140 : end do
141 10 : dqbar(1) = 0.5d0*dq(1)
142 10 : dqbar(nz) = 0.5d0*dq(nz-1) + dq(nz)
143 :
144 12085 : do k=2,nz_old-1
145 12085 : dqbar_old(k) = 0.5d0*(dq_old(k-1) + dq_old(k))
146 : end do
147 10 : dqbar_old(1) = 0.5d0*dq_old(1)
148 10 : dqbar_old(nz_old) = 0.5d0*dq_old(nz_old-1) + dq_old(nz_old)
149 :
150 12105 : do k=1,nz_old
151 12105 : xq_old_plus1(k) = xq_old(k)
152 : end do
153 : ! add point at true center so can interpolate xq_new > xq_old(nz_old)
154 10 : xq_old_plus1(nz_old+1) = 1
155 :
156 1981 : do k = 1, nzhi - nzlo + 1
157 1981 : xq_new(k) = xq(nzlo+k-1)
158 : end do
159 :
160 10 : xout_old(1) = xq_old(1)
161 12095 : do k=2,nz_old
162 12095 : xout_old(k) = xout_old(k-1) + dqbar_old(k-1)
163 : end do
164 :
165 10 : xout_new(1) = xq(1)
166 10623 : do k=2,nz
167 10623 : xout_new(k) = xout_new(k-1) + dqbar(k-1)
168 : end do
169 :
170 : if (dbg) write(*,*) 'call do_L'
171 : call do_L( &
172 : s, nz, nz_old, nzlo, nzhi, comes_from, &
173 : xh, xh_old, xq, xq_old_plus1, xq_new, &
174 10 : work, tmp1, tmp2, ierr)
175 0 : if (failed('do_L')) return
176 :
177 10 : if (s% RTI_flag) then
178 : if (dbg) write(*,*) 'call do_alpha_RTI'
179 : call do_alpha_RTI( &
180 : s, nz, nz_old, nzlo, nzhi, comes_from, &
181 : xh, xh_old, xq, xq_old_plus1, xq_new, &
182 0 : work, tmp1, tmp2, ierr)
183 0 : if (failed('do_alpha_RTI')) return
184 : call do_interp_pt_val( &
185 : s, nz, nz_old, nzlo, nzhi, s% dPdr_dRhodr_info, dPdr_dRhodr_info_old, &
186 0 : 0d0, xq, xq_old_plus1, xq_new, .true., work, tmp1, tmp2, ierr)
187 0 : if (failed('dPdr_dRhodr_info')) return
188 : end if
189 :
190 : if (dbg) write(*,*) 'call do_lnR_and_lnd'
191 : call do_lnR_and_lnd( &
192 : s, nz, nz_old, nzlo, nzhi, cell_type, comes_from, &
193 : xh, xh_old, xmstar, lnd_old, lnPgas_old, &
194 : dqbar, dqbar_old, old_r, old_m, old_rho, &
195 : dq, dq_old, xq, xq_old_plus1, density_new, work, &
196 10 : tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, ierr)
197 0 : if (failed('do_lnR_and_lnd')) return
198 :
199 10 : if (s% v_flag) then ! calculate new v to conserve kinetic energy
200 : if (dbg) write(*,*) 'call do_v'
201 : call do_v( &
202 : s, nz, nz_old, cell_type, comes_from, &
203 : xq_old, xq, dq_old, dq, xh, xh_old, &
204 0 : xout_old, xout_new, dqbar_old, dqbar, tmp1, ierr)
205 0 : if (failed('do_v')) return
206 : end if
207 :
208 10 : if (s% u_flag) then ! calculate new u to conserve kinetic energy
209 : if (dbg) write(*,*) 'call do_u'
210 : call do_u( &
211 : s, nz, nz_old, cell_type, comes_from, &
212 : xq_old, xq, dq_old, dq, xh, xh_old, &
213 0 : xout_old, xout_new, tmp1, ierr)
214 0 : if (failed('do_u')) return
215 : end if
216 :
217 10 : if (s% RSP2_flag) then ! calculate new etrb to conserve turbulent energy
218 : if (dbg) write(*,*) 'call do_etrb'
219 : call do_etrb( &
220 : s, nz, nz_old, cell_type, comes_from, &
221 : xq_old, xq, dq_old, dq, xh, xh_old, &
222 0 : xout_old, xout_new, tmp1, ierr)
223 0 : if (failed('do_etrb')) return
224 : if (dbg) write(*,*) 'call do_Hp_face'
225 : call do_Hp_face( &
226 : s, nz, nz_old, nzlo, nzhi, comes_from, &
227 : xh, xh_old, xq, xq_old_plus1, xq_new, &
228 0 : work, tmp1, tmp2, ierr)
229 0 : if (failed('do_Hp_face')) return
230 : end if
231 :
232 10 : if (s% rotation_flag) then
233 : call adjust_omega(s, nz, nz_old, comes_from, &
234 : xq_old, xq, dq_old, dq, xh, j_rot_old, &
235 0 : xout_old, xout_new, dqbar_old, dqbar, ierr)
236 0 : if (failed('adjust_omega')) return
237 0 : if (s% D_omega_flag) then
238 : call do_interp_pt_val( &
239 : s, nz, nz_old, nzlo, nzhi, s% D_omega, D_omega_old, &
240 0 : 0d0, xq, xq_old_plus1, xq_new, .true., work, tmp1, tmp2, ierr)
241 0 : if (failed('D_omega')) return
242 : end if
243 : end if
244 :
245 : call do_interp_pt_val( &
246 : s, nz, nz_old, nzlo, nzhi, s% mlt_vc, mlt_vc_old, &
247 20 : 0d0, xq, xq_old_plus1, xq_new, .true., work, tmp1, tmp2, ierr)
248 0 : if (failed('mlt_cv')) return
249 :
250 : call do_interp_pt_val( &
251 : s, nz, nz_old, nzlo, nzhi, s% D_mix, D_mix_old, &
252 20 : 0d0, xq, xq_old_plus1, xq_new, .true., work, tmp1, tmp2, ierr)
253 0 : if (failed('D_mix')) return
254 :
255 3457 : do k=nzlo_old,nzhi_old ! 1,nz_old !
256 :
257 : ! since we must adjust things to make the sum of xa's = 1,
258 : ! only do linear reconstruction.
259 31023 : do j=1,species
260 : call get1_lpp(k, species, nz_old, j, dq_old, xa_old, &
261 31023 : .false., xa_c0(:,j), xa_c1(:,j), xa_c2(:,j))
262 : end do
263 :
264 31023 : sumx = sum(xa_old(1:species,k))
265 31023 : do j=1,species
266 27576 : xa_c0(k,j) = xa_old(j,k)/sumx ! make sure that adds to 1
267 31023 : xa_c2(k,j) = 0 ! no curvature terms
268 : end do
269 :
270 : ! only reduce magnitude of slopes
271 : ! so don't risk producing values out of [0..1] range
272 31023 : if (sum(xa_c1(k,:)) > 0) then
273 14085 : j = maxloc(xa_c1(k,:), dim=1)
274 : else
275 16938 : j = minloc(xa_c1(k,:), dim=1)
276 : end if
277 3447 : xa_c1(k,j) = 0
278 31023 : xa_c1(k,j) = -sum(xa_c1(k,:))
279 : ! check for valid fractions at boundaries; set slopes to 0 if find a bad one.
280 31033 : do j=1,species
281 27576 : dxa = abs(xa_c1(k,j))*dq_old(k)/2
282 31023 : if (xa_c0(k,j) + dxa > 1 .or. xa_c0(k,j) - dxa < 0) then
283 0 : xa_c1(k,:) = 0
284 : exit
285 : end if
286 : end do
287 :
288 : end do
289 :
290 0 : if (failed('adjust_mesh nz_old parallel loop')) return
291 :
292 : if (dbg) write(*,*) 'do xa and lnT'
293 :
294 : total_internal_energy1 = &
295 12105 : dot_product(dq_old(1:nz_old), energy_old(1:nz_old))
296 :
297 10633 : do k = 1, nz
298 :
299 : op_err = 0
300 :
301 : ! calculate new abundances to conserve species
302 : call do_xa( &
303 : s, nz, nz_old, k, species, cell_type, comes_from, xa, xa_old, &
304 : xa_c0, xa_c1, xa_c2, xq, dq, xq_old, dq_old, &
305 10623 : .true., op_err)
306 10633 : if (op_err /= 0) then
307 0 : write(*,2) 'failed for do_xa', k
308 0 : stop
309 : write(message,*) 'do_xa for k', k
310 : ierr = op_err
311 : end if
312 :
313 : end do
314 :
315 : ! ensure the things we need to calculate m_grav and grav are up to date
316 10633 : do k = 1, nz
317 :
318 10623 : op_err = 0
319 :
320 : ! ensure that mass corrections are up to date
321 : call basic_composition_info( &
322 : species, s% chem_id, xa(1:species,k), s% X(k), s% Y(k), s% Z(k), &
323 : s% abar(k), s% zbar(k), s% z2bar(k), s% z53bar(k), &
324 10623 : s% ye(k), s% mass_correction(k), sumx)
325 :
326 : ! ensure that r is up to date
327 10633 : s% r(k) = exp(s% xh(s% i_lnr,k))
328 :
329 : end do
330 :
331 : ! needed because do1_lnT evaluates PE
332 10 : call set_m_grav_and_grav(s)
333 :
334 10633 : do k = 1, nz
335 :
336 : op_err = 0
337 :
338 : ! calculate new temperatures to conserve energy
339 : call do1_lnT( &
340 : s, nz_old, k, species, cell_type, comes_from, &
341 : xa, xh, xh_old, &
342 : xq, dq, xq_old, dq_old, eta_old, energy_old, lnT_old, &
343 : specific_PE_old, specific_KE_old, w_old, &
344 10623 : density_new, energy_new, op_err)
345 10623 : if (op_err /= 0) then
346 0 : write(*,2) 'failed for do1_lnT', k
347 0 : stop
348 : write(message,*) 'do1_lnT for k', k
349 : ierr = op_err
350 : end if
351 10633 : if (is_bad(energy_new(k)) .or. is_bad(dq(k))) then
352 0 : write(*,2) 'energy_new', k, energy_new(k)
353 0 : write(*,2) 'dq', k, dq(k)
354 0 : stop ''
355 : end if
356 :
357 : end do
358 :
359 0 : if (failed(message)) return
360 :
361 10633 : total_internal_energy2 = dot_product(dq(1:nz), energy_new(1:nz))
362 : err = abs(total_internal_energy1 - total_internal_energy2)/ &
363 10 : max(abs(total_internal_energy1),abs(total_internal_energy2),1d0)
364 10 : s% mesh_adjust_IE_conservation = err
365 :
366 10 : if (s% trace_mesh_adjust_error_in_conservation) then
367 :
368 0 : write(*,2) 'mesh adjust error in conservation of IE', &
369 0 : s% model_number, err, total_internal_energy2, total_internal_energy1
370 0 : if (err > 1d-8) then
371 0 : call show_errors
372 0 : write(*,*) 'err too large'
373 0 : call mesa_error(__FILE__,__LINE__,'mesh adjust')
374 : end if
375 :
376 : end if
377 :
378 : if (dbg) write(*,*) 'call check_species_conservation'
379 10 : call check_species_conservation(species,ierr)
380 10 : if (ierr /= 0) then
381 0 : write(*,*) 'failed in check_species_conservation'
382 0 : stop
383 : end if
384 :
385 10 : call dealloc
386 :
387 : contains
388 :
389 0 : subroutine show_errors
390 : integer :: k0, k0_from, k, k_from, k_outer
391 0 : real(dp) :: new_sum, old_sum
392 :
393 : include 'formats'
394 :
395 : k0 = 1
396 : k0_from = 1
397 0 : k_outer = 2
398 0 : write(*,'(A)')
399 : do
400 0 : if (cell_type(k_outer) == unchanged_type .and. k_outer /= nz) then
401 0 : k_outer = k_outer + 1
402 0 : cycle
403 : end if
404 0 : k0 = k_outer
405 0 : k0_from = comes_from(k_outer)
406 0 : do k = k0+1, nz
407 0 : if (cell_type(k) /= unchanged_type .and. k /= nz) cycle
408 0 : new_sum = dot_product(dq(k0:k),energy_new(k0:k))
409 0 : k_from = comes_from(k)
410 : old_sum = &
411 0 : dot_product(dq_old(k0_from:k_from),energy_old(k0_from:k_from))
412 0 : write(*,5) 'section err', k0, k, k0_from, k_from, &
413 0 : abs(old_sum - new_sum)/max(abs(old_sum),abs(new_sum),1d0), &
414 0 : old_sum, new_sum
415 0 : k_outer = k
416 0 : exit
417 : end do
418 0 : if (k_outer == nz) exit
419 0 : k_outer = k_outer + 1
420 : end do
421 0 : write(*,2) 'nz_old', nz_old
422 0 : write(*,2) 'nz', nz
423 0 : write(*,'(A)')
424 :
425 10 : end subroutine show_errors
426 :
427 10 : subroutine do_alloc(ierr)
428 : integer, intent(out) :: ierr
429 : integer :: sz
430 10 : sz = max(nz, nz_old) + 1
431 10 : call do_work_arrays(.true.,ierr)
432 : allocate( &
433 0 : dqbar(sz), dqbar_old(sz), new_r(sz), Vol_new(sz), xq_old_plus1(sz), &
434 0 : xout_old(sz), xout_new(sz), xq_new(sz), energy_new(sz), density_new(sz), &
435 0 : tmp1(sz), tmp2(sz), tmp3(sz), tmp4(sz), tmp5(sz), tmp6(sz), tmp7(sz), &
436 20 : xa_c0(sz,species), xa_c1(sz,species), xa_c2(sz,species))
437 10 : end subroutine do_alloc
438 :
439 10 : subroutine dealloc
440 10 : call do_work_arrays(.false.,ierr)
441 : end subroutine dealloc
442 :
443 20 : subroutine do_work_arrays(alloc_flag, ierr)
444 : use interp_1d_def
445 : use alloc
446 : logical, intent(in) :: alloc_flag
447 : integer, intent(out) :: ierr
448 : logical, parameter :: crit = .false.
449 : ierr = 0
450 : call work_array(s, alloc_flag, crit, &
451 20 : work, (nz_old+1)*pm_work_size, nz_alloc_extra, 'mesh_adjust', ierr)
452 20 : if (ierr /= 0) return
453 20 : end subroutine do_work_arrays
454 :
455 :
456 60 : logical function failed(msg)
457 : character (len=*) :: msg
458 60 : if (ierr == 0) then
459 60 : failed = .false.
460 : return
461 : end if
462 0 : failed = .true.
463 : if (dbg) write(*, *) 'mesh_revisions failed in ' // trim(msg)
464 : call dealloc
465 : return
466 20 : end function failed
467 :
468 :
469 10 : subroutine check_species_conservation(species,ierr)
470 : integer, intent(in) :: species
471 : integer, intent(out) :: ierr
472 : integer :: j, k, jbad
473 10 : real(dp) :: old_total, new_total
474 : logical :: okay
475 : include 'formats'
476 10 : ierr = 0
477 10 : okay = .true.
478 10 : jbad = -1
479 90 : do j=1,species
480 96840 : old_total = dot_product(xa_old(j,1:nz_old),dq_old(1:nz_old))
481 80 : if (old_total < 1d-9) cycle
482 85064 : new_total = dot_product(xa(j,1:nz),dq(1:nz))
483 90 : if (abs(new_total - old_total) > 1d-4) then ! check for major problems
484 0 : ierr = -1
485 0 : jbad = j
486 0 : okay = .false.
487 : if (dbg) then
488 : write(*,*) 'problem with conservation of species ' // &
489 : chem_isos% name(s% chem_id(j))
490 : write(*,1) 'new mass fraction', new_total
491 : write(*,1) 'old mass fraction', old_total
492 : write(*,1) 'new - old', new_total - old_total
493 : write(*,1) '(new - old)/old', (new_total - old_total) / old_total
494 : write(*,'(A)')
495 : end if
496 : end if
497 : end do
498 10 : if (okay) return
499 0 : ierr = -1
500 0 : write(*,'(A)')
501 0 : do j=1,species
502 0 : old_total = dot_product(xa_old(j,1:nz_old),dq_old(1:nz_old))
503 0 : if (old_total < 1d-9) cycle
504 0 : new_total = dot_product(xa(j,1:nz),dq(1:nz))
505 0 : write(*,2) 'new - old mass fraction ' // chem_isos% name(s% chem_id(j)), &
506 0 : j, new_total-old_total
507 : end do
508 0 : write(*,'(A)')
509 0 : j = jbad
510 0 : do k=2, nz
511 0 : if (comes_from(k) == comes_from(k-1)) cycle
512 0 : old_total = dot_product(xa_old(j,1:comes_from(k)-1),dq_old(1:comes_from(k)-1))
513 0 : if (old_total < 1d-9) cycle
514 0 : new_total = dot_product(xa(j,1:k-1),dq(1:k-1))
515 0 : write(*,2) 'partial new - old ' // chem_isos% name(s% chem_id(j)), k, &
516 0 : new_total-old_total, new_total, old_total
517 : end do
518 0 : write(*,'(A)')
519 0 : do k=415, nz
520 0 : write(*,'(a30,99i6)') 'cell_type(k)', k, cell_type(k), comes_from(k)
521 : end do
522 0 : write(*,'(A)')
523 0 : write(*,2) 'xq', 439, xq(439)
524 0 : write(*,2) 'xq_old', 429, xq_old(429)
525 0 : write(*,2) 'dq_old', 429, dq_old(429)
526 0 : write(*,2) 'dq', 439, dq(439)
527 0 : write(*,'(A)')
528 0 : write(*,2) 'xq', 424, xq(424)
529 0 : write(*,2) 'xq_old', 428, xq_old(428)
530 0 : write(*,2) 'dq_old', 428, dq_old(428)
531 0 : write(*,2) 'sum dq', 424, sum(dq(424:438))
532 0 : write(*,'(A)')
533 0 : write(*,2) 'xq_old + dq_old', 428, xq_old(428) + dq_old(428)
534 0 : write(*,2) 'xq_old', 429, xq_old(429)
535 0 : write(*,'(A)')
536 0 : write(*,2) 'xq + sum dq', 424, xq(424) + sum(dq(424:438))
537 0 : write(*,2) 'xq', 439, xq(439)
538 0 : write(*,'(A)')
539 0 : write(*,1) 'sum dq_old', sum(dq_old(1:nz_old))
540 :
541 0 : write(*,2) 'dq_old', 427, dq_old(427)
542 0 : write(*,2) 'sum new', 416, sum(dq(416:423))
543 0 : write(*,2) 'dq_old - sum new', 427, dq_old(427) - sum(dq(416:423))
544 0 : write(*,2) 'dq_old', 428, dq_old(428)
545 0 : write(*,2) 'sum new', 424, sum(dq(424:438))
546 0 : write(*,2) 'dq_old - sum new', 428, dq_old(428) - sum(dq(424:438))
547 : end subroutine check_species_conservation
548 :
549 :
550 : end subroutine do_mesh_adjust
551 :
552 :
553 0 : subroutine do_prune_mesh_surface( &
554 : s, nz, nz_old, xh_old, xa_old, &
555 : j_rot_old, i_rot_old, omega_old, D_omega_old, am_nu_rot_old, &
556 : mlt_vc_old, lnT_old, &
557 : dPdr_dRhodr_info_old, nu_ST_old, D_ST_old, D_DSI_old, D_SH_old, &
558 : D_SSI_old, D_ES_old, D_GSF_old, D_mix_old, &
559 : xh, xa, ierr)
560 : use auto_diff_support
561 : type (star_info), pointer :: s
562 : integer, intent(in) :: nz, nz_old
563 : real(dp), dimension(:), pointer :: &
564 : j_rot_old, omega_old, &
565 : D_omega_old, am_nu_rot_old, mlt_vc_old, lnT_old, &
566 : dPdr_dRhodr_info_old, nu_ST_old, D_ST_old, D_DSI_old, D_SH_old, &
567 : D_SSI_old, D_ES_old, D_GSF_old, D_mix_old
568 : type(auto_diff_real_star_order1), dimension(:), pointer :: i_rot_old
569 : real(dp), dimension(:,:), pointer :: xh_old, xa_old
570 : real(dp), dimension(:,:), pointer :: xh, xa
571 : integer, intent(out) :: ierr
572 :
573 : integer :: skip, k, k_old, j
574 :
575 : include 'formats'
576 :
577 0 : ierr = 0
578 0 : skip = nz_old - nz
579 0 : if (skip < 0) then
580 0 : write(*,3) 'bad nz prune_surface do_prune_mesh_surface', nz_old, nz
581 0 : ierr = -1
582 : return
583 : end if
584 :
585 0 : do k = 1, nz
586 0 : k_old = k + skip
587 0 : do j = 1, s% nvar_hydro
588 0 : xh(j,k) = xh_old(j,k_old)
589 : end do
590 0 : do j = 1, s% species
591 0 : xa(j,k) = xa_old(j,k_old)
592 : end do
593 : end do
594 :
595 0 : call prune1(s% lnT, lnT_old, skip)
596 0 : call prune1(s% D_mix, D_mix_old, skip)
597 :
598 0 : if (s% rotation_flag) then
599 0 : call prune1(s% j_rot, j_rot_old, skip)
600 0 : call prune1_ad(s% i_rot, i_rot_old, skip)
601 0 : call prune1(s% omega, omega_old, skip)
602 0 : call prune1(s% nu_ST, nu_ST_old, skip)
603 0 : call prune1(s% D_ST, D_ST_old, skip)
604 0 : call prune1(s% D_DSI, D_DSI_old, skip)
605 0 : call prune1(s% D_SH, D_SH_old, skip)
606 0 : call prune1(s% D_SSI, D_SSI_old, skip)
607 0 : call prune1(s% D_ES, D_ES_old, skip)
608 0 : call prune1(s% D_GSF, D_GSF_old, skip)
609 : end if
610 :
611 0 : if (s% D_omega_flag) then
612 0 : call prune1(s% D_omega, D_omega_old, skip)
613 : end if
614 :
615 0 : if (s% RTI_flag) then
616 0 : call prune1(s% dPdr_dRhodr_info, dPdr_dRhodr_info_old, skip)
617 : end if
618 :
619 : contains
620 :
621 0 : subroutine prune1(p,p_old,skip)
622 : real(dp), dimension(:), pointer :: p, p_old
623 : integer, intent(in) :: skip
624 : integer :: k
625 0 : do k=1,nz
626 0 : p(k) = p_old(k+skip)
627 : end do
628 0 : end subroutine prune1
629 :
630 0 : subroutine prune1_ad(p,p_old,skip)
631 : type(auto_diff_real_star_order1), dimension(:), pointer :: p, p_old
632 : integer, intent(in) :: skip
633 : integer :: k
634 0 : do k=1,nz
635 0 : p(k) = p_old(k+skip)
636 : end do
637 0 : end subroutine prune1_ad
638 :
639 : end subroutine do_prune_mesh_surface
640 :
641 :
642 : subroutine do_xh_pt_var( &
643 : s, i_var, nz, nz_old, nzlo, nzhi, comes_from, xh, xh_old, center_val, &
644 : xq, xq_old_plus1, xq_new, work, var_old_plus1, var_new, ierr)
645 : use interp_1d_def
646 : use interp_1d_lib
647 : type (star_info), pointer :: s
648 : integer, intent(in) :: i_var, nz, nz_old, nzlo, nzhi, comes_from(:)
649 : real(dp),intent(in) :: center_val
650 : real(dp), dimension(:,:), pointer :: xh, xh_old
651 : real(dp), dimension(:), pointer :: work
652 : real(dp), dimension(:) :: &
653 : xq, xq_old_plus1, var_old_plus1, var_new, xq_new
654 : integer, intent(out) :: ierr
655 :
656 : integer :: n, k
657 :
658 : include 'formats'
659 :
660 : ierr = 0
661 : n = nzhi - nzlo + 1
662 :
663 : do k=1,nz_old
664 : var_old_plus1(k) = xh_old(i_var,k)
665 : end do
666 : var_old_plus1(nz_old+1) = center_val
667 :
668 : call interpolate_vector( &
669 : nz_old+1, xq_old_plus1, n, xq_new, &
670 : var_old_plus1, var_new, interp_pm, nwork, work, &
671 : 'mesh_adjust do_xh_pt_var', ierr)
672 : if (ierr /= 0) then
673 : return
674 : end if
675 :
676 : do k=nzlo,nzhi
677 : xh(i_var,k) = max(0d0,var_new(k+1-nzlo))
678 : end do
679 :
680 : n = nzlo - 1
681 : if (n > 0) then
682 : do k=1,n
683 : xh(i_var,k) = xh_old(i_var,k)
684 : end do
685 : end if
686 :
687 : if (nzhi < nz) then
688 : n = nz - nzhi - 1 ! nz-n = nzhi+1
689 : do k=0,n
690 : xh(i_var,nz-k) = xh_old(i_var,nz_old-k)
691 : end do
692 : end if
693 :
694 : end subroutine do_xh_pt_var
695 :
696 :
697 10 : subroutine do_L( &
698 : s, nz, nz_old, nzlo, nzhi, comes_from, xh, xh_old, &
699 10 : xq, xq_old_plus1, xq_new, work, L_old_plus1, L_new, ierr)
700 : use interp_1d_def
701 : use interp_1d_lib
702 : type (star_info), pointer :: s
703 : integer, intent(in) :: nz, nz_old, nzlo, nzhi, comes_from(:)
704 : real(dp), dimension(:,:), pointer :: xh, xh_old
705 : real(dp), dimension(:), pointer :: work
706 : real(dp), dimension(:) :: &
707 : xq, xq_old_plus1, L_old_plus1, L_new, xq_new
708 : integer, intent(out) :: ierr
709 :
710 : integer :: n, i_lum, k
711 :
712 : include 'formats'
713 :
714 10 : ierr = 0
715 10 : i_lum = s% i_lum
716 10 : if (i_lum == 0) return
717 10 : n = nzhi - nzlo + 1
718 :
719 12105 : do k=1,nz_old
720 12105 : L_old_plus1(k) = xh_old(i_lum,k)
721 : end do
722 10 : L_old_plus1(nz_old+1) = s% L_center
723 :
724 : call interpolate_vector( &
725 : nz_old+1, xq_old_plus1, n, xq_new, &
726 : L_old_plus1, L_new, interp_pm, nwork, work, &
727 10 : 'mesh_adjust do_L', ierr)
728 10 : if (ierr /= 0) then
729 : return
730 :
731 : write(*,*) 'interpolate_vector failed in do_L for remesh'
732 : call mesa_error(__FILE__,__LINE__,'debug: mesh adjust: do_L')
733 : end if
734 :
735 1981 : do k=nzlo,nzhi
736 1981 : xh(i_lum,k) = L_new(k+1-nzlo)
737 : end do
738 :
739 10 : n = nzlo - 1
740 10 : if (n > 0) then
741 8430 : do k=1,n
742 8430 : xh(i_lum,k) = xh_old(i_lum,k)
743 : end do
744 : end if
745 :
746 10 : if (nzhi < nz) then
747 4 : n = nz - nzhi - 1 ! nz-n = nzhi+1
748 236 : do k=0,n
749 236 : xh(i_lum,nz-k) = xh_old(i_lum,nz_old-k)
750 : end do
751 : end if
752 :
753 10 : end subroutine do_L
754 :
755 :
756 0 : subroutine do_alpha_RTI( &
757 : s, nz, nz_old, nzlo, nzhi, comes_from, xh, xh_old, &
758 0 : xq, xq_old_plus1, xq_new, work, alpha_RTI_old_plus1, alpha_RTI_new, ierr)
759 10 : use interp_1d_def
760 : use interp_1d_lib
761 : type (star_info), pointer :: s
762 : integer, intent(in) :: nz, nz_old, nzlo, nzhi, comes_from(:)
763 : real(dp), dimension(:,:), pointer :: xh, xh_old
764 : real(dp), dimension(:), pointer :: work
765 : real(dp), dimension(:) :: &
766 : xq, xq_old_plus1, alpha_RTI_old_plus1, alpha_RTI_new, xq_new
767 : integer, intent(out) :: ierr
768 :
769 : integer :: n, i_alpha_RTI, k
770 :
771 : include 'formats'
772 :
773 : ierr = 0
774 0 : i_alpha_RTI = s% i_alpha_RTI
775 0 : n = nzhi - nzlo + 1
776 :
777 0 : do k=1,nz_old
778 0 : alpha_RTI_old_plus1(k) = xh_old(i_alpha_RTI,k)
779 : end do
780 0 : alpha_RTI_old_plus1(nz_old+1) = 0
781 :
782 : call interpolate_vector( &
783 : nz_old+1, xq_old_plus1, n, xq_new, &
784 : alpha_RTI_old_plus1, alpha_RTI_new, interp_pm, nwork, work, &
785 0 : 'mesh_adjust do_alpha_RTI', ierr)
786 0 : if (ierr /= 0) then
787 : return
788 : end if
789 :
790 0 : do k=nzlo,nzhi
791 0 : xh(i_alpha_RTI,k) = max(0d0,alpha_RTI_new(k+1-nzlo))
792 : end do
793 :
794 0 : n = nzlo - 1
795 0 : if (n > 0) then
796 0 : do k=1,n
797 0 : xh(i_alpha_RTI,k) = xh_old(i_alpha_RTI,k)
798 : end do
799 : end if
800 :
801 0 : if (nzhi < nz) then
802 0 : n = nz - nzhi - 1 ! nz-n = nzhi+1
803 0 : do k=0,n
804 0 : xh(i_alpha_RTI,nz-k) = xh_old(i_alpha_RTI,nz_old-k)
805 : end do
806 : end if
807 :
808 0 : end subroutine do_alpha_RTI
809 :
810 :
811 20 : subroutine do_interp_pt_val( &
812 : s, nz, nz_old, nzlo, nzhi, val, val_old, center_val, &
813 20 : xq, xq_old_plus1, xq_new, force_non_negative, &
814 40 : work, val_old_plus1, val_new, ierr)
815 0 : use interp_1d_def
816 : use interp_1d_lib
817 : type (star_info), pointer :: s
818 : integer, intent(in) :: nz, nz_old, nzlo, nzhi
819 : real(dp), dimension(:), pointer :: val, val_old
820 : real(dp), intent(in) :: center_val
821 : real(dp), dimension(:), pointer :: work
822 : real(dp), dimension(:) :: &
823 : xq, xq_old_plus1, xq_new, val_old_plus1, val_new
824 : logical, intent(in) :: force_non_negative
825 : integer, intent(out) :: ierr
826 : integer :: n, k
827 :
828 : include 'formats'
829 :
830 : ierr = 0
831 20 : n = nzhi - nzlo + 1
832 :
833 24210 : do k=1,nz_old
834 24210 : val_old_plus1(k) = val_old(k)
835 : end do
836 20 : val_old_plus1(nz_old+1) = center_val
837 :
838 : call interpolate_vector( &
839 : nz_old+1, xq_old_plus1, n, xq_new, &
840 : val_old_plus1, val_new, interp_pm, nwork, work, &
841 20 : 'mesh_adjust do_interp_pt_val', ierr)
842 20 : if (ierr /= 0) then
843 : return
844 : end if
845 :
846 3962 : do k=nzlo,nzhi
847 3962 : val(k) = val_new(k+1-nzlo)
848 : end do
849 :
850 20 : n = nzlo - 1
851 20 : if (n > 0) then
852 16860 : do k=1,n
853 16860 : val(k) = val_old(k)
854 : end do
855 : end if
856 :
857 20 : if (nzhi < nz) then
858 8 : n = nz - nzhi - 1 ! nz-n = nzhi+1
859 472 : do k=0,n
860 472 : val(nz-k) = val_old(nz_old-k)
861 : end do
862 : end if
863 :
864 20 : if (force_non_negative) then
865 3962 : do k=nzlo,nzhi
866 3962 : if (val(k) < 0) val(k) = 0
867 : end do
868 : end if
869 :
870 20 : end subroutine do_interp_pt_val
871 :
872 :
873 : subroutine do_interp_cell_val( &
874 : s, nz, nz_old, nzlo, nzhi, val_new_out, val_old, &
875 : xq, xq_old_plus1, dq, dq_old, work, val_old_plus1, val_new, ierr)
876 20 : use interp_1d_def
877 : use interp_1d_lib
878 : type (star_info), pointer :: s
879 : integer, intent(in) :: nz, nz_old, nzlo, nzhi
880 : real(dp), dimension(:), pointer :: val_new_out, val_old
881 : real(dp), dimension(:), pointer :: work
882 : real(dp), dimension(:) :: &
883 : xq, xq_old_plus1, dq, dq_old, val_old_plus1, val_new
884 : integer, intent(out) :: ierr
885 :
886 : real(dp), pointer, dimension(:) :: &
887 : mid_xq_new, mid_xq_old_plus1
888 : integer :: n, i, k
889 :
890 : ierr = 0
891 : n = nzhi - nzlo + 1
892 : call do_alloc(ierr)
893 : if (ierr /= 0) return
894 :
895 : do k=1,nz_old
896 : val_old_plus1(k) = val_old(k)
897 : mid_xq_old_plus1(k) = xq_old_plus1(k) + 0.5d0*dq_old(k)
898 : end do
899 : val_old_plus1(nz_old+1) = val_old_plus1(nz_old)
900 : mid_xq_old_plus1(nz_old+1) = 1
901 : do i=1,n
902 : mid_xq_new(i) = xq(nzlo+i-1) + 0.5d0*dq(nzlo+i-1)
903 : end do
904 :
905 : call interpolate_vector( &
906 : nz_old+1, mid_xq_old_plus1, n, mid_xq_new, &
907 : val_old_plus1, val_new, interp_pm, nwork, work, &
908 : 'mesh_adjust do_interp_cell_val', ierr)
909 : if (ierr /= 0) then
910 : call dealloc
911 : return
912 : end if
913 :
914 : do i=1,n
915 : val_new_out(nzlo+i-1) = val_new(i)
916 : end do
917 :
918 : n = nzlo - 1
919 : if (n > 0) then
920 : do i=1,n
921 : val_new_out(i) = val_old(i)
922 : end do
923 : end if
924 :
925 : if (nzhi < nz) then
926 : n = nz - nzhi - 1 ! nz-n = nzhi+1
927 : do i=0,n
928 : val_new_out(nz-i) = val_old(nz_old-i)
929 : end do
930 : end if
931 :
932 : call dealloc
933 :
934 : contains
935 :
936 : subroutine do_alloc(ierr)
937 : integer, intent(out) :: ierr
938 : call do_work_arrays(.true.,ierr)
939 : end subroutine do_alloc
940 :
941 : subroutine dealloc
942 : integer :: ierr
943 : call do_work_arrays(.false.,ierr)
944 : end subroutine dealloc
945 :
946 : subroutine do_work_arrays(alloc_flag, ierr)
947 : use alloc, only: work_array
948 : logical, intent(in) :: alloc_flag
949 : integer, intent(out) :: ierr
950 : logical, parameter :: crit = .false.
951 : ierr = 0
952 : call work_array(s, alloc_flag, crit, &
953 : mid_xq_old_plus1, nz_old+1, nz_alloc_extra, 'mesh_adjust', ierr)
954 : if (ierr /= 0) return
955 : call work_array(s, alloc_flag, crit, &
956 : mid_xq_new, n, nz_alloc_extra, 'mesh_adjust', ierr)
957 : if (ierr /= 0) return
958 : end subroutine do_work_arrays
959 :
960 : end subroutine do_interp_cell_val
961 :
962 :
963 10 : subroutine do_lnR_and_lnd( &
964 40 : s, nz, nz_old, nzlo, nzhi, cell_type, comes_from, &
965 10 : xh, xh_old, xmstar, lnd_old, lnPgas_old, &
966 40 : dqbar, dqbar_old, old_r, old_m, old_rho, &
967 20 : dq, dq_old, xq, xq_old_plus1, density_new, work, &
968 30 : Vol_old_plus1, Vol_new, new_r, Vol_init, &
969 20 : interp_Vol_new, interp_xq, density_init, ierr)
970 : use interp_1d_def
971 : use interp_1d_lib
972 : use num_lib
973 : type (star_info), pointer :: s
974 : integer, intent(in) :: nz, nz_old, nzlo, nzhi, comes_from(:)
975 : integer :: cell_type(:)
976 : real(dp), dimension(:,:), pointer :: xh, xh_old
977 : real(dp), intent(in) :: xmstar
978 : real(dp), dimension(:), pointer :: work
979 : real(dp), dimension(:) :: &
980 : lnd_old, lnPgas_old, dqbar, dqbar_old, old_r, old_m, old_rho, &
981 : xq, dq, dq_old, xq_old_plus1, density_new, &
982 : Vol_old_plus1, Vol_new, new_r, Vol_init, &
983 : interp_Vol_new, interp_xq, density_init
984 : integer, intent(out) :: ierr
985 :
986 : integer :: k, from_k, n, interp_lo, interp_hi, interp_n
987 10 : real(dp) :: Vol_min, Vol_max, cell_Vol, Vol_center, Vm1, V00, Vp1
988 :
989 : logical, parameter :: dbg = .false., trace_PE_residual = .false.
990 :
991 : include 'formats'
992 :
993 : ! NOTE: for interpolating volume, need to add point at center
994 :
995 10 : ierr = 0
996 :
997 10 : interp_lo = max(1, nzlo-1)
998 10 : interp_hi = min(nz, nzhi+1)
999 10 : interp_n = interp_hi - interp_lo + 1
1000 :
1001 12105 : do k=1,nz_old
1002 12105 : Vol_old_plus1(k) = four_thirds_pi*old_r(k)*old_r(k)*old_r(k)
1003 : end do
1004 10 : Vol_center = four_thirds_pi*s% R_center*s% R_center*s% R_center
1005 10 : Vol_old_plus1(nz_old+1) = Vol_center
1006 :
1007 : ! testing -- check for Vol_old_plus1 strictly decreasing
1008 12105 : do k = 2, nz_old+1
1009 12105 : if (Vol_old_plus1(k) >= Vol_old_plus1(k-1)) then
1010 0 : ierr = -1
1011 0 : if (.not. dbg) return
1012 : write(*,3) 'bad old vol', k, nz_old
1013 : write(*,1) 'Vol_old_plus1(k)', Vol_old_plus1(k)
1014 : write(*,1) 'Vol_old_plus1(k-1)', Vol_old_plus1(k-1)
1015 : write(*,'(A)')
1016 : call mesa_error(__FILE__,__LINE__,'debug: mesh adjust: do_lnR_and_lnd')
1017 : end if
1018 : end do
1019 :
1020 : ! testing -- check for q strictly decreasing
1021 10623 : do k = 2, nz
1022 10623 : if (xq(k) <= xq(k-1)) then
1023 0 : ierr = -1
1024 0 : if (.not. dbg) return
1025 :
1026 : write(*,3) 'bad xq', k, nz, xq(k), xq(k-1)
1027 : call mesa_error(__FILE__,__LINE__,'debug: mesh adjust: do_lnR_and_lnd')
1028 : end if
1029 : end do
1030 :
1031 1995 : do k=1,interp_n
1032 1995 : interp_xq(k) = xq(interp_lo+k-1)
1033 : end do
1034 :
1035 : call interpolate_vector( &
1036 : nz_old+1, xq_old_plus1, interp_n, interp_xq, Vol_old_plus1, &
1037 10 : interp_Vol_new, interp_pm, nwork, work, 'mesh_adjust do_lnR_and_lnd', ierr)
1038 10 : if (ierr /= 0) then
1039 : if (.not. dbg) return
1040 : write(*,*) 'failed in interpolate_vector'
1041 : call mesa_error(__FILE__,__LINE__,'debug: mesh_adjust')
1042 : end if
1043 :
1044 1995 : do k=1,interp_n
1045 1995 : Vol_new(interp_lo+k-1) = interp_Vol_new(k)
1046 : end do
1047 :
1048 10 : if (Vol_new(interp_lo+1) >= Vol_new(interp_lo)) then
1049 0 : Vol_new(interp_lo+1) = (Vol_new(interp_lo) + Vol_new(interp_lo+2))/2
1050 : if (dbg) write(*,2) 'fix Vol_new at lo+1', interp_lo+1, Vol_new(interp_lo+1)
1051 0 : if (Vol_new(interp_lo+1) >= Vol_new(interp_lo)) then
1052 0 : ierr = -1
1053 0 : if (.not. dbg) return
1054 : write(*,*) '(Vol_new(interp_lo+1) >= Vol_new(interp_lo))'
1055 : call mesa_error(__FILE__,__LINE__,'debug: mesh_adjust')
1056 : end if
1057 : end if
1058 :
1059 1975 : do k = interp_lo+1, interp_hi-1
1060 1975 : if (Vol_new(k+1) >= Vol_new(k) .or. Vol_new(k) >= Vol_new(k-1)) then
1061 : if (dbg) write(*,2) 'fix interpolated Vol_new', &
1062 : k, Vol_new(k+1), Vol_new(k), Vol_new(k-1)
1063 0 : Vol_min = minval(Vol_new(k-1:k+1))
1064 0 : Vol_max = maxval(Vol_new(k-1:k+1))
1065 0 : if (Vol_min == Vol_max .or. is_bad(Vol_min) .or. is_bad(Vol_max)) then
1066 0 : ierr = -1
1067 0 : if (s% stop_for_bad_nums) then
1068 0 : write(*,2) 'Vol_min', k, Vol_min
1069 0 : write(*,2) 'Vol_max', k, Vol_max
1070 0 : call mesa_error(__FILE__,__LINE__,'mesh_adjust')
1071 : end if
1072 0 : if (.not. dbg) return
1073 : write(*,1) 'Vol_min', Vol_min
1074 : write(*,1) 'Vol_max', Vol_max
1075 : call mesa_error(__FILE__,__LINE__,'debug: mesh_adjust')
1076 : end if
1077 0 : Vm1 = Vol_new(k-1)
1078 0 : V00 = Vol_new(k)
1079 0 : Vp1 = Vol_new(k+1)
1080 0 : Vol_new(k-1) = Vol_max
1081 0 : Vol_new(k) = (Vol_max + Vol_min)/2
1082 0 : Vol_new(k+1) = Vol_min
1083 : if (dbg) write(*,2) 'new Vol_new', &
1084 : k, Vol_new(k+1), Vol_new(k), Vol_new(k-1)
1085 0 : if (Vol_new(k+1) >= Vol_new(k) .or. Vol_new(k) >= Vol_new(k-1)) then
1086 0 : ierr = -1
1087 0 : if (.not. dbg) return
1088 : write(*,1) 'Vol_new(k-1)', Vol_new(k-1)
1089 : write(*,1) 'Vol_new(k)', Vol_new(k)
1090 : write(*,1) 'Vol_new(k+1)', Vol_new(k+1)
1091 : call mesa_error(__FILE__,__LINE__,'debug: do_lnR_and_lnd in mesh adjust: interpolation gave non-pos volume')
1092 : end if
1093 : end if
1094 : end do
1095 :
1096 10 : call set1_new_r(nzlo)
1097 1975 : do k = nzlo, min(nzhi,nz-1)
1098 1965 : if (ierr /= 0) cycle
1099 :
1100 1965 : call set1_new_r(k+1)
1101 :
1102 1965 : if (cell_type(k) == unchanged_type) then
1103 493 : call store_lnd_in_xh(s, k, lnd_old(comes_from(k)), xh)
1104 493 : density_new(k) = old_rho(comes_from(k))
1105 493 : cycle
1106 : end if
1107 :
1108 1472 : if (new_r(k) <= new_r(k+1)) then
1109 : if (dbg) then
1110 : write(*,*) 'do_lnR_and_lnd: (new_r(k) <= new_r(k+1))'
1111 : stop
1112 : end if
1113 0 : ierr = -1; cycle
1114 : end if
1115 :
1116 1472 : cell_Vol = Vol_new(k)-Vol_new(k+1)
1117 1472 : if (cell_Vol <= 0) then
1118 : if (dbg) then
1119 : write(*,2) 'do_lnR_and_lnd: cell_Vol <= 0', k
1120 : end if
1121 0 : ierr = -1; cycle
1122 : end if
1123 1472 : if (dq(k) <= 0) then
1124 : if (dbg) then
1125 : write(*,2) 'do_lnR_and_lnd: dq(k) <= 0', k
1126 : end if
1127 0 : ierr = -1; cycle
1128 : end if
1129 1472 : density_new(k) = xmstar*dq(k)/cell_Vol
1130 1482 : call store_rho_in_xh(s, k, density_new(k), xh)
1131 :
1132 : end do
1133 :
1134 10 : if (ierr /= 0) then
1135 : if (.not. dbg) return
1136 : call mesa_error(__FILE__,__LINE__,'debug: failed in mesh adjust do_lnR_and_lnd')
1137 : end if
1138 :
1139 10 : n = nzlo - 1
1140 10 : if (n > 0) then
1141 8430 : do k=1,n
1142 8420 : new_r(k) = old_r(k)
1143 8420 : density_new(k) = old_rho(k)
1144 8430 : Vol_new(k) = four_thirds_pi*new_r(k)*new_r(k)*new_r(k)
1145 : end do
1146 : end if
1147 :
1148 10 : if (nzhi < nz) then
1149 4 : n = nz - nzhi - 1 ! nz-n = nzhi+1
1150 236 : do k=0,n
1151 232 : new_r(nz-k) = old_r(nz_old-k)
1152 232 : density_new(nz-k) = old_rho(nz_old-k)
1153 236 : Vol_new(nz-k) = four_thirds_pi*new_r(nz-k)*new_r(nz-k)*new_r(nz-k)
1154 : end do
1155 : else ! nzhi == nz
1156 6 : density_new(nz) = xmstar*dq(nz)/(Vol_new(nz) - Vol_center)
1157 6 : new_r(nz) = pow(Vol_new(nz)/four_thirds_pi, one_third)
1158 :
1159 : if (dbg) then
1160 : write(*,2) 'old_rho(nz_old)', nz_old, old_rho(nz_old)
1161 : write(*,2) 'density_new(nz)', nz, density_new(nz)
1162 : end if
1163 :
1164 : end if
1165 :
1166 10633 : do k=1,nz
1167 10623 : Vol_init(k) = Vol_new(k)
1168 10633 : density_init(k) = density_new(k)
1169 : end do
1170 :
1171 10643 : do k=1,nz
1172 10623 : from_k = comes_from(k)
1173 10623 : if (new_r(k) == old_r(from_k)) then
1174 9187 : xh(s%i_lnR,k) = xh_old(s%i_lnR,from_k)
1175 : else
1176 1436 : call store_r_in_xh(s,k,new_r(k),xh)
1177 : end if
1178 10633 : if (density_new(k) == old_rho(from_k)) then
1179 9147 : xh(s%i_lnd,k) = xh_old(s%i_lnd,from_k)
1180 : else
1181 1476 : call store_rho_in_xh(s,k,density_new(k),xh)
1182 : end if
1183 : end do
1184 :
1185 :
1186 : contains
1187 :
1188 1975 : subroutine set1_new_r(k)
1189 : integer, intent(in) :: k
1190 : include 'formats'
1191 1975 : if (cell_type(k) == unchanged_type) then
1192 503 : new_r(k) = old_r(comes_from(k))
1193 : else
1194 1472 : new_r(k) = pow(Vol_new(k)/four_thirds_pi, one_third)
1195 : end if
1196 10 : end subroutine set1_new_r
1197 :
1198 :
1199 : end subroutine do_lnR_and_lnd
1200 :
1201 :
1202 10623 : subroutine do_xa( &
1203 10623 : s, nz, nz_old, k, species, cell_type, comes_from, &
1204 10623 : xa, xa_old, xa_c0, xa_c1, xa_c2, &
1205 : xq, dq, xq_old, dq_old, mesh_adjust_use_quadratic, ierr)
1206 : use chem_def, only: chem_isos
1207 : type (star_info), pointer :: s
1208 : integer, intent(in) :: nz, nz_old, species, k, cell_type(:), comes_from(:)
1209 : real(dp), dimension(:,:), pointer :: xa, xa_old
1210 : real(dp), dimension(:,:) :: xa_c0, xa_c1, xa_c2
1211 : real(dp), dimension(:), pointer :: xq, dq, xq_old, dq_old
1212 : logical, intent(in) :: mesh_adjust_use_quadratic
1213 : integer, intent(out) :: ierr
1214 :
1215 : integer :: j, jj, k_old, k_old_last, kdbg, order
1216 95607 : real(dp) :: xq_outer, cell_dq, xa_sum, total(species)
1217 : logical :: dbg_get_integral
1218 :
1219 : include 'formats'
1220 :
1221 10623 : ierr = 0
1222 :
1223 10623 : kdbg = -1074
1224 :
1225 10623 : if (mesh_adjust_use_quadratic) then
1226 10623 : order = 2
1227 : else
1228 0 : order = 1
1229 : end if
1230 :
1231 10623 : if (cell_type(k) == unchanged_type .or. &
1232 : cell_type(k) == revised_type) then
1233 82359 : do j=1,species
1234 82359 : xa(j,k) = xa_old(j,comes_from(k))
1235 : end do
1236 : return
1237 : end if
1238 :
1239 1472 : xq_outer = xq(k)
1240 1472 : if (k == nz) then
1241 0 : cell_dq = 1 - xq_outer
1242 : else
1243 1472 : cell_dq = dq(k)
1244 : end if
1245 :
1246 1472 : k_old = max(comes_from(k)-1,1)
1247 :
1248 : ! sum the old abundances between xq_outer and xq_inner
1249 1472 : dbg_get_integral = .false.
1250 13248 : total(:) = 0
1251 13248 : do j=1,species
1252 11776 : dbg_get_integral = (k == kdbg) .and. (j == 1) ! h1
1253 11776 : if (dbg_get_integral) write(*,2) trim(chem_isos% name(s% chem_id(j)))
1254 : call get_xq_integral( &
1255 : k_old, nz_old, xq_old, xq_outer, cell_dq, &
1256 : order, xa_c0(:,j), xa_c1(:,j), xa_c2(:,j), &
1257 13248 : total(j), dbg_get_integral, k_old_last, ierr)
1258 : end do
1259 :
1260 13248 : xa(:,k) = total(:)/cell_dq
1261 :
1262 1472 : if (k == kdbg) then
1263 0 : do j=1,species
1264 0 : write(*,2) 'new ' // trim(chem_isos% name(s% chem_id(j))), k, xa(j,k)
1265 : end do
1266 : end if
1267 :
1268 13248 : do j=1,species
1269 13248 : if (xa(j,k) > 1 + 1d-8 .or. xa(j,k) < -1d-8) then
1270 0 : ierr = -1
1271 0 : return
1272 :
1273 : do jj=1,species
1274 : write(*,1) 'xa ' // trim(chem_isos% name(s% chem_id(jj))), xa(jj,k)
1275 : end do
1276 : write(*,'(A)')
1277 : write(*,2) 'sum xa', k, sum(xa(:,k))
1278 : write(*,'(A)')
1279 : write(*,2) 'xa ' // trim(chem_isos% name(s% chem_id(j))), k, xa(j,k)
1280 : write(*,'(A)')
1281 : write(*,2) 'xq_outer', k, xq_outer
1282 : write(*,2) 'xq_inner', k, xq_outer + cell_dq
1283 : write(*,2) 'cell_dq', k, cell_dq
1284 : write(*,'(A)')
1285 : write(*,2) 'xq_old(k_old)', k_old, xq_old(k_old)
1286 : write(*,2) 'xq_inner(k_old)', k_old, xq_old(k_old)+dq_old(k_old)
1287 : write(*,2) 'dq_old(k_old)', k_old, dq_old(k_old)
1288 : write(*,'(A)')
1289 : write(*,2) 'xa_c0(k_old,j)', k_old, xa_c0(k_old,j)
1290 : write(*,2) 'xa_c1(k_old,j)', k_old, xa_c1(k_old,j)
1291 : write(*,2) 'xa_c2(k_old,j)', k_old, xa_c2(k_old,j)
1292 : write(*,'(A)')
1293 : write(*,2) 'old outer', k_old, xa_c0(k_old,j) + xa_c1(k_old,j)*dq_old(k_old)/2
1294 : write(*,2) 'old inner', k_old, xa_c0(k_old,j) - xa_c1(k_old,j)*dq_old(k_old)/2
1295 : write(*,'(A)')
1296 : call mesa_error(__FILE__,__LINE__,'debug: mesh adjust: do_xa')
1297 : end if
1298 : end do
1299 :
1300 13248 : xa_sum = sum(xa(:,k))
1301 : !write(*,1) 'xa_sum', xa_sum
1302 :
1303 1472 : if (is_bad(xa_sum)) then
1304 0 : ierr = -1
1305 0 : if (s% stop_for_bad_nums) then
1306 0 : write(*,2) 'xa_sum', k, xa_sum
1307 0 : call mesa_error(__FILE__,__LINE__,'mesh adjust: do_xa')
1308 : end if
1309 0 : return
1310 :
1311 : write(*,*) 'xa_sum', xa_sum
1312 : write(*,*) 'bug in revise mesh, do_xa bad num: k', k
1313 : call mesa_error(__FILE__,__LINE__,'debug: mesh adjust: do_xa')
1314 : end if
1315 :
1316 1472 : if (abs(1-xa_sum) > 1d-3) then
1317 0 : ierr = -1
1318 0 : return
1319 : end if
1320 :
1321 13248 : xa(:,k) = xa(:,k) / xa_sum
1322 :
1323 10623 : end subroutine do_xa
1324 :
1325 :
1326 10623 : subroutine do1_lnT( &
1327 : s, nz_old, k, &
1328 10623 : species, cell_type, comes_from, &
1329 : xa, xh, xh_old, &
1330 10623 : xq, dq, xq_old, dq_old, eta_old, energy_old, lnT_old, &
1331 10623 : specific_PE_old, specific_KE_old, w_old, &
1332 21246 : density_new, energy_new, ierr)
1333 10623 : use eos_def
1334 : use star_utils, only: set_rmid, cell_specific_PE, cell_specific_KE
1335 : type (star_info), pointer :: s
1336 : integer, intent(in) :: nz_old, k, species, cell_type(:), comes_from(:)
1337 : real(dp), dimension(:,:), pointer :: xa, xh, xh_old
1338 : real(dp), dimension(:) :: &
1339 : xq, dq, xq_old, dq_old, eta_old, energy_old, lnT_old, &
1340 : specific_PE_old, specific_KE_old, w_old, density_new, energy_new
1341 : integer, intent(out) :: ierr
1342 :
1343 : integer :: k_old
1344 : real(dp) :: &
1345 21246 : Rho, logRho, xq_outer, cell_dq, avg_energy, avg_PE, avg_KE, &
1346 21246 : new_PE, new_KE, max_delta_energy, delta_energy, revised_energy, &
1347 95607 : sum_lnT, avg_lnT, new_lnT, sum_energy, new_xa(species), &
1348 10623 : d_dlnR00, d_dlnRp1, d_dv00, d_dvp1
1349 :
1350 : include 'formats'
1351 :
1352 10623 : ierr = 0
1353 95607 : new_xa(:) = xa(:,k)
1354 10623 : k_old = comes_from(k)
1355 :
1356 10623 : if (cell_type(k) == unchanged_type) then
1357 9151 : xh(s%i_lnT, k) = xh_old(s%i_lnT, k_old)
1358 9151 : energy_new(k) = energy_old(k_old)
1359 9151 : if (is_bad(energy_old(k_old))) then
1360 0 : write(*,2) 'energy_old(k_old)', k_old, energy_old(k_old)
1361 0 : call mesa_error(__FILE__,__LINE__,'debug: mesh adjust: do1_lnT')
1362 : end if
1363 9151 : return
1364 : end if
1365 :
1366 1472 : xq_outer = xq(k)
1367 1472 : cell_dq = dq(k)
1368 :
1369 1472 : if (cell_type(k) == revised_type) then
1370 0 : avg_lnT = get_lnT_from_xh(s, k, xh_old)
1371 : else ! find average lnT between xq_outer and xq_inner
1372 : call get_old_value_integral( &
1373 : k, k_old, nz_old, xq_old, dq_old, xq_outer, cell_dq, &
1374 1472 : lnT_old, sum_lnT, dbg, ierr)
1375 1472 : if (ierr /= 0) then
1376 : if (dbg) write(*,*) 'get_old_value_integral lnT failed for do1_lnT'
1377 : if (.not. dbg) return
1378 : call mesa_error(__FILE__,__LINE__,'debug: mesh adjust: do1_lnT')
1379 : end if
1380 1472 : avg_lnT = sum_lnT/cell_dq
1381 : end if
1382 :
1383 1472 : if (is_bad(avg_lnT) .or. avg_lnT < 0 .or. avg_lnT > 100) then
1384 0 : ierr = -1
1385 0 : if (s% stop_for_bad_nums) then
1386 0 : write(*,2) 'avg_lnT', k, avg_lnT
1387 0 : call mesa_error(__FILE__,__LINE__,'mesh adjust: do1_lnT')
1388 : end if
1389 0 : return
1390 : end if
1391 :
1392 1472 : if (.not. s% mesh_adjust_get_T_from_E) then
1393 0 : call store_lnT_in_xh(s, k, avg_lnT, xh)
1394 0 : energy_new(k) = energy_old(k_old)
1395 0 : if (is_bad(energy_old(k_old))) then
1396 0 : write(*,2) 'energy_old(k_old)', k_old, energy_old(k_old)
1397 0 : call mesa_error(__FILE__,__LINE__,'debug: mesh adjust: do1_lnT')
1398 : end if
1399 0 : return
1400 : end if
1401 :
1402 1472 : if (eta_old(k_old) >= eta_limit) then
1403 0 : call store_lnT_in_xh(s, k, avg_lnT, xh)
1404 0 : energy_new(k) = energy_old(k_old)
1405 0 : if (is_bad(energy_old(k_old))) then
1406 0 : write(*,2) 'eta_old(k_old)', k_old, eta_old(k_old)
1407 0 : write(*,2) 'energy_old(k_old)', k_old, energy_old(k_old)
1408 0 : call mesa_error(__FILE__,__LINE__,'debug: mesh adjust: do1_lnT')
1409 : end if
1410 0 : return
1411 : end if
1412 :
1413 : if (dbg) write(*,2) 'eta_old(k_old)', k_old, eta_old(k_old)
1414 :
1415 1472 : if (cell_type(k) == revised_type) then
1416 0 : avg_energy = energy_old(k_old)
1417 : else ! find average internal energy between q_outer and q_inner
1418 : call get_old_value_integral( &
1419 : k, k_old, nz_old, xq_old, dq_old, xq_outer, cell_dq, &
1420 1472 : energy_old, sum_energy, dbg, ierr)
1421 1472 : if (ierr /= 0) then
1422 : if (dbg) write(*,*) 'get_old_value_integral failed for do1_lnT'
1423 : if (.not. dbg) return
1424 : call mesa_error(__FILE__,__LINE__,'debug: mesh adjust: energy_old do1_lnT')
1425 : end if
1426 1472 : avg_energy = sum_energy/cell_dq
1427 : end if
1428 :
1429 1472 : if (s% max_rel_delta_IE_for_mesh_total_energy_balance == 0d0) then
1430 :
1431 0 : energy_new(k) = avg_energy
1432 :
1433 : else
1434 :
1435 1472 : if (cell_type(k) == revised_type) then
1436 0 : avg_PE = specific_PE_old(k_old)
1437 : else ! find average potential energy between q_outer and q_inner
1438 : call get_old_value_integral( &
1439 : k, k_old, nz_old, xq_old, dq_old, xq_outer, cell_dq, &
1440 1472 : specific_PE_old, sum_energy, dbg, ierr)
1441 1472 : if (ierr /= 0) then
1442 : if (dbg) write(*,*) 'get_old_value_integral failed for do1_lnT'
1443 : if (.not. dbg) return
1444 : call mesa_error(__FILE__,__LINE__,'debug: mesh adjust: specific_PE_old do1_lnT')
1445 : end if
1446 1472 : avg_PE = sum_energy/cell_dq
1447 : end if
1448 :
1449 1472 : if (cell_type(k) == revised_type) then
1450 0 : avg_KE = specific_KE_old(k_old)
1451 : else ! find average kinetic energy between q_outer and q_inner
1452 : call get_old_value_integral( &
1453 : k, k_old, nz_old, xq_old, dq_old, xq_outer, cell_dq, &
1454 1472 : specific_KE_old, sum_energy, dbg, ierr)
1455 1472 : if (ierr /= 0) then
1456 : if (dbg) write(*,*) 'get_old_value_integral failed for do1_lnT'
1457 : if (.not. dbg) return
1458 : call mesa_error(__FILE__,__LINE__,'debug: mesh adjust: specific_KE_old do1_lnT')
1459 : end if
1460 1472 : avg_KE = sum_energy/cell_dq
1461 : end if
1462 :
1463 1472 : if (ierr /= 0) return
1464 1472 : new_PE = cell_specific_PE(s,k,d_dlnR00,d_dlnRp1)
1465 1472 : if (s% u_flag) then
1466 0 : s% u(k) = s% xh(s% i_u,k)
1467 1472 : else if (s% v_flag) then
1468 0 : s% v(k) = s% xh(s% i_v,k)
1469 0 : if (k < s% nz) s% v(k+1) = s% xh(s% i_v,k+1)
1470 : end if
1471 1472 : new_KE = cell_specific_KE(s,k,d_dv00,d_dvp1)
1472 :
1473 1472 : max_delta_energy = avg_energy*s% max_rel_delta_IE_for_mesh_total_energy_balance
1474 1472 : delta_energy = avg_PE + avg_KE - (new_PE + new_KE)
1475 1472 : if (abs(delta_energy) > max_delta_energy) then
1476 0 : delta_energy = sign(max_delta_energy,delta_energy)
1477 : end if
1478 1472 : energy_new(k) = avg_energy + delta_energy
1479 :
1480 1472 : if (energy_new(k) <= 0d0) then
1481 0 : write(*,2) 'energy_new(k) <= 0d0', k, energy_new(k), avg_energy
1482 0 : energy_new(k) = avg_energy
1483 : end if
1484 :
1485 : end if
1486 :
1487 : ! call eos to calculate lnT from new internal energy
1488 :
1489 1472 : Rho = density_new(k)
1490 1472 : logRho = log10(Rho)
1491 : call set_lnT_for_energy( &
1492 : s, k, &
1493 : s% net_iso(ih1), s% net_iso(ihe3), s% net_iso(ihe4), species, new_xa, &
1494 1472 : Rho, logRho, energy_new(k), avg_lnT, new_lnT, revised_energy, ierr)
1495 1472 : if (ierr /= 0) then
1496 : if (dbg) write(*,*) 'set_lnT_for_energy failed', k
1497 0 : new_lnT = avg_lnT
1498 0 : energy_new(k) = energy_old(k_old)
1499 0 : ierr = 0
1500 : end if
1501 :
1502 1472 : call store_lnT_in_xh(s, k, new_lnT, xh)
1503 :
1504 1472 : if (ierr /= 0) then
1505 0 : write(*,2) 'mesh_adjust do1_lnT ierr', ierr
1506 0 : call mesa_error(__FILE__,__LINE__,'do1_lnT')
1507 : end if
1508 :
1509 10623 : end subroutine do1_lnT
1510 :
1511 :
1512 5888 : subroutine get_old_value_integral( &
1513 5888 : k_new, k_old_in, nz_old, xq_old, dq_old, xq_outer, dq_range, &
1514 5888 : value_old, integral, dbg, ierr)
1515 : integer, intent(in) :: k_new, k_old_in, nz_old
1516 : real(dp), intent(in) :: xq_old(:), dq_old(:), xq_outer, dq_range
1517 : real(dp), intent(in), dimension(:) :: value_old
1518 : real(dp), intent(out) :: integral
1519 : logical, intent(in) :: dbg
1520 : integer, intent(out) :: ierr
1521 : real(dp), pointer :: p(:,:)
1522 : nullify(p)
1523 : call get_old_integral( &
1524 : k_new, k_old_in, nz_old, xq_old, dq_old, xq_outer, dq_range, &
1525 5888 : value_old, p, integral, dbg, ierr)
1526 10623 : end subroutine get_old_value_integral
1527 :
1528 :
1529 5888 : subroutine get_old_integral( &
1530 5888 : k_new, k_old_in, nz_old, xq_old, dq_old, xq_outer, dq_range, &
1531 5888 : value_old, xh_old, integral, dbg, ierr)
1532 : integer, intent(in) :: k_new, k_old_in, nz_old
1533 : real(dp), intent(in) :: xq_old(:), dq_old(:), xq_outer, dq_range
1534 : real(dp), intent(in) :: value_old(:)
1535 : real(dp), intent(in), pointer :: xh_old(:,:)
1536 : real(dp), intent(out) :: integral
1537 : logical, intent(in) :: dbg
1538 : integer, intent(out) :: ierr
1539 :
1540 : integer :: k, k_old
1541 5888 : real(dp) :: xq_inner, sum_dqs, old_xq_outer, old_xq_inner, &
1542 5888 : dq_overlap, val
1543 :
1544 : include 'formats'
1545 :
1546 5888 : if (dbg) write(*,*)
1547 :
1548 5888 : ierr = 0
1549 :
1550 5888 : k_old = k_old_in
1551 : ! move starting k_old if necessary
1552 0 : do
1553 5888 : if (k_old <= 1) exit
1554 5888 : if (xq_old(k_old) <= xq_outer) exit
1555 5888 : k_old = k_old - 1
1556 : end do
1557 :
1558 5888 : xq_inner = xq_outer + dq_range
1559 5888 : old_xq_inner = xq_old(k_old)
1560 5888 : sum_dqs = 0
1561 5888 : integral = 0d0
1562 :
1563 5888 : if (dbg) write(*,*)
1564 5888 : if (dbg) write(*,3) 'k_new k_old xq_outer xq_inner dq_range', &
1565 0 : k_new, k_old, xq_outer, xq_inner, dq_range
1566 :
1567 17664 : do k = k_old, nz_old
1568 :
1569 17664 : if (dq_range <= sum_dqs) exit
1570 11776 : old_xq_outer = old_xq_inner
1571 11776 : if (k == nz_old) then
1572 0 : old_xq_inner = 1
1573 : else
1574 11776 : old_xq_inner = xq_old(k+1)
1575 : end if
1576 :
1577 11776 : if (dbg) write(*,3) 'k_new k_old old_xq_outer old_xq_inner', &
1578 0 : k_new, k, old_xq_outer, old_xq_inner
1579 :
1580 11776 : val = value_old(k)
1581 :
1582 17664 : if (old_xq_inner <= xq_inner .and. old_xq_outer >= xq_outer) then
1583 :
1584 11776 : if (dbg) write(*,1) 'entire old cell is in new range'
1585 :
1586 11776 : sum_dqs = sum_dqs + dq_old(k)
1587 11776 : integral = integral + val*dq_old(k)
1588 :
1589 0 : else if (old_xq_inner >= xq_inner .and. old_xq_outer <= xq_outer) then
1590 :
1591 0 : if (dbg) write(*,1) 'entire new range is in this old cell'
1592 :
1593 0 : sum_dqs = dq_range
1594 0 : integral = val*dq_range
1595 :
1596 : else ! only use the part of old cell that is in new range
1597 :
1598 0 : if (xq_inner <= old_xq_inner) then
1599 :
1600 0 : if (dbg) write(*,1) 'last part of the new range'
1601 :
1602 0 : integral = integral + val*(dq_range - sum_dqs)
1603 0 : sum_dqs = dq_range
1604 :
1605 : else ! partial overlap -- general case
1606 :
1607 0 : dq_overlap = max(0d0, old_xq_inner - xq_outer)
1608 0 : sum_dqs = sum_dqs + dq_overlap
1609 0 : integral = integral + val*dq_overlap
1610 :
1611 0 : if (dbg) write(*,1) 'partial overlap'
1612 :
1613 : end if
1614 :
1615 : end if
1616 :
1617 : end do
1618 :
1619 5888 : if (dbg) write(*,2) 'integral/dq_range', &
1620 0 : k_new, integral/dq_range, integral, dq_range
1621 5888 : if (dbg) write(*,*)
1622 :
1623 5888 : end subroutine get_old_integral
1624 :
1625 :
1626 1472 : subroutine set_lnT_for_energy( &
1627 1472 : s, k, h1, he3, he4, species, xa, &
1628 : Rho, logRho, energy, lnT_guess, lnT, result_energy, ierr)
1629 : type (star_info), pointer :: s
1630 : integer, intent(in) :: k, h1, he3, he4, species
1631 : real(dp), intent(in) :: xa(species), Rho, logRho, energy, lnT_guess
1632 : real(dp), intent(out) :: lnT, result_energy
1633 : integer, intent(out) :: ierr
1634 : call set_lnT_for_energy_with_tol( &
1635 : s, k, h1, he3, he4, species, xa, 1d-11, &
1636 1472 : Rho, logRho, energy, lnT_guess, lnT, result_energy, ierr)
1637 0 : end subroutine set_lnT_for_energy
1638 :
1639 :
1640 1472 : subroutine set_lnT_for_energy_with_tol( &
1641 1472 : s, k, h1, he3, he4, species, xa, tol, &
1642 : Rho, logRho, energy, lnT_guess, lnT, result_energy, ierr)
1643 : use eos_support, only: solve_eos_given_DE
1644 : use eos_def, only: num_eos_basic_results, num_eos_d_dxa_results, i_lnE
1645 : use chem_lib, only: basic_composition_info
1646 : type (star_info), pointer :: s
1647 : integer, intent(in) :: k, h1, he3, he4, species
1648 : real(dp), intent(in) :: xa(species), tol, Rho, logRho, energy, lnT_guess
1649 : real(dp), intent(out) :: lnT, result_energy
1650 : integer, intent(out) :: ierr
1651 :
1652 : real(dp) :: &
1653 39744 : logT, res(num_eos_basic_results), &
1654 78016 : d_dlnd(num_eos_basic_results), d_dlnT(num_eos_basic_results), &
1655 38272 : d_dxa(num_eos_d_dxa_results, s% species), &
1656 1472 : logT_tol, logE_tol
1657 :
1658 : include 'formats'
1659 :
1660 1472 : ierr = 0
1661 :
1662 1472 : logT_tol = tol ! 1d-11
1663 1472 : logE_tol = tol ! 1d-11
1664 : call solve_eos_given_DE( &
1665 : s, k, xa(:), &
1666 : logRho, log10(energy), lnT_guess/ln10, &
1667 : logT_tol, logE_tol, &
1668 : logT, res, d_dlnd, d_dlnT, &
1669 : d_dxa, &
1670 1472 : ierr)
1671 1472 : lnT = logT*ln10
1672 :
1673 1472 : result_energy = exp(res(i_lnE))
1674 :
1675 1472 : if (ierr /= 0 .or. is_bad_num(lnT)) then
1676 0 : ierr = -1
1677 0 : if (s% stop_for_bad_nums) then
1678 0 : write(*,2) 'lnT', k, lnT
1679 0 : call mesa_error(__FILE__,__LINE__,'mesh adjust: do1_lnT')
1680 : end if
1681 : return
1682 : end if
1683 :
1684 : contains
1685 :
1686 : subroutine show
1687 : include 'formats'
1688 : write(*,'(A)')
1689 : write(*,*) 'set_lnT_for_energy ierr', ierr
1690 : write(*,*) 'k', k
1691 : write(*,'(A)')
1692 : write(*,1) 'lnT =', lnT
1693 : write(*,'(A)')
1694 : write(*,1) 'logRho =', logRho
1695 : write(*,1) 'logT_guess =', lnT_guess/ln10
1696 : write(*,1) 'energy =', energy
1697 : write(*,'(A)')
1698 : write(*,'(A)')
1699 1472 : end subroutine show
1700 :
1701 : end subroutine set_lnT_for_energy_with_tol
1702 :
1703 :
1704 : ! Stiriba, Youssef, Appl, Numer. Math. 45, 499-511. 2003.
1705 :
1706 : ! LPP-HARMOD -- local piecewise parabolic reconstruction
1707 :
1708 : ! interpolant is derived to conserve integral of v in cell k
1709 : ! interpolant slope at cell midpoint is harmonic mean of slopes between adjacent cells
1710 : ! where these slopes between cells are defined as the difference between cell averages
1711 : ! divided by the distance between cell midpoints.
1712 : ! interpolant curvature based on difference between the midpoint slope
1713 : ! and the smaller in magnitude of the slopes between adjacent cells.
1714 :
1715 : ! interpolant f(dq) = a + b*dq + (c/2)*dq^2, with dq = q - q_midpoint
1716 : ! c0 holds a's, c1 holds b's, and c2 holds c's.
1717 :
1718 :
1719 27576 : subroutine get1_lpp(k, ldv, nz, j, dq, v, quad, c0, c1, c2)
1720 : integer, intent(in) :: k, ldv, nz, j
1721 : real(dp), intent(in) :: dq(:) ! (nz)
1722 : real(dp), intent(in) :: v(:,:) ! (ldv,nz)
1723 : logical, intent(in) :: quad
1724 : real(dp), dimension(:) :: c0, c1, c2
1725 :
1726 27576 : real(dp) :: vbdy1, vbdy2, dqhalf, sm1, s00, sprod
1727 : real(dp), parameter :: rel_curvature_limit = 0.1d0
1728 :
1729 : logical :: dbg
1730 :
1731 : include 'formats'
1732 :
1733 27576 : if (k == 1 .or. k == nz) then
1734 48 : call set_const
1735 48 : return
1736 : end if
1737 :
1738 27528 : dbg = .false.
1739 : !dbg = (k == 30 .and. j == 3) ! .false.
1740 :
1741 27528 : sm1 = (v(j,k-1) - v(j,k)) / ((dq(k-1) + dq(k))/2)
1742 27528 : s00 = (v(j,k) - v(j,k+1)) / ((dq(k) + dq(k+1))/2)
1743 :
1744 27528 : sprod = sm1*s00
1745 27528 : if (sprod <= 0) then
1746 : ! at local min or max, so set slope and curvature to 0.
1747 8218 : call set_const
1748 8218 : return
1749 : end if
1750 :
1751 19310 : if (.not. quad) then
1752 19310 : c0(k) = v(j,k)
1753 19310 : c1(k) = (sm1 + s00)/2 ! use average to smooth abundance transitions
1754 19310 : c2(k) = 0 ! Yan Wang fixed this -- it was left out initially.
1755 : else
1756 0 : c1(k) = sprod*2/(s00 + sm1) ! harmonic mean slope
1757 0 : if (abs(sm1) <= abs(s00)) then
1758 0 : c2(k) = (sm1 - c1(k))/(2*dq(k))
1759 : else
1760 0 : c2(k) = (c1(k) - s00)/(2*dq(k))
1761 : end if
1762 0 : c0(k) = v(j,k) - c2(k)*dq(k)*dq(k)/24
1763 : end if
1764 :
1765 : ! check values at edges for monotonicity
1766 19310 : dqhalf = dq(k)/2
1767 19310 : vbdy1 = c0(k) + c1(k)*dqhalf + c2(k)/2*dqhalf*dqhalf ! value at face(k)
1768 19310 : vbdy2 = c0(k) - c1(k)*dqhalf + c2(k)/2*dqhalf*dqhalf ! value at face(k+1)
1769 19310 : if ((v(j,k-1) - vbdy1)*(vbdy1 - v(j,k)) < 0 .or. &
1770 : (v(j,k) - vbdy2)*(vbdy2 - v(j,k+1)) < 0) then
1771 : if (dbg) then
1772 : write(*,*) 'non-monotonic'
1773 : write(*,2) 'v(j,k-1)', k-1, v(j,k-1)
1774 : write(*,2) 'vbdy1', k, vbdy1
1775 : write(*,2) 'v(j,k)', k, v(j,k)
1776 : write(*,2) 'vbdy2', k, vbdy2
1777 : write(*,2) 'v(j,k+1)', k+1, v(j,k+1)
1778 : write(*,'(A)')
1779 : write(*,2) 'v(j,k-1) - vbdy1', k, v(j,k-1) - vbdy1
1780 : write(*,2) 'vbdy1 - v(j,k+1)', k, vbdy1 - v(j,k+1)
1781 : write(*,'(A)')
1782 : write(*,2) 'v(j,k-1) - vbdy2', k, v(j,k-1) - vbdy2
1783 : write(*,2) 'vbdy2 - v(j,k+1)', k, vbdy2 - v(j,k+1)
1784 : write(*,'(A)')
1785 : write(*,2) 'sm1', k, sm1
1786 : write(*,2) 's00', k, s00
1787 : write(*,'(A)')
1788 : call mesa_error(__FILE__,__LINE__,'debug: get1_lpp')
1789 : end if
1790 160 : c2(k) = 0
1791 160 : if (abs(sm1) <= abs(s00)) then
1792 84 : c1(k) = sm1
1793 : else
1794 76 : c1(k) = s00
1795 : end if
1796 : end if
1797 :
1798 : contains
1799 :
1800 8266 : subroutine set_const
1801 8266 : c0(k) = v(j,k)
1802 8266 : c1(k) = 0
1803 8266 : if (quad) c2(k) = 0
1804 8266 : end subroutine set_const
1805 :
1806 : end subroutine get1_lpp
1807 :
1808 :
1809 11776 : subroutine get_xq_integral( &
1810 11776 : k_old_in, nz_old, xq_old, xq_outer, dq, &
1811 11776 : order, c0, c1, c2, integral, dbg, k_old_last, ierr)
1812 : ! integrate val(j,:) from xq_inner to xq_outer, with xq_inner = xq_outer + dq
1813 : integer, intent(in) :: k_old_in, nz_old
1814 : real(dp), intent(in) :: xq_old(:), xq_outer, dq
1815 : integer, intent(in) :: order ! 0, 1, 2
1816 : real(dp), intent(in), dimension(:) :: c0, c1, c2 ! coefficients
1817 : real(dp), intent(out) :: integral
1818 : logical, intent(in) :: dbg
1819 : integer, intent(out) :: k_old_last, ierr
1820 :
1821 : integer :: k, k_old
1822 11776 : real(dp) :: a, b, c, old_xq_inner, old_xq_outer, xq_inner, &
1823 11776 : xq_overlap_outer, xq_overlap_inner, dq1, sum_dqs, old_xq_mid, &
1824 11776 : v_overlap_outer, v_overlap_inner, dq_outer, dq_inner, avg_value
1825 :
1826 : include 'formats'
1827 :
1828 11776 : if (dbg) write(*,*)
1829 :
1830 11776 : ierr = 0
1831 11776 : k_old = k_old_in
1832 : ! move starting k_old if necessary
1833 0 : do
1834 11776 : if (k_old <= 1) exit
1835 11776 : if (xq_old(k_old) <= xq_outer) exit
1836 11776 : k_old = k_old - 1
1837 : end do
1838 11776 : xq_inner = xq_outer + dq
1839 11776 : old_xq_inner = xq_old(k_old)
1840 11776 : sum_dqs = 0
1841 11776 : integral = 0d0
1842 :
1843 35328 : do k = k_old, nz_old
1844 35328 : if (dq <= sum_dqs) exit
1845 35328 : old_xq_outer = old_xq_inner
1846 35328 : if (k == nz_old) then
1847 0 : old_xq_inner = 1
1848 : else
1849 35328 : old_xq_inner = xq_old(k+1)
1850 : end if
1851 35328 : xq_overlap_outer = max(xq_outer, old_xq_outer)
1852 35328 : xq_overlap_inner = min(xq_inner, old_xq_inner)
1853 :
1854 35328 : if (dbg) then
1855 0 : write(*,2) 'xq_overlap_outer', k, xq_overlap_outer
1856 0 : write(*,2) 'xq_outer', k, xq_outer
1857 0 : write(*,2) 'old_xq_outer', k, old_xq_outer
1858 0 : write(*,2) 'xq_overlap_inner', k, xq_overlap_inner
1859 0 : write(*,2) 'xq_inner', k, xq_inner
1860 0 : write(*,2) 'old_xq_inner', k, old_xq_inner
1861 : end if
1862 :
1863 35328 : if (sum_dqs == 0 .and. xq_overlap_outer == xq_outer .and. &
1864 : xq_overlap_inner == xq_inner) then
1865 : ! fully contained
1866 0 : xq_inner = xq_outer + dq
1867 0 : xq_overlap_inner = xq_inner
1868 0 : dq1 = dq
1869 35328 : else if (old_xq_inner >= xq_inner) then ! this is the last one
1870 11776 : dq1 = dq - sum_dqs
1871 : else
1872 23552 : dq1 = max(0d0, xq_overlap_inner-xq_overlap_outer)
1873 : end if
1874 35328 : sum_dqs = sum_dqs + dq1
1875 : ! interpolant f(dq) = a + b*dq + (c/2)*dq^2, with dq = q - q_midpoint
1876 35328 : a = c0(k)
1877 35328 : b = c1(k)
1878 35328 : if (order == 2) then
1879 35328 : c = c2(k)
1880 : else
1881 0 : c = 0
1882 : end if
1883 :
1884 35328 : if (dq1 == 0 .or. (b==0 .and. c==0)) then
1885 17203 : avg_value = a
1886 : else
1887 18125 : old_xq_mid = (old_xq_outer + old_xq_inner)/2
1888 18125 : dq_outer = old_xq_mid - xq_overlap_outer
1889 18125 : dq_inner = old_xq_mid - xq_overlap_inner
1890 :
1891 18125 : if (order == 0) then
1892 0 : avg_value = a
1893 18125 : else if (order == 1 .or. c == 0) then
1894 : ! use slope to estimate average value in the region being used
1895 18125 : if (dbg) write(*,*) 'use slope to estimate average value'
1896 18125 : v_overlap_outer = a + dq_outer*b
1897 18125 : v_overlap_inner = a + dq_inner*b
1898 18125 : avg_value = (v_overlap_outer + v_overlap_inner)/2
1899 : else ! use quadratic reconstruction
1900 0 : if (dbg) write(*,*) 'use quadratic reconstruction'
1901 : avg_value = &
1902 : a + b*(dq_inner + dq_outer)/2 + &
1903 0 : c*(dq_inner*dq_inner + dq_inner*dq_outer + dq_outer*dq_outer)/6
1904 : end if
1905 : end if
1906 35328 : integral = integral + dq1*avg_value
1907 35328 : if (dbg) then
1908 0 : write(*,2) 'a', k, a
1909 0 : write(*,2) 'b', k, b
1910 0 : write(*,2) 'c', k, c
1911 0 : write(*,2) 'dq1', k, dq1
1912 0 : write(*,2) 'avg_value', k, avg_value
1913 0 : write(*,2) 'integral', k, integral
1914 0 : write(*,'(A)')
1915 : end if
1916 35328 : k_old_last = k
1917 35328 : if (old_xq_inner >= xq_inner) exit ! this is the last one
1918 :
1919 : end do
1920 :
1921 11776 : if (dbg) write(*,1) 'integral/dq', integral/dq, integral, dq
1922 :
1923 11776 : end subroutine get_xq_integral
1924 :
1925 :
1926 0 : subroutine adjust_omega(s, nz, nz_old, comes_from, &
1927 0 : old_xq, new_xq, old_dq, new_dq, xh, old_j_rot, &
1928 0 : xout_old, xout_new, old_dqbar, new_dqbar, ierr)
1929 : use alloc
1930 : type (star_info), pointer :: s
1931 : integer, intent(in) :: nz, nz_old
1932 : integer, dimension(:) :: comes_from
1933 : real(dp), dimension(:) :: &
1934 : old_xq, new_xq, old_dq, new_dq, old_j_rot, &
1935 : xout_old, xout_new, old_dqbar, new_dqbar
1936 : real(dp), intent(in) :: xh(:,:)
1937 : integer, intent(out) :: ierr
1938 : integer :: k, op_err
1939 : include 'formats'
1940 0 : ierr = 0
1941 :
1942 0 : !$OMP PARALLEL DO PRIVATE(k, op_err) SCHEDULE(dynamic,2)
1943 : do k = 1, nz
1944 : op_err = 0
1945 : call adjust1_omega(s, k, nz, nz_old, comes_from, &
1946 : xout_old, xout_new, old_dqbar, new_dqbar, old_j_rot, xh, op_err)
1947 : if (op_err /= 0) ierr = op_err
1948 : end do
1949 : !$OMP END PARALLEL DO
1950 :
1951 0 : end subroutine adjust_omega
1952 :
1953 :
1954 0 : subroutine adjust1_omega(s, k, nz, nz_old, comes_from, &
1955 0 : xout_old, xout_new, old_dqbar, new_dqbar, old_j_rot, xh, ierr)
1956 0 : use hydro_rotation, only: w_div_w_roche_jrot, update1_i_rot_from_xh
1957 : ! set new value for s% omega(k)
1958 : type (star_info), pointer :: s
1959 : integer, intent(in) :: k, nz, nz_old
1960 : integer, dimension(:) :: comes_from
1961 : real(dp), dimension(:), intent(in) :: &
1962 : xout_old, xout_new, old_dqbar, new_dqbar, old_j_rot
1963 : real(dp), intent(in) :: xh(:,:)
1964 : integer, intent(out) :: ierr
1965 :
1966 0 : real(dp) :: xq_outer, xq_inner, j_tot, xq0, xq1, new_point_dqbar, dq_sum, dq, r00
1967 : integer :: kk, k_outer
1968 :
1969 : integer, parameter :: k_dbg = -1
1970 :
1971 : include 'formats'
1972 :
1973 0 : ierr = 0
1974 0 : xq_outer = xout_new(k)
1975 0 : new_point_dqbar = new_dqbar(k)
1976 0 : if (k < nz) then
1977 0 : xq_inner = xq_outer + new_point_dqbar
1978 : else
1979 0 : xq_inner = 1d0
1980 : end if
1981 :
1982 0 : if (k == k_dbg) then
1983 0 : write(*,2) 'xq_outer', k, xq_outer
1984 0 : write(*,2) 'xq_inner', k, xq_inner
1985 0 : write(*,2) 'new_point_dqbar', k, new_point_dqbar
1986 : end if
1987 :
1988 0 : dq_sum = 0d0
1989 0 : j_tot = 0
1990 0 : if (xq_outer >= xout_old(nz_old)) then
1991 : ! new contained entirely in old center zone
1992 0 : k_outer = nz_old
1993 0 : if (k == k_dbg) &
1994 0 : write(*,2) 'new contained in old center', &
1995 0 : k_outer, xout_old(k_outer)
1996 0 : else if (k == 1) then
1997 0 : k_outer = 1
1998 : else
1999 0 : k_outer = comes_from(k-1)
2000 : end if
2001 :
2002 0 : do kk = k_outer, nz_old ! loop until reach m_inner
2003 :
2004 0 : if (kk == nz_old) then
2005 0 : xq1 = 1d0
2006 : else
2007 0 : xq1 = xout_old(kk+1)
2008 : end if
2009 0 : if (xq1 <= xq_outer) cycle
2010 :
2011 0 : xq0 = xout_old(kk)
2012 0 : if (xq0 >= xq_inner) then
2013 0 : if (dq_sum < new_point_dqbar .and. kk > 1) then
2014 : ! need to add a bit more from the previous source
2015 0 : dq = new_point_dqbar - dq_sum
2016 0 : dq_sum = new_point_dqbar
2017 0 : j_tot = j_tot + old_j_rot(kk-1)*dq
2018 : end if
2019 : exit
2020 : end if
2021 :
2022 0 : if (xq1 < xq_outer) then
2023 0 : ierr = -1
2024 0 : return
2025 : end if
2026 :
2027 0 : if (xq0 >= xq_outer .and. xq1 <= xq_inner) then ! entire old kk is in new k
2028 :
2029 0 : dq = old_dqbar(kk)
2030 0 : dq_sum = dq_sum + dq
2031 :
2032 0 : if (dq_sum > new_point_dqbar) then
2033 : ! dq too large -- numerical roundoff problems
2034 0 : dq = dq - (new_point_dqbar - dq_sum)
2035 0 : dq_sum = new_point_dqbar
2036 : end if
2037 :
2038 0 : j_tot = j_tot + old_j_rot(kk)*dq
2039 :
2040 0 : else if (xq0 <= xq_outer .and. xq1 >= xq_inner) then ! entire new k is in old kk
2041 :
2042 0 : dq = new_dqbar(k)
2043 0 : dq_sum = dq_sum + dq
2044 0 : j_tot = j_tot + old_j_rot(kk)*dq
2045 :
2046 : else ! only use the part of old kk that is in new k
2047 :
2048 0 : if (k == k_dbg) then
2049 0 : write(*,*) 'only use the part of old kk that is in new k', xq_inner <= xq1
2050 0 : write(*,1) 'xq_outer', xq_outer
2051 0 : write(*,1) 'xq_inner', xq_inner
2052 0 : write(*,1) 'xq0', xq0
2053 0 : write(*,1) 'xq1', xq1
2054 0 : write(*,1) 'dq_sum', dq_sum
2055 0 : write(*,1) 'new_point_dqbar', new_point_dqbar
2056 0 : write(*,1) 'new_point_dqbar - dq_sum', new_point_dqbar - dq_sum
2057 : end if
2058 :
2059 0 : if (xq_inner <= xq1) then ! this is the last part of new k
2060 :
2061 0 : if (k == k_dbg) write(*,3) 'this is the last part of new k', k, kk
2062 :
2063 0 : dq = new_point_dqbar - dq_sum
2064 0 : dq_sum = new_point_dqbar
2065 :
2066 : else ! we avoid this case if possible because of numerical roundoff
2067 :
2068 0 : if (k == k_dbg) write(*,3) 'we avoid this case if possible', k, kk
2069 :
2070 0 : dq = max(0d0, xq1 - xq_outer)
2071 0 : if (dq_sum + dq > new_point_dqbar) dq = new_point_dqbar - dq_sum
2072 0 : dq_sum = dq_sum + dq
2073 :
2074 : end if
2075 :
2076 0 : j_tot = j_tot + old_j_rot(kk)*dq
2077 :
2078 0 : if (dq <= 0) then
2079 0 : ierr = -1
2080 0 : return
2081 : end if
2082 :
2083 : end if
2084 :
2085 0 : if (dq_sum >= new_point_dqbar) then
2086 0 : if (k == k_dbg) then
2087 0 : write(*,2) 'exit for k', k
2088 0 : write(*,2) 'dq_sum', kk, dq_sum
2089 0 : write(*,2) 'new_point_dqbar', kk, new_point_dqbar
2090 : end if
2091 : exit
2092 : end if
2093 :
2094 : end do
2095 :
2096 0 : s% j_rot(k) = j_tot/dq_sum
2097 0 : r00 = get_r_from_xh(s,k)
2098 : s% w_div_w_crit_roche(k) = &
2099 : w_div_w_roche_jrot(r00,s% m(k),s% j_rot(k),s% cgrav(k), &
2100 0 : s% w_div_wcrit_max, s% w_div_wcrit_max2, s% w_div_wc_flag)
2101 0 : call update1_i_rot_from_xh(s, k)
2102 0 : s% omega(k) = s% j_rot(k)/s% i_rot(k)% val
2103 :
2104 0 : if (k_dbg == k) then
2105 0 : write(*,2) 's% omega(k)', k, s% omega(k)
2106 0 : write(*,2) 's% j_rot(k)', k, s% j_rot(k)
2107 0 : write(*,2) 's% i_rot(k)', k, s% i_rot(k)
2108 0 : if (s% model_number > 1925) call mesa_error(__FILE__,__LINE__,'debugging: adjust1_omega')
2109 : end if
2110 :
2111 0 : end subroutine adjust1_omega
2112 :
2113 :
2114 : ! like adjust_omega. conserve kinetic energy
2115 0 : subroutine do_v( &
2116 0 : s, nz, nz_old, cell_type, comes_from, &
2117 0 : old_xq, new_xq, old_dq, new_dq, xh, xh_old, &
2118 0 : xout_old, xout_new, old_dqbar, new_dqbar, old_ke, ierr)
2119 0 : use alloc
2120 : type (star_info), pointer :: s
2121 : integer, intent(in) :: nz, nz_old
2122 : integer, dimension(:) :: cell_type, comes_from
2123 : real(dp), dimension(:) :: &
2124 : xout_old, xout_new, old_dqbar, new_dqbar, &
2125 : old_xq, new_xq, old_dq, new_dq, old_ke
2126 : real(dp), dimension(:,:) :: xh, xh_old
2127 : integer, intent(out) :: ierr
2128 :
2129 : integer :: k, op_err, i_v
2130 0 : real(dp) :: old_ke_tot, new_ke_tot, xmstar, err
2131 :
2132 : include 'formats'
2133 0 : ierr = 0
2134 0 : i_v = s% i_v
2135 0 : xmstar = s% xmstar
2136 :
2137 0 : old_ke_tot = 0d0
2138 0 : do k=1,nz_old ! skip common factor 1/2 xmstar in ke
2139 0 : old_ke(k) = old_dqbar(k)*xh_old(i_v,k)*xh_old(i_v,k)
2140 0 : old_ke_tot = old_ke_tot + old_ke(k)
2141 : end do
2142 :
2143 0 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
2144 : do k = 1, nz
2145 : op_err = 0
2146 : call adjust1_v( &
2147 : s, k, nz, nz_old, cell_type, comes_from, xout_old, xout_new, &
2148 : old_dqbar, new_dqbar, old_ke, i_v, xh, xh_old, op_err)
2149 : if (op_err /= 0) ierr = op_err
2150 : end do
2151 : !$OMP END PARALLEL DO
2152 0 : if (ierr /= 0) then
2153 : return
2154 : end if
2155 :
2156 0 : new_ke_tot = 0
2157 0 : do k=1,nz
2158 0 : new_ke_tot = new_ke_tot + new_dqbar(k)*xh(i_v,k)*xh(i_v,k)
2159 : end do
2160 :
2161 0 : if (abs(old_ke_tot - new_ke_tot) > 1d-7*new_ke_tot) then
2162 0 : ierr = -1
2163 0 : s% retry_message = 'failed in mesh_adjust do_v'
2164 0 : if (s% report_ierr) write(*, *) s% retry_message
2165 : !call mesa_error(__FILE__,__LINE__,'do_v')
2166 0 : return
2167 : end if
2168 :
2169 0 : err = abs(old_ke_tot - new_ke_tot)/max(new_ke_tot,old_ke_tot,1d0)
2170 0 : s% mesh_adjust_KE_conservation = err
2171 :
2172 0 : if (s% trace_mesh_adjust_error_in_conservation) then
2173 0 : write(*,2) 'mesh adjust error in conservation of KE', s% model_number, &
2174 0 : err, new_ke_tot, old_ke_tot
2175 0 : if (err > 1d-10) then
2176 0 : write(*,*) 'err too large'
2177 0 : call mesa_error(__FILE__,__LINE__,'do_v')
2178 : end if
2179 : end if
2180 :
2181 0 : end subroutine do_v
2182 :
2183 :
2184 0 : subroutine adjust1_v( &
2185 0 : s, k, nz, nz_old, cell_type, comes_from, xout_old, xout_new, &
2186 0 : old_dqbar, new_dqbar, old_ke, i_v, xh, xh_old, ierr)
2187 : ! set new value for s% v(k) to conserve kinetic energy
2188 : type (star_info), pointer :: s
2189 : integer, intent(in) :: k, nz, nz_old, i_v
2190 : integer, dimension(:) :: cell_type, comes_from
2191 : real(dp), dimension(:), intent(in) :: &
2192 : xout_old, xout_new, old_dqbar, new_dqbar, old_ke
2193 : real(dp), dimension(:,:) :: xh, xh_old
2194 : integer, intent(out) :: ierr
2195 :
2196 0 : real(dp) :: xq_outer, xq_inner, ke_sum, &
2197 0 : xq0, xq1, new_point_dqbar, dq_sum, dq
2198 : integer :: kk, k_outer
2199 :
2200 : integer, parameter :: k_dbg = -1
2201 :
2202 : include 'formats'
2203 :
2204 0 : ierr = 0
2205 :
2206 0 : if (cell_type(k) == unchanged_type .or. &
2207 : cell_type(k) == revised_type) then
2208 0 : if (k == 1) then
2209 0 : xh(i_v,k) = xh_old(i_v,comes_from(k))
2210 0 : return
2211 : end if
2212 0 : if (cell_type(k-1) == unchanged_type) then
2213 0 : xh(i_v,k) = xh_old(i_v,comes_from(k))
2214 0 : return
2215 : end if
2216 : end if
2217 :
2218 0 : xq_outer = xout_new(k)
2219 0 : new_point_dqbar = new_dqbar(k)
2220 0 : if (k < nz) then
2221 0 : xq_inner = xq_outer + new_point_dqbar
2222 : else
2223 0 : xq_inner = 1d0
2224 : end if
2225 :
2226 0 : if (k == k_dbg) then
2227 0 : write(*,2) 'xq_outer', k, xq_outer
2228 0 : write(*,2) 'xq_inner', k, xq_inner
2229 0 : write(*,2) 'new_point_dqbar', k, new_point_dqbar
2230 : end if
2231 :
2232 0 : dq_sum = 0d0
2233 0 : ke_sum = 0
2234 0 : if (xq_outer >= xout_old(nz_old)) then
2235 : ! new contained entirely in old center zone
2236 0 : k_outer = nz_old
2237 0 : if (k == k_dbg) &
2238 0 : write(*,2) 'new contained in old center', &
2239 0 : k_outer, xout_old(k_outer)
2240 0 : else if (k == 1) then
2241 0 : k_outer = 1
2242 : else
2243 0 : k_outer = comes_from(k-1)
2244 : end if
2245 :
2246 0 : do kk = k_outer, nz_old ! loop until reach xq_inner
2247 :
2248 0 : if (kk == nz_old) then
2249 0 : xq1 = 1d0
2250 : else
2251 0 : xq1 = xout_old(kk+1)
2252 : end if
2253 0 : if (xq1 <= xq_outer) cycle
2254 :
2255 0 : xq0 = xout_old(kk)
2256 0 : if (xq0 >= xq_inner) then
2257 0 : if (dq_sum < new_point_dqbar .and. kk > 1) then
2258 : ! need to add a bit more from the previous source
2259 0 : dq = new_point_dqbar - dq_sum
2260 0 : dq_sum = new_point_dqbar
2261 0 : ke_sum = ke_sum + old_ke(kk-1)*dq/old_dqbar(kk-1)
2262 :
2263 0 : if (k == k_dbg) &
2264 0 : write(*,3) 'new k contains some of old kk-1', &
2265 0 : k, kk, old_ke(kk-1)*dq, dq_sum
2266 :
2267 : end if
2268 : exit
2269 : end if
2270 :
2271 0 : if (xq1 < xq_outer) then
2272 0 : ierr = -1
2273 0 : return
2274 : end if
2275 :
2276 0 : if (xq0 >= xq_outer .and. xq1 <= xq_inner) then ! entire old kk is in new k
2277 :
2278 0 : dq = old_dqbar(kk)
2279 0 : dq_sum = dq_sum + dq
2280 :
2281 0 : if (dq_sum > new_point_dqbar) then
2282 : ! dq too large -- numerical roundoff problems
2283 0 : dq = dq - (new_point_dqbar - dq_sum)
2284 0 : dq_sum = new_point_dqbar
2285 : end if
2286 :
2287 0 : ke_sum = ke_sum + old_ke(kk)*dq/old_dqbar(kk)
2288 :
2289 0 : if (k == k_dbg) &
2290 0 : write(*,3) 'new k contains all of old kk', &
2291 0 : k, kk, old_ke(kk)*dq, ke_sum
2292 :
2293 0 : else if (xq0 <= xq_outer .and. xq1 >= xq_inner) then ! entire new k is in old kk
2294 :
2295 0 : dq = new_dqbar(k)
2296 0 : dq_sum = dq_sum + dq
2297 0 : ke_sum = ke_sum + old_ke(kk)*dq/old_dqbar(kk)
2298 :
2299 0 : if (k == k_dbg) &
2300 0 : write(*,3) 'all new k is in old kk', &
2301 0 : k, kk, old_ke(kk)*dq, ke_sum
2302 :
2303 : else ! only use the part of old kk that is in new k
2304 :
2305 0 : if (k == k_dbg) then
2306 0 : write(*,*) 'only use the part of old kk that is in new k', xq_inner <= xq1
2307 0 : write(*,1) 'xq_outer', xq_outer
2308 0 : write(*,1) 'xq_inner', xq_inner
2309 0 : write(*,1) 'xq0', xq0
2310 0 : write(*,1) 'xq1', xq1
2311 0 : write(*,1) 'dq_sum', dq_sum
2312 0 : write(*,1) 'new_point_dqbar', new_point_dqbar
2313 0 : write(*,1) 'new_point_dqbar - dq_sum', new_point_dqbar - dq_sum
2314 : end if
2315 :
2316 0 : if (xq_inner <= xq1) then ! this is the last part of new k
2317 :
2318 0 : dq = new_point_dqbar - dq_sum
2319 0 : dq_sum = new_point_dqbar
2320 :
2321 : else ! we avoid this case if possible because of numerical roundoff
2322 :
2323 0 : if (k == k_dbg) write(*,3) 'we avoid this case if possible', k, kk
2324 :
2325 0 : dq = max(0d0, xq1 - xq_outer)
2326 0 : if (dq_sum + dq > new_point_dqbar) dq = new_point_dqbar - dq_sum
2327 0 : dq_sum = dq_sum + dq
2328 :
2329 : end if
2330 :
2331 0 : if (k == k_dbg) then
2332 0 : write(*,3) 'new k use only part of old kk', k, kk
2333 0 : write(*,2) 'dq_sum', k, dq_sum
2334 0 : write(*,2) 'dq', k, dq
2335 0 : write(*,2) 'old_ke(kk)', kk, old_ke(kk)
2336 0 : write(*,2) 'old ke_sum', k, ke_sum
2337 0 : write(*,2) 'new ke_sum', k, ke_sum + old_ke(kk)*dq
2338 : end if
2339 :
2340 0 : ke_sum = ke_sum + old_ke(kk)*dq/old_dqbar(kk)
2341 :
2342 0 : if (dq <= 0) then
2343 0 : ierr = -1
2344 : !return
2345 0 : write(*,*) 'dq <= 0', dq
2346 0 : call mesa_error(__FILE__,__LINE__,'debugging: adjust1_v')
2347 : end if
2348 :
2349 : end if
2350 :
2351 0 : if (dq_sum >= new_point_dqbar) then
2352 : exit
2353 : end if
2354 :
2355 : end do
2356 :
2357 0 : xh(i_v,k) = sqrt(ke_sum/new_point_dqbar) ! we have skipped the 1/2 xmstar factor
2358 0 : if (xh_old(i_v,comes_from(k)) < 0d0) xh(i_v,k) = -xh(i_v,k)
2359 :
2360 0 : if (k == k_dbg) then
2361 0 : !$OMP critical (adjust1_v_dbg)
2362 0 : write(*,2) 'xh(i_v,k) new_dqbar', k, xh(i_v,k), new_dqbar(k)
2363 0 : write(*,2) 'xh_old(i_v,comes_from(k)) old_dqbar', &
2364 0 : comes_from(k), xh_old(i_v,comes_from(k)), old_dqbar(comes_from(k))
2365 0 : if (k == k_dbg) call mesa_error(__FILE__,__LINE__,'adjust1_v')
2366 : !$OMP end critical (adjust1_v_dbg)
2367 : !stop
2368 : end if
2369 :
2370 0 : end subroutine adjust1_v
2371 :
2372 :
2373 0 : subroutine do_u( &
2374 0 : s, nz, nz_old, cell_type, comes_from, &
2375 0 : old_xq, new_xq, old_dq, new_dq, xh, xh_old, &
2376 0 : xout_old, xout_new, old_ke, ierr)
2377 : use alloc
2378 : type (star_info), pointer :: s
2379 : integer, intent(in) :: nz, nz_old
2380 : integer, dimension(:) :: cell_type, comes_from
2381 : real(dp), dimension(:) :: &
2382 : xout_old, xout_new, old_xq, new_xq, old_dq, new_dq, old_ke
2383 : real(dp), dimension(:,:) :: xh, xh_old
2384 : integer, intent(out) :: ierr
2385 :
2386 : integer :: k, op_err, i_u
2387 0 : real(dp) :: old_ke_tot, new_ke_tot, xmstar, err
2388 :
2389 : include 'formats'
2390 0 : ierr = 0
2391 0 : i_u = s% i_u
2392 0 : xmstar = s% xmstar
2393 :
2394 0 : old_ke_tot = 0d0
2395 0 : do k=1,nz_old ! skip common factor 1/2 xmstar in ke
2396 0 : old_ke(k) = old_dq(k)*xh_old(i_u,k)*xh_old(i_u,k)
2397 0 : old_ke_tot = old_ke_tot + old_ke(k)
2398 : end do
2399 :
2400 0 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
2401 : do k = 1, nz
2402 : op_err = 0
2403 : call adjust1_u( &
2404 : s, k, nz, nz_old, cell_type, comes_from, xout_old, xout_new, &
2405 : old_dq, new_dq, old_ke, i_u, xh, xh_old, op_err)
2406 : if (op_err /= 0) ierr = op_err
2407 : end do
2408 : !$OMP END PARALLEL DO
2409 0 : if (ierr /= 0) then
2410 : return
2411 : end if
2412 :
2413 0 : new_ke_tot = 0
2414 0 : do k=1,nz
2415 0 : new_ke_tot = new_ke_tot + new_dq(k)*xh(i_u,k)*xh(i_u,k)
2416 : end do
2417 :
2418 0 : err = abs(old_ke_tot - new_ke_tot)/max(new_ke_tot,old_ke_tot,1d0)
2419 0 : s% mesh_adjust_KE_conservation = err
2420 :
2421 0 : if (s% trace_mesh_adjust_error_in_conservation) then
2422 0 : write(*,2) 'mesh adjust error in conservation of KE', s% model_number, &
2423 0 : err, new_ke_tot, old_ke_tot
2424 0 : if (err > 1d-10) then
2425 0 : write(*,*) 'err too large'
2426 0 : call mesa_error(__FILE__,__LINE__,'do_u')
2427 : end if
2428 : end if
2429 :
2430 0 : end subroutine do_u
2431 :
2432 :
2433 0 : subroutine adjust1_u( &
2434 0 : s, k, nz, nz_old, cell_type, comes_from, xout_old, xout_new, &
2435 0 : old_dq, new_dq, old_ke, i_u, xh, xh_old, ierr)
2436 : ! set new value for s% u(k) to conserve kinetic energy
2437 : type (star_info), pointer :: s
2438 : integer, intent(in) :: k, nz, nz_old, i_u
2439 : integer, dimension(:) :: cell_type, comes_from
2440 : real(dp), dimension(:), intent(in) :: &
2441 : xout_old, xout_new, old_dq, new_dq, old_ke
2442 : real(dp), dimension(:,:) :: xh, xh_old
2443 : integer, intent(out) :: ierr
2444 :
2445 0 : real(dp) :: xq_outer, xq_inner, ke_sum, &
2446 0 : xq0, xq1, new_cell_dq, dq_sum, dq
2447 : integer :: kk, k_outer
2448 :
2449 : integer, parameter :: k_dbg = -1
2450 :
2451 : include 'formats'
2452 :
2453 0 : ierr = 0
2454 :
2455 0 : if (cell_type(k) == unchanged_type .or. &
2456 : cell_type(k) == revised_type) then
2457 0 : if (k == 1) then
2458 0 : xh(i_u,k) = xh_old(i_u,comes_from(k))
2459 0 : return
2460 : end if
2461 0 : if (cell_type(k-1) == unchanged_type) then
2462 0 : xh(i_u,k) = xh_old(i_u,comes_from(k))
2463 0 : return
2464 : end if
2465 : end if
2466 :
2467 0 : xq_outer = xout_new(k)
2468 0 : new_cell_dq = new_dq(k)
2469 0 : if (k < nz) then
2470 0 : xq_inner = xq_outer + new_cell_dq
2471 : else
2472 0 : xq_inner = 1d0
2473 : end if
2474 :
2475 0 : if (k == k_dbg) then
2476 0 : write(*,2) 'xq_outer', k, xq_outer
2477 0 : write(*,2) 'xq_inner', k, xq_inner
2478 0 : write(*,2) 'new_cell_dq', k, new_cell_dq
2479 : end if
2480 :
2481 0 : dq_sum = 0d0
2482 0 : ke_sum = 0
2483 0 : if (xq_outer >= xout_old(nz_old)) then
2484 : ! new contained entirely in old center zone
2485 0 : k_outer = nz_old
2486 0 : if (k == k_dbg) &
2487 0 : write(*,2) 'new contained in old center', &
2488 0 : k_outer, xout_old(k_outer)
2489 0 : else if (k == 1) then
2490 0 : k_outer = 1
2491 : else
2492 0 : k_outer = comes_from(k-1)
2493 : end if
2494 :
2495 0 : do kk = k_outer, nz_old ! loop until reach xq_inner
2496 :
2497 0 : if (kk == nz_old) then
2498 0 : xq1 = 1d0
2499 : else
2500 0 : xq1 = xout_old(kk+1)
2501 : end if
2502 0 : if (xq1 <= xq_outer) cycle
2503 :
2504 0 : if (xq1 < xq_outer) then
2505 0 : ierr = -1
2506 0 : return
2507 : end if
2508 :
2509 0 : xq0 = xout_old(kk)
2510 0 : if (xq0 >= xq_outer .and. xq1 <= xq_inner) then ! entire old kk is in new k
2511 :
2512 0 : dq = old_dq(kk)
2513 0 : dq_sum = dq_sum + dq
2514 :
2515 0 : if (dq_sum > new_cell_dq) then
2516 : ! dq too large -- numerical roundoff problems
2517 0 : dq = dq - (new_cell_dq - dq_sum)
2518 0 : dq_sum = new_cell_dq
2519 : end if
2520 :
2521 0 : ke_sum = ke_sum + old_ke(kk)*dq/old_dq(kk)
2522 :
2523 0 : if (k == k_dbg) &
2524 0 : write(*,3) 'new k contains all of old kk', &
2525 0 : k, kk, old_ke(kk)*dq, ke_sum
2526 :
2527 0 : else if (xq0 <= xq_outer .and. xq1 >= xq_inner) then ! entire new k is in old kk
2528 :
2529 0 : dq = new_dq(k)
2530 0 : dq_sum = dq_sum + dq
2531 0 : ke_sum = ke_sum + old_ke(kk)*dq/old_dq(kk)
2532 :
2533 0 : if (k == k_dbg) &
2534 0 : write(*,3) 'all new k is in old kk', &
2535 0 : k, kk, old_ke(kk)*dq, ke_sum
2536 :
2537 : else ! only use the part of old kk that is in new k
2538 :
2539 0 : if (k == k_dbg) then
2540 0 : write(*,*) 'only use the part of old kk that is in new k', xq_inner <= xq1
2541 0 : write(*,1) 'xq_outer', xq_outer
2542 0 : write(*,1) 'xq_inner', xq_inner
2543 0 : write(*,1) 'xq0', xq0
2544 0 : write(*,1) 'xq1', xq1
2545 0 : write(*,1) 'dq_sum', dq_sum
2546 0 : write(*,1) 'new_cell_dq', new_cell_dq
2547 0 : write(*,1) 'new_cell_dq - dq_sum', new_cell_dq - dq_sum
2548 : end if
2549 :
2550 0 : if (xq_inner <= xq1) then ! this is the last part of new k
2551 :
2552 0 : dq = new_cell_dq - dq_sum
2553 0 : dq_sum = new_cell_dq
2554 :
2555 : else ! we avoid this case if possible because of numerical roundoff
2556 :
2557 0 : if (k == k_dbg) write(*,3) 'we avoid this case if possible', k, kk
2558 :
2559 0 : dq = max(0d0, xq1 - xq_outer)
2560 0 : if (dq_sum + dq > new_cell_dq) dq = new_cell_dq - dq_sum
2561 0 : dq_sum = dq_sum + dq
2562 :
2563 : end if
2564 :
2565 0 : if (k == k_dbg) then
2566 0 : write(*,3) 'new k use only part of old kk', k, kk
2567 0 : write(*,2) 'dq_sum', k, dq_sum
2568 0 : write(*,2) 'dq', k, dq
2569 0 : write(*,2) 'old_ke(kk)', kk, old_ke(kk)
2570 0 : write(*,2) 'old ke_sum', k, ke_sum
2571 0 : write(*,2) 'new ke_sum', k, ke_sum + old_ke(kk)*dq
2572 : end if
2573 :
2574 0 : ke_sum = ke_sum + old_ke(kk)*dq/old_dq(kk)
2575 :
2576 0 : if (dq <= 0) then
2577 0 : ierr = -1
2578 : !return
2579 0 : write(*,*) 'dq <= 0', dq
2580 0 : call mesa_error(__FILE__,__LINE__,'debugging: adjust1_u')
2581 : end if
2582 :
2583 : end if
2584 :
2585 0 : if (dq_sum >= new_cell_dq) then
2586 : exit
2587 : end if
2588 :
2589 : end do
2590 :
2591 0 : xh(i_u,k) = sqrt(ke_sum/new_cell_dq) ! we have skipped the 1/2 xmstar factor
2592 0 : if (xh_old(i_u,comes_from(k)) < 0d0) xh(i_u,k) = -xh(i_u,k)
2593 :
2594 0 : if (k == k_dbg) then
2595 0 : !$OMP critical (adjust1_u_dbg)
2596 0 : write(*,2) 'xh(i_u,k) new_dq', k, xh(i_u,k), new_dq(k)
2597 0 : write(*,2) 'xh_old(i_u,comes_from(k)) old_dq', &
2598 0 : comes_from(k), xh_old(i_u,comes_from(k)), old_dq(comes_from(k))
2599 0 : if (k == k_dbg) call mesa_error(__FILE__,__LINE__,'adjust1_u')
2600 : !$OMP end critical (adjust1_u_dbg)
2601 : !stop
2602 : end if
2603 :
2604 0 : end subroutine adjust1_u
2605 :
2606 :
2607 0 : subroutine do_Hp_face( &
2608 : s, nz, nz_old, nzlo, nzhi, comes_from, xh, xh_old, &
2609 0 : xq, xq_old_plus1, xq_new, work, Hp_face_old_plus1, Hp_face_new, ierr)
2610 : use interp_1d_def
2611 : use interp_1d_lib
2612 : type (star_info), pointer :: s
2613 : integer, intent(in) :: nz, nz_old, nzlo, nzhi, comes_from(:)
2614 : real(dp), dimension(:,:), pointer :: xh, xh_old
2615 : real(dp), dimension(:), pointer :: work
2616 : real(dp), dimension(:) :: &
2617 : xq, xq_old_plus1, Hp_face_old_plus1, Hp_face_new, xq_new
2618 : integer, intent(out) :: ierr
2619 :
2620 : integer :: n, i_Hp, k
2621 :
2622 : include 'formats'
2623 :
2624 0 : ierr = 0
2625 0 : i_Hp = s% i_Hp
2626 0 : if (i_Hp == 0) return
2627 0 : n = nzhi - nzlo + 1
2628 :
2629 0 : do k=1,nz_old
2630 0 : Hp_face_old_plus1(k) = xh_old(i_Hp,k)
2631 : end do
2632 0 : Hp_face_old_plus1(nz_old+1) = Hp_face_old_plus1(nz_old)
2633 :
2634 : call interpolate_vector( &
2635 : nz_old+1, xq_old_plus1, n, xq_new, &
2636 : Hp_face_old_plus1, Hp_face_new, interp_pm, nwork, work, &
2637 0 : 'mesh_adjust do_Hp_face', ierr)
2638 0 : if (ierr /= 0) then
2639 : return
2640 : write(*,*) 'interpolate_vector failed in do_Hp_face for remesh'
2641 : call mesa_error(__FILE__,__LINE__,'debug: mesh adjust: do_Hp_face')
2642 : end if
2643 :
2644 0 : do k=nzlo,nzhi
2645 0 : xh(i_Hp,k) = Hp_face_new(k+1-nzlo)
2646 : end do
2647 :
2648 0 : n = nzlo - 1
2649 0 : if (n > 0) then
2650 0 : do k=1,n
2651 0 : xh(i_Hp,k) = xh_old(i_Hp,k)
2652 : end do
2653 : end if
2654 :
2655 0 : if (nzhi < nz) then
2656 0 : n = nz - nzhi - 1 ! nz-n = nzhi+1
2657 0 : do k=0,n
2658 0 : xh(i_Hp,nz-k) = xh_old(i_Hp,nz_old-k)
2659 : end do
2660 : end if
2661 :
2662 0 : end subroutine do_Hp_face
2663 :
2664 :
2665 0 : subroutine do_etrb( & ! same logic as do_u
2666 0 : s, nz, nz_old, cell_type, comes_from, &
2667 0 : old_xq, new_xq, old_dq, new_dq, xh, xh_old, &
2668 0 : xout_old, xout_new, old_eturb, ierr)
2669 0 : use alloc
2670 : type (star_info), pointer :: s
2671 : integer, intent(in) :: nz, nz_old
2672 : integer, dimension(:) :: cell_type, comes_from
2673 : real(dp), dimension(:) :: &
2674 : xout_old, xout_new, old_xq, new_xq, old_dq, new_dq, old_eturb
2675 : real(dp), dimension(:,:) :: xh, xh_old
2676 : integer, intent(out) :: ierr
2677 :
2678 : integer :: k, op_err, i_w
2679 0 : real(dp) :: old_eturb_tot, new_eturb_tot, xmstar, err
2680 :
2681 : include 'formats'
2682 0 : ierr = 0
2683 0 : i_w = s% i_w
2684 0 : xmstar = s% xmstar
2685 :
2686 0 : old_eturb_tot = 0d0
2687 0 : do k=1,nz_old
2688 0 : old_eturb(k) = old_dq(k)*pow2(xh_old(i_w,k))
2689 0 : old_eturb_tot = old_eturb_tot + old_eturb(k)
2690 : end do
2691 :
2692 0 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
2693 : do k = 1, nz
2694 : op_err = 0
2695 : call adjust1_etrb( &
2696 : s, k, nz, nz_old, cell_type, comes_from, xout_old, xout_new, &
2697 : old_dq, new_dq, old_eturb, i_w, xh, xh_old, op_err)
2698 : if (op_err /= 0) ierr = op_err
2699 : end do
2700 : !$OMP END PARALLEL DO
2701 0 : if (ierr /= 0) then
2702 : return
2703 : end if
2704 :
2705 0 : new_eturb_tot = 0
2706 0 : do k=1,nz
2707 0 : new_eturb_tot = new_eturb_tot + new_dq(k)*pow2(xh(i_w,k))
2708 : end do
2709 :
2710 0 : err = abs(old_eturb_tot - new_eturb_tot)/max(new_eturb_tot,old_eturb_tot,1d0)
2711 0 : s% mesh_adjust_Eturb_conservation = err
2712 :
2713 0 : if (s% trace_mesh_adjust_error_in_conservation) then
2714 0 : write(*,2) 'mesh adjust error in conservation of turbulent energy', s% model_number, &
2715 0 : err, new_eturb_tot, old_eturb_tot
2716 0 : if (err > 1d-10) then
2717 0 : write(*,*) 'err too large'
2718 0 : call mesa_error(__FILE__,__LINE__,'do_etrb')
2719 : end if
2720 : end if
2721 :
2722 0 : end subroutine do_etrb
2723 :
2724 :
2725 0 : subroutine adjust1_etrb( &
2726 0 : s, k, nz, nz_old, cell_type, comes_from, xout_old, xout_new, &
2727 0 : old_dq, new_dq, old_eturb, i_w, xh, xh_old, ierr)
2728 : ! set new value for s% w(k) to conserve turbulent energy
2729 : type (star_info), pointer :: s
2730 : integer, intent(in) :: k, nz, nz_old, i_w
2731 : integer, dimension(:) :: cell_type, comes_from
2732 : real(dp), dimension(:), intent(in) :: &
2733 : xout_old, xout_new, old_dq, new_dq, old_eturb
2734 : real(dp), dimension(:,:) :: xh, xh_old
2735 : integer, intent(out) :: ierr
2736 :
2737 0 : real(dp) :: xq_outer, xq_inner, eturb_sum, &
2738 0 : xq0, xq1, new_cell_dq, dq_sum, dq
2739 : integer :: kk, k_outer
2740 :
2741 : integer, parameter :: k_dbg = -1
2742 :
2743 : include 'formats'
2744 :
2745 0 : ierr = 0
2746 :
2747 0 : if (cell_type(k) == unchanged_type .or. &
2748 : cell_type(k) == revised_type) then
2749 : ! just copy the old value
2750 0 : if (k == 1) then
2751 0 : xh(i_w,k) = xh_old(i_w,comes_from(k))
2752 0 : return
2753 : end if
2754 0 : if (cell_type(k-1) == unchanged_type) then
2755 0 : xh(i_w,k) = xh_old(i_w,comes_from(k))
2756 0 : return
2757 : end if
2758 : end if
2759 :
2760 0 : xq_outer = xout_new(k)
2761 0 : new_cell_dq = new_dq(k)
2762 0 : if (k < nz) then
2763 0 : xq_inner = xq_outer + new_cell_dq
2764 : else
2765 0 : xq_inner = 1d0
2766 : end if
2767 :
2768 0 : if (k == k_dbg) then
2769 0 : write(*,2) 'xq_outer', k, xq_outer
2770 0 : write(*,2) 'xq_inner', k, xq_inner
2771 0 : write(*,2) 'new_cell_dq', k, new_cell_dq
2772 : end if
2773 :
2774 0 : dq_sum = 0d0
2775 0 : eturb_sum = 0
2776 0 : if (xq_outer >= xout_old(nz_old)) then
2777 : ! new contained entirely in old center zone
2778 0 : k_outer = nz_old
2779 0 : if (k == k_dbg) &
2780 0 : write(*,2) 'new contained in old center', &
2781 0 : k_outer, xout_old(k_outer)
2782 0 : else if (k == 1) then
2783 0 : k_outer = 1
2784 : else
2785 0 : k_outer = comes_from(k-1)
2786 : end if
2787 :
2788 0 : do kk = k_outer, nz_old ! loop until reach xq_inner
2789 :
2790 0 : if (kk == nz_old) then
2791 0 : xq1 = 1d0
2792 : else
2793 0 : xq1 = xout_old(kk+1)
2794 : end if
2795 0 : if (xq1 <= xq_outer) cycle
2796 :
2797 0 : if (xq1 < xq_outer) then
2798 0 : ierr = -1
2799 0 : return
2800 : end if
2801 :
2802 0 : xq0 = xout_old(kk)
2803 0 : if (xq0 >= xq_outer .and. xq1 <= xq_inner) then ! entire old kk is in new k
2804 :
2805 0 : dq = old_dq(kk)
2806 0 : dq_sum = dq_sum + dq
2807 :
2808 0 : if (dq_sum > new_cell_dq) then
2809 : ! dq too large -- numerical roundoff problems
2810 0 : dq = dq - (new_cell_dq - dq_sum)
2811 0 : dq_sum = new_cell_dq
2812 : end if
2813 :
2814 0 : eturb_sum = eturb_sum + old_eturb(kk)*dq/old_dq(kk)
2815 :
2816 0 : if (k == k_dbg) &
2817 0 : write(*,3) 'new k contains all of old kk', &
2818 0 : k, kk, old_eturb(kk)*dq, eturb_sum
2819 :
2820 0 : else if (xq0 <= xq_outer .and. xq1 >= xq_inner) then ! entire new k is in old kk
2821 :
2822 0 : dq = new_dq(k)
2823 0 : dq_sum = dq_sum + dq
2824 0 : eturb_sum = eturb_sum + old_eturb(kk)*dq/old_dq(kk)
2825 :
2826 0 : if (k == k_dbg) &
2827 0 : write(*,3) 'all new k is in old kk', &
2828 0 : k, kk, old_eturb(kk)*dq, eturb_sum
2829 :
2830 : else ! only use the part of old kk that is in new k
2831 :
2832 0 : if (k == k_dbg) then
2833 0 : write(*,*) 'only use the part of old kk that is in new k', xq_inner <= xq1
2834 0 : write(*,1) 'xq_outer', xq_outer
2835 0 : write(*,1) 'xq_inner', xq_inner
2836 0 : write(*,1) 'xq0', xq0
2837 0 : write(*,1) 'xq1', xq1
2838 0 : write(*,1) 'dq_sum', dq_sum
2839 0 : write(*,1) 'new_cell_dq', new_cell_dq
2840 0 : write(*,1) 'new_cell_dq - dq_sum', new_cell_dq - dq_sum
2841 : end if
2842 :
2843 0 : if (xq_inner <= xq1) then ! this is the last part of new k
2844 :
2845 0 : dq = new_cell_dq - dq_sum
2846 0 : dq_sum = new_cell_dq
2847 :
2848 : else ! we avoid this case if possible because of numerical roundoff
2849 :
2850 0 : if (k == k_dbg) write(*,3) 'we avoid this case if possible', k, kk
2851 :
2852 0 : dq = max(0d0, xq1 - xq_outer)
2853 0 : if (dq_sum + dq > new_cell_dq) dq = new_cell_dq - dq_sum
2854 0 : dq_sum = dq_sum + dq
2855 :
2856 : end if
2857 :
2858 0 : if (k == k_dbg) then
2859 0 : write(*,3) 'new k use only part of old kk', k, kk
2860 0 : write(*,2) 'dq_sum', k, dq_sum
2861 0 : write(*,2) 'dq', k, dq
2862 0 : write(*,2) 'old_eturb(kk)', kk, old_eturb(kk)
2863 0 : write(*,2) 'old eturb_sum', k, eturb_sum
2864 0 : write(*,2) 'new eturb_sum', k, eturb_sum + old_eturb(kk)*dq
2865 : end if
2866 :
2867 0 : eturb_sum = eturb_sum + old_eturb(kk)*dq/old_dq(kk)
2868 :
2869 0 : if (dq <= 0) then
2870 0 : ierr = -1
2871 : !return
2872 0 : write(*,*) 'dq <= 0', dq
2873 0 : call mesa_error(__FILE__,__LINE__,'debugging: adjust1_etrb')
2874 : end if
2875 :
2876 : end if
2877 :
2878 0 : if (dq_sum >= new_cell_dq) then
2879 : exit
2880 : end if
2881 :
2882 : end do
2883 :
2884 0 : xh(i_w,k) = sqrt(max(0d0,eturb_sum/new_cell_dq))
2885 :
2886 0 : end subroutine adjust1_etrb
2887 :
2888 : end module mesh_adjust
|
