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