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 : real(dp) :: dxa, xmstar, mstar, sumx, &
72 : 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, omega_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 14130 : j = maxloc(xa_c1(k,:), dim=1)
274 : else
275 16893 : 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 : 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 : write(message,*) 'do1_lnT for k', k
348 0 : ierr = op_err
349 0 : call mesa_error(__FILE__,__LINE__,'mesh adjust: do1_lnT')
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 : call mesa_error(__FILE__,__LINE__,'mesh adjust: bad energy_new or dq')
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 : 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 0 : 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 260 : 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 : if (ierr /= 0) return
453 : 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 : 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 : 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 0 : 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 : 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 : 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 0 : 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 : 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 : 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 0 : 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 : 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 : 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 20 : s, nz, nz_old, nzlo, nzhi, cell_type, comes_from, &
966 10 : xh, xh_old, xmstar, lnd_old, lnPgas_old, &
967 10 : 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 : 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 : 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 : 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 10633 : do k=1,nz
1173 10623 : from_k = comes_from(k)
1174 10623 : if (new_r(k) == old_r(from_k)) then
1175 9197 : xh(s%i_lnR,k) = xh_old(s%i_lnR,from_k)
1176 : else
1177 1426 : 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 9146 : xh(s%i_lnd,k) = xh_old(s%i_lnd,from_k)
1181 : else
1182 1477 : 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 1975 : 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 10623 : 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 : 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 : 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 : 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 : 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 : Rho, logRho, xq_outer, cell_dq, avg_energy, avg_PE, avg_KE, &
1347 : new_PE, new_KE, max_delta_energy, delta_energy, revised_energy, &
1348 10623 : sum_lnT, avg_lnT, new_lnT, sum_energy, new_xa(species), &
1349 : 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 : 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 : 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 : 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 : 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 : 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 5888 : 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 : real(dp) :: xq_inner, sum_dqs, old_xq_outer, old_xq_inner, &
1543 : 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 1472 : 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 : logT, res(num_eos_basic_results), &
1655 : d_dlnd(num_eos_basic_results), d_dlnT(num_eos_basic_results), &
1656 1472 : d_dxa(num_eos_d_dxa_results, s% species), &
1657 : 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 : 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 : 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 155 : c2(k) = 0
1792 155 : if (abs(sm1) <= abs(s00)) then
1793 81 : c1(k) = sm1
1794 : else
1795 74 : 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 : real(dp) :: a, b, c, old_xq_inner, old_xq_outer, xq_inner, &
1824 : xq_overlap_outer, xq_overlap_inner, dq1, sum_dqs, old_xq_mid, &
1825 : 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 17260 : avg_value = a
1887 : else
1888 18068 : old_xq_mid = (old_xq_outer + old_xq_inner)/2
1889 18068 : dq_outer = old_xq_mid - xq_overlap_outer
1890 18068 : dq_inner = old_xq_mid - xq_overlap_inner
1891 :
1892 18068 : if (order == 0) then
1893 0 : avg_value = a
1894 18068 : else if (order == 1 .or. c == 0) then
1895 : ! use slope to estimate average value in the region being used
1896 18068 : if (dbg) write(*,*) 'use slope to estimate average value'
1897 18068 : v_overlap_outer = a + dq_outer*b
1898 18068 : v_overlap_inner = a + dq_inner*b
1899 18068 : 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, old_omega, &
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, old_omega, &
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, &
1948 : old_j_rot, old_omega, xh, op_err)
1949 : if (op_err /= 0) ierr = op_err
1950 : end do
1951 : !$OMP END PARALLEL DO
1952 :
1953 0 : end subroutine adjust_omega
1954 :
1955 :
1956 0 : subroutine adjust1_omega(s, k, nz, nz_old, comes_from, &
1957 0 : xout_old, xout_new, old_dqbar, new_dqbar, &
1958 0 : old_j_rot, old_omega, xh, ierr)
1959 : use hydro_rotation, only: w_div_w_roche_jrot, update1_i_rot_from_xh
1960 : ! set new value for s% omega(k)
1961 : type (star_info), pointer :: s
1962 : integer, intent(in) :: k, nz, nz_old
1963 : integer, dimension(:) :: comes_from
1964 : real(dp), dimension(:), intent(in) :: &
1965 : xout_old, xout_new, old_dqbar, new_dqbar, old_j_rot, old_omega
1966 : real(dp), intent(in) :: xh(:,:)
1967 : integer, intent(out) :: ierr
1968 :
1969 : real(dp) :: xq_outer, xq_inner, j_tot, omega_tot, &
1970 : xq0, xq1, new_point_dqbar, dq_sum, dq, r00
1971 : integer :: kk, k_outer
1972 :
1973 : integer, parameter :: k_dbg = -1
1974 :
1975 : include 'formats'
1976 :
1977 0 : ierr = 0
1978 0 : xq_outer = xout_new(k)
1979 0 : new_point_dqbar = new_dqbar(k)
1980 0 : if (k < nz) then
1981 0 : xq_inner = xq_outer + new_point_dqbar
1982 : else
1983 0 : xq_inner = 1d0
1984 : end if
1985 :
1986 0 : if (k == k_dbg) then
1987 0 : write(*,2) 'xq_outer', k, xq_outer
1988 0 : write(*,2) 'xq_inner', k, xq_inner
1989 0 : write(*,2) 'new_point_dqbar', k, new_point_dqbar
1990 : end if
1991 :
1992 0 : dq_sum = 0d0
1993 0 : j_tot = 0d0
1994 0 : omega_tot = 0d0
1995 0 : if (xq_outer >= xout_old(nz_old)) then
1996 : ! new contained entirely in old center zone
1997 0 : k_outer = nz_old
1998 0 : if (k == k_dbg) &
1999 0 : write(*,2) 'new contained in old center', &
2000 0 : k_outer, xout_old(k_outer)
2001 0 : else if (k == 1) then
2002 0 : k_outer = 1
2003 : else
2004 0 : k_outer = comes_from(k-1)
2005 : end if
2006 :
2007 0 : do kk = k_outer, nz_old ! loop until reach m_inner
2008 :
2009 0 : if (kk == nz_old) then
2010 0 : xq1 = 1d0
2011 : else
2012 0 : xq1 = xout_old(kk+1)
2013 : end if
2014 0 : if (xq1 <= xq_outer) cycle
2015 :
2016 0 : xq0 = xout_old(kk)
2017 0 : if (xq0 >= xq_inner) then
2018 0 : if (dq_sum < new_point_dqbar .and. kk > 1) then
2019 : ! need to add a bit more from the previous source
2020 0 : dq = new_point_dqbar - dq_sum
2021 0 : dq_sum = new_point_dqbar
2022 0 : j_tot = j_tot + old_j_rot(kk-1)*dq
2023 0 : omega_tot = omega_tot + old_omega(kk-1)*dq
2024 : end if
2025 : exit
2026 : end if
2027 :
2028 : if (xq1 < xq_outer) then
2029 : ierr = -1
2030 0 : return
2031 : end if
2032 :
2033 0 : if (xq0 >= xq_outer .and. xq1 <= xq_inner) then ! entire old kk is in new k
2034 :
2035 0 : dq = old_dqbar(kk)
2036 0 : dq_sum = dq_sum + dq
2037 :
2038 0 : if (dq_sum > new_point_dqbar) then
2039 : ! dq too large -- numerical roundoff problems
2040 0 : dq = dq - (new_point_dqbar - dq_sum)
2041 0 : dq_sum = new_point_dqbar
2042 : end if
2043 :
2044 0 : j_tot = j_tot + old_j_rot(kk)*dq
2045 0 : omega_tot = omega_tot + old_omega(kk)*dq
2046 :
2047 0 : else if (xq0 <= xq_outer .and. xq1 >= xq_inner) then ! entire new k is in old kk
2048 :
2049 0 : dq = new_dqbar(k)
2050 0 : dq_sum = dq_sum + dq
2051 0 : j_tot = j_tot + old_j_rot(kk)*dq
2052 0 : omega_tot = omega_tot + old_omega(kk)*dq
2053 :
2054 : else ! only use the part of old kk that is in new k
2055 :
2056 0 : if (k == k_dbg) then
2057 0 : write(*,*) 'only use the part of old kk that is in new k', xq_inner <= xq1
2058 0 : write(*,1) 'xq_outer', xq_outer
2059 0 : write(*,1) 'xq_inner', xq_inner
2060 0 : write(*,1) 'xq0', xq0
2061 0 : write(*,1) 'xq1', xq1
2062 0 : write(*,1) 'dq_sum', dq_sum
2063 0 : write(*,1) 'new_point_dqbar', new_point_dqbar
2064 0 : write(*,1) 'new_point_dqbar - dq_sum', new_point_dqbar - dq_sum
2065 : end if
2066 :
2067 0 : if (xq_inner <= xq1) then ! this is the last part of new k
2068 :
2069 0 : if (k == k_dbg) write(*,3) 'this is the last part of new k', k, kk
2070 :
2071 0 : dq = new_point_dqbar - dq_sum
2072 0 : dq_sum = new_point_dqbar
2073 :
2074 : else ! we avoid this case if possible because of numerical roundoff
2075 :
2076 0 : if (k == k_dbg) write(*,3) 'we avoid this case if possible', k, kk
2077 :
2078 0 : dq = max(0d0, xq1 - xq_outer)
2079 0 : if (dq_sum + dq > new_point_dqbar) dq = new_point_dqbar - dq_sum
2080 0 : dq_sum = dq_sum + dq
2081 :
2082 : end if
2083 :
2084 0 : j_tot = j_tot + old_j_rot(kk)*dq
2085 0 : omega_tot = omega_tot + old_omega(kk)*dq
2086 :
2087 0 : if (dq <= 0) then
2088 0 : ierr = -1
2089 : return
2090 : end if
2091 :
2092 : end if
2093 :
2094 0 : if (dq_sum >= new_point_dqbar) then
2095 0 : if (k == k_dbg) then
2096 0 : write(*,2) 'exit for k', k
2097 0 : write(*,2) 'dq_sum', kk, dq_sum
2098 0 : write(*,2) 'new_point_dqbar', kk, new_point_dqbar
2099 : end if
2100 : exit
2101 : end if
2102 :
2103 : end do
2104 :
2105 0 : s% j_rot(k) = j_tot/dq_sum
2106 : ! set an omega seed before evaluating a rotation-dependent moment of inertia
2107 0 : s% omega(k) = omega_tot/dq_sum
2108 0 : r00 = get_r_from_xh(s,k)
2109 : s% w_div_w_crit_roche(k) = &
2110 : w_div_w_roche_jrot(r00,s% m(k),s% j_rot(k),s% cgrav(k), &
2111 0 : s% w_div_wcrit_max, s% w_div_wcrit_max2, s% w_div_wc_flag)
2112 0 : call update1_i_rot_from_xh(s, k)
2113 0 : s% omega(k) = s% j_rot(k)/s% i_rot(k)% val
2114 :
2115 0 : if (k_dbg == k) then
2116 0 : write(*,2) 's% omega(k)', k, s% omega(k)
2117 0 : write(*,2) 's% j_rot(k)', k, s% j_rot(k)
2118 0 : write(*,2) 's% i_rot(k)', k, s% i_rot(k)
2119 0 : if (s% model_number > 1925) call mesa_error(__FILE__,__LINE__,'debugging: adjust1_omega')
2120 : end if
2121 :
2122 : end subroutine adjust1_omega
2123 :
2124 :
2125 : ! like adjust_omega. conserve kinetic energy
2126 0 : subroutine do_v( &
2127 0 : s, nz, nz_old, cell_type, comes_from, &
2128 0 : old_xq, new_xq, old_dq, new_dq, xh, xh_old, &
2129 0 : xout_old, xout_new, old_dqbar, new_dqbar, old_ke, ierr)
2130 : use alloc
2131 : type (star_info), pointer :: s
2132 : integer, intent(in) :: nz, nz_old
2133 : integer, dimension(:) :: cell_type, comes_from
2134 : real(dp), dimension(:) :: &
2135 : xout_old, xout_new, old_dqbar, new_dqbar, &
2136 : old_xq, new_xq, old_dq, new_dq, old_ke
2137 : real(dp), dimension(:,:) :: xh, xh_old
2138 : integer, intent(out) :: ierr
2139 :
2140 : integer :: k, op_err, i_v
2141 : real(dp) :: old_ke_tot, new_ke_tot, xmstar, err
2142 :
2143 : include 'formats'
2144 0 : ierr = 0
2145 0 : i_v = s% i_v
2146 0 : xmstar = s% xmstar
2147 :
2148 0 : old_ke_tot = 0d0
2149 0 : do k=1,nz_old ! skip common factor 1/2 xmstar in ke
2150 0 : old_ke(k) = old_dqbar(k)*xh_old(i_v,k)*xh_old(i_v,k)
2151 0 : old_ke_tot = old_ke_tot + old_ke(k)
2152 : end do
2153 :
2154 0 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
2155 : do k = 1, nz
2156 : op_err = 0
2157 : call adjust1_v( &
2158 : s, k, nz, nz_old, cell_type, comes_from, xout_old, xout_new, &
2159 : old_dqbar, new_dqbar, old_ke, i_v, xh, xh_old, op_err)
2160 : if (op_err /= 0) ierr = op_err
2161 : end do
2162 : !$OMP END PARALLEL DO
2163 0 : if (ierr /= 0) then
2164 0 : return
2165 : end if
2166 :
2167 0 : new_ke_tot = 0
2168 0 : do k=1,nz
2169 0 : new_ke_tot = new_ke_tot + new_dqbar(k)*xh(i_v,k)*xh(i_v,k)
2170 : end do
2171 :
2172 0 : if (abs(old_ke_tot - new_ke_tot) > 1d-7*new_ke_tot) then
2173 0 : ierr = -1
2174 0 : s% retry_message = 'failed in mesh_adjust do_v'
2175 0 : if (s% report_ierr) write(*, *) s% retry_message
2176 : !call mesa_error(__FILE__,__LINE__,'do_v')
2177 0 : return
2178 : end if
2179 :
2180 0 : err = abs(old_ke_tot - new_ke_tot)/max(new_ke_tot,old_ke_tot,1d0)
2181 0 : s% mesh_adjust_KE_conservation = err
2182 :
2183 0 : if (s% trace_mesh_adjust_error_in_conservation) then
2184 0 : write(*,2) 'mesh adjust error in conservation of KE', s% model_number, &
2185 0 : err, new_ke_tot, old_ke_tot
2186 0 : if (err > 1d-10) then
2187 0 : write(*,*) 'err too large'
2188 0 : call mesa_error(__FILE__,__LINE__,'do_v')
2189 : end if
2190 : end if
2191 :
2192 : end subroutine do_v
2193 :
2194 :
2195 0 : subroutine adjust1_v( &
2196 0 : s, k, nz, nz_old, cell_type, comes_from, xout_old, xout_new, &
2197 0 : old_dqbar, new_dqbar, old_ke, i_v, xh, xh_old, ierr)
2198 : ! set new value for s% v(k) to conserve kinetic energy
2199 : type (star_info), pointer :: s
2200 : integer, intent(in) :: k, nz, nz_old, i_v
2201 : integer, dimension(:) :: cell_type, comes_from
2202 : real(dp), dimension(:), intent(in) :: &
2203 : xout_old, xout_new, old_dqbar, new_dqbar, old_ke
2204 : real(dp), dimension(:,:) :: xh, xh_old
2205 : integer, intent(out) :: ierr
2206 :
2207 : real(dp) :: xq_outer, xq_inner, ke_sum, &
2208 : xq0, xq1, new_point_dqbar, dq_sum, dq
2209 : integer :: kk, k_outer
2210 :
2211 : integer, parameter :: k_dbg = -1
2212 :
2213 : include 'formats'
2214 :
2215 0 : ierr = 0
2216 :
2217 0 : if (cell_type(k) == unchanged_type .or. &
2218 : cell_type(k) == revised_type) then
2219 0 : if (k == 1) then
2220 0 : xh(i_v,k) = xh_old(i_v,comes_from(k))
2221 0 : return
2222 : end if
2223 0 : if (cell_type(k-1) == unchanged_type) then
2224 0 : xh(i_v,k) = xh_old(i_v,comes_from(k))
2225 0 : return
2226 : end if
2227 : end if
2228 :
2229 0 : xq_outer = xout_new(k)
2230 0 : new_point_dqbar = new_dqbar(k)
2231 0 : if (k < nz) then
2232 0 : xq_inner = xq_outer + new_point_dqbar
2233 : else
2234 0 : xq_inner = 1d0
2235 : end if
2236 :
2237 0 : if (k == k_dbg) then
2238 0 : write(*,2) 'xq_outer', k, xq_outer
2239 0 : write(*,2) 'xq_inner', k, xq_inner
2240 0 : write(*,2) 'new_point_dqbar', k, new_point_dqbar
2241 : end if
2242 :
2243 0 : dq_sum = 0d0
2244 0 : ke_sum = 0
2245 0 : if (xq_outer >= xout_old(nz_old)) then
2246 : ! new contained entirely in old center zone
2247 0 : k_outer = nz_old
2248 0 : if (k == k_dbg) &
2249 0 : write(*,2) 'new contained in old center', &
2250 0 : k_outer, xout_old(k_outer)
2251 0 : else if (k == 1) then
2252 0 : k_outer = 1
2253 : else
2254 0 : k_outer = comes_from(k-1)
2255 : end if
2256 :
2257 0 : do kk = k_outer, nz_old ! loop until reach xq_inner
2258 :
2259 0 : if (kk == nz_old) then
2260 0 : xq1 = 1d0
2261 : else
2262 0 : xq1 = xout_old(kk+1)
2263 : end if
2264 0 : if (xq1 <= xq_outer) cycle
2265 :
2266 0 : xq0 = xout_old(kk)
2267 0 : if (xq0 >= xq_inner) then
2268 0 : if (dq_sum < new_point_dqbar .and. kk > 1) then
2269 : ! need to add a bit more from the previous source
2270 0 : dq = new_point_dqbar - dq_sum
2271 0 : dq_sum = new_point_dqbar
2272 0 : ke_sum = ke_sum + old_ke(kk-1)*dq/old_dqbar(kk-1)
2273 :
2274 0 : if (k == k_dbg) &
2275 0 : write(*,3) 'new k contains some of old kk-1', &
2276 0 : k, kk, old_ke(kk-1)*dq, dq_sum
2277 :
2278 : end if
2279 : exit
2280 : end if
2281 :
2282 : if (xq1 < xq_outer) then
2283 : ierr = -1
2284 : return
2285 : end if
2286 :
2287 0 : if (xq0 >= xq_outer .and. xq1 <= xq_inner) then ! entire old kk is in new k
2288 :
2289 0 : dq = old_dqbar(kk)
2290 0 : dq_sum = dq_sum + dq
2291 :
2292 0 : if (dq_sum > new_point_dqbar) then
2293 : ! dq too large -- numerical roundoff problems
2294 0 : dq = dq - (new_point_dqbar - dq_sum)
2295 0 : dq_sum = new_point_dqbar
2296 : end if
2297 :
2298 0 : ke_sum = ke_sum + old_ke(kk)*dq/old_dqbar(kk)
2299 :
2300 0 : if (k == k_dbg) &
2301 0 : write(*,3) 'new k contains all of old kk', &
2302 0 : k, kk, old_ke(kk)*dq, ke_sum
2303 :
2304 0 : else if (xq0 <= xq_outer .and. xq1 >= xq_inner) then ! entire new k is in old kk
2305 :
2306 0 : dq = new_dqbar(k)
2307 0 : dq_sum = dq_sum + dq
2308 0 : ke_sum = ke_sum + old_ke(kk)*dq/old_dqbar(kk)
2309 :
2310 0 : if (k == k_dbg) &
2311 0 : write(*,3) 'all new k is in old kk', &
2312 0 : k, kk, old_ke(kk)*dq, ke_sum
2313 :
2314 : else ! only use the part of old kk that is in new k
2315 :
2316 0 : if (k == k_dbg) then
2317 0 : write(*,*) 'only use the part of old kk that is in new k', xq_inner <= xq1
2318 0 : write(*,1) 'xq_outer', xq_outer
2319 0 : write(*,1) 'xq_inner', xq_inner
2320 0 : write(*,1) 'xq0', xq0
2321 0 : write(*,1) 'xq1', xq1
2322 0 : write(*,1) 'dq_sum', dq_sum
2323 0 : write(*,1) 'new_point_dqbar', new_point_dqbar
2324 0 : write(*,1) 'new_point_dqbar - dq_sum', new_point_dqbar - dq_sum
2325 : end if
2326 :
2327 0 : if (xq_inner <= xq1) then ! this is the last part of new k
2328 :
2329 0 : dq = new_point_dqbar - dq_sum
2330 0 : dq_sum = new_point_dqbar
2331 :
2332 : else ! we avoid this case if possible because of numerical roundoff
2333 :
2334 0 : if (k == k_dbg) write(*,3) 'we avoid this case if possible', k, kk
2335 :
2336 0 : dq = max(0d0, xq1 - xq_outer)
2337 0 : if (dq_sum + dq > new_point_dqbar) dq = new_point_dqbar - dq_sum
2338 0 : dq_sum = dq_sum + dq
2339 :
2340 : end if
2341 :
2342 0 : if (k == k_dbg) then
2343 0 : write(*,3) 'new k use only part of old kk', k, kk
2344 0 : write(*,2) 'dq_sum', k, dq_sum
2345 0 : write(*,2) 'dq', k, dq
2346 0 : write(*,2) 'old_ke(kk)', kk, old_ke(kk)
2347 0 : write(*,2) 'old ke_sum', k, ke_sum
2348 0 : write(*,2) 'new ke_sum', k, ke_sum + old_ke(kk)*dq
2349 : end if
2350 :
2351 0 : ke_sum = ke_sum + old_ke(kk)*dq/old_dqbar(kk)
2352 :
2353 0 : if (dq <= 0) then
2354 0 : ierr = -1
2355 : !return
2356 0 : write(*,*) 'dq <= 0', dq
2357 0 : call mesa_error(__FILE__,__LINE__,'debugging: adjust1_v')
2358 : end if
2359 :
2360 : end if
2361 :
2362 0 : if (dq_sum >= new_point_dqbar) then
2363 : exit
2364 : end if
2365 :
2366 : end do
2367 :
2368 0 : xh(i_v,k) = sqrt(ke_sum/new_point_dqbar) ! we have skipped the 1/2 xmstar factor
2369 0 : if (xh_old(i_v,comes_from(k)) < 0d0) xh(i_v,k) = -xh(i_v,k)
2370 :
2371 0 : if (k == k_dbg) then
2372 0 : !$OMP critical (adjust1_v_dbg)
2373 0 : write(*,2) 'xh(i_v,k) new_dqbar', k, xh(i_v,k), new_dqbar(k)
2374 0 : write(*,2) 'xh_old(i_v,comes_from(k)) old_dqbar', &
2375 0 : comes_from(k), xh_old(i_v,comes_from(k)), old_dqbar(comes_from(k))
2376 0 : if (k == k_dbg) call mesa_error(__FILE__,__LINE__,'adjust1_v')
2377 : !$OMP end critical (adjust1_v_dbg)
2378 : !stop
2379 : end if
2380 :
2381 : end subroutine adjust1_v
2382 :
2383 :
2384 0 : subroutine do_u( &
2385 0 : s, nz, nz_old, cell_type, comes_from, &
2386 0 : old_xq, new_xq, old_dq, new_dq, xh, xh_old, &
2387 0 : xout_old, xout_new, old_ke, ierr)
2388 : use alloc
2389 : type (star_info), pointer :: s
2390 : integer, intent(in) :: nz, nz_old
2391 : integer, dimension(:) :: cell_type, comes_from
2392 : real(dp), dimension(:) :: &
2393 : xout_old, xout_new, old_xq, new_xq, old_dq, new_dq, old_ke
2394 : real(dp), dimension(:,:) :: xh, xh_old
2395 : integer, intent(out) :: ierr
2396 :
2397 : integer :: k, op_err, i_u
2398 : real(dp) :: old_ke_tot, new_ke_tot, xmstar, err
2399 :
2400 : include 'formats'
2401 0 : ierr = 0
2402 0 : i_u = s% i_u
2403 0 : xmstar = s% xmstar
2404 :
2405 0 : old_ke_tot = 0d0
2406 0 : do k=1,nz_old ! skip common factor 1/2 xmstar in ke
2407 0 : old_ke(k) = old_dq(k)*xh_old(i_u,k)*xh_old(i_u,k)
2408 0 : old_ke_tot = old_ke_tot + old_ke(k)
2409 : end do
2410 :
2411 0 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
2412 : do k = 1, nz
2413 : op_err = 0
2414 : call adjust1_u( &
2415 : s, k, nz, nz_old, cell_type, comes_from, xout_old, xout_new, &
2416 : old_dq, new_dq, old_ke, i_u, xh, xh_old, op_err)
2417 : if (op_err /= 0) ierr = op_err
2418 : end do
2419 : !$OMP END PARALLEL DO
2420 0 : if (ierr /= 0) then
2421 0 : return
2422 : end if
2423 :
2424 0 : new_ke_tot = 0
2425 0 : do k=1,nz
2426 0 : new_ke_tot = new_ke_tot + new_dq(k)*xh(i_u,k)*xh(i_u,k)
2427 : end do
2428 :
2429 0 : err = abs(old_ke_tot - new_ke_tot)/max(new_ke_tot,old_ke_tot,1d0)
2430 0 : s% mesh_adjust_KE_conservation = err
2431 :
2432 0 : if (s% trace_mesh_adjust_error_in_conservation) then
2433 0 : write(*,2) 'mesh adjust error in conservation of KE', s% model_number, &
2434 0 : err, new_ke_tot, old_ke_tot
2435 0 : if (err > 1d-10) then
2436 0 : write(*,*) 'err too large'
2437 0 : call mesa_error(__FILE__,__LINE__,'do_u')
2438 : end if
2439 : end if
2440 :
2441 : end subroutine do_u
2442 :
2443 :
2444 0 : subroutine adjust1_u( &
2445 0 : s, k, nz, nz_old, cell_type, comes_from, xout_old, xout_new, &
2446 0 : old_dq, new_dq, old_ke, i_u, xh, xh_old, ierr)
2447 : ! set new value for s% u(k) to conserve kinetic energy
2448 : type (star_info), pointer :: s
2449 : integer, intent(in) :: k, nz, nz_old, i_u
2450 : integer, dimension(:) :: cell_type, comes_from
2451 : real(dp), dimension(:), intent(in) :: &
2452 : xout_old, xout_new, old_dq, new_dq, old_ke
2453 : real(dp), dimension(:,:) :: xh, xh_old
2454 : integer, intent(out) :: ierr
2455 :
2456 : real(dp) :: xq_outer, xq_inner, ke_sum, &
2457 : xq0, xq1, new_cell_dq, dq_sum, dq
2458 : integer :: kk, k_outer
2459 :
2460 : integer, parameter :: k_dbg = -1
2461 :
2462 : include 'formats'
2463 :
2464 0 : ierr = 0
2465 :
2466 0 : if (cell_type(k) == unchanged_type .or. &
2467 : cell_type(k) == revised_type) then
2468 0 : if (k == 1) then
2469 0 : xh(i_u,k) = xh_old(i_u,comes_from(k))
2470 0 : return
2471 : end if
2472 0 : if (cell_type(k-1) == unchanged_type) then
2473 0 : xh(i_u,k) = xh_old(i_u,comes_from(k))
2474 0 : return
2475 : end if
2476 : end if
2477 :
2478 0 : xq_outer = xout_new(k)
2479 0 : new_cell_dq = new_dq(k)
2480 0 : if (k < nz) then
2481 0 : xq_inner = xq_outer + new_cell_dq
2482 : else
2483 0 : xq_inner = 1d0
2484 : end if
2485 :
2486 0 : if (k == k_dbg) then
2487 0 : write(*,2) 'xq_outer', k, xq_outer
2488 0 : write(*,2) 'xq_inner', k, xq_inner
2489 0 : write(*,2) 'new_cell_dq', k, new_cell_dq
2490 : end if
2491 :
2492 0 : dq_sum = 0d0
2493 0 : ke_sum = 0
2494 0 : if (xq_outer >= xout_old(nz_old)) then
2495 : ! new contained entirely in old center zone
2496 0 : k_outer = nz_old
2497 0 : if (k == k_dbg) &
2498 0 : write(*,2) 'new contained in old center', &
2499 0 : k_outer, xout_old(k_outer)
2500 0 : else if (k == 1) then
2501 0 : k_outer = 1
2502 : else
2503 0 : k_outer = comes_from(k-1)
2504 : end if
2505 :
2506 0 : do kk = k_outer, nz_old ! loop until reach xq_inner
2507 :
2508 0 : if (kk == nz_old) then
2509 0 : xq1 = 1d0
2510 : else
2511 0 : xq1 = xout_old(kk+1)
2512 : end if
2513 0 : if (xq1 <= xq_outer) cycle
2514 :
2515 : if (xq1 < xq_outer) then
2516 : ierr = -1
2517 : return
2518 : end if
2519 :
2520 0 : xq0 = xout_old(kk)
2521 0 : if (xq0 >= xq_outer .and. xq1 <= xq_inner) then ! entire old kk is in new k
2522 :
2523 0 : dq = old_dq(kk)
2524 0 : dq_sum = dq_sum + dq
2525 :
2526 0 : if (dq_sum > new_cell_dq) then
2527 : ! dq too large -- numerical roundoff problems
2528 0 : dq = dq - (new_cell_dq - dq_sum)
2529 0 : dq_sum = new_cell_dq
2530 : end if
2531 :
2532 0 : ke_sum = ke_sum + old_ke(kk)*dq/old_dq(kk)
2533 :
2534 0 : if (k == k_dbg) &
2535 0 : write(*,3) 'new k contains all of old kk', &
2536 0 : k, kk, old_ke(kk)*dq, ke_sum
2537 :
2538 0 : else if (xq0 <= xq_outer .and. xq1 >= xq_inner) then ! entire new k is in old kk
2539 :
2540 0 : dq = new_dq(k)
2541 0 : dq_sum = dq_sum + dq
2542 0 : ke_sum = ke_sum + old_ke(kk)*dq/old_dq(kk)
2543 :
2544 0 : if (k == k_dbg) &
2545 0 : write(*,3) 'all new k is in old kk', &
2546 0 : k, kk, old_ke(kk)*dq, ke_sum
2547 :
2548 : else ! only use the part of old kk that is in new k
2549 :
2550 0 : if (k == k_dbg) then
2551 0 : write(*,*) 'only use the part of old kk that is in new k', xq_inner <= xq1
2552 0 : write(*,1) 'xq_outer', xq_outer
2553 0 : write(*,1) 'xq_inner', xq_inner
2554 0 : write(*,1) 'xq0', xq0
2555 0 : write(*,1) 'xq1', xq1
2556 0 : write(*,1) 'dq_sum', dq_sum
2557 0 : write(*,1) 'new_cell_dq', new_cell_dq
2558 0 : write(*,1) 'new_cell_dq - dq_sum', new_cell_dq - dq_sum
2559 : end if
2560 :
2561 0 : if (xq_inner <= xq1) then ! this is the last part of new k
2562 :
2563 0 : dq = new_cell_dq - dq_sum
2564 0 : dq_sum = new_cell_dq
2565 :
2566 : else ! we avoid this case if possible because of numerical roundoff
2567 :
2568 0 : if (k == k_dbg) write(*,3) 'we avoid this case if possible', k, kk
2569 :
2570 0 : dq = max(0d0, xq1 - xq_outer)
2571 0 : if (dq_sum + dq > new_cell_dq) dq = new_cell_dq - dq_sum
2572 0 : dq_sum = dq_sum + dq
2573 :
2574 : end if
2575 :
2576 0 : if (k == k_dbg) then
2577 0 : write(*,3) 'new k use only part of old kk', k, kk
2578 0 : write(*,2) 'dq_sum', k, dq_sum
2579 0 : write(*,2) 'dq', k, dq
2580 0 : write(*,2) 'old_ke(kk)', kk, old_ke(kk)
2581 0 : write(*,2) 'old ke_sum', k, ke_sum
2582 0 : write(*,2) 'new ke_sum', k, ke_sum + old_ke(kk)*dq
2583 : end if
2584 :
2585 0 : ke_sum = ke_sum + old_ke(kk)*dq/old_dq(kk)
2586 :
2587 0 : if (dq <= 0) then
2588 0 : ierr = -1
2589 : !return
2590 0 : write(*,*) 'dq <= 0', dq
2591 0 : call mesa_error(__FILE__,__LINE__,'debugging: adjust1_u')
2592 : end if
2593 :
2594 : end if
2595 :
2596 0 : if (dq_sum >= new_cell_dq) then
2597 : exit
2598 : end if
2599 :
2600 : end do
2601 :
2602 0 : xh(i_u,k) = sqrt(ke_sum/new_cell_dq) ! we have skipped the 1/2 xmstar factor
2603 0 : if (xh_old(i_u,comes_from(k)) < 0d0) xh(i_u,k) = -xh(i_u,k)
2604 :
2605 0 : if (k == k_dbg) then
2606 0 : !$OMP critical (adjust1_u_dbg)
2607 0 : write(*,2) 'xh(i_u,k) new_dq', k, xh(i_u,k), new_dq(k)
2608 0 : write(*,2) 'xh_old(i_u,comes_from(k)) old_dq', &
2609 0 : comes_from(k), xh_old(i_u,comes_from(k)), old_dq(comes_from(k))
2610 0 : if (k == k_dbg) call mesa_error(__FILE__,__LINE__,'adjust1_u')
2611 : !$OMP end critical (adjust1_u_dbg)
2612 : !stop
2613 : end if
2614 :
2615 : end subroutine adjust1_u
2616 :
2617 :
2618 0 : subroutine do_Hp_face( &
2619 : s, nz, nz_old, nzlo, nzhi, comes_from, xh, xh_old, &
2620 0 : xq, xq_old_plus1, xq_new, work, Hp_face_old_plus1, Hp_face_new, ierr)
2621 : use interp_1d_def
2622 : use interp_1d_lib
2623 : type (star_info), pointer :: s
2624 : integer, intent(in) :: nz, nz_old, nzlo, nzhi, comes_from(:)
2625 : real(dp), dimension(:,:), pointer :: xh, xh_old
2626 : real(dp), dimension(:), pointer :: work
2627 : real(dp), dimension(:) :: &
2628 : xq, xq_old_plus1, Hp_face_old_plus1, Hp_face_new, xq_new
2629 : integer, intent(out) :: ierr
2630 :
2631 : integer :: n, i_Hp, k
2632 :
2633 : include 'formats'
2634 :
2635 0 : ierr = 0
2636 0 : i_Hp = s% i_Hp
2637 0 : if (i_Hp == 0) return
2638 0 : n = nzhi - nzlo + 1
2639 :
2640 0 : do k=1,nz_old
2641 0 : Hp_face_old_plus1(k) = xh_old(i_Hp,k)
2642 : end do
2643 0 : Hp_face_old_plus1(nz_old+1) = Hp_face_old_plus1(nz_old)
2644 :
2645 : call interpolate_vector( &
2646 : nz_old+1, xq_old_plus1, n, xq_new, &
2647 : Hp_face_old_plus1, Hp_face_new, interp_pm, nwork, work, &
2648 0 : 'mesh_adjust do_Hp_face', ierr)
2649 0 : if (ierr /= 0) then
2650 : return
2651 : write(*,*) 'interpolate_vector failed in do_Hp_face for remesh'
2652 : call mesa_error(__FILE__,__LINE__,'debug: mesh adjust: do_Hp_face')
2653 : end if
2654 :
2655 0 : do k=nzlo,nzhi
2656 0 : xh(i_Hp,k) = Hp_face_new(k+1-nzlo)
2657 : end do
2658 :
2659 0 : n = nzlo - 1
2660 0 : if (n > 0) then
2661 0 : do k=1,n
2662 0 : xh(i_Hp,k) = xh_old(i_Hp,k)
2663 : end do
2664 : end if
2665 :
2666 0 : if (nzhi < nz) then
2667 0 : n = nz - nzhi - 1 ! nz-n = nzhi+1
2668 0 : do k=0,n
2669 0 : xh(i_Hp,nz-k) = xh_old(i_Hp,nz_old-k)
2670 : end do
2671 : end if
2672 :
2673 : end subroutine do_Hp_face
2674 :
2675 :
2676 0 : subroutine do_etrb( & ! same logic as do_u
2677 0 : s, nz, nz_old, cell_type, comes_from, &
2678 0 : old_xq, new_xq, old_dq, new_dq, xh, xh_old, &
2679 0 : xout_old, xout_new, old_eturb, ierr)
2680 : use alloc
2681 : type (star_info), pointer :: s
2682 : integer, intent(in) :: nz, nz_old
2683 : integer, dimension(:) :: cell_type, comes_from
2684 : real(dp), dimension(:) :: &
2685 : xout_old, xout_new, old_xq, new_xq, old_dq, new_dq, old_eturb
2686 : real(dp), dimension(:,:) :: xh, xh_old
2687 : integer, intent(out) :: ierr
2688 :
2689 : integer :: k, op_err, i_w
2690 : real(dp) :: old_eturb_tot, new_eturb_tot, xmstar, err
2691 :
2692 : include 'formats'
2693 0 : ierr = 0
2694 0 : i_w = s% i_w
2695 0 : xmstar = s% xmstar
2696 :
2697 0 : old_eturb_tot = 0d0
2698 0 : do k=1,nz_old
2699 0 : old_eturb(k) = old_dq(k)*pow2(xh_old(i_w,k))
2700 0 : old_eturb_tot = old_eturb_tot + old_eturb(k)
2701 : end do
2702 :
2703 0 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
2704 : do k = 1, nz
2705 : op_err = 0
2706 : call adjust1_etrb( &
2707 : s, k, nz, nz_old, cell_type, comes_from, xout_old, xout_new, &
2708 : old_dq, new_dq, old_eturb, i_w, xh, xh_old, op_err)
2709 : if (op_err /= 0) ierr = op_err
2710 : end do
2711 : !$OMP END PARALLEL DO
2712 0 : if (ierr /= 0) then
2713 0 : return
2714 : end if
2715 :
2716 0 : new_eturb_tot = 0
2717 0 : do k=1,nz
2718 0 : new_eturb_tot = new_eturb_tot + new_dq(k)*pow2(xh(i_w,k))
2719 : end do
2720 :
2721 0 : err = abs(old_eturb_tot - new_eturb_tot)/max(new_eturb_tot,old_eturb_tot,1d0)
2722 0 : s% mesh_adjust_Eturb_conservation = err
2723 :
2724 0 : if (s% trace_mesh_adjust_error_in_conservation) then
2725 0 : write(*,2) 'mesh adjust error in conservation of turbulent energy', s% model_number, &
2726 0 : err, new_eturb_tot, old_eturb_tot
2727 0 : if (err > 1d-10) then
2728 0 : write(*,*) 'err too large'
2729 0 : call mesa_error(__FILE__,__LINE__,'do_etrb')
2730 : end if
2731 : end if
2732 :
2733 : end subroutine do_etrb
2734 :
2735 :
2736 0 : subroutine adjust1_etrb( &
2737 0 : s, k, nz, nz_old, cell_type, comes_from, xout_old, xout_new, &
2738 0 : old_dq, new_dq, old_eturb, i_w, xh, xh_old, ierr)
2739 : ! set new value for s% w(k) to conserve turbulent energy
2740 : type (star_info), pointer :: s
2741 : integer, intent(in) :: k, nz, nz_old, i_w
2742 : integer, dimension(:) :: cell_type, comes_from
2743 : real(dp), dimension(:), intent(in) :: &
2744 : xout_old, xout_new, old_dq, new_dq, old_eturb
2745 : real(dp), dimension(:,:) :: xh, xh_old
2746 : integer, intent(out) :: ierr
2747 :
2748 : real(dp) :: xq_outer, xq_inner, eturb_sum, &
2749 : xq0, xq1, new_cell_dq, dq_sum, dq
2750 : integer :: kk, k_outer
2751 :
2752 : integer, parameter :: k_dbg = -1
2753 :
2754 : include 'formats'
2755 :
2756 0 : ierr = 0
2757 :
2758 0 : if (cell_type(k) == unchanged_type .or. &
2759 : cell_type(k) == revised_type) then
2760 : ! just copy the old value
2761 0 : if (k == 1) then
2762 0 : xh(i_w,k) = xh_old(i_w,comes_from(k))
2763 0 : return
2764 : end if
2765 0 : if (cell_type(k-1) == unchanged_type) then
2766 0 : xh(i_w,k) = xh_old(i_w,comes_from(k))
2767 0 : return
2768 : end if
2769 : end if
2770 :
2771 0 : xq_outer = xout_new(k)
2772 0 : new_cell_dq = new_dq(k)
2773 0 : if (k < nz) then
2774 0 : xq_inner = xq_outer + new_cell_dq
2775 : else
2776 0 : xq_inner = 1d0
2777 : end if
2778 :
2779 0 : if (k == k_dbg) then
2780 0 : write(*,2) 'xq_outer', k, xq_outer
2781 0 : write(*,2) 'xq_inner', k, xq_inner
2782 0 : write(*,2) 'new_cell_dq', k, new_cell_dq
2783 : end if
2784 :
2785 0 : dq_sum = 0d0
2786 0 : eturb_sum = 0
2787 0 : if (xq_outer >= xout_old(nz_old)) then
2788 : ! new contained entirely in old center zone
2789 0 : k_outer = nz_old
2790 0 : if (k == k_dbg) &
2791 0 : write(*,2) 'new contained in old center', &
2792 0 : k_outer, xout_old(k_outer)
2793 0 : else if (k == 1) then
2794 0 : k_outer = 1
2795 : else
2796 0 : k_outer = comes_from(k-1)
2797 : end if
2798 :
2799 0 : do kk = k_outer, nz_old ! loop until reach xq_inner
2800 :
2801 0 : if (kk == nz_old) then
2802 0 : xq1 = 1d0
2803 : else
2804 0 : xq1 = xout_old(kk+1)
2805 : end if
2806 0 : if (xq1 <= xq_outer) cycle
2807 :
2808 : if (xq1 < xq_outer) then
2809 : ierr = -1
2810 : return
2811 : end if
2812 :
2813 0 : xq0 = xout_old(kk)
2814 0 : if (xq0 >= xq_outer .and. xq1 <= xq_inner) then ! entire old kk is in new k
2815 :
2816 0 : dq = old_dq(kk)
2817 0 : dq_sum = dq_sum + dq
2818 :
2819 0 : if (dq_sum > new_cell_dq) then
2820 : ! dq too large -- numerical roundoff problems
2821 0 : dq = dq - (new_cell_dq - dq_sum)
2822 0 : dq_sum = new_cell_dq
2823 : end if
2824 :
2825 0 : eturb_sum = eturb_sum + old_eturb(kk)*dq/old_dq(kk)
2826 :
2827 0 : if (k == k_dbg) &
2828 0 : write(*,3) 'new k contains all of old kk', &
2829 0 : k, kk, old_eturb(kk)*dq, eturb_sum
2830 :
2831 0 : else if (xq0 <= xq_outer .and. xq1 >= xq_inner) then ! entire new k is in old kk
2832 :
2833 0 : dq = new_dq(k)
2834 0 : dq_sum = dq_sum + dq
2835 0 : eturb_sum = eturb_sum + old_eturb(kk)*dq/old_dq(kk)
2836 :
2837 0 : if (k == k_dbg) &
2838 0 : write(*,3) 'all new k is in old kk', &
2839 0 : k, kk, old_eturb(kk)*dq, eturb_sum
2840 :
2841 : else ! only use the part of old kk that is in new k
2842 :
2843 0 : if (k == k_dbg) then
2844 0 : write(*,*) 'only use the part of old kk that is in new k', xq_inner <= xq1
2845 0 : write(*,1) 'xq_outer', xq_outer
2846 0 : write(*,1) 'xq_inner', xq_inner
2847 0 : write(*,1) 'xq0', xq0
2848 0 : write(*,1) 'xq1', xq1
2849 0 : write(*,1) 'dq_sum', dq_sum
2850 0 : write(*,1) 'new_cell_dq', new_cell_dq
2851 0 : write(*,1) 'new_cell_dq - dq_sum', new_cell_dq - dq_sum
2852 : end if
2853 :
2854 0 : if (xq_inner <= xq1) then ! this is the last part of new k
2855 :
2856 0 : dq = new_cell_dq - dq_sum
2857 0 : dq_sum = new_cell_dq
2858 :
2859 : else ! we avoid this case if possible because of numerical roundoff
2860 :
2861 0 : if (k == k_dbg) write(*,3) 'we avoid this case if possible', k, kk
2862 :
2863 0 : dq = max(0d0, xq1 - xq_outer)
2864 0 : if (dq_sum + dq > new_cell_dq) dq = new_cell_dq - dq_sum
2865 0 : dq_sum = dq_sum + dq
2866 :
2867 : end if
2868 :
2869 0 : if (k == k_dbg) then
2870 0 : write(*,3) 'new k use only part of old kk', k, kk
2871 0 : write(*,2) 'dq_sum', k, dq_sum
2872 0 : write(*,2) 'dq', k, dq
2873 0 : write(*,2) 'old_eturb(kk)', kk, old_eturb(kk)
2874 0 : write(*,2) 'old eturb_sum', k, eturb_sum
2875 0 : write(*,2) 'new eturb_sum', k, eturb_sum + old_eturb(kk)*dq
2876 : end if
2877 :
2878 0 : eturb_sum = eturb_sum + old_eturb(kk)*dq/old_dq(kk)
2879 :
2880 0 : if (dq <= 0) then
2881 0 : ierr = -1
2882 : !return
2883 0 : write(*,*) 'dq <= 0', dq
2884 0 : call mesa_error(__FILE__,__LINE__,'debugging: adjust1_etrb')
2885 : end if
2886 :
2887 : end if
2888 :
2889 0 : if (dq_sum >= new_cell_dq) then
2890 : exit
2891 : end if
2892 :
2893 : end do
2894 :
2895 0 : xh(i_w,k) = sqrt(max(0d0,eturb_sum/new_cell_dq))
2896 :
2897 : end subroutine adjust1_etrb
2898 :
2899 : end module mesh_adjust
|
