Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2021 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 turb_support
21 :
22 : use star_private_def
23 : use const_def, only: dp, crad, no_mixing
24 : use num_lib
25 : use utils_lib
26 : use auto_diff_support
27 : use star_utils
28 : use turb
29 :
30 : implicit none
31 :
32 : private
33 : public :: get_gradT
34 : public :: do1_mlt_eval
35 : public :: Get_results
36 :
37 : contains
38 :
39 : !> Determines if it is safe (physically) to use TDC instead of MLT.
40 : !!
41 : !! Currently we only know we have to fall back to MLT in cells that get touched
42 : !! by adjust_mass, because there the convection speeds at the start of the
43 : !! step can be badly out of whack.
44 : !!
45 : !! @param s star pointer
46 : !! @param k face index
47 : !! @param fallback False if we can use TDC, True if we can fall back to MLT.
48 65210 : logical function check_if_must_fall_back_to_MLT(s, k) result(fallback)
49 : type (star_info), pointer :: s
50 : integer, intent(in) :: k
51 :
52 65210 : fallback = .false.
53 0 : if (abs(s%mstar_dot) > 1d-99 .and. k < s% k_const_mass) then
54 : fallback = .true.
55 : end if
56 : end function check_if_must_fall_back_to_MLT
57 :
58 0 : subroutine get_gradT(s, MLT_option, & ! used to create models
59 : r, L, T, P, opacity, rho, chiRho, chiT, Cp, gradr, grada, scale_height, &
60 : iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, &
61 : mixing_type, gradT, Y_face, conv_vel, D, Gamma, ierr)
62 : type (star_info), pointer :: s
63 : character (len=*), intent(in) :: MLT_option
64 : real(dp), intent(in) :: &
65 : r, L, T, P, opacity, rho, chiRho, chiT, Cp, gradr, grada, scale_height, &
66 : XH1, cgrav, m, gradL_composition_term, mixing_length_alpha
67 : integer, intent(in) :: iso
68 : real(dp), intent(out) :: gradT, Y_face, conv_vel, D, Gamma
69 : integer, intent(out) :: mixing_type, ierr
70 : type(auto_diff_real_star_order1) :: &
71 : gradr_ad, grada_ad, scale_height_ad, gradT_ad, Y_face_ad, mlt_vc_ad, D_ad, &
72 : Gamma_ad, r_ad, L_ad, T_ad, P_ad, opacity_ad, rho_ad, dV_ad, chiRho_ad, chiT_ad, Cp_ad
73 0 : ierr = 0
74 0 : r_ad = r
75 0 : L_ad = L
76 0 : T_ad = T
77 0 : P_ad = P
78 0 : opacity_ad = opacity
79 0 : rho_ad = rho
80 0 : dV_ad = 0d0
81 0 : chiRho_ad = chiRho
82 0 : chiT_ad = chiT
83 0 : Cp_ad = Cp
84 0 : gradr_ad = gradr
85 0 : grada_ad = grada
86 0 : scale_height_ad = scale_height
87 : call Get_results(s, 0, MLT_option, &
88 : r_ad, L_ad, T_ad, P_ad, opacity_ad, rho_ad, dV_ad, chiRho_ad, chiT_ad, Cp_ad, &
89 : gradr_ad, grada_ad, scale_height_ad, &
90 : iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, &
91 : s% alpha_semiconvection, s% thermohaline_coeff, &
92 0 : mixing_type, gradT_ad, Y_face_ad, mlt_vc_ad, D_ad, Gamma_ad, ierr)
93 0 : gradT = gradT_ad%val
94 0 : Y_face = Y_face_ad%val
95 0 : conv_vel = mlt_vc_ad%val
96 0 : D = D_ad%val
97 0 : Gamma = Gamma_ad%val
98 0 : end subroutine get_gradT
99 :
100 :
101 79994 : subroutine do1_mlt_eval( &
102 : s, k, MLT_option, gradL_composition_term, &
103 : gradr_in, grada, scale_height, mixing_length_alpha, &
104 : mixing_type, gradT, Y_face, mlt_vc, D, Gamma, ierr)
105 : use chem_def, only: ih1
106 : use starspots, only: starspot_tweak_gradr
107 : type (star_info), pointer :: s
108 : integer, intent(in) :: k
109 : character (len=*), intent(in) :: MLT_option
110 : type(auto_diff_real_star_order1), intent(in) :: gradr_in, grada, scale_height
111 : real(dp), intent(in) :: gradL_composition_term, mixing_length_alpha
112 : integer, intent(out) :: mixing_type
113 : type(auto_diff_real_star_order1), intent(out) :: &
114 : gradT, Y_face, mlt_vc, D, Gamma
115 : integer, intent(out) :: ierr
116 :
117 : real(dp) :: cgrav, m, XH1
118 : integer :: iso
119 : type(auto_diff_real_star_order1) :: gradr, r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp
120 : include 'formats'
121 79994 : ierr = 0
122 :
123 79994 : gradr = gradr_in
124 :
125 79994 : cgrav = s% cgrav(k)
126 79994 : m = s% m_grav(k)
127 79994 : L = wrap_L_00(s,k)
128 79994 : T = get_T_face(s,k)
129 79994 : P = get_Peos_face(s,k)
130 79994 : r = wrap_r_00(s,k)
131 79994 : opacity = get_kap_face(s,k)
132 79994 : rho = get_Rho_face(s,k)
133 79994 : dV = 1d0/rho - 1d0/s% rho_start(k)
134 79994 : chiRho = get_ChiRho_face(s,k)
135 79994 : chiT = get_ChiT_face(s,k)
136 79994 : Cp = get_Cp_face(s,k)
137 79994 : iso = s% dominant_iso_for_thermohaline(k)
138 79994 : XH1 = s% xa(s% net_iso(ih1),k)
139 :
140 79994 : if (s% use_other_mlt_results) then
141 : call s% other_mlt_results(s% id, k, MLT_option, &
142 : r, L, T, P, opacity, rho, chiRho, chiT, Cp, gradr, grada, scale_height, &
143 : iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, &
144 : s% alpha_semiconvection, s% thermohaline_coeff, &
145 0 : mixing_type, gradT, Y_face, mlt_vc, D, Gamma, ierr)
146 : else
147 : ! starspot YREC routine
148 79994 : if (s% do_starspots) then
149 : !dV = 0d0 ! dV = 1/rho - 1/rho_start and we assume rho = rho_start.
150 0 : call starspot_tweak_gradr(s, P, gradr_in, gradr)
151 : end if
152 : call Get_results(s, k, MLT_option, &
153 : r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, &
154 : iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, &
155 : s% alpha_semiconvection, s% thermohaline_coeff, &
156 79994 : mixing_type, gradT, Y_face, mlt_vc, D, Gamma, ierr)
157 : end if
158 :
159 79994 : end subroutine do1_mlt_eval
160 :
161 :
162 79994 : subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg
163 : r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, &
164 : iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, &
165 : alpha_semiconvection, thermohaline_coeff, &
166 : mixing_type, gradT, Y_face, conv_vel, D, Gamma, ierr)
167 79994 : use star_utils
168 : type (star_info), pointer :: s
169 : integer, intent(in) :: k
170 : character (len=*), intent(in) :: MLT_option
171 : type(auto_diff_real_star_order1), intent(in) :: &
172 : r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height
173 : integer, intent(in) :: iso
174 : real(dp), intent(in) :: &
175 : XH1, cgrav, m, gradL_composition_term, &
176 : mixing_length_alpha, alpha_semiconvection, thermohaline_coeff
177 : integer, intent(out) :: mixing_type
178 : type(auto_diff_real_star_order1), intent(out) :: gradT, Y_face, conv_vel, D, Gamma
179 : integer, intent(out) :: ierr
180 :
181 : type(auto_diff_real_star_order1) :: Pr, Pg, grav, Lambda, gradL, beta
182 79994 : real(dp) :: conv_vel_start, scale, max_conv_vel
183 :
184 : ! these are used by use_superad_reduction
185 79994 : real(dp) :: Gamma_limit, scale_value1, scale_value2, diff_grads_limit, reduction_limit, lambda_limit
186 : type(auto_diff_real_star_order1) :: Lrad_div_Ledd, Gamma_inv_threshold, Gamma_factor, alfa0, &
187 : diff_grads_factor, Gamma_term, exp_limit, grad_scale, gradr_scaled
188 :
189 : logical :: test_partials, using_TDC
190 : logical, parameter :: report = .false.
191 : include 'formats'
192 :
193 : ! Pre-calculate some things.
194 79994 : Pr = crad*pow4(T)/3d0
195 79994 : Pg = P - Pr
196 79994 : beta = Pg / P
197 79994 : Lambda = mixing_length_alpha*scale_height
198 79994 : grav = cgrav*m/pow2(r)
199 79994 : max_conv_vel = 1d99
200 79994 : if (s% use_Ledoux_criterion) then
201 0 : gradL = grada + gradL_composition_term ! Ledoux temperature gradient
202 : else
203 : gradL = grada
204 : end if
205 :
206 : ! maximum convection velocity.
207 79994 : if (k>=1) then
208 79994 : if (s% q(k) <= s% max_conv_vel_div_csound_maxq) then
209 0 : max_conv_vel = s% csound_face(k) * s% max_conv_vel_div_csound
210 : else
211 : max_conv_vel = 1d99
212 : end if
213 : end if
214 :
215 :
216 : ! Initialize with no mixing
217 79994 : mixing_type = no_mixing
218 79994 : gradT = gradr
219 79994 : Y_face = gradT - gradL
220 79994 : conv_vel = 0d0
221 79994 : D = 0d0
222 79994 : Gamma = 0d0
223 79994 : if (k /= 0) s% superad_reduction_factor(k) = 1d0
224 :
225 : ! Bail if we asked for no mixing, or if parameters are bad.
226 : if (MLT_option == 'none' .or. beta < 1d-10 .or. mixing_length_alpha <= 0d0 .or. &
227 : opacity%val < 1d-10 .or. P%val < 1d-20 .or. T%val < 1d-10 .or. Rho%val < 1d-20 &
228 79994 : .or. m < 1d-10 .or. r%val < 1d-10 .or. cgrav < 1d-10) return
229 :
230 : !test_partials = (k == s% solver_test_partials_k)
231 79994 : test_partials = .false.
232 79994 : ierr = 0
233 79994 : if (k > 0) then
234 79994 : s% tdc_num_iters(k) = 0
235 : end if
236 :
237 : if (report) then
238 : write(*,'(A)')
239 : write(*,4) 'enter Get_results k slvr_itr model gradr grada scale_height ' // trim(MLT_option), &
240 : k, s% solver_iter, s% model_number, gradr%val, grada%val, scale_height%val
241 : end if
242 :
243 :
244 : ! check if this particular k can be done with TDC
245 79994 : using_TDC = .false.
246 79994 : if (s% MLT_option == 'TDC') using_TDC = .true.
247 79994 : if (.not. s% have_mlt_vc) using_TDC = .false.
248 79994 : if (k <= 0 .or. s%dt <= 0d0) using_TDC = .false.
249 140276 : if (using_TDC) using_TDC = .not. check_if_must_fall_back_to_MLT(s, k)
250 :
251 79994 : if (k >= 1) then
252 79994 : s% dvc_dt_TDC(k) = 0d0
253 : end if
254 79994 : if (using_TDC) then
255 : if (report) write(*,3) 'call set_TDC', k, s% solver_iter
256 65210 : if (s% okay_to_set_mlt_vc) then
257 42492 : conv_vel_start = s% mlt_vc_old(k)
258 : else
259 22718 : conv_vel_start = s% mlt_vc(k)
260 : end if
261 :
262 : ! Set scale for judging the TDC luminosity equation Q(Y)=0.
263 : ! Q has units of a luminosity, so the scale should be a luminosity.
264 65210 : if (s% solver_iter == 0) then
265 39866797 : scale = max(abs(s% L(k)), 1d-3*maxval(s% L(1:s% nz)))
266 : else
267 34776621 : scale = max(abs(s% L_start(k)), 1d-3*maxval(s% L_start(1:s% nz)))
268 : end if
269 :
270 : call set_TDC(&
271 : conv_vel_start, mixing_length_alpha, &
272 : s% alpha_TDC_DAMP, s%alpha_TDC_DAMPR, s%alpha_TDC_PtdVdt, s%dt, cgrav, m, report, &
273 : mixing_type, scale, chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, &
274 65210 : scale_height, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), max_conv_vel, ierr)
275 65210 : s% dvc_dt_TDC(k) = (conv_vel%val - conv_vel_start) / s%dt
276 :
277 65210 : if (ierr /= 0) then
278 0 : if (s% report_ierr) write(*,*) 'ierr from set_TDC'
279 0 : return
280 : end if
281 :
282 : ! Experimental method to lower superadiabaticity. Call TDC again with an artificially reduced
283 : ! gradr if the resulting gradT would lead to the radiative luminosity approaching the Eddington
284 : ! limit, or when a density inversion is expected to happen.
285 : ! This is meant as an implicit alternative to okay_to_reduce_gradT_excess
286 65210 : if (s% use_superad_reduction) then
287 0 : call set_superad_reduction
288 0 : if (Gamma_factor > 1d0) then
289 : call set_TDC(&
290 : conv_vel_start, mixing_length_alpha, &
291 : s% alpha_TDC_DAMP, s%alpha_TDC_DAMPR, s%alpha_TDC_PtdVdt, s%dt, cgrav, m, report, &
292 : mixing_type, scale, chiT, chiRho, gradr_scaled, r, P, T, rho, dV, Cp, opacity, &
293 0 : scale_height, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), max_conv_vel, ierr)
294 0 : s% dvc_dt_TDC(k) = (conv_vel%val - conv_vel_start) / s%dt
295 0 : if (ierr /= 0) then
296 0 : if (s% report_ierr) write(*,*) 'ierr from set_TDC when using superad_reduction'
297 0 : return
298 : end if
299 : end if
300 : end if
301 :
302 14784 : else if (gradr > gradL) then
303 : if (report) write(*,3) 'call set_MLT', k, s% solver_iter
304 : call set_MLT(MLT_option, mixing_length_alpha, s% Henyey_MLT_nu_param, s% Henyey_MLT_y_param, &
305 : chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, &
306 : gradr, grada, gradL, &
307 2729 : Gamma, gradT, Y_face, conv_vel, D, mixing_type, max_conv_vel, ierr)
308 :
309 :
310 2729 : if (ierr /= 0) then
311 0 : if (s% report_ierr) write(*,*) 'ierr from set_MLT'
312 0 : return
313 : end if
314 :
315 : ! Experimental method to lower superadiabaticity. Call MLT again with an artificially reduced
316 : ! gradr if the resulting gradT would lead to the radiative luminosity approaching the Eddington
317 : ! limit, or when a density inversion is expected to happen.
318 : ! This is meant as an implicit alternative to okay_to_reduce_gradT_excess
319 2729 : if (s% use_superad_reduction) then
320 0 : call set_superad_reduction
321 0 : if (Gamma_factor > 1d0) then
322 : call set_MLT(MLT_option, mixing_length_alpha, s% Henyey_MLT_nu_param, s% Henyey_MLT_y_param, &
323 : chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, &
324 : gradr_scaled, grada, gradL, &
325 0 : Gamma, gradT, Y_face, conv_vel, D, mixing_type, max_conv_vel, ierr)
326 :
327 0 : if (ierr /= 0) then
328 0 : if (s% report_ierr) write(*,*) 'ierr from set_MLT when using superad_reduction'
329 0 : return
330 : end if
331 : end if
332 : end if
333 : end if
334 :
335 : ! If we're not convecting, try thermohaline and semiconvection.
336 79994 : if (mixing_type == no_mixing) then
337 68285 : if (gradL_composition_term < 0) then
338 : if (report) write(*,3) 'call set_thermohaline', k, s% solver_iter
339 : call set_thermohaline(s%thermohaline_option, Lambda, grada, gradr, T, opacity, rho, Cp, gradL_composition_term, &
340 : iso, XH1, thermohaline_coeff, &
341 0 : D, gradT, Y_face, conv_vel, mixing_type, ierr)
342 0 : if (ierr /= 0) then
343 0 : if (s% report_ierr) write(*,*) 'ierr from set_thermohaline'
344 0 : return
345 : end if
346 68285 : else if (gradr > grada) then
347 : if (report) write(*,3) 'call set_semiconvection', k, s% solver_iter
348 : call set_semiconvection(L, Lambda, m, T, P, Pr, beta, opacity, rho, alpha_semiconvection, &
349 : s% semiconvection_option, cgrav, Cp, gradr, grada, gradL, &
350 : gradL_composition_term, &
351 0 : gradT, Y_face, conv_vel, D, mixing_type, ierr)
352 0 : if (ierr /= 0) then
353 0 : if (s% report_ierr) write(*,*) 'ierr from set_semiconvection'
354 0 : return
355 : end if
356 : end if
357 : end if
358 :
359 : ! If there's too-little mixing to bother, or we hit a bad value, fall back on no mixing.
360 148279 : if (D%val < s% remove_small_D_limit .or. is_bad(D%val)) then
361 : if (report) write(*,2) 'D < s% remove_small_D_limit', k, D%val, s% remove_small_D_limit
362 68285 : mixing_type = no_mixing
363 68285 : gradT = gradr
364 68285 : Y_face = gradT - gradL
365 68285 : conv_vel = 0d0
366 68285 : D = 0d0
367 68285 : Gamma = 0d0
368 : end if
369 :
370 : contains
371 :
372 0 : subroutine set_superad_reduction()
373 0 : Gamma_limit = s% superad_reduction_Gamma_limit
374 0 : scale_value1 = s% superad_reduction_Gamma_limit_scale
375 0 : scale_value2 = s% superad_reduction_Gamma_inv_scale
376 0 : diff_grads_limit = s% superad_reduction_diff_grads_limit
377 0 : reduction_limit = s% superad_reduction_limit
378 0 : Lrad_div_Ledd = 4d0*crad/3d0*pow4(T)/P*gradT
379 0 : Gamma_inv_threshold = 4d0*(1d0-beta)/(4d0-3*beta)
380 :
381 0 : Gamma_factor = 1d0
382 0 : if (gradT > gradL) then
383 0 : if (Lrad_div_Ledd > Gamma_limit .or. Lrad_div_Ledd > Gamma_inv_threshold) then
384 0 : alfa0 = (gradT-gradL)/diff_grads_limit
385 0 : if (alfa0 < 1d0) then
386 0 : diff_grads_factor = -alfa0*alfa0*alfa0*(-10d0 + alfa0*(15d0 - 6d0*alfa0))
387 : else
388 0 : diff_grads_factor = 1d0
389 : end if
390 :
391 0 : Gamma_term = 0d0
392 : !if (Lrad_div_Ledd > Gamma_limit) then
393 : ! Gamma_term = Gamma_term + scale_value1*pow2(Lrad_div_Ledd/Gamma_limit-1d0)
394 : !end if
395 : !if (Lrad_div_Ledd% val > Gamma_inv_threshold) then
396 : ! Gamma_term = Gamma_term + scale_value2*pow2(Lrad_div_Ledd/Gamma_inv_threshold-1d0)
397 : !end if
398 0 : if (Lrad_div_Ledd > Gamma_limit) then
399 0 : alfa0 = Lrad_div_Ledd/Gamma_limit-1d0
400 0 : if (alfa0 < 1d0) then
401 0 : Gamma_term = Gamma_term + scale_value1*(0.5d0*alfa0*alfa0)
402 : else
403 0 : Gamma_term = Gamma_term + scale_value1*(alfa0-0.5d0)
404 : end if
405 : !Gamma_term = Gamma_term + scale_value1*pow2(Lrad_div_Ledd/Gamma_limit-1d0)
406 : end if
407 0 : if (Lrad_div_Ledd% val > Gamma_inv_threshold) then
408 0 : alfa0 = Lrad_div_Ledd/Gamma_inv_threshold-1d0
409 0 : if (alfa0 < 1d0) then
410 0 : Gamma_term = Gamma_term + scale_value1*(0.5d0*alfa0*alfa0)
411 : else
412 0 : Gamma_term = Gamma_term + scale_value1*(alfa0-0.5d0)
413 : end if
414 : !Gamma_term = Gamma_term + scale_value2*pow2(Lrad_div_Ledd/Gamma_inv_threshold-1d0)
415 : end if
416 :
417 0 : if (Gamma_term > 0d0) then
418 0 : Gamma_factor = Gamma_term/pow(beta,0.5d0)*diff_grads_factor
419 0 : Gamma_factor = Gamma_factor + 1d0
420 0 : if (reduction_limit > 1d0) then
421 0 : lambda_limit = 2d0/(reduction_limit-1d0)
422 0 : exp_limit = exp(-lambda_limit*(Gamma_factor-1d0))
423 0 : Gamma_factor = 2d0*(reduction_limit-1d0)*(1d0/(1d0+exp_limit)-0.5d0)+1d0
424 : end if
425 : end if
426 : end if
427 : end if
428 0 : if (k /= 0) s% superad_reduction_factor(k) = Gamma_factor% val
429 0 : if (Gamma_factor > 1d0) then
430 0 : grad_scale = (gradr-gradL)/(Gamma_factor*gradr) + gradL/gradr
431 0 : gradr_scaled = grad_scale*gradr
432 : end if
433 79994 : end subroutine set_superad_reduction
434 : end subroutine Get_results
435 :
436 : end module turb_support
|
