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 adjust_mass
21 :
22 : use star_private_def
23 : use const_def, only: dp, i8, ln10, msun, msun, secyer, one_third, four_thirds_pi
24 : use utils_lib
25 :
26 : implicit none
27 :
28 : private
29 : public :: do_adjust_mass
30 : public :: compute_prev_mesh_dm
31 :
32 : logical, parameter :: dbg_adjm = .false.
33 :
34 : contains
35 :
36 11 : real(dp) function compute_delta_m(s) result(delta_m)
37 : type (star_info), pointer :: s
38 : include 'formats'
39 11 : delta_m = s% dt*s% mstar_dot
40 :
41 11 : if (s% super_eddington_wind_mdot /= 0 .and. s% super_eddington_wind_mdot > -s% mstar_dot) then
42 : ! NOTE: super_eddington should not be greater than max_wind
43 0 : if (s% max_wind > 0 .and. &
44 : s% max_wind*Msun/secyer < s% super_eddington_wind_mdot) then
45 0 : s% mstar_dot = -s% max_wind*Msun/secyer
46 : else
47 0 : s% mstar_dot = -s% super_eddington_wind_mdot
48 : end if
49 : else if (delta_m == 0 &
50 : .or. (delta_m < 0 .and. s% star_mass <= s% min_star_mass_for_loss) &
51 11 : .or. (delta_m > 0 .and. s% max_star_mass_for_gain > 0 &
52 : .and. s% star_mass >= s% max_star_mass_for_gain)) then
53 11 : s% mstar_dot = 0d0
54 : end if
55 :
56 11 : delta_m = s% dt*s% mstar_dot
57 :
58 11 : if (is_bad(s% dt)) then
59 0 : write(*,1) 's% dt', s% dt
60 0 : call mesa_error(__FILE__,__LINE__,'compute_delta_m')
61 : end if
62 :
63 11 : if (is_bad(s% mstar_dot)) then
64 0 : write(*,1) 's% mstar_dot', s% mstar_dot
65 0 : call mesa_error(__FILE__,__LINE__,'compute_delta_m')
66 : end if
67 :
68 11 : if (is_bad(delta_m)) then
69 0 : write(*,1) 'delta_m', delta_m
70 0 : call mesa_error(__FILE__,__LINE__,'compute_delta_m')
71 : end if
72 :
73 11 : end function compute_delta_m
74 :
75 11 : subroutine save_for_eps_mdot(s)
76 : use star_utils, only: eval_deltaM_total_energy_integrals
77 : type (star_info), pointer :: s
78 :
79 : integer :: j
80 :
81 : call eval_deltaM_total_energy_integrals( &
82 : s, 1, s% nz, s% mstar, .true., &
83 : s% total_energy_profile_before_adjust_mass, &
84 : s% total_internal_energy_before_adjust_mass, &
85 : s% total_gravitational_energy_before_adjust_mass, &
86 : s% total_radial_kinetic_energy_before_adjust_mass, &
87 : s% total_rotational_kinetic_energy_before_adjust_mass, &
88 : s% total_turbulent_energy_before_adjust_mass, &
89 11 : s% total_energy_before_adjust_mass)
90 :
91 13098 : do j=1,s%nz
92 13098 : s% dm_before_adjust_mass(j) = s% dm(j)
93 : end do
94 :
95 11 : end subroutine save_for_eps_mdot
96 :
97 :
98 : ! Reconstruct previous mesh based on what the mass differences should have been.
99 : ! This is to deal with the changes in cell mass being small (relatively speaking) at depth.
100 : ! That results in differences at the limit of double precision, so we need to go to quad precision.
101 : ! However dm and prev_mesh_dm are only known to double precision, whereas the difference between them
102 : ! is known to quad precision. So we take dm as given and reconstruct prev_mesh_dm using the difference.
103 0 : subroutine compute_prev_mesh_dm(s, prev_mesh_dm, dm, change_in_dm)
104 : ! Inputs
105 : type (star_info), pointer :: s
106 :
107 : ! Intermediates
108 : integer :: j, nz
109 0 : real(dp) :: change_sum
110 :
111 : ! Outputs
112 : real(qp), dimension(:), intent(out) :: prev_mesh_dm, dm, change_in_dm
113 :
114 :
115 0 : nz = s% nz
116 :
117 0 : change_sum = 0
118 0 : do j=1,nz
119 0 : dm(j) = s%dm(j)
120 0 : prev_mesh_dm(j) = s%dm_before_adjust_mass(j)
121 0 : if (j < s%k_below_const_q) then
122 0 : change_in_dm(j) = s%adjust_mass_outer_frac_sub1 * prev_mesh_dm(j)
123 0 : else if (j < s%k_const_mass) then
124 0 : change_in_dm(j) = s%adjust_mass_mid_frac_sub1 * prev_mesh_dm(j)
125 : else
126 0 : change_in_dm(j) = s%adjust_mass_inner_frac_sub1 * prev_mesh_dm(j)
127 : end if
128 0 : change_sum = change_sum + change_in_dm(j)
129 0 : prev_mesh_dm(j) = dm(j) - change_in_dm(j)
130 : end do
131 :
132 11 : end subroutine compute_prev_mesh_dm
133 :
134 : ! Updates the radii of cells after adjust mass. Note that the radius appears in many
135 : ! places in s. Currently those are: r, r_start, rmid, rmid_start, R2 (r^2), and lnR (log(r)).
136 : ! If more of these are added they should also be updated here.
137 0 : subroutine update_radius(s)
138 : use star_utils, only: store_r_in_xh, get_r_and_lnR_from_xh
139 : ! Inputs
140 : type (star_info), pointer :: s
141 :
142 : ! Intermediates
143 : integer :: j, nz
144 0 : real(dp) :: r_new, vol00, volp1, cell_vol
145 :
146 0 : nz = s%nz
147 :
148 0 : vol00 = four_thirds_pi*s% R_center*s% R_center*s% R_center
149 0 : do j=nz,1,-1
150 :
151 0 : volp1 = vol00
152 0 : cell_vol = s% dm(j)/s% rho(j)
153 0 : vol00 = volp1 + cell_vol
154 :
155 0 : r_new = exp(log(vol00/four_thirds_pi)*one_third)
156 :
157 0 : s%r(j) = r_new
158 0 : s%r_start(j) = s%r(j)
159 0 : s%rmid(j) = r_new
160 0 : s%rmid_start(j) = r_new
161 0 : call store_r_in_xh(s, j, r_new)
162 0 : call get_r_and_lnR_from_xh(s, j, s% r(j), s% lnR(j))
163 0 : s% xh_start(s% i_lnR,j) = s% xh(s% i_lnR,j)
164 :
165 : end do
166 :
167 0 : end subroutine update_radius
168 :
169 11 : subroutine do_adjust_mass(s, species, ierr)
170 0 : use adjust_xyz, only: get_xa_for_accretion
171 : use star_utils, only: report_xa_bad_nums, &
172 : start_time, update_time
173 : use hydro_rotation, only: set_omega
174 : use chem_def
175 :
176 : type (star_info), pointer :: s
177 : integer, intent(in) :: species
178 : integer, intent(out) :: ierr
179 :
180 : real(dp) :: &
181 22 : dt, delta_m, old_mstar, new_mstar, total, &
182 22 : frac, env_mass, mmax, new_xmstar, old_xmstar, removed, &
183 33 : q_for_just_added, xq_for_CpT_absMdot_div_L, sum_dq, dm, sumx
184 :
185 198 : real(dp), dimension(species) :: xaccrete, mtot_init
186 : real(dp), dimension(:), allocatable :: &
187 11 : rxm_old, rxm_new, old_cell_mass, new_cell_mass, &
188 11 : oldloc, newloc, oldval, newval, xm_old, xm_new, &
189 11 : xq_center_old, xq_center_new
190 11 : real(dp), dimension(:,:), allocatable :: xa_old
191 11 : real(dp), pointer :: work(:)
192 :
193 : integer :: j, k, l, m, k_const_mass, nz, k_below_just_added, k_newval
194 : integer(i8) :: time0
195 11 : real(dp) :: partial_xa_mass, region_total_mass
196 11 : real(dp) :: starting_j_rot_surf
197 :
198 : logical, parameter :: dbg = .false.
199 :
200 : include 'formats'
201 :
202 : if (dbg) write(*,*) 'do_adjust_mass'
203 :
204 11 : call save_for_eps_mdot(s)
205 :
206 11 : ierr = 0
207 :
208 11 : nz = s% nz
209 11 : dt = s% dt
210 :
211 :
212 11 : s% k_const_mass = 1
213 :
214 11 : s% k_for_test_CpT_absMdot_div_L = 1
215 :
216 11 : q_for_just_added = 1d0
217 11 : s% k_below_just_added = 1
218 :
219 : ! NOTE: don't assume that vars are all set at this point.
220 : ! exceptions: omega, i_rot, j_rot, and total_angular_momentum have been set.
221 : ! store surface j_rot, to later adjust angular momentum lost
222 11 : if (s% rotation_flag) then
223 0 : starting_j_rot_surf = s% j_rot(1)
224 : end if
225 :
226 11 : delta_m = compute_delta_m(s)
227 :
228 11 : if (delta_m == 0) then
229 : return
230 : end if
231 :
232 0 : if (s% doing_timing) call start_time(s, time0, total)
233 :
234 0 : s% need_to_setvars = .true.
235 :
236 0 : old_mstar = s% mstar
237 0 : old_xmstar = s% xmstar
238 :
239 0 : new_mstar = old_mstar + delta_m
240 0 : new_xmstar = old_xmstar + delta_m
241 :
242 0 : if (is_bad(new_xmstar)) then
243 0 : write(*,1) 'new_xmstar', new_xmstar
244 0 : call mesa_error(__FILE__,__LINE__,'do_adjust_mass')
245 : end if
246 :
247 : if (delta_m > 0 .and. s% max_star_mass_for_gain > 0 &
248 0 : .and. new_mstar > Msun*s% max_star_mass_for_gain) then
249 0 : new_mstar = Msun*s% max_star_mass_for_gain
250 0 : delta_m = new_mstar - old_mstar
251 0 : s% mstar_dot = delta_m/dt
252 0 : else if (delta_m < 0 .and. new_mstar < Msun*s% min_star_mass_for_loss) then
253 0 : new_mstar = Msun*s% min_star_mass_for_loss
254 0 : delta_m = new_mstar - old_mstar
255 0 : s% mstar_dot = delta_m/dt
256 : end if
257 :
258 0 : frac = old_xmstar/new_xmstar
259 0 : new_xmstar = old_xmstar/frac
260 0 : if (new_xmstar <= 0) then
261 0 : ierr = -1
262 0 : return
263 : end if
264 0 : s% xmstar = new_xmstar
265 0 : s% mstar = s% xmstar + s% M_center
266 :
267 : if (dbg_adjm) then
268 : env_mass = old_mstar - s% he_core_mass*Msun
269 : write(*,'(a40,f26.16)') 'env_mass/old_mstar', env_mass/old_mstar
270 : write(*,'(A)')
271 : write(*,1) 'delta_m/old_mstar', delta_m/old_mstar
272 : write(*,1) 's% he_core_mass*Msun', s% he_core_mass*Msun
273 : write(*,1) 'env_mass', env_mass
274 : write(*,1) 'delta_m/env_mass', delta_m/env_mass
275 : write(*,1) 'log10(abs(delta_m/env_mass))', safe_log10(abs(delta_m/env_mass))
276 : write(*,'(A)')
277 : end if
278 :
279 0 : call do_alloc(ierr)
280 0 : if (ierr /= 0) return
281 :
282 0 : do k=1,nz
283 0 : old_cell_mass(k) = old_xmstar*s% dq(k)
284 : end do
285 0 : xm_old(1) = 0
286 0 : xq_center_old(1) = 0.5d0*s% dq(1)
287 0 : do k=2,nz
288 0 : xm_old(k) = xm_old(k-1) + old_cell_mass(k-1)
289 0 : xq_center_old(k) = xq_center_old(k-1) + 0.5d0*(s% dq(k-1) + s% dq(k))
290 : end do
291 :
292 0 : s% xa_removed(1:species) = 0
293 0 : s% angular_momentum_removed = 0
294 0 : s% angular_momentum_added = 0
295 0 : if (delta_m < 0) then
296 0 : s% angular_momentum_removed = angular_momentum_removed(ierr)
297 0 : if (ierr /= 0) then
298 0 : call dealloc
299 0 : return
300 : end if
301 0 : removed = -delta_m
302 0 : do k=1,nz
303 0 : if (xm_old(k) >= removed) then
304 0 : do j=1,species
305 0 : s% xa_removed(j) = s% xa_removed(j)/removed
306 : end do
307 : exit
308 : end if
309 0 : do j=1,species
310 0 : dm = min(old_cell_mass(k), removed-xm_old(k))
311 0 : s% xa_removed(j) = s% xa_removed(j) + s% xa(j,k)*dm
312 : end do
313 : end do
314 : end if
315 : !s% angular_momentum_added is set after j_rot is adjusted in outer layers
316 :
317 : call revise_q_and_dq( &
318 0 : s, nz, old_xmstar, new_xmstar, delta_m, k_const_mass, ierr)
319 0 : if (ierr /= 0) then
320 0 : s% retry_message = 'revise_q_and_dq failed in adjust mass'
321 0 : if (s% report_ierr) write(*,*) s% retry_message
322 0 : call dealloc
323 0 : return
324 : end if
325 :
326 0 : mtot_init(1:species) = 0
327 0 : do k=1,k_const_mass ! save total mass by species down to where constant
328 0 : sumx = 0
329 0 : do j=1,species
330 0 : s% xa(j,k) = max(0d0, min(1d0, s% xa(j,k)))
331 0 : sumx = sumx + s% xa(j,k)
332 : end do
333 0 : do j=1,species
334 0 : s% xa(j,k) = s% xa(j,k)/sumx
335 0 : mtot_init(j) = mtot_init(j) + old_cell_mass(k)*s% xa(j,k)
336 : end do
337 : end do
338 :
339 0 : do k=1,nz
340 0 : new_cell_mass(k) = new_xmstar*s% dq(k)
341 : end do
342 :
343 0 : xm_new(1) = 0
344 0 : xq_center_new(1) = 0.5d0*s% dq(1)
345 0 : do k=2,nz
346 0 : xm_new(k) = xm_new(k-1) + new_cell_mass(k-1)
347 0 : xq_center_new(k) = xq_center_new(k-1) + 0.5d0*(s% dq(k-1) + s% dq(k))
348 : end do
349 :
350 0 : xq_for_CpT_absMdot_div_L = abs(delta_m)/new_xmstar
351 0 : s% k_for_test_CpT_absMdot_div_L = nz
352 0 : sum_dq = 0
353 0 : do k = 1, nz
354 0 : if (sum_dq >= xq_for_CpT_absMdot_div_L) then
355 0 : s% k_for_test_CpT_absMdot_div_L = k; exit
356 : end if
357 0 : sum_dq = sum_dq + s% dq(k)
358 : end do
359 :
360 0 : if (delta_m > 0) then
361 0 : q_for_just_added = old_xmstar/new_xmstar
362 0 : s% k_below_just_added = nz
363 0 : do k = 1, nz
364 0 : if (s% q(k) <= q_for_just_added) then
365 0 : s% k_below_just_added = k; exit
366 : end if
367 : end do
368 : end if
369 :
370 0 : k_below_just_added = s% k_below_just_added
371 0 : k_newval = 1
372 :
373 0 : do k=1,nz
374 0 : do j=1,species
375 0 : xa_old(j,k) = s% xa(j,k)
376 : end do
377 : end do
378 :
379 0 : if (delta_m < 0) then
380 0 : xaccrete(1:species) = 0 ! xaccrete not used when removing mass
381 : else ! set xaccrete for composition of added material
382 0 : if (s% accrete_same_as_surface) then
383 0 : do j=1,species
384 0 : xaccrete(j) = xa_old(j,1)
385 : end do
386 : else
387 0 : call get_xa_for_accretion(s, xaccrete, ierr)
388 0 : if (ierr /= 0) then
389 0 : s% retry_message = 'get_xa_for_accretion failed in adjust mass'
390 0 : if (s% report_ierr) write(*, *) s% retry_message
391 0 : call dealloc
392 0 : return
393 : end if
394 : end if
395 : end if
396 :
397 0 : mmax = max(old_mstar, new_mstar)
398 :
399 : ! rxm_old and rxm_new are for interpolating by mass
400 : ! but instead of using mass coord, we want to use external mass
401 : ! in order to get better accuracy near the surface.
402 : ! to simplify this, the zero point is the same for both rxm_old and rxm_new.
403 : ! that makes rxm_new = rxm_old for k >= k_const_mass.
404 0 : if (delta_m < 0) then
405 0 : rxm_old(1) = 0
406 0 : rxm_new(1) = -delta_m ! note that rxm_new(1) > 0 since delta_m < 0
407 : else
408 0 : rxm_old(1) = delta_m
409 0 : rxm_new(1) = 0
410 : end if
411 0 : do k = 2, nz
412 0 : rxm_old(k) = rxm_old(k-1) + old_cell_mass(k-1)
413 0 : if (k >= k_const_mass) then
414 0 : rxm_new(k) = rxm_old(k)
415 0 : new_cell_mass(k) = old_cell_mass(k)
416 : else
417 0 : rxm_new(k) = rxm_new(k-1) + new_cell_mass(k-1)
418 : end if
419 : end do
420 :
421 : call set_xa(s, nz, k_const_mass, species, xa_old, xaccrete, &
422 0 : rxm_old, rxm_new, mmax, old_cell_mass, new_cell_mass, ierr)
423 0 : if (ierr /= 0) then
424 0 : s% retry_message = 'set_xa failed in adjust mass'
425 0 : if (s% report_ierr) write(*,*) s% retry_message
426 0 : call dealloc
427 0 : return
428 : end if
429 :
430 0 : if (s% rotation_flag) then
431 : ! update j_rot if have extra angular momentum change
432 : call set_omega_adjust_mass( &
433 : s, nz, k_const_mass, k_below_just_added, delta_m, &
434 : rxm_old, rxm_new, mmax, old_cell_mass, new_cell_mass, &
435 : ! reuse these work arrays
436 : oldloc, newloc, oldval, newval, xm_old, xm_new, &
437 0 : ierr)
438 0 : if (ierr /= 0) then
439 0 : s% retry_message = 'set_omega_adjust_mass failed in adjust mass'
440 0 : if (s% report_ierr) write(*,*) s% retry_message
441 0 : call dealloc
442 0 : return
443 : end if
444 0 : if (s% do_adjust_J_lost) then
445 0 : call adjust_J_lost(s, k_below_just_added, starting_j_rot_surf, ierr)
446 0 : if (ierr /= 0) then
447 0 : s% retry_message = 'do_adjust_J_lost failed in adjust mass'
448 0 : if (s% report_ierr) write(*,*) s% retry_message
449 0 : call dealloc
450 0 : return
451 : end if
452 : end if
453 0 : if (s% D_omega_flag) then
454 : call set_D_omega( &
455 : s, nz, k_const_mass, k_newval, &
456 : rxm_old, rxm_new, delta_m, old_xmstar, new_xmstar, &
457 0 : s% D_omega, oldloc, newloc, oldval, newval, work, ierr)
458 0 : if (ierr /= 0) then
459 0 : s% retry_message = 'set_D_omega failed in adjust mass'
460 0 : if (s% report_ierr) write(*,*) s% retry_message
461 0 : call dealloc
462 0 : return
463 : end if
464 : end if
465 0 : if (s% am_nu_rot_flag) then
466 : call set_D_omega( & ! reuse the set_D_omega routine to set am_nu_rot
467 : s, nz, k_const_mass, k_newval, &
468 : rxm_old, rxm_new, delta_m, old_xmstar, new_xmstar, &
469 0 : s% am_nu_rot, oldloc, newloc, oldval, newval, work, ierr)
470 0 : if (ierr /= 0) then
471 0 : s% retry_message = 'set_am_nu_rot failed in adjust mass'
472 0 : if (s% report_ierr) write(*,*) s% retry_message
473 0 : call dealloc
474 0 : return
475 : end if
476 : end if
477 : end if
478 :
479 : ! soften outer xa (Pablo -- this should be the very last step)
480 0 : if (s% smooth_outer_xa_big > 0.0d0 .and. s% smooth_outer_xa_small > 0.0d0) then
481 0 : m = 1
482 0 : do k = 1, s% nz
483 0 : if (s% m(1) - s% m(k) > s% smooth_outer_xa_big * s% m(1)) exit
484 0 : region_total_mass = 0
485 0 : m = k
486 0 : do l = k, s% nz
487 0 : if (s% m(k) - s% m(l) > s% smooth_outer_xa_small * s% m(1) * &
488 : (1-(s% m(1) - s% m(k))/(s% smooth_outer_xa_big * s% m(1)))) exit
489 0 : m = l
490 0 : region_total_mass = region_total_mass + s% dm(l)
491 : end do
492 0 : if (m == k) cycle
493 0 : do j=1,s% species
494 : partial_xa_mass = 0
495 0 : do l = k, m
496 0 : partial_xa_mass = partial_xa_mass + s% dm(l) * s% xa(j,l)
497 : end do
498 0 : do l = k, m
499 0 : s% xa(j,l) = partial_xa_mass / region_total_mass
500 : end do
501 : end do
502 : end do
503 : ! fix small errors to ensure xa's add up to unity
504 0 : do k=1,m
505 0 : frac = 1d0/sum(s% xa(1:s% species,k))
506 0 : do j=1,s% species
507 0 : s% xa(j,k) = s% xa(j,k)*frac
508 : end do
509 : end do
510 : end if
511 :
512 :
513 0 : call update_radius(s)
514 :
515 :
516 0 : call dealloc
517 :
518 :
519 0 : if (s% doing_timing) call update_time(s, time0, total, s% time_adjust_mass)
520 :
521 : if (dbg_adjm) call mesa_error(__FILE__,__LINE__,'debugging: do_adjust_mass')
522 11 : if (dbg) write(*,*) 'do_adjust_mass return'
523 :
524 : contains
525 :
526 :
527 0 : real(dp) function angular_momentum_removed(ierr) result(J)
528 : ! when call this, s% j_rot is still for old mass
529 : integer, intent(out) :: ierr
530 : integer :: k
531 0 : real(dp) :: dmm1, dm00, dm, dm_sum, dm_lost
532 : include 'formats'
533 0 : ierr = 0
534 0 : J = 0
535 0 : if (.not. s% rotation_flag) return
536 0 : dm00 = 0
537 0 : dm_sum = 0
538 0 : dm_lost = -delta_m
539 0 : do k = 1, nz
540 0 : dmm1 = dm00
541 0 : dm00 = old_cell_mass(k)
542 0 : dm = 0.5d0*(dmm1+dm00)
543 0 : if (dm_sum + dm > dm_lost) then
544 0 : dm = dm_lost - dm_sum
545 0 : dm_sum = dm_lost
546 : else
547 : dm_sum = dm_sum + dm
548 : end if
549 0 : J = J + dm*s% j_rot(k)
550 0 : if (dm_sum == dm_lost) exit
551 : end do
552 11 : end function angular_momentum_removed
553 :
554 :
555 : real(dp) function eval_total_angular_momentum(s,cell_mass,nz_last) result(J)
556 : type (star_info), pointer :: s
557 : real(dp) :: cell_mass(:)
558 : integer, intent(in) :: nz_last
559 : integer :: k
560 : real(dp) :: dmm1, dm00, dm
561 : include 'formats'
562 : J = 0
563 : if (.not. s% rotation_flag) return
564 : dm00 = 0
565 : do k = 1, nz_last
566 : dmm1 = dm00
567 : dm00 = cell_mass(k)
568 : if (k == s% nz) then
569 : dm = 0.5d0*dmm1+dm00
570 : else if (k == nz_last) then
571 : dm = 0.5d0*dmm1
572 : else
573 : dm = 0.5d0*(dmm1+dm00)
574 : end if
575 : J = J + dm*s% j_rot(k)
576 : end do
577 : end function eval_total_angular_momentum
578 :
579 0 : subroutine do_alloc(ierr)
580 : integer, intent(out) :: ierr
581 0 : allocate(rxm_old(nz), rxm_new(nz), old_cell_mass(nz), new_cell_mass(nz), &
582 0 : xa_old(species,nz), oldloc(nz), newloc(nz), oldval(nz), newval(nz), &
583 0 : xm_old(nz), xm_new(nz), xq_center_old(nz), xq_center_new(nz))
584 0 : call do_work_arrays(.true.,ierr)
585 0 : end subroutine do_alloc
586 :
587 0 : subroutine dealloc
588 : integer :: ierr
589 0 : call do_work_arrays(.false.,ierr)
590 : end subroutine dealloc
591 :
592 0 : subroutine do_work_arrays(alloc_flag, ierr)
593 : use interp_1d_def
594 : use alloc, only: work_array
595 : logical, intent(in) :: alloc_flag
596 : integer, intent(out) :: ierr
597 : logical, parameter :: crit = .false.
598 : ierr = 0
599 : call work_array(s, alloc_flag, crit, &
600 0 : work, nz*pm_work_size, nz_alloc_extra, 'adjust_mass work', ierr)
601 0 : if (ierr /= 0) return
602 0 : end subroutine do_work_arrays
603 :
604 : end subroutine do_adjust_mass
605 :
606 :
607 0 : subroutine revise_q_and_dq( &
608 : s, nz, old_xmstar, new_xmstar, delta_m, k_const_mass, ierr)
609 : use star_utils, only: get_lnT_from_xh
610 : type (star_info), pointer :: s
611 : integer, intent(in) :: nz
612 : real(dp), intent(in) :: old_xmstar, new_xmstar, delta_m
613 : integer, intent(out) :: k_const_mass, ierr
614 :
615 : integer :: k, kA, kB, j00, jp1
616 0 : real(dp) :: lnTlim_A, lnTlim_B, sumdq, &
617 0 : mold_o_mnew, lnTmax, lnT_A, lnT_B, &
618 0 : xqA, xqB_old, xqB_new, qfrac, frac, dqacc
619 0 : real(dp) :: xq(nz)
620 0 : real(qp) :: qfrac_qp, frac_qp, mold_o_mnew_qp, q1, q2
621 0 : real(qp) :: adjust_mass_outer_frac, adjust_mass_mid_frac, adjust_mass_inner_frac
622 : integer, parameter :: min_kA = 5
623 : logical :: dbg, flag
624 : logical :: okay_to_move_kB_inward
625 :
626 : include 'formats'
627 :
628 0 : ierr = 0
629 0 : dbg = .false.
630 0 : flag = .false.
631 :
632 0 : if (is_bad(old_xmstar)) then
633 0 : write(*,1) 'old_xmstar', old_xmstar
634 0 : call mesa_error(__FILE__,__LINE__,'revise_q_and_dq')
635 : end if
636 :
637 0 : if (is_bad(new_xmstar)) then
638 0 : write(*,1) 'new_xmstar', new_xmstar
639 0 : call mesa_error(__FILE__,__LINE__,'revise_q_and_dq')
640 : end if
641 :
642 0 : if (is_bad(delta_m)) then
643 0 : write(*,1) 'delta_m', delta_m
644 0 : call mesa_error(__FILE__,__LINE__,'revise_q_and_dq')
645 : end if
646 :
647 :
648 0 : okay_to_move_kB_inward = .false.
649 :
650 0 : lnTlim_A = ln10*s% max_logT_for_k_below_const_q
651 0 : lnTlim_B = ln10*s% max_logT_for_k_const_mass
652 :
653 0 : q1 = old_xmstar
654 0 : q2 = new_xmstar
655 0 : mold_o_mnew_qp = q1/q2
656 0 : mold_o_mnew = mold_o_mnew_qp
657 0 : s% max_q_for_k_below_const_q = mold_o_mnew
658 0 : dqacc = delta_m / new_xmstar
659 :
660 : ! might have drifted away from summing to 1 in mesh adjustment, fix up
661 0 : sumdq = 0
662 0 : do k=1,nz
663 0 : sumdq = sumdq + s% dq(k)
664 : end do
665 0 : frac = 1.0d0/sumdq
666 0 : if (is_bad(frac)) then
667 0 : write(*,1) 'frac for initial renorm', frac
668 0 : call mesa_error(__FILE__,__LINE__,'revise_q_and_dq')
669 : end if
670 0 : do k = 1, nz
671 0 : if (is_bad(s% dq(k))) then
672 0 : write(*,2) 'bad dq input', s% dq(k)
673 0 : call mesa_error(__FILE__,__LINE__,'revise_q_and_dq')
674 : end if
675 0 : s% dq(k) = s% dq(k) * frac
676 : end do
677 :
678 : ! will work in xq := 1-q to reduce truncation errors
679 0 : xq(1)=0
680 0 : do k = 2, nz
681 0 : xq(k) = xq(k-1) + s% dq(k-1)
682 : end do
683 :
684 0 : k = maxloc(s% xh(s% i_lnT,1:nz),dim=1)
685 0 : lnTmax = get_lnT_from_xh(s, k)
686 0 : lnT_A = min(lnTmax, lnTlim_A)
687 :
688 0 : if (is_bad(s% max_q_for_k_below_const_q)) then
689 0 : write(*,*) 's% max_q_for_k_below_const_q', s% max_q_for_k_below_const_q
690 0 : call mesa_error(__FILE__,__LINE__,'revise_q_and_dq')
691 : end if
692 0 : if (is_bad(s% min_q_for_k_below_const_q)) then
693 0 : write(*,*) 's% min_q_for_k_below_const_q', s% min_q_for_k_below_const_q
694 0 : call mesa_error(__FILE__,__LINE__,'revise_q_and_dq')
695 : end if
696 :
697 0 : kA = min_kA
698 0 : do k = min_kA, nz-1
699 0 : kA = k
700 0 : if ( (1-xq(k)) > s% max_q_for_k_below_const_q) cycle
701 0 : if ( 1.0d0-xq(k) <= s% min_q_for_k_below_const_q) exit
702 0 : if (get_lnT_from_xh(s, k) >= lnT_A) exit
703 : end do
704 0 : xqA = xq(kA)
705 :
706 0 : lnT_B = min(lnTmax, lnTlim_B)
707 :
708 0 : if (is_bad(s% max_q_for_k_const_mass)) then
709 0 : write(*,*) 's% max_q_for_k_const_mass', s% max_q_for_k_const_mass
710 0 : call mesa_error(__FILE__,__LINE__,'revise_q_and_dq')
711 : end if
712 0 : if (is_bad(s% min_q_for_k_const_mass)) then
713 0 : write(*,*) 's% min_q_for_k_const_mass', s% min_q_for_k_const_mass
714 0 : call mesa_error(__FILE__,__LINE__,'revise_q_and_dq')
715 : end if
716 :
717 0 : kB = kA+1
718 0 : do k = kB, nz
719 0 : kB = k
720 0 : if ( (1-xq(k)) > s% max_q_for_k_const_mass) cycle
721 0 : if ( (1-xq(k)) <= s% min_q_for_k_const_mass) exit
722 0 : if (get_lnT_from_xh(s, k) >= lnT_B) exit
723 : end do
724 :
725 0 : xqB_old = xq(kB)
726 :
727 : if (dbg_adjm) then
728 : write(*,*) 'kA', kA
729 : write(*,1) 'xqA', xqA
730 : write(*,*) 'kB', kB
731 : write(*,1) 'xqB_old', xqB_old
732 : write(*,'(A)')
733 : write(*,1) 'xqA-xqB_old', xqA-xqB_old
734 : write(*,'(A)')
735 : end if
736 :
737 0 : xqB_new = dqacc + xqB_old*mold_o_mnew ! in order to keep m interior to kB constant
738 : ! same as qB_new = qB_old*mold/mnew, but without the truncation problems
739 :
740 0 : do ! make sure qfrac is not too far from 1 by moving kB inward
741 0 : q1 = xqB_new - xqA
742 0 : q2 = max(1d-99,xqB_old - xqA)
743 0 : qfrac_qp = q1/q2
744 0 : qfrac = qfrac_qp
745 0 : if (kB == nz) exit
746 0 : if (kB-kA > 10) then
747 0 : if (qfrac > 0.67d0 .and. qfrac < 1.5d0) exit
748 0 : if (qfrac > 0.50d0 .and. qfrac < 2.0d0) then
749 0 : j00 = maxloc(s% xa(:,kB),dim=1) ! most abundant species at kB
750 0 : jp1 = maxloc(s% xa(:,kB+1),dim=1) ! most abundant species at kB+1
751 0 : if (j00 /= jp1) then ! change in composition.
752 : if (dbg) write(*,*) 'change in composition. back up kB.'
753 0 : kB = max(1,kB-5)
754 0 : exit
755 : end if
756 : end if
757 : end if
758 0 : kB = kB+1
759 0 : xqB_old = xq(kB)
760 0 : xqB_new = dqacc + xqB_old*mold_o_mnew
761 : end do
762 :
763 : ! set new dq's
764 : ! s% dq(1:kA-1) unchanged
765 0 : s% dq(kA:kB-1) = s% dq(kA:kB-1)*qfrac
766 0 : frac_qp = mold_o_mnew_qp
767 0 : frac = frac_qp
768 0 : if (is_bad(frac)) then
769 0 : write(*,1) 'frac for kA:kB-1', frac
770 0 : call mesa_error(__FILE__,__LINE__,'revise_q_and_dq')
771 : end if
772 0 : s% dq(kB:nz) = s% dq(kB:nz)*frac
773 :
774 0 : adjust_mass_outer_frac = 1d0
775 0 : adjust_mass_mid_frac = qfrac_qp
776 0 : adjust_mass_inner_frac = frac_qp
777 :
778 : if (dbg_adjm) then
779 : write(*,1) 'frac for kb:nz', frac
780 : write(*,1) 'qfrac for kA:kB-1', qfrac
781 : write(*,1) 'revise_q_and_dq sum dqs', sum(s% dq(1:nz))
782 : write(*,2) 'qfrac region', kB, qfrac, s% q(kB), s% lnT(kB)/ln10
783 : write(*,2) 'frac region', kA, frac, s% q(kA), s% lnT(kA)/ln10
784 : write(*,2) 'kA', kA
785 : write(*,2) 'kB', kB
786 : write(*,2) 'nz', nz
787 : write(*,'(A)')
788 : call mesa_error(__FILE__,__LINE__,'adjust_mass')
789 : end if
790 :
791 : ! renorm
792 0 : sumdq = 0
793 0 : do k=1,nz
794 0 : sumdq = sumdq + s% dq(k)
795 : end do
796 0 : q1 = 1d0
797 0 : q2 = sumdq
798 0 : frac_qp = q1/q2
799 0 : frac = frac_qp
800 0 : if (is_bad(frac)) then
801 0 : write(*,1) 'frac for renorm', frac
802 0 : call mesa_error(__FILE__,__LINE__,'revise_q_and_dq')
803 : end if
804 0 : do k = 1, nz
805 0 : s% dq(k) = s% dq(k) * frac
806 : end do
807 :
808 0 : s% adjust_mass_outer_frac_sub1 = frac_qp*adjust_mass_outer_frac*new_xmstar / old_xmstar - 1.0_qp
809 0 : s% adjust_mass_mid_frac_sub1 = frac_qp*adjust_mass_mid_frac*new_xmstar / old_xmstar - 1.0_qp
810 0 : s% adjust_mass_inner_frac_sub1 = frac_qp*adjust_mass_inner_frac*new_xmstar / old_xmstar - 1.0_qp
811 :
812 : ! set q's to match new dq's
813 0 : s% q(1) = 1d0
814 0 : sumdq = s% dq(1)
815 0 : do k = 2, nz-1
816 0 : s% q(k) = 1d0 - sumdq
817 0 : sumdq = sumdq + s% dq(k)
818 : end do
819 0 : s% q(nz) = 1d0 - sumdq
820 :
821 0 : s% dq(nz) = s% q(nz)
822 :
823 : ! update dm's for new dq's
824 0 : do k=1,nz
825 0 : s% dm(k) = s% dq(k) * new_xmstar
826 : end do
827 :
828 0 : s% k_below_const_q = kA
829 0 : s% k_const_mass = kB
830 0 : k_const_mass = kB
831 :
832 0 : end subroutine revise_q_and_dq
833 :
834 :
835 0 : subroutine set_xa( &
836 0 : s, nz, k_const_mass, species, xa_old, xaccrete, &
837 0 : old_cell_xbdy, new_cell_xbdy, mmax, old_cell_mass, new_cell_mass, ierr)
838 : ! set new values for s% xa(:,:)
839 : type (star_info), pointer :: s
840 : integer, intent(in) :: nz, k_const_mass, species
841 : real(dp), intent(in) :: mmax
842 : real(dp), intent(in) :: xa_old(:, :), xaccrete(:)
843 : real(dp), dimension(:), intent(in) :: &
844 : old_cell_xbdy, new_cell_xbdy, old_cell_mass, new_cell_mass ! (nz)
845 : integer, intent(out) :: ierr
846 : integer :: k, j, op_err
847 : real(dp), parameter :: max_sum_abs = 10d0
848 : real(dp), parameter :: xsum_tol = 1d-2
849 : include 'formats'
850 0 : ierr = 0
851 : if (dbg_adjm) &
852 : write(*,2) 'set_xa: k_const_mass', k_const_mass
853 0 : if (k_const_mass < nz) then
854 : ! for k >= k_const_mass have m_new(k) = m_old(k),
855 : ! so no change in xa_new(:,k) for k > k_const_mass
856 0 : do k=k_const_mass+1,nz
857 0 : do j=1,species
858 0 : s% xa(j,k) = xa_old(j,k)
859 : end do
860 : end do
861 : end if
862 0 : !$OMP PARALLEL DO PRIVATE(k, op_err) SCHEDULE(dynamic,2)
863 : do k = 1, k_const_mass
864 : op_err = 0
865 : call set1_xa(s, k, nz, species, xa_old, xaccrete, &
866 : old_cell_xbdy, new_cell_xbdy, mmax, old_cell_mass, new_cell_mass, op_err)
867 : if (op_err /= 0) then
868 : if (s% report_ierr) write(*,2) 'set1_xa error', k
869 : ierr = op_err
870 : end if
871 : end do
872 : !$OMP END PARALLEL DO
873 :
874 0 : end subroutine set_xa
875 :
876 :
877 0 : subroutine set1_xa(s, k, nz, species, xa_old, xaccrete, &
878 0 : old_cell_xbdy, new_cell_xbdy, mmax, old_cell_mass, new_cell_mass, ierr)
879 : ! set new values for s% xa(:,k)
880 : use num_lib, only: binary_search
881 : type (star_info), pointer :: s
882 : integer, intent(in) :: k, nz, species
883 : real(dp), intent(in) :: mmax
884 : real(dp), intent(in) :: xa_old(:,:), xaccrete(:)
885 : real(dp), dimension(:), intent(in) :: &
886 : old_cell_xbdy, new_cell_xbdy, old_cell_mass, new_cell_mass
887 : integer, intent(out) :: ierr
888 :
889 0 : real(dp) :: msum(species), xm_inner, sumx, &
890 0 : xm0, xm1, new_cell_dm, dm_sum, dm, xm_outer
891 : integer :: kk, k_outer, j
892 :
893 : integer, parameter :: k_dbg = -1
894 : logical, parameter :: xa_dbg = .false.
895 :
896 : logical, parameter :: do_not_mix_accretion = .false.
897 :
898 : logical :: partially_accreted_cell
899 :
900 : include 'formats'
901 :
902 0 : ierr = 0
903 0 : msum(:) = -1 ! for testing
904 :
905 0 : xm_outer = new_cell_xbdy(k)
906 0 : if (k == nz) then
907 0 : new_cell_dm = mmax - xm_outer - s% M_center
908 : else
909 0 : new_cell_dm = new_cell_mass(k)
910 : end if
911 0 : xm_inner = xm_outer + new_cell_dm
912 :
913 0 : dm_sum = 0d0
914 :
915 0 : partially_accreted_cell = .false.
916 0 : if (xm_outer < old_cell_xbdy(1)) then ! there is some accreted material in new cell
917 0 : if (do_not_mix_accretion .or. xm_inner <= old_cell_xbdy(1)) then
918 : ! new cell is entirely accreted material
919 : !write(*,2) 'new cell is entirely accreted material', k, new_cell_dm
920 0 : do j=1,species
921 0 : s% xa(j,k) = xaccrete(j)
922 : end do
923 : return
924 : end if
925 0 : dm = min(new_cell_dm, old_cell_xbdy(1) - xm_outer)
926 0 : dm_sum = dm
927 0 : do j=1,species
928 0 : msum(j) = xaccrete(j)*dm
929 : end do
930 0 : partially_accreted_cell = .true.
931 0 : xm_outer = old_cell_xbdy(1)
932 0 : k_outer = 1
933 : else ! new cell entirely composed of old material
934 0 : msum(:) = 0
935 0 : if (xm_outer >= old_cell_xbdy(nz)) then
936 : ! new cell contained entirely in old cell nz
937 : k_outer = nz
938 : else
939 : ! binary search for k_outer such that
940 : ! xm_outer >= old_cell_xbdy(k_outer)
941 : ! and old_cell_xbdy(k_outer+1) > xm_outer
942 0 : k_outer = binary_search(nz, old_cell_xbdy, 0, xm_outer)
943 :
944 : ! check
945 0 : if (k_outer <= 0 .or. k_outer > nz) then
946 :
947 0 : ierr = -1
948 0 : if (.not. xa_dbg) return
949 :
950 : write(*,2) 'k', k
951 : write(*,2) 'k_outer', k_outer
952 : write(*,1) 'xm_outer', xm_outer
953 : write(*,2) 'old_cell_xbdy(1)', 1, old_cell_xbdy(1)
954 : write(*,2) 'old_cell_xbdy(nz)', nz, old_cell_xbdy(nz)
955 : call mesa_error(__FILE__,__LINE__,'debugging: set1_xa')
956 : end if
957 :
958 0 : if (xm_outer < old_cell_xbdy(k_outer)) then
959 :
960 0 : ierr = -1
961 0 : if (.not. xa_dbg) return
962 :
963 : write(*,*) 'k', k
964 : write(*,*) 'k_outer', k_outer
965 : write(*,1) 'xm_outer', xm_outer
966 : write(*,1) 'old_cell_xbdy(k_outer)', old_cell_xbdy(k_outer)
967 : write(*,*) '(xm_outer < old_cell_xbdy(k_outer))'
968 : call mesa_error(__FILE__,__LINE__,'debugging: set1_xa')
969 : end if
970 :
971 0 : if (k_outer < nz) then
972 0 : if (old_cell_xbdy(k_outer+1) <= xm_outer) then
973 :
974 0 : ierr = -1
975 0 : if (.not. xa_dbg) return
976 :
977 : write(*,*) 'k', k
978 : write(*,*) 'k_outer', k_outer
979 : write(*,1) 'xm_outer', xm_outer
980 : write(*,1) 'old_cell_xbdy(k_outer+1)', old_cell_xbdy(k_outer+1)
981 : write(*,*) '(old_cell_xbdy(k_outer+1) <= xm_outer)'
982 : call mesa_error(__FILE__,__LINE__,'debugging: set1_xa')
983 : end if
984 : end if
985 :
986 : end if
987 : end if
988 :
989 0 : if (k == -1) then
990 0 : ierr = -1
991 0 : if (.not. xa_dbg) return
992 :
993 : write(*,2) 'nz', nz
994 : write(*,2) 'k_outer', k_outer
995 : write(*,1) 'xm_outer', xm_outer
996 : write(*,1) 'xm_inner', xm_inner
997 : end if
998 :
999 0 : do kk = k_outer, nz ! loop until reach m_inner
1000 0 : xm0 = old_cell_xbdy(kk)
1001 :
1002 0 : if (xm0 >= xm_inner) then
1003 0 : if (dm_sum < new_cell_dm .and. kk > 1) then
1004 : ! need to add a bit more from the previous source cell
1005 0 : dm = new_cell_dm - dm_sum
1006 0 : dm_sum = new_cell_dm
1007 0 : do j=1,species
1008 0 : msum(j) = msum(j) + xa_old(j,kk-1)*dm
1009 : end do
1010 : end if
1011 : exit
1012 : end if
1013 :
1014 0 : if (kk == nz) then
1015 0 : xm1 = mmax - s% M_center
1016 : else
1017 0 : xm1 = old_cell_xbdy(kk+1)
1018 : end if
1019 :
1020 0 : if (xm1 < xm_outer) then
1021 0 : ierr = -1
1022 0 : if (.not. xa_dbg) return
1023 : write(*,'(A)')
1024 : write(*,*) 'k', k
1025 : write(*,*) 'kk', kk
1026 : write(*,1) 'xm1', xm1
1027 : write(*,1) 'xm_outer', xm_outer
1028 : write(*,*) 'xm1 < xm_outer'
1029 : call mesa_error(__FILE__,__LINE__,'debugging: set1_xa')
1030 : end if
1031 :
1032 0 : if (xm0 >= xm_outer .and. xm1 <= xm_inner) then
1033 : ! entire old cell kk is in new cell k
1034 :
1035 0 : dm = old_cell_mass(kk)
1036 0 : dm_sum = dm_sum + dm
1037 :
1038 0 : if (dm_sum > new_cell_dm) then
1039 : ! dm too large -- numerical roundoff problems
1040 0 : dm = dm - (new_cell_dm - dm_sum)
1041 0 : dm_sum = new_cell_dm
1042 : end if
1043 :
1044 0 : do j=1,species
1045 0 : msum(j) = msum(j) + xa_old(j,kk)*dm
1046 : end do
1047 :
1048 0 : else if (xm0 <= xm_outer .and. xm1 >= xm_inner) then
1049 : ! entire new cell k is in old cell kk
1050 : ! except if it is a partially accreted cell for which xm_outer has been reset
1051 0 : if (partially_accreted_cell) then
1052 0 : dm = xm_inner-xm_outer
1053 : else
1054 0 : dm = new_cell_mass(k)
1055 : end if
1056 0 : dm_sum = dm_sum + dm
1057 0 : do j=1,species
1058 0 : msum(j) = msum(j) + xa_old(j,kk)*dm
1059 : end do
1060 :
1061 : else ! only use the part of old cell kk that is in new cell k
1062 :
1063 0 : if (xm_inner <= xm1) then ! this is the last part of new cell k
1064 :
1065 0 : dm = new_cell_dm - dm_sum
1066 0 : dm_sum = new_cell_dm
1067 :
1068 : else ! notice that we avoid this case if possible because of numerical roundoff
1069 :
1070 0 : dm = max(0d0, xm1 - xm_outer)
1071 0 : if (dm_sum + dm > new_cell_dm) dm = new_cell_dm - dm_sum
1072 0 : dm_sum = dm_sum + dm
1073 :
1074 : end if
1075 :
1076 0 : do j=1,species
1077 0 : msum(j) = msum(j) + xa_old(j,kk)*dm
1078 : end do
1079 :
1080 0 : if (dm <= 0) then
1081 0 : ierr = -1
1082 0 : if (.not. xa_dbg) return
1083 : write(*,*) 'dm <= 0', dm
1084 : call mesa_error(__FILE__,__LINE__,'debugging: set1_xa')
1085 : end if
1086 :
1087 : end if
1088 :
1089 0 : if (dm_sum >= new_cell_dm) then
1090 : exit
1091 : end if
1092 :
1093 : end do
1094 :
1095 : ! revise and renormalize
1096 0 : do j=1,species
1097 0 : s% xa(j,k) = msum(j) / new_cell_mass(k)
1098 : end do
1099 0 : sumx = sum(s% xa(:,k))
1100 :
1101 0 : sumx = sum(s% xa(:,k))
1102 0 : if (sumx <= 0d0) then
1103 0 : ierr = -1
1104 0 : return
1105 : end if
1106 :
1107 0 : do j=1,species
1108 0 : s% xa(j,k) = s% xa(j,k)/sumx
1109 : end do
1110 :
1111 0 : end subroutine set1_xa
1112 :
1113 :
1114 0 : subroutine set_omega_adjust_mass( &
1115 : s, nz, k_const_mass, k_below_just_added, delta_m, &
1116 0 : old_cell_xbdy, new_cell_xbdy, mmax, old_cell_mass, new_cell_mass, &
1117 : ! work arrays (nz)
1118 0 : old_xout, new_xout, old_dmbar, new_dmbar, old_j_rot, extra_work, &
1119 : ierr)
1120 0 : use star_utils, only: total_angular_momentum
1121 : type (star_info), pointer :: s
1122 : integer, intent(in) :: nz, k_const_mass, k_below_just_added
1123 : real(dp), intent(in) :: mmax, delta_m
1124 : real(dp), dimension(:), intent(in) :: &
1125 : old_cell_xbdy, new_cell_xbdy, old_cell_mass, new_cell_mass ! (nz)
1126 : real(dp), dimension(:) :: &
1127 : old_xout, new_xout, old_dmbar, new_dmbar, old_j_rot, extra_work
1128 : integer, intent(out) :: ierr
1129 :
1130 : integer :: k, k0, op_err
1131 : logical :: okay
1132 0 : real(dp) :: old_j_tot, goal_total_added, actual_total_added, &
1133 0 : goal_total, bdy_j, bdy_total, inner_total, outer_total
1134 :
1135 : include 'formats'
1136 :
1137 0 : ierr = 0
1138 :
1139 0 : old_xout(1) = old_cell_xbdy(1)
1140 0 : new_xout(1) = new_cell_xbdy(1)
1141 0 : old_dmbar(1) = old_cell_mass(1)/2
1142 0 : new_dmbar(1) = new_cell_mass(1)/2
1143 0 : old_j_rot(1) = s% j_rot(1)
1144 0 : do k=2,nz
1145 0 : old_xout(k) = old_xout(k-1) + old_dmbar(k-1)
1146 0 : new_xout(k) = new_xout(k-1) + new_dmbar(k-1)
1147 0 : old_dmbar(k) = (old_cell_mass(k-1) + old_cell_mass(k))/2
1148 0 : new_dmbar(k) = (new_cell_mass(k-1) + new_cell_mass(k))/2
1149 0 : old_j_rot(k) = s% j_rot(k)
1150 : end do
1151 0 : old_dmbar(nz) = old_cell_mass(nz-1)/2 + old_cell_mass(nz)
1152 0 : new_dmbar(nz) = new_cell_mass(nz-1)/2 + new_cell_mass(nz)
1153 :
1154 0 : old_j_tot = dot_product(s% j_rot(1:nz),old_dmbar(1:nz))
1155 :
1156 0 : if (k_below_just_added == 1) then
1157 : k0 = 1
1158 : else
1159 0 : k0 = k_below_just_added + 1
1160 : end if
1161 0 : !$OMP PARALLEL DO PRIVATE(k, op_err) SCHEDULE(dynamic,2)
1162 : do k = 1, k_const_mass
1163 : if (k < k0) then
1164 : cycle
1165 : end if
1166 : op_err = 0
1167 : call set1_omega(s, k, k_below_just_added, nz, &
1168 : old_xout, new_xout, mmax, old_dmbar, new_dmbar, old_j_rot, op_err)
1169 : if (op_err /= 0) ierr = op_err
1170 : end do
1171 : !$OMP END PARALLEL DO
1172 :
1173 0 : if (k_below_just_added > 1) then
1174 : ! set omega in cells with newly added material
1175 0 : if (s% use_accreted_material_j) then
1176 0 : actual_total_added = 0d0
1177 0 : do k=1,k_below_just_added-2 ! remaining 2 done below
1178 0 : s% j_rot(k) = s% accreted_material_j
1179 0 : call set1_irot(s, k, k_below_just_added, .true.)
1180 0 : s% omega(k) = s% j_rot(k)/s% i_rot(k)% val
1181 0 : actual_total_added = actual_total_added + s% j_rot(k)*new_dmbar(k)
1182 : end do
1183 : k = k_below_just_added
1184 0 : goal_total_added = delta_m*s% accreted_material_j
1185 0 : goal_total = old_j_tot + goal_total_added
1186 0 : inner_total = dot_product(s% j_rot(k+1:nz),new_dmbar(k+1:nz))
1187 0 : outer_total = dot_product(s% j_rot(1:k-2),new_dmbar(1:k-2))
1188 0 : bdy_total = goal_total - (inner_total + outer_total)
1189 0 : bdy_j = bdy_total/sum(new_dmbar(k-1:k))
1190 0 : if (bdy_j <= 0) then
1191 0 : ierr = -1
1192 : return
1193 : end if
1194 0 : do k=k_below_just_added-1,k_below_just_added
1195 0 : s% j_rot(k) = bdy_j
1196 0 : call set1_irot(s, k, k_below_just_added, .true.)
1197 0 : s% omega(k) = s% j_rot(k)/s% i_rot(k)% val
1198 : end do
1199 : else ! use old surface omega in all the new material
1200 0 : do k=1,k_below_just_added-1
1201 0 : s% omega(k) = s% omega(k_below_just_added)
1202 0 : call set1_irot(s, k, k_below_just_added, .false.)
1203 0 : s% j_rot(k) = s% omega(k)*s% i_rot(k)% val
1204 : end do
1205 : end if
1206 : end if
1207 :
1208 0 : okay = .true.
1209 0 : do k=1,nz
1210 0 : if (is_bad(s% omega(k)) .or. &
1211 0 : abs(s% omega(k)) > 1d50) then
1212 0 : if (s% report_ierr) then
1213 0 : !$OMP critical (star_adjust_mass_omega)
1214 0 : write(*,2) 's% omega(k)', k, s% omega(k)
1215 : !$OMP end critical (star_adjust_mass_omega)
1216 : end if
1217 0 : if (s% stop_for_bad_nums) then
1218 0 : write(*,2) 's% omega(k)', k, s% omega(k)
1219 0 : call mesa_error(__FILE__,__LINE__,'set_omega_adjust_mass')
1220 : end if
1221 : okay = .false.
1222 : end if
1223 : end do
1224 0 : if (.not. okay) then
1225 0 : write(*,2) 'model_number', s% model_number
1226 0 : call mesa_error(__FILE__,__LINE__,'set_omega_adjust_mass')
1227 : end if
1228 :
1229 0 : end subroutine set_omega_adjust_mass
1230 :
1231 :
1232 0 : subroutine set1_irot(s, k, k_below_just_added, jrot_known) ! using lnR_for_d_dt_const_m
1233 0 : use hydro_rotation, only: eval_i_rot, w_div_w_roche_jrot, w_div_w_roche_omega
1234 : use star_utils, only: get_r_from_xh
1235 : type (star_info), pointer :: s
1236 : integer, intent(in) :: k, k_below_just_added
1237 : logical, intent(in) :: jrot_known
1238 :
1239 : real(dp) :: r00
1240 0 : real(dp) :: w_div_wcrit_roche
1241 :
1242 0 : r00 = get_r_from_xh(s,k)
1243 : ! The moment of inertia depends
1244 : ! on the ratio of rotational frequency to its critical value.
1245 : ! This ratio is computed in two different ways depending on whether
1246 : ! omega or j_rot is known.
1247 0 : if (jrot_known) then
1248 : w_div_wcrit_roche = w_div_w_roche_jrot(r00,s% m(k),s% j_rot(k),s% cgrav(k), &
1249 0 : s% w_div_wcrit_max, s% w_div_wcrit_max2, s% w_div_wc_flag)
1250 : else
1251 : w_div_wcrit_roche = w_div_w_roche_omega(r00,s% m(k),s% omega(k),s% cgrav(k), &
1252 0 : s% w_div_wcrit_max, s% w_div_wcrit_max2, s% w_div_wc_flag)
1253 : end if
1254 :
1255 0 : call eval_i_rot(s, k, r00, w_div_wcrit_roche, s% i_rot(k))
1256 :
1257 0 : end subroutine set1_irot
1258 :
1259 :
1260 : ! this works like set1_xa except shifted to cell edge instead of cell center
1261 0 : subroutine set1_omega(s, k, k_below_just_added, nz, &
1262 0 : old_xout, new_xout, mmax, old_dmbar, new_dmbar, old_j_rot, ierr)
1263 : ! set new value for s% omega(k)
1264 0 : use num_lib, only: binary_search
1265 : type (star_info), pointer :: s
1266 : integer, intent(in) :: k, k_below_just_added, nz
1267 : real(dp), intent(in) :: mmax
1268 : real(dp), dimension(:), intent(in) :: &
1269 : old_xout, new_xout, old_dmbar, new_dmbar, old_j_rot
1270 : integer, intent(out) :: ierr
1271 :
1272 0 : real(dp) :: xm_outer, xm_inner, j_tot, xm0, xm1, new_point_dmbar, &
1273 0 : dm_sum, dm
1274 : integer :: kk, k_outer
1275 :
1276 : integer, parameter :: k_dbg = -1
1277 :
1278 : include 'formats'
1279 :
1280 0 : ierr = 0
1281 :
1282 0 : xm_outer = new_xout(k)
1283 0 : if (k == nz) then
1284 0 : new_point_dmbar = mmax - xm_outer - s% M_center
1285 : else
1286 0 : new_point_dmbar = new_dmbar(k)
1287 : end if
1288 0 : xm_inner = xm_outer + new_point_dmbar
1289 :
1290 0 : if (k == k_dbg) then
1291 0 : write(*,2) 'xm_outer', k, xm_outer
1292 0 : write(*,2) 'xm_inner', k, xm_inner
1293 0 : write(*,2) 'new_point_dmbar', k, new_point_dmbar
1294 : end if
1295 :
1296 : !write(*,*)
1297 : !write(*,2) 'xm_outer', k, xm_outer
1298 :
1299 0 : dm_sum = 0d0
1300 :
1301 0 : if (xm_outer < old_xout(1)) then ! there is some accreted material in new
1302 0 : if (xm_inner <= old_xout(1)) then
1303 : ! new is entirely accreted material
1304 : !write(*,2) 'new is entirely accreted material', k, new_point_dmbar
1305 0 : s% omega(k) = 0
1306 0 : return
1307 : end if
1308 0 : dm = min(new_point_dmbar, old_xout(1) - xm_outer)
1309 0 : dm_sum = dm
1310 0 : j_tot = 0
1311 0 : xm_outer = old_xout(1)
1312 0 : k_outer = 1
1313 : else ! new entirely composed of old material
1314 0 : if (k == k_dbg) write(*,*) 'new entirely composed of old material'
1315 0 : j_tot = 0
1316 0 : if (xm_outer >= old_xout(nz)) then
1317 : ! new contained entirely in old nz
1318 0 : k_outer = nz
1319 : else
1320 : ! binary search for k_outer such that
1321 : ! xm_outer >= old_xout(k_outer)
1322 : ! and old_xout(k_outer+1) > xm_outer
1323 0 : k_outer = binary_search(nz, old_xout, 0, xm_outer)
1324 :
1325 0 : if (k == k_dbg) &
1326 0 : write(*,2) 'k_outer', k_outer, old_xout(k_outer), old_xout(k_outer+1)
1327 :
1328 : ! check
1329 0 : if (k_outer <= 0 .or. k_outer > nz) then
1330 0 : ierr = -1
1331 0 : return
1332 : end if
1333 :
1334 0 : if (xm_outer < old_xout(k_outer)) then
1335 0 : ierr = -1
1336 0 : return
1337 : end if
1338 :
1339 0 : if (k_outer < nz) then
1340 0 : if (old_xout(k_outer+1) <= xm_outer) then
1341 0 : ierr = -1
1342 0 : return
1343 : end if
1344 : end if
1345 :
1346 : end if
1347 : end if
1348 :
1349 0 : if (k == -1) then
1350 0 : ierr = -1
1351 0 : return
1352 : end if
1353 :
1354 0 : do kk = k_outer, nz ! loop until reach xm_inner
1355 0 : xm0 = old_xout(kk)
1356 :
1357 : if (k == k_dbg) write(*,2) 'kk', kk, old_xout(kk), old_xout(kk+1)
1358 :
1359 0 : if (xm0 >= xm_inner) then
1360 0 : if (dm_sum < new_point_dmbar .and. kk > 1) then
1361 : ! need to add a bit more from the previous source
1362 0 : dm = new_point_dmbar - dm_sum
1363 0 : dm_sum = new_point_dmbar
1364 0 : j_tot = j_tot + old_j_rot(kk-1)*dm
1365 : end if
1366 : exit
1367 : end if
1368 :
1369 0 : if (kk == nz) then
1370 0 : xm1 = mmax - s% M_center
1371 : else
1372 0 : xm1 = old_xout(kk+1)
1373 : end if
1374 :
1375 0 : if (xm1 < xm_outer) then
1376 0 : ierr = -1
1377 0 : return
1378 : end if
1379 :
1380 0 : if (xm0 >= xm_outer .and. xm1 <= xm_inner) then ! entire old kk is in new k
1381 :
1382 0 : dm = old_dmbar(kk)
1383 0 : dm_sum = dm_sum + dm
1384 :
1385 0 : if (dm_sum > new_point_dmbar) then
1386 : ! dm too large -- numerical roundoff problems
1387 0 : dm = dm - (new_point_dmbar - dm_sum)
1388 0 : dm_sum = new_point_dmbar
1389 : end if
1390 :
1391 0 : j_tot = j_tot + old_j_rot(kk)*dm
1392 :
1393 0 : else if (xm0 <= xm_outer .and. xm1 >= xm_inner) then ! entire new k is in old kk
1394 :
1395 0 : dm = new_dmbar(k)
1396 0 : dm_sum = dm_sum + dm
1397 0 : j_tot = j_tot + old_j_rot(kk)*dm
1398 :
1399 : else ! only use the part of old kk that is in new k
1400 :
1401 : if (k == k_dbg) then
1402 : write(*,*) 'only use the part of old kk that is in new k', xm_inner <= xm1
1403 : write(*,1) 'xm_outer', xm_outer
1404 : write(*,1) 'xm_inner', xm_inner
1405 : write(*,1) 'xm0', xm0
1406 : write(*,1) 'xm1', xm1
1407 : write(*,1) 'dm_sum', dm_sum
1408 : write(*,1) 'new_point_dmbar', new_point_dmbar
1409 : write(*,1) 'new_point_dmbar - dm_sum', new_point_dmbar - dm_sum
1410 : end if
1411 :
1412 0 : if (xm_inner <= xm1) then ! this is the last part of new k
1413 :
1414 : if (k == k_dbg) write(*,3) 'this is the last part of new k', k, kk
1415 :
1416 0 : dm = new_point_dmbar - dm_sum
1417 0 : dm_sum = new_point_dmbar
1418 :
1419 : else
1420 : ! notice that we avoid this case if possible because of numerical roundoff
1421 :
1422 : if (k == k_dbg) write(*,3) 'we avoid this case if possible', k, kk
1423 :
1424 0 : dm = max(0d0, xm1 - xm_outer)
1425 0 : if (dm_sum + dm > new_point_dmbar) dm = new_point_dmbar - dm_sum
1426 0 : dm_sum = dm_sum + dm
1427 :
1428 : end if
1429 :
1430 0 : j_tot = j_tot + old_j_rot(kk)*dm
1431 :
1432 0 : if (dm <= 0) then
1433 0 : ierr = -1
1434 0 : return
1435 : end if
1436 :
1437 : end if
1438 :
1439 0 : if (dm_sum >= new_point_dmbar) then
1440 : if (k == k_dbg) then
1441 : write(*,2) 'exit for k', k
1442 : write(*,2) 'dm_sum', kk, dm_sum
1443 : write(*,2) 'new_point_dmbar', kk, new_point_dmbar
1444 : end if
1445 0 : dm_sum = new_point_dmbar
1446 0 : exit
1447 : end if
1448 :
1449 : end do
1450 :
1451 0 : if (dm_sum /= new_point_dmbar) then
1452 0 : if (k == nz) then
1453 0 : j_tot = j_tot + old_j_rot(nz)*(new_point_dmbar - dm_sum)
1454 : else
1455 0 : ierr = -1
1456 0 : return
1457 : end if
1458 : end if
1459 :
1460 0 : s% j_rot(k) = j_tot/new_point_dmbar
1461 0 : call set1_irot(s, k, k_below_just_added, .true.)
1462 0 : s% omega(k) = s% j_rot(k)/s% i_rot(k)% val
1463 :
1464 : if (k_dbg == k) then
1465 : write(*,2) 's% omega(k)', k, s% omega(k)
1466 : write(*,2) 's% j_rot(k)', k, s% j_rot(k)
1467 : write(*,2) 's% i_rot(k)% val', k, s% i_rot(k)% val
1468 : call mesa_error(__FILE__,__LINE__,'debugging: set1_omega')
1469 : end if
1470 :
1471 0 : end subroutine set1_omega
1472 :
1473 : ! before mix, remove actual_J_lost - s% angular_momentum_removed
1474 : ! then set s% angular_momentum_removed to actual_J_lost
1475 :
1476 0 : subroutine adjust_J_lost(s, k_below_just_added, starting_j_rot_surf, ierr)
1477 : type (star_info), pointer :: s
1478 : integer, intent(in) :: k_below_just_added
1479 : real(dp), intent(in) :: starting_j_rot_surf
1480 : integer, intent(out) :: ierr
1481 :
1482 0 : real(dp) :: dmm1, dm00, dm, J, mass_lost, j_for_mass_loss, &
1483 0 : actual_J_lost, delta_J, frac, qlast, jnew, J_removed
1484 :
1485 : integer :: k, last_k
1486 :
1487 : include 'formats'
1488 :
1489 0 : ierr = 0
1490 :
1491 0 : mass_lost = s% mstar_old - s% mstar
1492 0 : if (mass_lost <= 0) then
1493 0 : s% adjust_J_q = 1d0 ! nothing to do
1494 0 : return
1495 : end if
1496 :
1497 : ! can use a different j to account for things like wind braking
1498 0 : if (s% use_other_j_for_adjust_J_lost) then
1499 0 : call s% other_j_for_adjust_J_lost(s% id, starting_j_rot_surf, j_for_mass_loss, ierr)
1500 0 : if (ierr /= 0) then
1501 0 : write(*,*) "Error in other_j_for_adjust_J_lost"
1502 0 : return
1503 : end if
1504 : else
1505 0 : j_for_mass_loss = starting_j_rot_surf
1506 : end if
1507 :
1508 : actual_J_lost = &
1509 : s% adjust_J_fraction*mass_lost*j_for_mass_loss + &
1510 0 : (1d0 - s% adjust_J_fraction)*s% angular_momentum_removed
1511 0 : delta_J = actual_J_lost - s% angular_momentum_removed
1512 :
1513 0 : if (delta_J == 0d0) then
1514 : ! this means the model is not rotating, but with rotation flag on, nothing to do
1515 0 : s% adjust_J_q = 1d0
1516 0 : return
1517 : end if
1518 :
1519 0 : dm00 = 0
1520 0 : J = 0
1521 0 : last_k = s% nz
1522 0 : do k = 1, s% nz
1523 0 : dmm1 = dm00
1524 0 : dm00 = s% dm(k)
1525 0 : dm = 0.5d0*(dmm1+dm00)
1526 0 : J = J + dm*s% j_rot(k)
1527 : if (J < s% min_J_div_delta_J*abs(delta_J) &
1528 0 : .or. s% q(k) > s% min_q_for_adjust_J_lost) cycle
1529 : last_k = k
1530 0 : exit
1531 : end do
1532 0 : if (last_k == s% nz) then
1533 0 : ierr = 1
1534 0 : write(*,*) "Error in adjust_J_lost: last_q = nz"
1535 0 : return
1536 : end if
1537 :
1538 : ! adjust in a manner that distributes the shear
1539 : ! at each layer above last_k the specific angular momentum is scaled by a factor
1540 : ! (frac(q-qlast)+1), where frac is computed such that the total change in angular momentum
1541 : ! matches delta_J
1542 : ! NOTE: this could in principle cause counter-rotation depending on the choice of
1543 : ! parameters, but since j_rot normally rises steeply towards the surface its unlikely.
1544 0 : dm00 = 0
1545 0 : frac = 0
1546 0 : qlast = s% q(last_k)
1547 0 : s% adjust_J_q = qlast
1548 0 : do k = 1, last_k-1
1549 0 : dmm1 = dm00
1550 0 : dm00 = s% dm(k)
1551 0 : dm = 0.5d0*(dmm1+dm00)
1552 0 : frac = frac + (s% q(k) - qlast)*s% j_rot(k)*dm
1553 : end do
1554 0 : frac = -delta_J/frac
1555 :
1556 0 : dm00 = 0
1557 0 : J_removed = 0d0
1558 0 : do k = 1, last_k-1
1559 0 : dmm1 = dm00
1560 0 : dm00 = s% dm(k)
1561 0 : dm = 0.5d0*(dmm1+dm00)
1562 0 : jnew = (frac*(s% q(k)-qlast)+1)*s% j_rot(k)
1563 0 : J_removed = J_removed + dm*(s% j_rot(k) - jnew)
1564 0 : s% j_rot(k) = jnew
1565 0 : call set1_irot(s, k, k_below_just_added, .true.)
1566 0 : s% omega(k) = s% j_rot(k)/s% i_rot(k)% val
1567 : end do
1568 :
1569 0 : s% angular_momentum_removed = actual_J_lost
1570 :
1571 0 : end subroutine adjust_J_lost
1572 :
1573 :
1574 0 : subroutine set_D_omega( &
1575 : s, nz, k_const_mass, k_newval, &
1576 0 : rxm_old, rxm_new, delta_m, old_xmstar, new_xmstar, &
1577 0 : D_omega, oldloc, newloc, oldval, newval, work, ierr)
1578 : use interp_1d_lib
1579 : use interp_1d_def
1580 : type (star_info), pointer :: s
1581 : integer, intent(in) :: nz, k_const_mass, k_newval
1582 : real(dp), dimension(:), intent(in) :: rxm_old, rxm_new ! (nz)
1583 : real(dp), intent(in) :: delta_m, old_xmstar, new_xmstar
1584 : real(dp), dimension(:) :: &
1585 : D_omega, oldloc, newloc, oldval, newval
1586 :
1587 : real(dp), pointer :: work(:)
1588 :
1589 : integer, intent(out) :: ierr
1590 :
1591 : integer :: n, nwork, k
1592 : logical :: dbg
1593 :
1594 : include 'formats'
1595 :
1596 0 : ierr = 0
1597 :
1598 0 : dbg = .false.
1599 0 : n = nz ! k_const_mass
1600 0 : nwork = pm_work_size
1601 :
1602 0 : oldloc(1) = 0
1603 0 : do k=2,n
1604 0 : oldloc(k) = rxm_old(k)
1605 : end do
1606 0 : do k=1,n
1607 0 : newloc(k) = rxm_new(k)
1608 0 : oldval(k) = D_omega(k)
1609 : end do
1610 :
1611 : call interpolate_vector( &
1612 : n, oldloc, n, newloc, oldval, newval, interp_pm, nwork, work, &
1613 0 : 'adjust_mass set_D_omega', ierr)
1614 0 : if (ierr /= 0) return
1615 0 : do k=1,k_newval-1
1616 0 : D_omega(k) = 0d0
1617 : end do
1618 0 : do k=k_newval,n
1619 0 : D_omega(k) = newval(k)
1620 : end do
1621 :
1622 0 : end subroutine set_D_omega
1623 :
1624 : end module adjust_mass
|
