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 : ! Routine eval_fp_ft for computing rotation corrections to the stellar structure equations.
21 : ! Following the method of Kippenhahn & Thomas, 1970; Endal & Sofia 1976, as implemented in
22 : ! Paxton et al., 2019, using the updated fits from Fabry, Marchant & Sana, 2022
23 :
24 : module hydro_rotation
25 :
26 : use const_def, only: pi, pi4, ln10, two_thirds, one_third, one_sixth
27 : use star_utils, only: get_r_from_xh
28 : use star_private_def
29 :
30 : implicit none
31 :
32 : private
33 : public :: log_term_power
34 : public :: a
35 : public :: c
36 : public :: w_div_w_roche_jrot
37 : public :: w_div_w_roche_omega
38 : public :: set_rotation_info
39 : public :: compute_j_fluxes_and_extra_jdot
40 : public :: set_surf_avg_rotation_info
41 : public :: use_xh_to_update_i_rot_and_j_rot
42 : public :: set_i_rot_from_omega_and_j_rot
43 : public :: use_xh_to_update_i_rot
44 : public :: update1_i_rot_from_xh
45 : public :: get_rotation_sigmas
46 : public :: set_omega
47 : public :: set_uniform_omega
48 : public :: set_i_rot
49 : public :: eval_i_rot
50 :
51 : real(dp), parameter :: log_term_power = 5.626d0
52 :
53 : contains
54 :
55 : ! compute w_div_w_roche for a known angular frequency omega, rpsi, and Mpsi
56 0 : real(dp) function w_div_w_roche_omega(rpsi, Mphi, omega, cgrav, max_w, max_w2, w_div_wc_flag) result(w_roche)
57 : real(dp), intent(in) :: rpsi, Mphi, omega, cgrav, max_w, max_w2
58 : logical, intent(in) :: w_div_wc_flag
59 0 : real(dp) :: wr, wr_high, wr_low, dimless_rpsi, new_dimless_rpsi, rpsi_lim1, rpsi_lim2
60 :
61 0 : if (omega == 0d0) then
62 0 : w_roche = 0d0
63 : return
64 : end if
65 :
66 0 : dimless_rpsi = rpsi * pow(abs(omega), two_thirds) / pow(cgrav * Mphi, one_third)
67 0 : if (.not. w_div_wc_flag) then
68 : ! verify if w_div_w_roche is not above max, otherwise limit it to that
69 0 : wr = max_w
70 0 : new_dimless_rpsi = pow(wr, two_thirds) * rpsi_from_re_factor(pow2(wr), pow4(wr), pow6(wr))
71 0 : if (dimless_rpsi > new_dimless_rpsi) then
72 0 : w_roche = wr
73 : return
74 : end if
75 : else
76 : ! smoothly cap to max_w2 to get a continuous function
77 : ! nothing is done when we are below max_w, but between max_w and max_w2 we smoothly
78 : ! produce an asymptote that would result in w_div_wc=max_w2 for jrot->infinity
79 0 : wr = max_w
80 0 : rpsi_lim1 = pow(wr, two_thirds) * rpsi_from_re_factor(pow2(wr), pow4(wr), pow6(wr))
81 :
82 0 : wr = max_w2
83 0 : rpsi_lim2 = pow(wr, two_thirds) * rpsi_from_re_factor(pow2(wr), pow4(wr), pow6(wr))
84 :
85 0 : if (abs(dimless_rpsi) > rpsi_lim1) then
86 0 : dimless_rpsi = sigmoid(abs(dimless_rpsi), rpsi_lim1, rpsi_lim2)
87 : end if
88 : end if
89 :
90 : ! otherwise, bisect result
91 0 : wr_high = wr
92 0 : wr_low = 0
93 0 : do while (wr_high - wr_low > 1d-6)
94 0 : wr = 0.5d0 * (wr_high + wr_low)
95 0 : new_dimless_rpsi = pow(wr, two_thirds) * rpsi_from_re_factor(pow2(wr), pow4(wr), pow6(wr))
96 0 : if (dimless_rpsi > new_dimless_rpsi) then
97 : wr_low = wr
98 : else
99 0 : wr_high = wr
100 : end if
101 : end do
102 0 : w_roche = 0.5d0*(wr_high+wr_low)
103 :
104 0 : if (omega < 0d0) then
105 0 : w_roche = -w_roche
106 : end if
107 :
108 : end function w_div_w_roche_omega
109 :
110 : ! compute w_div_w_roche for a known specific angular momentum jrot, rpsi, and Mphi
111 0 : real(dp) function w_div_w_roche_jrot(rpsi, Mphi, jrot, cgrav, max_w, max_w2, w_div_wc_flag) result(w_roche)
112 : real(dp), intent(in) :: rpsi,Mphi,jrot,cgrav, max_w, max_w2
113 : logical, intent(in) :: w_div_wc_flag
114 0 : real(dp) :: wr, wr_high, wr_low, dimless_factor, new_dimless_factor
115 0 : real(dp) :: w2, w4, w6, w_log_term, jr_lim1, jr_lim2
116 :
117 0 : if (jrot == 0d0) then
118 0 : w_roche = 0d0
119 : return
120 : end if
121 :
122 0 : dimless_factor = abs(jrot) / sqrt(cgrav * Mphi * rpsi)
123 :
124 0 : if (.not. w_div_wc_flag) then
125 : ! verify if w_div_w_roche is not above max, otherwise limit it to that
126 0 : wr = max_w
127 0 : w2 = pow2(wr)
128 0 : w4 = pow4(wr)
129 0 : w6 = pow6(wr)
130 0 : w_log_term = log(1d0 - pow(wr, log_term_power))
131 0 : new_dimless_factor = two_thirds * wr * C(w2, w4, w6, w_log_term) / A(w4, w6, w_log_term)
132 0 : if (dimless_factor > new_dimless_factor) then
133 0 : w_roche = wr
134 : return
135 : end if
136 : else
137 : ! smoothly cap to max_w2 to get a continuous function
138 : ! nothing is done when we are below max_w, but between max_w and max_w2 we smoothly
139 : ! produce an asymptote that would result in w_div_wc=max_w2 for jrot->infinity
140 0 : wr = max_w
141 0 : w4 = pow4(wr)
142 0 : w6 = pow6(wr)
143 0 : w_log_term = log(1 - pow(wr, log_term_power))
144 0 : jr_lim1 = two_thirds * wr * C(pow2(wr), w4, w6, w_log_term) / A(w4, w6, w_log_term)
145 :
146 0 : wr = max_w2
147 0 : w4 = pow4(wr)
148 0 : w6 = pow6(wr)
149 0 : w_log_term = log(1 - pow(wr, log_term_power))
150 0 : jr_lim2 = two_thirds * wr * C(pow2(wr), w4, w6, w_log_term) / A(w4, w6, w_log_term)
151 :
152 0 : if (abs(dimless_factor) > jr_lim1) then
153 0 : dimless_factor = sigmoid(abs(dimless_factor), jr_lim1, jr_lim2)
154 : end if
155 : end if
156 :
157 : ! otherwise, bisect result
158 0 : wr_high = wr
159 0 : wr_low = 0
160 0 : do while (wr_high - wr_low > 1d-6)
161 0 : wr = 0.5d0 * (wr_high + wr_low)
162 0 : w2 = pow2(wr)
163 0 : w4 = pow4(wr)
164 0 : w6 = pow6(wr)
165 0 : w_log_term = log(1d0 - pow(wr, log_term_power))
166 0 : new_dimless_factor = two_thirds * wr * C(w2, w4, w6, w_log_term) / A(w4, w6, w_log_term)
167 0 : if (dimless_factor > new_dimless_factor) then
168 : wr_low = wr
169 : else
170 0 : wr_high = wr
171 : end if
172 : end do
173 0 : w_roche = 0.5d0 * (wr_high + wr_low)
174 :
175 0 : if (jrot < 0d0) then
176 0 : w_roche = -w_roche
177 : end if
178 :
179 : end function w_div_w_roche_jrot
180 :
181 0 : subroutine set_i_rot(s, skip_w_div_w_crit_roche)
182 : type (star_info), pointer :: s
183 : logical, intent(in) :: skip_w_div_w_crit_roche
184 : integer :: k
185 : include 'formats'
186 :
187 0 : !$OMP PARALLEL DO PRIVATE(k) SCHEDULE(dynamic,2)
188 : do k=1,s% nz
189 : if (.not. skip_w_div_w_crit_roche) then
190 : s% w_div_w_crit_roche(k) = &
191 : w_div_w_roche_jrot(s% r(k),s% m(k),s% j_rot(k),s% cgrav(k), &
192 : s% w_div_wcrit_max, s% w_div_wcrit_max2, s% w_div_wc_flag)
193 : end if
194 : call eval_i_rot(s, k, s% r(k), s% w_div_w_crit_roche(k), s% i_rot(k))
195 : end do
196 : !$OMP END PARALLEL DO
197 :
198 0 : end subroutine set_i_rot
199 :
200 0 : subroutine set_i_rot_from_omega_and_j_rot(s)
201 : use auto_diff_support
202 : type (star_info), pointer :: s
203 : integer :: k
204 : include 'formats'
205 0 : do k=1,s% nz
206 0 : if (s% omega(k) /= 0d0) then
207 : ! we can directly compute i_rot using j_rot and omega
208 0 : s% i_rot(k) = 0d0
209 0 : s% i_rot(k)% val = s% j_rot(k)/s% omega(k)
210 : else
211 0 : call update1_i_rot_from_xh(s, k)
212 : end if
213 : end do
214 0 : end subroutine set_i_rot_from_omega_and_j_rot
215 :
216 0 : subroutine set_j_rot(s)
217 : type (star_info), pointer :: s
218 : integer :: k
219 : include 'formats'
220 0 : do k=1,s% nz
221 0 : s% j_rot(k) = s% i_rot(k)% val*s% omega(k)
222 : end do
223 0 : end subroutine set_j_rot
224 :
225 0 : subroutine set_omega(s, str)
226 : type (star_info), pointer :: s
227 : character (len=*) :: str
228 : integer :: k
229 : include 'formats'
230 0 : do k=1,s% nz
231 0 : s% omega(k) = s% j_rot(k)/s% i_rot(k)% val
232 : end do
233 0 : end subroutine set_omega
234 :
235 : subroutine check_omega(s, str)
236 : type (star_info), pointer :: s
237 : character (len=*) :: str
238 : integer :: k
239 : logical :: okay
240 : include 'formats'
241 : okay = .true.
242 : do k=1,s% nz
243 : if (abs(s% omega(k) - s% j_rot(k)/s% i_rot(k)% val) > 1d-14) then
244 : write(*,2) 'omega error', k, s% omega(k) - s% j_rot(k)/s% i_rot(k)% val
245 : okay = .false.
246 : exit
247 : end if
248 : end do
249 : if (okay) return
250 : write(*,*) trim(str)
251 : call mesa_error(__FILE__,__LINE__,'check_omega')
252 : end subroutine check_omega
253 :
254 0 : subroutine update1_i_rot_from_xh(s, k)
255 : type (star_info), pointer :: s
256 : integer, intent(in) :: k
257 : real(dp) :: r00
258 : include 'formats'
259 :
260 0 : r00 = get_r_from_xh(s,k)
261 :
262 0 : call eval_i_rot(s, k, r00, s% w_div_w_crit_roche(k), s% i_rot(k))
263 0 : end subroutine update1_i_rot_from_xh
264 :
265 0 : subroutine use_xh_to_update_i_rot(s)
266 : type (star_info), pointer :: s
267 : integer :: k
268 0 : do k=1,s% nz
269 0 : if (s% j_rot(k) /= 0d0) then
270 : s% w_div_w_crit_roche(k) = &
271 : w_div_w_roche_jrot(get_r_from_xh(s,k),s% m(k),s% j_rot(k),s% cgrav(k), &
272 0 : s% w_div_wcrit_max, s% w_div_wcrit_max2, s% w_div_wc_flag)
273 : else
274 0 : s% w_div_w_crit_roche(k) = 0d0
275 : end if
276 : end do
277 0 : do k=1,s% nz
278 0 : call update1_i_rot_from_xh(s,k)
279 : end do
280 0 : end subroutine use_xh_to_update_i_rot
281 :
282 0 : subroutine use_xh_to_update_i_rot_and_j_rot(s)
283 : type (star_info), pointer :: s
284 : integer :: k
285 0 : do k=1,s% nz
286 : s% w_div_w_crit_roche(k) = &
287 : w_div_w_roche_omega(get_r_from_xh(s,k),s% m(k),s% omega(k),s% cgrav(k), &
288 0 : s% w_div_wcrit_max, s% w_div_wcrit_max2, s% w_div_wc_flag)
289 : end do
290 0 : do k=1,s% nz
291 0 : call update1_i_rot_from_xh(s,k)
292 : end do
293 0 : call set_j_rot(s)
294 0 : end subroutine use_xh_to_update_i_rot_and_j_rot
295 :
296 0 : subroutine get_rotation_sigmas(s, nzlo, nzhi, dt, ierr)
297 : type (star_info), pointer :: s
298 : integer, intent(in) :: nzlo, nzhi
299 : real(dp), intent(in) :: dt
300 : integer, intent(out) :: ierr
301 :
302 : integer :: k, nz
303 : real(dp), allocatable :: am_nu(:), am_sig(:)
304 :
305 : include 'formats'
306 :
307 0 : ierr = 0
308 0 : nz = s% nz
309 :
310 0 : allocate(am_nu(nz), am_sig(nz))
311 :
312 0 : call get1_am_sig(s, nzlo, nzhi, s% am_nu_j, s% am_sig_j, dt, ierr)
313 0 : if (ierr /= 0) then
314 0 : if (s% report_ierr) write(*,1) 'failed in get_rotation_sigmas'
315 0 : return
316 : end if
317 :
318 0 : do k=1,nz
319 0 : am_nu(k) = s% am_nu_j(k) + s% am_nu_omega(k)
320 : end do
321 : ! do it this way so apply limit to sum; sum is used as diffusion coeff for omega
322 0 : call get1_am_sig(s, nzlo, nzhi, am_nu, am_sig, dt, ierr)
323 0 : if (ierr /= 0) then
324 0 : if (s% report_ierr) write(*,1) 'failed in get_rotation_sigmas'
325 0 : return
326 : end if
327 :
328 0 : do k=1,nz
329 0 : s% am_sig_omega(k) = max(0d0, am_sig(k) - s% am_sig_j(k))
330 : end do
331 :
332 0 : end subroutine get_rotation_sigmas
333 :
334 0 : subroutine get1_am_sig(s, nzlo, nzhi, am_nu, am_sig, dt, ierr)
335 : type (star_info), pointer :: s
336 : integer, intent(in) :: nzlo, nzhi
337 : real(dp), intent(in) :: dt
338 : real(dp), dimension(:) :: am_nu, am_sig
339 : integer, intent(out) :: ierr
340 :
341 : integer :: k, nz, nz2
342 0 : real(dp) :: r, D, am_nu_E00, am_nu_Ep1, dmbar, s1, &
343 0 : sig_term_limit, xmstar, siglim
344 :
345 : include 'formats'
346 :
347 0 : ierr = 0
348 0 : xmstar = s% xmstar
349 0 : sig_term_limit = s% am_sig_term_limit
350 0 : nz = s% nz
351 : ! note: am_sig is cell centered, so combine adjacent am_nu face values.
352 0 : am_nu_E00 = 0; am_nu_Ep1 = 0
353 0 : nz2 = nzhi
354 0 : if (nzhi == nz) then
355 0 : k = nz
356 0 : D = am_nu_E00
357 0 : r = 0.5d0*s% r(k)
358 0 : s1 = pi4*r*r*s% rho(k)
359 0 : am_sig(k) = s1*s1*D/s% dm(k)
360 0 : nz2 = nz-1
361 : end if
362 0 : do k = nzlo, nz2
363 0 : am_nu_E00 = max(0d0, am_nu(k))
364 0 : am_nu_Ep1 = max(0d0, am_nu(k+1))
365 : ! Meynet, Maeder, & Mowlavi, A&A 416, 1023-1036, 2004, eqn 51 with f = 1/2.
366 0 : D = 2*(am_nu_E00*am_nu_Ep1)/max(1d-99, am_nu_E00 + am_nu_Ep1)
367 0 : r = 0.5d0*(s% r(k) + s% r(k+1)) ! consistent with f = 1/2
368 0 : s1 = pi4*r*r*s% rho(k)
369 0 : am_sig(k) = s1*s1*D/s% dm(k)
370 : end do
371 :
372 : ! can get numerical problems unless limit am_sig
373 : ! adjust am_sig to make sure am_sig*dt/dmbar is < allowed limit
374 0 : do k = nzlo, nzhi
375 0 : if (k < nz) then
376 0 : dmbar = xmstar*min(s% dq(k),s% dq(k+1))
377 0 : siglim = sig_term_limit*dmbar/dt
378 : else
379 0 : dmbar = xmstar*s% dq(k)
380 0 : siglim = sig_term_limit*dmbar/dt
381 : end if
382 0 : if (am_sig(k) > siglim) then
383 0 : am_sig(k) = siglim
384 : end if
385 : end do
386 :
387 0 : end subroutine get1_am_sig
388 :
389 0 : subroutine set_uniform_omega(id, omega, ierr)
390 : use auto_diff_support
391 : integer, intent(in) :: id
392 : real(dp), intent(in) :: omega
393 : integer, intent(out) :: ierr
394 : type (star_info), pointer :: s
395 : integer :: k, nz
396 : include 'formats'
397 0 : call get_star_ptr(id, s, ierr)
398 0 : if (ierr /= 0) return
399 0 : nz = s% nz
400 0 : do k=1, nz
401 0 : s% omega(k) = omega
402 0 : s% fp_rot(k) = 0d0
403 0 : s% ft_rot(k) = 0d0
404 : end do
405 0 : call use_xh_to_update_i_rot_and_j_rot(s)
406 0 : s% w_div_w_crit_roche(1:nz) = 0d0
407 0 : s% r_polar(1:nz) = 0d0
408 0 : s% r_equatorial(1:nz) = 0d0
409 0 : s% am_nu_rot(1:nz) = 0d0
410 0 : s% am_nu_non_rot(1:nz) = 0d0
411 0 : s% am_nu_omega(1:nz) = 0d0
412 0 : s% am_nu_j(1:nz) = 0d0
413 0 : s% am_sig_omega(1:nz) = 0d0
414 0 : s% am_sig_j(1:nz) = 0d0
415 0 : s% domega_dlnR(1:nz) = 0d0
416 0 : s% richardson_number(1:nz) = 0d0
417 0 : s% D_mix_non_rotation(1:nz) = 0d0
418 0 : s% D_visc(1:nz) = 0d0
419 0 : s% D_DSI(1:nz) = 0d0
420 0 : s% D_SH(1:nz) = 0d0
421 0 : s% D_SSI(1:nz) = 0d0
422 0 : s% D_ES(1:nz) = 0d0
423 0 : s% D_GSF(1:nz) = 0d0
424 0 : s% D_ST(1:nz) = 0d0
425 0 : s% nu_ST(1:nz) = 0d0
426 0 : s% omega_shear(1:nz) = 0d0
427 0 : s% dynamo_B_r(1:nz) = 0d0
428 0 : s% dynamo_B_phi(1:nz) = 0d0
429 0 : call set_rotation_info(s, .false., ierr)
430 0 : if (ierr /= 0) return
431 0 : s% need_to_setvars = .true.
432 0 : end subroutine set_uniform_omega
433 :
434 0 : subroutine set_rotation_info(s, skip_w_div_w_crit_roche, ierr)
435 : type (star_info), pointer :: s
436 : logical, intent(in) :: skip_w_div_w_crit_roche
437 : integer, intent(out) :: ierr
438 : include 'formats'
439 0 : ierr = 0
440 :
441 0 : if (.not. s% rotation_flag) return
442 :
443 0 : call set_i_rot(s, skip_w_div_w_crit_roche)
444 0 : call set_omega(s, 'set_rotation_info')
445 :
446 0 : if (.not. s% use_other_eval_fp_ft) then
447 : call eval_fp_ft( &
448 : s% id, s% nz, s% m, s% r, s% rho, s% omega, s% ft_rot, s% fp_rot, &
449 0 : s% r_polar, s% r_equatorial, s% report_ierr, ierr)
450 : else
451 : call s% other_eval_fp_ft( &
452 : s% id, s% nz, s% m, s% r, s% rho, s% omega, s% ft_rot, s% fp_rot, &
453 0 : s% r_polar, s% r_equatorial, s% report_ierr, ierr)
454 : end if
455 0 : if (ierr /= 0) then
456 0 : write(*,*) 'failed in eval_fp_ft'
457 : end if
458 0 : end subroutine set_rotation_info
459 :
460 33 : subroutine set_surf_avg_rotation_info(s)
461 : use star_utils, only: get_Lrad_div_Ledd
462 : type (star_info), pointer :: s
463 : real(dp) :: &
464 33 : dm, dmsum, omega_sum, omega_crit_sum, omega_div_omega_crit_sum, &
465 33 : v_rot_sum, v_crit_sum, v_div_v_crit_sum, Lrad_div_Ledd_sum, &
466 33 : gamma_factor, omega_crit, omega, kap_sum, &
467 33 : j_rot_sum, j_rot, v_rot, v_crit, Lrad_div_Ledd, dtau, tau, &
468 33 : cgrav, kap, mmid, Lmid, rmid, logT_sum, logRho_sum
469 : integer :: k, ierr
470 : logical, parameter :: dbg = .false.
471 : include 'formats'
472 :
473 33 : if (.not. s% rotation_flag) then
474 33 : s% omega_avg_surf = 0
475 33 : s% omega_crit_avg_surf = 0
476 33 : s% w_div_w_crit_avg_surf = 0
477 33 : s% j_rot_avg_surf = 0
478 33 : s% v_rot_avg_surf = 0
479 33 : s% v_crit_avg_surf = 0
480 33 : s% v_div_v_crit_avg_surf = 0
481 33 : s% Lrad_div_Ledd_avg_surf = 0
482 33 : s% opacity_avg_surf = 0
483 33 : s% logT_avg_surf = 0
484 33 : s% logRho_avg_surf = 0
485 : return
486 : end if
487 :
488 : ierr = 0
489 0 : call set_rotation_info(s,.true.,ierr)
490 0 : if (ierr /= 0) then
491 0 : write(*,*) 'got ierr from call set_rotation_info in set_surf_avg_rotation_info'
492 0 : write(*,*) 'just ignore it'
493 : end if
494 :
495 0 : tau = s% tau_factor*s% tau_base
496 0 : dmsum = 0d0
497 0 : Lrad_div_Ledd_sum = 0d0
498 0 : rmid = 0d0
499 :
500 0 : do k = 1, s% nz - 1
501 0 : kap = s% opacity(k)
502 0 : rmid = s% rmid(k)
503 0 : mmid = 0.5d0*(s% m_grav(k) + s% m_grav(k+1))
504 0 : Lmid = 0.5d0*(s% L(k) + s% L(k+1))
505 0 : cgrav = 0.5d0*(s% cgrav(k) + s% cgrav(k+1))
506 0 : dm = s% dm(k)
507 0 : dtau = dm*kap/(pi4*rmid*rmid)
508 :
509 0 : if (tau + dtau <= s% surf_avg_tau_min) then
510 : tau = tau + dtau
511 : cycle
512 : end if
513 :
514 : ! check for partial contribution from cell
515 : ! the tau < s% surf_avg_tau is meant for the case in which the surface tau is set
516 : ! equal or larger to surf_avg_tau. In that case we just use the values of the surface cell.
517 0 : if (tau < s% surf_avg_tau) then
518 0 : if (tau < s% surf_avg_tau_min) then ! only use part of this cell
519 0 : dm = dm*(tau + dtau - s% surf_avg_tau_min)/dtau
520 0 : else if (tau + dtau > s% surf_avg_tau) then ! only use part of this cell
521 0 : dm = dm*(s% surf_avg_tau - tau)/dtau
522 : !write(*,2) 'tau limit', k, (s% surf_avg_tau - tau)/dtau
523 : end if
524 : end if
525 0 : dmsum = dmsum + dm
526 0 : Lrad_div_Ledd = get_Lrad_div_Ledd(s,k)
527 0 : Lrad_div_Ledd_sum = Lrad_div_Ledd_sum + dm*Lrad_div_Ledd
528 0 : tau = tau + dtau
529 0 : if (tau >= s% surf_avg_tau) exit
530 : end do
531 :
532 0 : s% Lrad_div_Ledd_avg_surf = Lrad_div_Ledd_sum/dmsum
533 0 : gamma_factor = 1d0 - min(s% Lrad_div_Ledd_avg_surf, 0.9999d0)
534 :
535 0 : tau = s% tau_factor*s% tau_base
536 0 : dmsum = 0
537 0 : j_rot_sum = 0
538 0 : omega_sum = 0
539 0 : omega_crit_sum = 0
540 0 : omega_div_omega_crit_sum = 0
541 0 : v_rot_sum = 0
542 0 : v_crit_sum = 0
543 0 : v_div_v_crit_sum = 0
544 0 : kap_sum = 0
545 0 : logT_sum = 0
546 0 : logRho_sum = 0
547 :
548 0 : do k = 1, s% nz - 1
549 :
550 0 : kap = s% opacity(k)
551 : ! TODO: better explain
552 : ! Use equatorial radius
553 0 : rmid = 0.5d0*(s% r_equatorial(k) + s% r_equatorial(k+1))
554 0 : dm = s% dm(k)
555 0 : dtau = dm*kap/(pi4*rmid*rmid)
556 :
557 0 : if (tau + dtau <= s% surf_avg_tau_min) then
558 : tau = tau + dtau
559 : cycle
560 : end if
561 :
562 : ! check for partial contribution from cell
563 : ! the tau < s% surf_avg_tau is meant for the case in which the surface tau is set
564 : ! equal or larger to surf_avg_tau. In this case we just use the values of the surface cell.
565 0 : if (tau < s% surf_avg_tau) then
566 0 : if (tau < s% surf_avg_tau_min) then ! only use part of this cell
567 0 : dm = dm*(tau + dtau - s% surf_avg_tau_min)/dtau
568 0 : else if (tau + dtau > s% surf_avg_tau) then ! only use part of this cell
569 0 : dm = dm*(s% surf_avg_tau - tau)/dtau
570 : end if
571 : end if
572 :
573 0 : dmsum = dmsum + dm
574 0 : cgrav = 0.5d0*(s% cgrav(k) + s% cgrav(k+1))
575 0 : mmid = 0.5d0*(s% m_grav(k) + s% m_grav(k+1))
576 0 : omega = 0.5d0*(s% omega(k) + s% omega(k+1))
577 0 : j_rot = 0.5d0*(s% j_rot(k) + s% j_rot(k+1))
578 :
579 0 : kap_sum = kap_sum + dm*kap
580 0 : j_rot_sum = j_rot_sum + dm*j_rot
581 :
582 0 : omega_crit = sqrt(gamma_factor*cgrav*mmid/pow3(rmid))
583 0 : omega_div_omega_crit_sum = omega_div_omega_crit_sum + dm*abs(omega/omega_crit)
584 :
585 0 : v_rot = omega*rmid
586 0 : v_crit = omega_crit*rmid
587 0 : omega_sum = omega_sum + dm*omega
588 0 : omega_crit_sum = omega_crit_sum + dm*omega_crit
589 0 : v_rot_sum = v_rot_sum + dm*v_rot
590 0 : v_crit_sum = v_crit_sum + dm*v_crit
591 0 : v_div_v_crit_sum = v_div_v_crit_sum + dm*abs(v_rot/v_crit)
592 0 : logT_sum = logT_sum + dm*s% lnT(k)/ln10
593 0 : logRho_sum = logRho_sum + dm*s% lnd(k)/ln10
594 0 : kap_sum = kap_sum + dm*kap
595 0 : tau = tau + dtau
596 0 : if (tau >= s% surf_avg_tau) exit
597 :
598 : end do
599 :
600 0 : s% logT_avg_surf = logT_sum/dmsum
601 0 : s% logRho_avg_surf = logRho_sum/dmsum
602 0 : s% opacity_avg_surf = kap_sum/dmsum
603 0 : s% j_rot_avg_surf = j_rot_sum/dmsum
604 0 : s% omega_avg_surf = omega_sum/dmsum
605 0 : s% omega_crit_avg_surf = omega_crit_sum/dmsum
606 0 : s% w_div_w_crit_avg_surf = omega_div_omega_crit_sum/dmsum
607 0 : s% v_rot_avg_surf = v_rot_sum/dmsum
608 0 : s% v_crit_avg_surf = v_crit_sum/dmsum
609 0 : s% v_div_v_crit_avg_surf = v_div_v_crit_sum/dmsum
610 :
611 33 : end subroutine set_surf_avg_rotation_info
612 :
613 : ! Input variables:
614 : ! N Number of meshpoints used by the model (arrays are this size)
615 : ! XM Mass coordinate [gram]
616 : ! R Radius coordinate [cm]
617 : ! RHO Density [gram/cm^3]
618 : ! AW Angular velocity [rad/sec]
619 : ! Output variables:
620 : ! Correction factor FT at each meshpoint
621 : ! Correction factor FP at each meshpoint
622 : ! r_polar, r_equatorial at each meshpoint
623 0 : subroutine eval_fp_ft(id, nz, xm, r, rho, aw, ft, fp, r_polar, r_equatorial, report_ierr, ierr)
624 33 : use num_lib
625 : use auto_diff_support
626 : integer, intent(in) :: id
627 : integer, intent(in) :: nz
628 : real(dp), intent(in) :: aw(:), r(:), rho(:), xm(:) ! (nz)
629 : type(auto_diff_real_star_order1), intent(out) :: ft(:), fp(:) ! (nz)
630 : real(dp), intent(inout) :: r_polar(:), r_equatorial(:) ! (nz)
631 : logical, intent(in) :: report_ierr
632 : integer, intent(out) :: ierr
633 :
634 : type (star_info), pointer :: s
635 : integer :: j
636 :
637 : logical :: dbg
638 :
639 : type (auto_diff_real_1var_order1) :: A_omega, w, w2, w4, w6, w_log_term, fp_temp, ft_temp
640 :
641 : include 'formats'
642 :
643 : ierr = 0
644 0 : call star_ptr(id, s, ierr)
645 0 : if (ierr /= 0) return
646 :
647 0 : dbg = .false. ! (s% model_number >= 5)
648 :
649 0 : !$OMP PARALLEL DO PRIVATE(j, A_omega, fp_temp, ft_temp, w, w2, w4, w6, w_log_term) SCHEDULE(dynamic,2)
650 : do j=1, s% nz
651 : ! Compute fp, ft, re and rp using fits to the Roche geometry of a single star.
652 : ! by this point in the code, w_div_w_crit_roche is set
653 : w = abs(s% w_div_w_crit_roche(j))
654 : w% d1val1 = 1d0
655 :
656 : w2 = pow2(w)
657 : w4 = pow4(w)
658 : w6 = pow6(w)
659 : w_log_term = log(1d0 - pow(w, log_term_power))
660 : ! cannot use real function below because need derivatives
661 : A_omega = 1d0 + 0.3293d0 * w4 - 0.4926d0 * w6 - 0.5560d0 * w_log_term
662 :
663 : ! fits for fp, ft; Fabry+2022, Eqs. A.10, A.11
664 : fp_temp = (1d0 - two_thirds * w2 - 0.2133d0 * w4 - 0.1068d0 * w6) / A_omega
665 : ft_temp = (1d0 - 0.07955d0 * w4 - 0.2322d0 * w6) / A_omega
666 :
667 : ! re and rp can be derived analytically from w_div_wcrit
668 : r_equatorial(j) = r(j) * re_from_rpsi_factor(w2% val, w4% val, w6% val)
669 : r_polar(j) = r_equatorial(j) / (1d0 + 0.5d0 * w2% val)
670 :
671 : ! Be sure they are consistent with r_Psi
672 : r_equatorial(j) = max(r_equatorial(j), r(j))
673 : r_polar(j) = min(r_polar(j), r(j))
674 :
675 : fp(j) = 0d0
676 : ft(j) = 0d0
677 : fp(j)% val = fp_temp% val
678 : ft(j)% val = ft_temp% val
679 : if (s% w_div_wc_flag) then
680 : fp(j)% d1Array(i_w_div_wc_00) = fp_temp% d1val1
681 : ft(j)% d1Array(i_w_div_wc_00) = ft_temp% d1val1
682 : end if
683 : !if (j == s% solver_test_partials_k) then
684 : ! s% solver_test_partials_val = fp(j)
685 : ! s% solver_test_partials_var = s% i_w_div_wc
686 : ! s% solver_test_partials_dval_dx = s% dfp_rot_dw_div_wc(j)
687 : !end if
688 : end do
689 : !$OMP END PARALLEL DO
690 :
691 0 : end subroutine eval_fp_ft
692 :
693 0 : subroutine eval_i_rot(s, k, r00, w_div_w_crit_roche, i_rot)
694 0 : use auto_diff_support
695 : type (star_info), pointer :: s
696 : integer, intent(in) :: k ! just for debugging
697 : real(dp), intent(in) :: r00, w_div_w_crit_roche
698 : type(auto_diff_real_star_order1), intent(out) :: i_rot
699 :
700 : type(auto_diff_real_2var_order1) :: ir, r, re, w, w2, w4, w6, w_log_term
701 :
702 : include 'formats'
703 :
704 0 : i_rot = 0d0
705 0 : if (s% use_other_eval_i_rot) then
706 0 : call s% other_eval_i_rot(s% id, k, r00, w_div_w_crit_roche, i_rot)
707 0 : else if (s% simple_i_rot_flag) then
708 0 : i_rot = two_thirds * r00 * r00
709 0 : i_rot% d1Array(i_lnR_00) = 2d0 * i_rot% val
710 0 : i_rot% d1Array(i_w_div_wc_00) = 0d0
711 : else
712 0 : w = abs(w_div_w_crit_roche)
713 0 : w% d1val1 = 1d0
714 :
715 0 : r = r00
716 0 : r% d1val2 = r00 ! Makes the independent variable lnR
717 :
718 0 : w2 = pow2(w)
719 0 : w4 = pow4(w)
720 0 : w6 = pow6(w)
721 0 : w_log_term = log(1d0 - pow(w, log_term_power))
722 : ! cannot use real functions below, since need derivatives
723 0 : re = r * (1d0 + one_sixth * w2 - 0.005124d0 * w4 + 0.06562d0 * w6)
724 : ! Compute i_rot following Fabry+2022, Eq. A.9
725 : ir = two_thirds * pow2(re) * (1d0 + w2 / 5d0 + 0.4140d0 * w4 - 0.8650d0 * w6 - 0.8370d0 * w_log_term) &
726 0 : / (1d0 + 0.3293d0 * w4 - 0.4926d0 * w6 - 0.5560d0 * w_log_term)
727 :
728 0 : i_rot = 0d0
729 0 : i_rot = ir% val
730 0 : i_rot% d1Array(i_w_div_wc_00) = ir% d1val1
731 0 : i_rot% d1Array(i_lnR_00) = ir% d1val2
732 : end if
733 :
734 0 : end subroutine eval_i_rot
735 :
736 0 : subroutine compute_j_fluxes_and_extra_jdot(id, ierr)
737 0 : use auto_diff_support
738 : integer, intent(in) :: id
739 : integer, intent(out) :: ierr
740 : type (star_info), pointer :: s
741 : type(auto_diff_real_star_order1) :: omega00, omegap1, r00, rp1, i_rot00, i_rotp1, rho00, part1
742 :
743 : integer :: k
744 : logical :: dbg
745 :
746 0 : real(dp) :: pi2_div4
747 : ierr = 0
748 0 : dbg = .false.
749 :
750 0 : call star_ptr(id, s, ierr)
751 0 : if (ierr /= 0) return
752 0 : if (s% am_D_mix_factor==0d0) then
753 0 : s% am_nu_omega(:) = 0d0
754 : end if
755 :
756 0 : s% extra_jdot(:) = 0d0
757 0 : if (s% use_other_torque) then
758 0 : call s% other_torque(s% id, ierr)
759 0 : if (ierr /= 0) then
760 0 : if (s% report_ierr .or. dbg) &
761 0 : write(*, *) 'solve_omega_mix: other_torque returned ierr', ierr
762 0 : return
763 : end if
764 : end if
765 0 : if (associated(s% binary_other_torque)) then
766 0 : call s% binary_other_torque(s% id, ierr)
767 0 : if (ierr /= 0) then
768 0 : if (s% report_ierr .or. dbg) &
769 0 : write(*, *) 'solve_omega_mix: binary_other_torque returned ierr', ierr
770 0 : return
771 : end if
772 : end if
773 :
774 0 : pi2_div4 = pow2(pi)/4d0
775 0 : do k=1, s% nz-1
776 :
777 0 : r00 = wrap_r_00(s, k)
778 0 : rp1 = wrap_r_p1(s, k)
779 0 : rho00 = wrap_d_00(s, k)
780 :
781 0 : omega00 = wrap_omega_00(s, k)
782 0 : omegap1 = wrap_omega_p1(s, k)
783 :
784 0 : i_rot00 = s% i_rot(k)
785 0 : i_rotp1 = shift_p1(s% i_rot(k+1))
786 :
787 0 : part1 = -pi2_div4*pow4(r00+rp1)*pow2(rho00)*(i_rot00+i_rotp1)*(s% am_nu_omega(k)+s% am_nu_omega(k+1))
788 0 : s% j_flux(k) = part1*(omega00-omegap1)/s% dm(k)
789 :
790 : !! this is to test partials
791 : !if (k==188) then
792 : ! part1 = (omega00-omegap1)/s% dm(k)
793 : ! s% solver_test_partials_val = part1% val
794 : ! s% solver_test_partials_var = s% i_lnR !s% i_w_div_wc
795 : ! s% solver_test_partials_dval_dx = part1% d1Array(i_lnR_00) !part1% d1Array(i_w_div_wc_00)
796 : !end if
797 :
798 : end do
799 0 : s% j_flux(s% nz) = 0d0
800 0 : end subroutine compute_j_fluxes_and_extra_jdot
801 :
802 0 : real(dp) function rpsi_from_re_factor(o2, o4, o6)
803 : real(dp), intent(in) :: o2
804 : real(dp), intent(in) :: o4
805 : real(dp), intent(in) :: o6
806 : ! Fabry+2022, Eq. A.2
807 0 : rpsi_from_re_factor = 1d0 - one_sixth * o2 + 0.02025d0 * o4 - 0.03870d0 * o6
808 0 : end function rpsi_from_re_factor
809 :
810 : real(dp) function re_from_rpsi_factor(o2, o4, o6)
811 : real(dp), intent(in) :: o2
812 : real(dp), intent(in) :: o4
813 : real(dp), intent(in) :: o6
814 : ! Fabry+2022, Eq. A.3
815 : re_from_rpsi_factor = 1d0 + one_sixth * o2 - 0.005124d0 * o4 + 0.06562d0 * o6
816 : end function re_from_rpsi_factor
817 :
818 0 : real(dp) function A(o4, o6, o_log_term)
819 : real(dp), intent(in) :: o4
820 : real(dp), intent(in) :: o6
821 : real(dp), intent(in) :: o_log_term
822 : ! Fabry+2022, Eq. A.6
823 0 : A = 1d0 + 0.3293d0 * o4 - 0.4926d0 * o6 - 0.5560d0 * o_log_term
824 0 : end function A
825 :
826 : real(dp) function B(o2, o4, o6, o_log_term)
827 : real(dp), intent(in) :: o2
828 : real(dp), intent(in) :: o4
829 : real(dp), intent(in) :: o6
830 : real(dp), intent(in) :: o_log_term
831 : ! Fabry+2022, Eq. A.7
832 : B = 1d0 + o2 / 5d0 + 0.4140d0 * o4 - 0.8650d0 * o6 - 0.8370d0 * o_log_term
833 : end function B
834 :
835 0 : real(dp) function C(o2, o4, o6, o_log_term)
836 : real(dp), intent(in) :: o2
837 : real(dp), intent(in) :: o4
838 : real(dp), intent(in) :: o6
839 : real(dp), intent(in) :: o_log_term
840 : ! Fabry+2022, Eq. A.12
841 0 : C = 1d0 + 17d0 / 60d0 * o2 + 0.4010d0 * o4 - 0.8606d0 * o6 - 0.9305d0 * o_log_term
842 0 : end function C
843 :
844 0 : real(dp) function sigmoid(x, xmax1, xmax2)
845 : real(dp), intent(in) :: x, xmax1, xmax2
846 : ! smoothly maps abs(x) = [xmax1, infty] to sigmoid = [xmax1, xmax2]
847 0 : sigmoid = 2d0 * (xmax2 - xmax1) / (1d0 + exp(-2d0 * (abs(x) - xmax1) / (xmax2 - xmax1))) - xmax2 + 2d0 * xmax1
848 0 : end function sigmoid
849 :
850 : end module hydro_rotation
|
