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 (MESA V), using the updated fits from Fabry, Marchant & Sana, 2022, A&A 661:A123.
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 eval_i_rot(id, r00, w_div_w_crit_roche, i_rot)
182 : use auto_diff_support
183 :
184 : integer, intent(in) :: id
185 : real(dp), intent(in) :: r00,w_div_w_crit_roche
186 : type (auto_diff_real_star_order1), intent(out) :: i_rot
187 :
188 : type (star_info), pointer :: s
189 : type (auto_diff_real_2var_order1) :: ir, r, re, w, w2, w4, w6, lg_one_sub_w4, B, A
190 : integer :: ierr
191 :
192 : include 'formats'
193 :
194 0 : call star_ptr(id, s, ierr)
195 0 : if (ierr /= 0) return
196 :
197 0 : i_rot = 0d0
198 0 : if (s% simple_i_rot_flag) then
199 0 : i_rot = two_thirds*r00*r00
200 0 : i_rot% d1Array(i_lnR_00) = 2*i_rot% val
201 0 : i_rot% d1Array(i_w_div_wc_00) = 0d0
202 : else
203 : ! Compute i_rot following Paxton et al. 2019 (ApJs, 243, 10)
204 0 : w = w_div_w_crit_roche
205 0 : w%d1val1 = 1d0
206 :
207 0 : r = r00
208 0 : r%d1val2 = r00 ! Makes the independent variable lnR
209 :
210 0 : w2 = pow2(w)
211 0 : w4 = pow4(w)
212 0 : w6 = pow6(w)
213 0 : lg_one_sub_w4 = log(1d0-w4)
214 0 : re = r*(1d0+w2/6d0-0.0002507d0*w4+0.06075d0*w6)
215 0 : B = (1d0+w2/5d0-0.2735d0*w4-0.4327d0*w6-3d0/2d0*0.5583d0*lg_one_sub_w4)
216 0 : A = (1d0-0.1076d0*w4-0.2336d0*w6-0.5583d0*lg_one_sub_w4)
217 :
218 0 : ir = two_thirds*pow2(re)*B/A
219 :
220 0 : i_rot = 0d0
221 0 : i_rot = ir%val
222 0 : i_rot% d1Array(i_w_div_wc_00) = ir%d1val1
223 0 : i_rot% d1Array(i_lnR_00) = ir%d1val2
224 : end if
225 :
226 0 : end subroutine eval_i_rot
227 :
228 :
229 0 : subroutine set_i_rot(s, skip_w_div_w_crit_roche)
230 : type (star_info), pointer :: s
231 : logical, intent(in) :: skip_w_div_w_crit_roche
232 : integer :: k
233 : include 'formats'
234 :
235 0 : !$OMP PARALLEL DO PRIVATE(k) SCHEDULE(dynamic,2)
236 : do k=1,s% nz
237 : if (.not. skip_w_div_w_crit_roche) then
238 : s% w_div_w_crit_roche(k) = &
239 : w_div_w_roche_jrot(s% r(k), s% m(k), s% j_rot(k), s% cgrav(k), &
240 : s% w_div_wcrit_max, s% w_div_wcrit_max2, s% w_div_wc_flag)
241 : end if
242 : call set1_i_rot(s, k, s% r(k))
243 : end do
244 : !$OMP END PARALLEL DO
245 0 : end subroutine set_i_rot
246 :
247 0 : subroutine set_i_rot_from_omega_and_j_rot(s)
248 : use auto_diff_support
249 : type (star_info), pointer :: s
250 : integer :: k
251 : include 'formats'
252 :
253 0 : !$OMP PARALLEL DO PRIVATE(k) SCHEDULE(dynamic, 2)
254 : do k=1,s% nz
255 : if (s% omega(k) /= 0d0) then
256 : ! we can directly compute i_rot using j_rot and omega
257 : s% i_rot(k) = 0d0
258 : s% i_rot(k)% val = s% j_rot(k)/s% omega(k)
259 : else
260 : call update1_i_rot_from_xh(s, k)
261 : end if
262 : end do
263 : !$OMP END PARALLEL DO
264 0 : end subroutine set_i_rot_from_omega_and_j_rot
265 :
266 0 : subroutine set_j_rot(s)
267 : type (star_info), pointer :: s
268 : integer :: k
269 : include 'formats'
270 0 : !$OMP PARALLEL DO PRIVATE(k) SCHEDULE(dynamic, 2)
271 : do k=1,s% nz
272 : s% j_rot(k) = s% i_rot(k)% val*s% omega(k)
273 : end do
274 : !$OMP END PARALLEL DO
275 0 : end subroutine set_j_rot
276 :
277 0 : subroutine set_omega(s, str)
278 : type (star_info), pointer :: s
279 : character (len=*) :: str
280 : integer :: k
281 : include 'formats'
282 0 : !$OMP PARALLEL DO PRIVATE(k) SCHEDULE(dynamic, 2)
283 : do k=1,s% nz
284 : s% omega(k) = s% j_rot(k)/s% i_rot(k)% val
285 : end do
286 : !$OMP END PARALLEL DO
287 0 : end subroutine set_omega
288 :
289 : subroutine check_omega(s, str)
290 : type (star_info), pointer :: s
291 : character (len=*) :: str
292 : integer :: k
293 : logical :: okay
294 : include 'formats'
295 : okay = .true.
296 : do k=1,s% nz
297 : if (abs(s% omega(k) - s% j_rot(k)/s% i_rot(k)% val) > 1d-14) then
298 : write(*,2) 'omega error', k, s% omega(k) - s% j_rot(k)/s% i_rot(k)% val
299 : okay = .false.
300 : exit
301 : end if
302 : end do
303 : if (okay) return
304 : write(*,*) trim(str)
305 : call mesa_error(__FILE__,__LINE__,'check_omega')
306 : end subroutine check_omega
307 :
308 0 : subroutine update1_i_rot_from_xh(s, k)
309 : type (star_info), pointer :: s
310 : integer, intent(in) :: k
311 : real(dp) :: r00
312 : include 'formats'
313 0 : r00 = get_r_from_xh(s, k)
314 0 : call set1_i_rot(s, k, r00)
315 0 : end subroutine update1_i_rot_from_xh
316 :
317 0 : subroutine set1_i_rot(s, k, r)
318 : use auto_diff
319 : type (star_info), pointer :: s
320 : real(dp), intent(in) :: r
321 : integer, intent(in) :: k
322 :
323 : type (auto_diff_real_star_order1) :: i_rot_single, i_rot_tidal
324 :
325 0 : if (s% use_other_eval_i_rot) then
326 0 : call s% other_eval_i_rot(s% id, k, r, s% w_div_w_crit_roche(k), s% i_rot(k))
327 : else
328 0 : call eval_i_rot(s% id, r, s% w_div_w_crit_roche(k), i_rot_single)
329 0 : if (associated(s% binary_other_irot)) then
330 0 : call s% binary_other_irot(s% id, r, i_rot_tidal)
331 0 : call blend_tidal_values(s, k, i_rot_single, i_rot_tidal, s% omega(k), s% i_rot(k))
332 : else
333 0 : s% i_rot(k) = i_rot_single
334 : end if
335 : end if
336 0 : end subroutine set1_i_rot
337 :
338 0 : subroutine blend_tidal_values(s, k, single, tidal, omega, final)
339 0 : use auto_diff
340 :
341 : type (star_info), pointer :: s
342 : integer, intent(in) :: k
343 : real(dp), intent(in) :: omega
344 : type (auto_diff_real_star_order1), intent(in) :: single, tidal
345 : type (auto_diff_real_star_order1), intent(out) :: final
346 :
347 0 : real(dp) :: tidal_frac
348 : integer :: ierr
349 :
350 : include 'formats'
351 :
352 : ! do blending with tidal case according to synchronicity
353 0 : if (associated(s% binary_deformation_switch_fraction)) then
354 0 : call s% binary_deformation_switch_fraction(s% id, k, omega, tidal_frac, ierr) ! tells how much of the tidal value we should use
355 0 : if (ierr /= 0) then
356 0 : if (s% report_ierr) write(*, 1) "failed in blend_tidal_values"
357 : end if
358 : end if
359 0 : final = tidal_frac * tidal + (1d0 - tidal_frac) * single
360 0 : end subroutine blend_tidal_values
361 :
362 0 : subroutine use_xh_to_update_i_rot(s)
363 : type (star_info), pointer :: s
364 : integer :: k
365 0 : real(dp) :: r00
366 :
367 0 : !$OMP PARALLEL DO PRIVATE(k, r00) SCHEDULE(dynamic,2)
368 : do k=1,s% nz
369 : if (s% j_rot(k) /= 0d0) then
370 : r00 = get_r_from_xh(s,k)
371 : s% w_div_w_crit_roche(k) = &
372 : w_div_w_roche_jrot(r00, s% m(k), s% j_rot(k), s% cgrav(k), &
373 : s% w_div_wcrit_max, s% w_div_wcrit_max2, s% w_div_wc_flag)
374 : else
375 : s% w_div_w_crit_roche(k) = 0d0
376 : end if
377 : call update1_i_rot_from_xh(s, k)
378 : end do
379 : !$OMP END PARALLEL DO
380 0 : end subroutine use_xh_to_update_i_rot
381 :
382 0 : subroutine use_xh_to_update_i_rot_and_j_rot(s)
383 : type (star_info), pointer :: s
384 : integer :: k
385 0 : real(dp) :: r00
386 :
387 0 : !$OMP PARALLEL DO PRIVATE(k, r00) SCHEDULE(dynamic,2)
388 : do k=1,s% nz
389 : r00 = get_r_from_xh(s,k)
390 : s% w_div_w_crit_roche(k) = &
391 : w_div_w_roche_omega(r00, s% m(k), s% omega(k), s% cgrav(k), &
392 : s% w_div_wcrit_max, s% w_div_wcrit_max2, s% w_div_wc_flag)
393 : call update1_i_rot_from_xh(s,k)
394 : end do
395 : !$OMP END PARALLEL DO
396 0 : call set_j_rot(s)
397 0 : end subroutine use_xh_to_update_i_rot_and_j_rot
398 :
399 0 : subroutine get_rotation_sigmas(s, nzlo, nzhi, dt, ierr)
400 : type (star_info), pointer :: s
401 : integer, intent(in) :: nzlo, nzhi
402 : real(dp), intent(in) :: dt
403 : integer, intent(out) :: ierr
404 :
405 : integer :: k, nz
406 : real(dp), allocatable :: am_nu(:), am_sig(:)
407 :
408 : include 'formats'
409 :
410 0 : ierr = 0
411 0 : nz = s% nz
412 :
413 0 : allocate(am_nu(nz), am_sig(nz))
414 :
415 0 : call get1_am_sig(s, nzlo, nzhi, s% am_nu_j, s% am_sig_j, dt, ierr)
416 0 : if (ierr /= 0) then
417 0 : if (s% report_ierr) write(*,1) 'failed in get_rotation_sigmas'
418 0 : return
419 : end if
420 :
421 0 : do k=1,nz
422 0 : am_nu(k) = s% am_nu_j(k) + s% am_nu_omega(k)
423 : end do
424 : ! do it this way so apply limit to sum; sum is used as diffusion coeff for omega
425 0 : call get1_am_sig(s, nzlo, nzhi, am_nu, am_sig, dt, ierr)
426 0 : if (ierr /= 0) then
427 0 : if (s% report_ierr) write(*,1) 'failed in get_rotation_sigmas'
428 0 : return
429 : end if
430 :
431 0 : do k=1,nz
432 0 : s% am_sig_omega(k) = max(0d0, am_sig(k) - s% am_sig_j(k))
433 : end do
434 :
435 0 : end subroutine get_rotation_sigmas
436 :
437 0 : subroutine get1_am_sig(s, nzlo, nzhi, am_nu, am_sig, dt, ierr)
438 : type (star_info), pointer :: s
439 : integer, intent(in) :: nzlo, nzhi
440 : real(dp), intent(in) :: dt
441 : real(dp), dimension(:) :: am_nu, am_sig
442 : integer, intent(out) :: ierr
443 :
444 : integer :: k, nz, nz2
445 0 : real(dp) :: r, D, am_nu_E00, am_nu_Ep1, dmbar, s1, &
446 0 : sig_term_limit, xmstar, siglim
447 :
448 : include 'formats'
449 :
450 0 : ierr = 0
451 0 : xmstar = s% xmstar
452 0 : sig_term_limit = s% am_sig_term_limit
453 0 : nz = s% nz
454 : ! note: am_sig is cell centered, so combine adjacent am_nu face values.
455 0 : am_nu_E00 = 0; am_nu_Ep1 = 0
456 0 : nz2 = nzhi
457 0 : if (nzhi == nz) then
458 0 : k = nz
459 0 : D = am_nu_E00
460 0 : r = 0.5d0*s% r(k)
461 0 : s1 = pi4*r*r*s% rho(k)
462 0 : am_sig(k) = s1*s1*D/s% dm(k)
463 0 : nz2 = nz-1
464 : end if
465 0 : do k = nzlo, nz2
466 0 : am_nu_E00 = max(0d0, am_nu(k))
467 0 : am_nu_Ep1 = max(0d0, am_nu(k+1))
468 : ! Meynet, Maeder, & Mowlavi, A&A 416, 1023-1036, 2004, eqn 51 with f = 1/2.
469 0 : D = 2*(am_nu_E00*am_nu_Ep1)/max(1d-99, am_nu_E00 + am_nu_Ep1)
470 0 : r = 0.5d0*(s% r(k) + s% r(k+1)) ! consistent with f = 1/2
471 0 : s1 = pi4*r*r*s% rho(k)
472 0 : am_sig(k) = s1*s1*D/s% dm(k)
473 : end do
474 :
475 : ! can get numerical problems unless limit am_sig
476 : ! adjust am_sig to make sure am_sig*dt/dmbar is < allowed limit
477 0 : do k = nzlo, nzhi
478 0 : if (k < nz) then
479 0 : dmbar = xmstar*min(s% dq(k),s% dq(k+1))
480 0 : siglim = sig_term_limit*dmbar/dt
481 : else
482 0 : dmbar = xmstar*s% dq(k)
483 0 : siglim = sig_term_limit*dmbar/dt
484 : end if
485 0 : if (am_sig(k) > siglim) then
486 0 : am_sig(k) = siglim
487 : end if
488 : end do
489 :
490 0 : end subroutine get1_am_sig
491 :
492 0 : subroutine set_uniform_omega(id, omega, ierr)
493 : use auto_diff_support
494 : integer, intent(in) :: id
495 : real(dp), intent(in) :: omega
496 : integer, intent(out) :: ierr
497 : type (star_info), pointer :: s
498 : integer :: k, nz
499 : include 'formats'
500 0 : call get_star_ptr(id, s, ierr)
501 0 : if (ierr /= 0) return
502 0 : nz = s% nz
503 0 : do k=1, nz
504 0 : s% omega(k) = omega
505 0 : s% fp_rot(k) = 0d0
506 0 : s% ft_rot(k) = 0d0
507 : end do
508 0 : call use_xh_to_update_i_rot_and_j_rot(s)
509 0 : s% w_div_w_crit_roche(1:nz) = 0d0
510 0 : s% r_polar(1:nz) = 0d0
511 0 : s% r_equatorial(1:nz) = 0d0
512 0 : s% am_nu_rot(1:nz) = 0d0
513 0 : s% am_nu_non_rot(1:nz) = 0d0
514 0 : s% am_nu_omega(1:nz) = 0d0
515 0 : s% am_nu_j(1:nz) = 0d0
516 0 : s% am_sig_omega(1:nz) = 0d0
517 0 : s% am_sig_j(1:nz) = 0d0
518 0 : s% domega_dlnR(1:nz) = 0d0
519 0 : s% richardson_number(1:nz) = 0d0
520 0 : s% D_mix_non_rotation(1:nz) = 0d0
521 0 : s% D_visc(1:nz) = 0d0
522 0 : s% D_DSI(1:nz) = 0d0
523 0 : s% D_SH(1:nz) = 0d0
524 0 : s% D_SSI(1:nz) = 0d0
525 0 : s% D_ES(1:nz) = 0d0
526 0 : s% D_GSF(1:nz) = 0d0
527 0 : s% D_ST(1:nz) = 0d0
528 0 : s% nu_ST(1:nz) = 0d0
529 0 : s% omega_shear(1:nz) = 0d0
530 0 : s% dynamo_B_r(1:nz) = 0d0
531 0 : s% dynamo_B_phi(1:nz) = 0d0
532 0 : call set_rotation_info(s, .false., ierr)
533 0 : if (ierr /= 0) return
534 0 : s% need_to_setvars = .true.
535 0 : end subroutine set_uniform_omega
536 :
537 0 : subroutine set_rotation_info(s, skip_w_div_w_crit_roche, ierr)
538 0 : use auto_diff
539 : type (star_info), pointer :: s
540 : logical, intent(in) :: skip_w_div_w_crit_roche
541 : integer, intent(out) :: ierr
542 :
543 : type (auto_diff_real_star_order1) :: fp_single, fp_tidal, ft_single, ft_tidal
544 : integer :: k
545 : include 'formats'
546 0 : ierr = 0
547 :
548 0 : if (.not. s% rotation_flag) return
549 :
550 0 : call set_i_rot(s, skip_w_div_w_crit_roche)
551 0 : call set_omega(s, 'set_rotation_info')
552 :
553 0 : !$OMP PARALLEL DO PRIVATE(k, fp_single, ft_single, fp_tidal, ft_tidal, ierr) SCHEDULE(dynamic,2)
554 : do k=1, s% nz
555 : if (s% use_other_eval_fp_ft) then
556 : call s% other_eval_fp_ft( &
557 : s% id, k, s% m(k), s% r(k), s% rho(k), s% omega(k), s% fp_rot(k), s% ft_rot(k), &
558 : s% r_polar(k), s% r_equatorial(k), s% report_ierr, ierr)
559 : else
560 : call eval_fp_ft( &
561 : s% id, s% r(k), s% w_div_w_crit_roche(k), fp_single, ft_single, &
562 : s% r_polar(k), s% r_equatorial(k), s% report_ierr, ierr)
563 : !if (k == s% solver_test_partials_k) then
564 : ! s% solver_test_partials_val = fp_single
565 : ! s% solver_test_partials_var = s% i_w_div_wc
566 : ! s% solver_test_partials_dval_dx = s% dfp_rot_dw_div_wc(k)
567 : !end if
568 : if (associated(s% binary_other_fp_ft)) then
569 : call s% binary_other_fp_ft(&
570 : s% id, s% r(k), fp_tidal, ft_tidal, s% r_polar(k), s% r_equatorial(k), &
571 : s% report_ierr, ierr)
572 : call blend_tidal_values(s, k, fp_single, fp_tidal, s% omega(k), s% fp_rot(k))
573 : call blend_tidal_values(s, k, ft_single, ft_tidal, s% omega(k), s% ft_rot(k))
574 : else
575 : s% fp_rot(k) = fp_single
576 : s% ft_rot(k) = ft_single
577 : end if
578 :
579 : if (ierr /= 0) then
580 : write(*,1) 'failed in eval_fp_ft'
581 : end if
582 :
583 : end if
584 : end do
585 : !$OMP END PARALLEL DO
586 0 : end subroutine set_rotation_info
587 :
588 33 : subroutine set_surf_avg_rotation_info(s)
589 0 : use star_utils, only: get_Lrad_div_Ledd
590 : type (star_info), pointer :: s
591 : real(dp) :: &
592 33 : dm, dmsum, omega_sum, omega_crit_sum, omega_div_omega_crit_sum, &
593 33 : v_rot_sum, v_crit_sum, v_div_v_crit_sum, Lrad_div_Ledd_sum, &
594 33 : gamma_factor, omega_crit, omega, kap_sum, &
595 33 : j_rot_sum, j_rot, v_rot, v_crit, Lrad_div_Ledd, dtau, tau, &
596 33 : cgrav, kap, mmid, Lmid, rmid, logT_sum, logRho_sum
597 : integer :: k, ierr
598 : logical, parameter :: dbg = .false.
599 : include 'formats'
600 :
601 33 : if (.not. s% rotation_flag) then
602 33 : s% omega_avg_surf = 0
603 33 : s% omega_crit_avg_surf = 0
604 33 : s% w_div_w_crit_avg_surf = 0
605 33 : s% j_rot_avg_surf = 0
606 33 : s% v_rot_avg_surf = 0
607 33 : s% v_crit_avg_surf = 0
608 33 : s% v_div_v_crit_avg_surf = 0
609 33 : s% Lrad_div_Ledd_avg_surf = 0
610 33 : s% opacity_avg_surf = 0
611 33 : s% logT_avg_surf = 0
612 33 : s% logRho_avg_surf = 0
613 : return
614 : end if
615 :
616 : ierr = 0
617 0 : call set_rotation_info(s,.true.,ierr)
618 0 : if (ierr /= 0) then
619 0 : write(*,*) 'got ierr from call set_rotation_info in set_surf_avg_rotation_info'
620 0 : write(*,*) 'just ignore it'
621 : end if
622 :
623 0 : tau = s% tau_factor*s% tau_base
624 0 : dmsum = 0d0
625 0 : Lrad_div_Ledd_sum = 0d0
626 0 : rmid = 0d0
627 :
628 0 : do k = 1, s% nz - 1
629 0 : kap = s% opacity(k)
630 0 : rmid = s% rmid(k)
631 0 : mmid = 0.5d0*(s% m_grav(k) + s% m_grav(k+1))
632 0 : Lmid = 0.5d0*(s% L(k) + s% L(k+1))
633 0 : cgrav = 0.5d0*(s% cgrav(k) + s% cgrav(k+1))
634 0 : dm = s% dm(k)
635 0 : dtau = dm*kap/(pi4*rmid*rmid)
636 :
637 0 : if (tau + dtau <= s% surf_avg_tau_min) then
638 : tau = tau + dtau
639 : cycle
640 : end if
641 :
642 : ! check for partial contribution from cell
643 : ! the tau < s% surf_avg_tau is meant for the case in which the surface tau is set
644 : ! equal or larger to surf_avg_tau. In that case we just use the values of the surface cell.
645 0 : if (tau < s% surf_avg_tau) then
646 0 : if (tau < s% surf_avg_tau_min) then ! only use part of this cell
647 0 : dm = dm*(tau + dtau - s% surf_avg_tau_min)/dtau
648 0 : else if (tau + dtau > s% surf_avg_tau) then ! only use part of this cell
649 0 : dm = dm*(s% surf_avg_tau - tau)/dtau
650 : !write(*,2) 'tau limit', k, (s% surf_avg_tau - tau)/dtau
651 : end if
652 : end if
653 0 : dmsum = dmsum + dm
654 0 : Lrad_div_Ledd = get_Lrad_div_Ledd(s,k)
655 0 : Lrad_div_Ledd_sum = Lrad_div_Ledd_sum + dm*Lrad_div_Ledd
656 0 : tau = tau + dtau
657 0 : if (tau >= s% surf_avg_tau) exit
658 : end do
659 :
660 0 : s% Lrad_div_Ledd_avg_surf = Lrad_div_Ledd_sum/dmsum
661 0 : gamma_factor = 1d0 - min(s% Lrad_div_Ledd_avg_surf, 0.9999d0)
662 :
663 0 : tau = s% tau_factor*s% tau_base
664 0 : dmsum = 0
665 0 : j_rot_sum = 0
666 0 : omega_sum = 0
667 0 : omega_crit_sum = 0
668 0 : omega_div_omega_crit_sum = 0
669 0 : v_rot_sum = 0
670 0 : v_crit_sum = 0
671 0 : v_div_v_crit_sum = 0
672 0 : kap_sum = 0
673 0 : logT_sum = 0
674 0 : logRho_sum = 0
675 :
676 0 : do k = 1, s% nz - 1
677 :
678 0 : kap = s% opacity(k)
679 : ! TODO: better explain
680 : ! Use equatorial radius
681 0 : rmid = 0.5d0*(s% r_equatorial(k) + s% r_equatorial(k+1))
682 0 : dm = s% dm(k)
683 0 : dtau = dm*kap/(pi4*rmid*rmid)
684 :
685 0 : if (tau + dtau <= s% surf_avg_tau_min) then
686 : tau = tau + dtau
687 : cycle
688 : end if
689 :
690 : ! check for partial contribution from cell
691 : ! the tau < s% surf_avg_tau is meant for the case in which the surface tau is set
692 : ! equal or larger to surf_avg_tau. In this case we just use the values of the surface cell.
693 0 : if (tau < s% surf_avg_tau) then
694 0 : if (tau < s% surf_avg_tau_min) then ! only use part of this cell
695 0 : dm = dm*(tau + dtau - s% surf_avg_tau_min)/dtau
696 0 : else if (tau + dtau > s% surf_avg_tau) then ! only use part of this cell
697 0 : dm = dm*(s% surf_avg_tau - tau)/dtau
698 : end if
699 : end if
700 :
701 0 : dmsum = dmsum + dm
702 0 : cgrav = 0.5d0*(s% cgrav(k) + s% cgrav(k+1))
703 0 : mmid = 0.5d0*(s% m_grav(k) + s% m_grav(k+1))
704 0 : omega = 0.5d0*(s% omega(k) + s% omega(k+1))
705 0 : j_rot = 0.5d0*(s% j_rot(k) + s% j_rot(k+1))
706 :
707 0 : kap_sum = kap_sum + dm*kap
708 0 : j_rot_sum = j_rot_sum + dm*j_rot
709 :
710 0 : omega_crit = sqrt(gamma_factor*cgrav*mmid/pow3(rmid))
711 0 : omega_div_omega_crit_sum = omega_div_omega_crit_sum + dm*abs(omega/omega_crit)
712 :
713 0 : v_rot = omega*rmid
714 0 : v_crit = omega_crit*rmid
715 0 : omega_sum = omega_sum + dm*omega
716 0 : omega_crit_sum = omega_crit_sum + dm*omega_crit
717 0 : v_rot_sum = v_rot_sum + dm*v_rot
718 0 : v_crit_sum = v_crit_sum + dm*v_crit
719 0 : v_div_v_crit_sum = v_div_v_crit_sum + dm*abs(v_rot/v_crit)
720 0 : logT_sum = logT_sum + dm*s% lnT(k)/ln10
721 0 : logRho_sum = logRho_sum + dm*s% lnd(k)/ln10
722 0 : kap_sum = kap_sum + dm*kap
723 0 : tau = tau + dtau
724 0 : if (tau >= s% surf_avg_tau) exit
725 :
726 : end do
727 :
728 0 : s% logT_avg_surf = logT_sum/dmsum
729 0 : s% logRho_avg_surf = logRho_sum/dmsum
730 0 : s% opacity_avg_surf = kap_sum/dmsum
731 0 : s% j_rot_avg_surf = j_rot_sum/dmsum
732 0 : s% omega_avg_surf = omega_sum/dmsum
733 0 : s% omega_crit_avg_surf = omega_crit_sum/dmsum
734 0 : s% w_div_w_crit_avg_surf = omega_div_omega_crit_sum/dmsum
735 0 : s% v_rot_avg_surf = v_rot_sum/dmsum
736 0 : s% v_crit_avg_surf = v_crit_sum/dmsum
737 0 : s% v_div_v_crit_avg_surf = v_div_v_crit_sum/dmsum
738 :
739 33 : end subroutine set_surf_avg_rotation_info
740 :
741 : ! Input variables:
742 : ! r Radius coordinate [cm]
743 : ! aw fractional critical angular velocity Omega/Omega_crit
744 : ! Output variables:
745 : ! Correction factor Fp at each meshpoint
746 : ! Correction factor Ft at each meshpoint
747 : ! r_polar, r_equatorial at each meshpoint
748 0 : subroutine eval_fp_ft(id, r, aw, fp, ft, r_polar, r_equatorial, report_ierr, ierr)
749 33 : use num_lib
750 : use auto_diff_support
751 : integer, intent(in) :: id
752 : real(dp), intent(in) :: r, aw
753 : type(auto_diff_real_star_order1), intent(out) :: fp, ft
754 : real(dp), intent(inout) :: r_polar, r_equatorial
755 : logical, intent(in) :: report_ierr
756 : integer, intent(out) :: ierr
757 :
758 : type (star_info), pointer :: s
759 : integer :: j
760 :
761 : logical :: dbg
762 :
763 : type(auto_diff_real_1var_order1) :: A_omega, fp_numerator, ft_numerator, &
764 : w, w2, w3, w4, w5, w6, w_log_term, fp_temp, ft_temp
765 :
766 : include 'formats'
767 :
768 : ierr = 0
769 0 : call star_ptr(id, s, ierr)
770 0 : if (ierr /= 0) return
771 :
772 0 : dbg = .false. ! (s% model_number >= 5)
773 :
774 : !Compute fp, ft, re and rp using fits to the Roche geometry of a single star.
775 0 : w = abs(aw) ! single rotation deformation is symmetric to prograde/retrograde
776 0 : w% d1val1 = 1d0
777 :
778 0 : w2 = pow2(w)
779 0 : w4 = pow4(w)
780 0 : w6 = pow6(w)
781 0 : w_log_term = log(1d0 - pow(w, log_term_power))
782 : ! cannot use real function below because these are auto_diff variables
783 0 : A_omega = 1d0 + 0.3293d0 * w4 - 0.4926d0 * w6 - 0.5560d0 * w_log_term
784 :
785 : ! fits for fp, ft; Fabry+2022, Eqs. A.10, A.11
786 0 : fp_temp = (1d0 - two_thirds * w2 - 0.2133d0 * w4 - 0.1068d0 * w6) / A_omega
787 0 : ft_temp = (1d0 - 0.07955d0 * w4 - 0.2322d0 * w6) / A_omega
788 :
789 : ! re and rp can be derived analytically from w_div_wcrit
790 0 : r_equatorial = r * re_from_rpsi_factor(w2% val, w4% val, w6% val)
791 0 : r_polar = r_equatorial / (1d0 + 0.5d0 * w2% val)
792 :
793 : ! Be sure they are consistent with r_Psi
794 0 : r_equatorial = max(r_equatorial, r)
795 0 : r_polar = min(r_polar, r)
796 :
797 0 : fp = 0d0
798 0 : ft = 0d0
799 0 : fp% val = fp_temp% val
800 0 : ft% val = ft_temp% val
801 0 : if (s% w_div_wc_flag) then
802 0 : fp% d1Array(i_w_div_wc_00) = fp_temp% d1val1
803 0 : ft% d1Array(i_w_div_wc_00) = ft_temp% d1val1
804 : end if
805 :
806 0 : end subroutine eval_fp_ft
807 :
808 0 : subroutine compute_j_fluxes_and_extra_jdot(id, ierr)
809 0 : use auto_diff_support
810 : integer, intent(in) :: id
811 : integer, intent(out) :: ierr
812 : type (star_info), pointer :: s
813 : type(auto_diff_real_star_order1) :: omega00, omegap1, r00, rp1, i_rot00, i_rotp1, rho00, part1
814 :
815 : integer :: k
816 : logical :: dbg
817 :
818 0 : real(dp) :: pi2_div4
819 : ierr = 0
820 0 : dbg = .false.
821 :
822 0 : call star_ptr(id, s, ierr)
823 0 : if (ierr /= 0) return
824 0 : if (s% am_D_mix_factor==0d0) then
825 0 : s% am_nu_omega(:) = 0d0
826 : end if
827 :
828 0 : s% extra_jdot(:) = 0d0
829 0 : if (s% use_other_torque) then
830 0 : call s% other_torque(s% id, ierr)
831 0 : if (ierr /= 0) then
832 0 : if (s% report_ierr .or. dbg) &
833 0 : write(*, *) 'solve_omega_mix: other_torque returned ierr', ierr
834 0 : return
835 : end if
836 : end if
837 0 : if (associated(s% binary_other_torque)) then
838 0 : call s% binary_other_torque(s% id, ierr)
839 0 : if (ierr /= 0) then
840 0 : if (s% report_ierr .or. dbg) &
841 0 : write(*, *) 'solve_omega_mix: binary_other_torque returned ierr', ierr
842 0 : return
843 : end if
844 : end if
845 :
846 0 : pi2_div4 = pow2(pi)/4d0
847 0 : do k=1, s% nz-1
848 :
849 0 : r00 = wrap_r_00(s, k)
850 0 : rp1 = wrap_r_p1(s, k)
851 0 : rho00 = wrap_d_00(s, k)
852 :
853 0 : omega00 = wrap_omega_00(s, k)
854 0 : omegap1 = wrap_omega_p1(s, k)
855 :
856 0 : i_rot00 = s% i_rot(k)
857 0 : i_rotp1 = shift_p1(s% i_rot(k+1))
858 :
859 0 : part1 = -pi2_div4*pow4(r00+rp1)*pow2(rho00)*(i_rot00+i_rotp1)*(s% am_nu_omega(k)+s% am_nu_omega(k+1))
860 0 : s% j_flux(k) = part1*(omega00-omegap1)/s% dm(k)
861 :
862 : !! this is to test partials
863 : !if (k==188) then
864 : ! part1 = (omega00-omegap1)/s% dm(k)
865 : ! s% solver_test_partials_val = part1% val
866 : ! s% solver_test_partials_var = s% i_lnR !s% i_w_div_wc
867 : ! s% solver_test_partials_dval_dx = part1% d1Array(i_lnR_00) !part1% d1Array(i_w_div_wc_00)
868 : !end if
869 :
870 : end do
871 0 : s% j_flux(s% nz) = 0d0
872 0 : end subroutine compute_j_fluxes_and_extra_jdot
873 :
874 0 : real(dp) function rpsi_from_re_factor(o2, o4, o6)
875 : real(dp), intent(in) :: o2
876 : real(dp), intent(in) :: o4
877 : real(dp), intent(in) :: o6
878 : ! Fabry+2022, Eq. A.2
879 0 : rpsi_from_re_factor = 1d0 - one_sixth * o2 + 0.02025d0 * o4 - 0.03870d0 * o6
880 0 : end function rpsi_from_re_factor
881 :
882 0 : real(dp) function re_from_rpsi_factor(o2, o4, o6)
883 : real(dp), intent(in) :: o2
884 : real(dp), intent(in) :: o4
885 : real(dp), intent(in) :: o6
886 : ! Fabry+2022, Eq. A.3
887 0 : re_from_rpsi_factor = 1d0 + one_sixth * o2 - 0.005124d0 * o4 + 0.06562d0 * o6
888 : end function re_from_rpsi_factor
889 :
890 0 : real(dp) function A(o4, o6, o_log_term)
891 : real(dp), intent(in) :: o4
892 : real(dp), intent(in) :: o6
893 : real(dp), intent(in) :: o_log_term
894 : ! Fabry+2022, Eq. A.6
895 0 : A = 1d0 + 0.3293d0 * o4 - 0.4926d0 * o6 - 0.5560d0 * o_log_term
896 0 : end function A
897 :
898 : real(dp) function B(o2, o4, o6, o_log_term)
899 : real(dp), intent(in) :: o2
900 : real(dp), intent(in) :: o4
901 : real(dp), intent(in) :: o6
902 : real(dp), intent(in) :: o_log_term
903 : ! Fabry+2022, Eq. A.7
904 : B = 1d0 + o2 / 5d0 + 0.4140d0 * o4 - 0.8650d0 * o6 - 0.8370d0 * o_log_term
905 : end function B
906 :
907 0 : real(dp) function C(o2, o4, o6, o_log_term)
908 : real(dp), intent(in) :: o2
909 : real(dp), intent(in) :: o4
910 : real(dp), intent(in) :: o6
911 : real(dp), intent(in) :: o_log_term
912 : ! Fabry+2022, Eq. A.12
913 0 : C = 1d0 + 17d0 / 60d0 * o2 + 0.4010d0 * o4 - 0.8606d0 * o6 - 0.9305d0 * o_log_term
914 0 : end function C
915 :
916 0 : real(dp) function sigmoid(x, xmax1, xmax2)
917 : real(dp), intent(in) :: x, xmax1, xmax2
918 : ! smoothly maps abs(x) = [xmax1, infty] to sigmoid = [xmax1, xmax2]
919 0 : sigmoid = 2d0 * (xmax2 - xmax1) / (1d0 + exp(-2d0 * (abs(x) - xmax1) / (xmax2 - xmax1))) - xmax2 + 2d0 * xmax1
920 0 : end function sigmoid
921 :
922 : end module hydro_rotation
|
