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