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 turb
21 : use const_def, only: dp, pi, boltz_sigma, sqrt_2_div_3, crad, clight, no_mixing, convective_mixing, thermohaline_mixing
22 : use num_lib
23 : use utils_lib
24 : use auto_diff
25 :
26 : implicit none
27 :
28 : private
29 : public :: set_thermohaline
30 : public :: set_mlt
31 : public :: set_tdc
32 : public :: set_semiconvection
33 :
34 : contains
35 :
36 : !> Computes the diffusivity of thermohaline mixing when the
37 : !! thermal gradient is stable and the composition gradient is unstable.
38 : !!
39 : !! @param thermohaline_option A string specifying which thermohaline prescription to use.
40 : !! @param grada Adiabatic gradient dlnT/dlnP
41 : !! @param gradr Radiative temperature gradient dlnT/dlnP, equals the actual gradient because there's no convection
42 : !! @param T Temperature
43 : !! @param opacity opacity
44 : !! @param rho Density
45 : !! @param Cp Heat capacity at constant pressure
46 : !! @param gradL_composition_term dlnMu/dlnP where Mu is the mean molecular weight.
47 : !! @param iso The index of the species that drives thermohaline mixing.
48 : !! @param XH1 Mass fraction of H1.
49 : !! @param thermohaline_coeff Free parameter multiplying the thermohaline diffusivity.
50 : !! @param D_thrm Output, diffusivity.
51 : !! @param ierr Output, error index.
52 0 : subroutine set_thermohaline(thermohaline_option, Lambda, grada, gradr, T, opacity, rho, Cp, gradL_composition_term, &
53 : iso, XH1, thermohaline_coeff, &
54 : D, gradT, Y_face, conv_vel, mixing_type, ierr)
55 : use thermohaline
56 : character(len=*), intent(in) :: thermohaline_option
57 : type(auto_diff_real_star_order1), intent(in) :: Lambda, grada, gradr, T, opacity, rho, Cp
58 : real(dp), intent(in) :: gradL_composition_term, XH1, thermohaline_coeff
59 : integer, intent(in) :: iso
60 :
61 : type(auto_diff_real_star_order1), intent(out) :: gradT, Y_face, conv_vel, D
62 : integer, intent(out) :: mixing_type, ierr
63 :
64 : real(dp) :: D_thrm
65 :
66 : call get_D_thermohaline(&
67 : thermohaline_option, grada%val, gradr%val, T%val, opacity%val, rho%val, &
68 : Cp%val, gradL_composition_term, &
69 0 : iso, XH1, thermohaline_coeff, D_thrm, ierr)
70 :
71 0 : D = D_thrm
72 0 : gradT = gradr
73 0 : Y_face = gradT - grada
74 0 : conv_vel = 3d0*D/Lambda
75 0 : mixing_type = thermohaline_mixing
76 0 : end subroutine set_thermohaline
77 :
78 : !> Computes the outputs of time-dependent convection theory following the model specified in
79 : !! Radek Smolec's thesis [https://users.camk.edu.pl/smolec/phd_smolec.pdf], which in turn
80 : !! follows the model of Kuhfuss 1986.
81 : !!
82 : !! Internally this solves the equation L = L_conv + L_rad.
83 : !!
84 : !! @param conv_vel_start The convection speed at the start of the step.
85 : !! @param mixing_length_alpha The mixing length parameter.
86 : !! @param alpha_TDC_DAMP TDC turbulent damping parameter
87 : !! @param alpha_TDC_DAMPR TDC radiative damping parameter
88 : !! @param alpha_TDC_PtdVdt TDC coefficient on P_turb*dV/dt. Physically should probably be 1.
89 : !! @param The time-step (s).
90 : !! @param cgrav gravitational constant (erg*cm/g^2).
91 : !! @param m Mass inside the face (g).
92 : !! @param report Write debug output if true, not if false.
93 : !! @param mixing_type Set to semiconvective if convection operates (output).
94 : !! @param scale The scale for computing residuals to the luminosity equation (erg/s).
95 : !! @param chiT dlnP/dlnT|rho
96 : !! @param chiRho dlnP/dlnRho|T
97 : !! @param gradr Radiative temperature gradient.
98 : !! @param r radial coordinate of the face (cm).
99 : !! @param P pressure (erg/cm^3).
100 : !! @param T temperature (K).
101 : !! @param rho density (g/cm^3).
102 : !! @param dV The change in specific volume of the face (cm^3/g) since the start of the step.
103 : !! @param Cp Specific heat at constant pressure (erg/g/K).
104 : !! @param opacity opacity (cm^2/g).
105 : !! @param scale_height The pressure scale-height (cm).
106 : !! @param gradL The Ledoux temperature gradient dlnT/dlnP
107 : !! @param grada The adiabatic temperature gradient dlnT/dlnP|s
108 : !! @param conv_vel The convection speed (cm/s).
109 : !! @param D The chemical diffusion coefficient (cm^2/s).
110 : !! @param Y_face The superadiabaticity (dlnT/dlnP - grada, output).
111 : !! @param gradT The temperature gradient dlnT/dlnP (output).
112 : !! @param tdc_num_iters Number of iterations taken in the TDC solver.
113 : !! @param ierr Tracks errors (output).
114 65242 : subroutine set_TDC( &
115 : conv_vel_start, mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, report, &
116 : mixing_type, scale, chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, &
117 : scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, max_conv_vel, &
118 : Eq_div_w, grav, include_mlt_corr_to_TDC, alpha_TDC_C, alpha_TDC_S, ierr)
119 0 : use tdc
120 : use tdc_support
121 : real(dp), intent(in) :: conv_vel_start, mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt
122 : real(dp), intent(in) :: dt, cgrav, m, scale, max_conv_vel, alpha_TDC_c, alpha_TDC_s
123 : type(auto_diff_real_star_order1), intent(in) :: &
124 : chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, scale_height, gradL, grada, Eq_div_w, grav
125 : logical, intent(in) :: report, include_mlt_corr_to_TDC
126 : type(auto_diff_real_star_order1),intent(out) :: conv_vel, Y_face, gradT, D
127 : integer, intent(out) :: tdc_num_iters, mixing_type, ierr
128 : type(tdc_info) :: info
129 : type(auto_diff_real_star_order1) :: L, Lambda, Gamma
130 : real(dp), parameter :: alpha_c = (1d0/2d0)*sqrt_2_div_3
131 : real(dp), parameter :: lower_bound_Z = -1d2
132 : real(dp), parameter :: upper_bound_Z = 1d2
133 : real(dp), parameter :: eps = 1d-2 ! Threshold in logY for separating multiple solutions.
134 : type(auto_diff_real_tdc) :: Zub, Zlb
135 : include 'formats'
136 :
137 : ! Do a call to MLT
138 : !grav = cgrav * m / pow2(r)
139 65242 : L = 64d0 * pi * boltz_sigma * pow4(T) * grav * pow2(r) * gradr / (3d0 * P * opacity)
140 65242 : Lambda = mixing_length_alpha * scale_height
141 : call set_MLT('Cox', mixing_length_alpha, 0d0, 0d0, &
142 : chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, &
143 : gradr, grada, gradL, &
144 130484 : Gamma, gradT, Y_face, conv_vel, D, mixing_type,1d99, ierr)
145 :
146 :
147 : ! Pack TDC info
148 65242 : info%report = report
149 65242 : info%include_mlt_corr_to_TDC = include_mlt_corr_to_TDC
150 65242 : info%mixing_length_alpha = mixing_length_alpha
151 65242 : info%alpha_TDC_DAMP = alpha_TDC_DAMP
152 65242 : info%alpha_TDC_DAMPR = alpha_TDC_DAMPR
153 65242 : info%alpha_TDC_PtdVdt = alpha_TDC_PtdVdt
154 65242 : info%dt = dt
155 65242 : info%L = convert(L)
156 65242 : info%gradL = convert(gradL)
157 65242 : info%grada = convert(grada)
158 65242 : info%c0 = convert(alpha_TDC_C * mixing_length_alpha * alpha_c * rho * T * Cp * 4d0 * pi * pow2(r))
159 65242 : info%L0 = convert((16d0*pi*crad*clight/3d0)*cgrav*m*pow4(T)/(P*opacity)) ! assumes QHSE for dP/dm, needs correction for if s% make_mlt_hydrodynamic = .false.
160 65242 : info%A0 = conv_vel_start/sqrt_2_div_3
161 65242 : info%T = T
162 65242 : info%rho = rho
163 65242 : info%dV = dV
164 65242 : info%Cp = Cp
165 65242 : info%kap = opacity
166 65242 : info%Hp = scale_height
167 65242 : info%Gamma = Gamma
168 65242 : info%Eq_div_w = Eq_div_w
169 65242 : info%alpha_TDC_c = alpha_TDC_C
170 65242 : info%alpha_TDC_s = alpha_TDC_S
171 :
172 : ! Get solution
173 65242 : Zub = upper_bound_Z
174 65242 : Zlb = lower_bound_Z
175 65242 : call get_TDC_solution(info, scale, Zlb, Zub, conv_vel, Y_face, tdc_num_iters, ierr)
176 :
177 : ! Cap conv_vel at max_conv_vel_div_csound*cs
178 65242 : if (conv_vel%val > max_conv_vel) then
179 0 : conv_vel = max_conv_vel
180 : ! if max_conv_vel = csound,
181 : ! L = L0 * (gradL + Y) + c0 * Af * Y_env
182 : ! L = L0 * (gradL + Y) + c0 * sqrt_2_div_3 * csound * (Gamma / (1 + Gamma)) * Y
183 : ! L - L0 * gradL = Y * (L0 + c0 * sqrt_2_div_3 * csound * (Gamma / (1 + Gamma)))
184 0 : if (include_mlt_corr_to_TDC) then
185 : Y_face = unconvert(info%L - info%L0 * info%gradL) / &
186 0 : (unconvert(info%L0) + unconvert(info%c0) * sqrt_2_div_3 * max_conv_vel * (info%Gamma / (1d0 + info%Gamma)))
187 : else
188 : Y_face = unconvert(info%L - info%L0 * info%gradL) / &
189 0 : (unconvert(info%L0) + unconvert(info%c0) * sqrt_2_div_3 * max_conv_vel)
190 : end if
191 : end if
192 :
193 : ! Unpack output
194 65242 : gradT = Y_face + gradL
195 65242 : D = conv_vel*scale_height*mixing_length_alpha/3d0 ! diffusion coefficient [cm^2/sec]
196 65242 : if (conv_vel > 0d0) then
197 8997 : mixing_type = convective_mixing
198 : else
199 56245 : mixing_type = no_mixing
200 : end if
201 65242 : end subroutine set_TDC
202 :
203 : !> Calculates the outputs of semiconvective mixing theory.
204 : !!
205 : !! @param L Luminosity across a face (erg/s).
206 : !! @param Lambda The mixing length (cm).
207 : !! @param m Mass inside the face (g).
208 : !! @param T temperature (K).
209 : !! @param P pressure (erg/cm^3).
210 : !! @param Pr radiation pressure (erg/cm^3).
211 : !! @param beta ratio of gas pressure to radiation pressure.
212 : !! @param opacity opacity (cm^2/g).
213 : !! @param rho density (g/cm^3).
214 : !! @param alpha_semiconvection The semiconvective alpha parameter.
215 : !! @param semiconvection_option A string specifying which semiconvection theory to use. Currently supported are 'Langer_85 mixing; gradT = gradr' and 'Langer_85'.
216 : !! @param cgrav gravitational constant (erg*cm/g^2).
217 : !! @param Cp Specific heat at constant pressure (erg/g/K).
218 : !! @param gradr The radiative temperature gradient dlnT/dlnP_{rad}
219 : !! @param grada The adiabatic temperature gradient dlnT/dlnP|s
220 : !! @param gradL The Ledoux temperature gradient dlnT/dlnP
221 : !! @param gradL_composition_term The contribution of composition gradients to the Ledoux temperature gradient.
222 : !! @param gradT The temperature gradient dlnT/dlnP (output).
223 : !! @param Y_face The superadiabaticity (dlnT/dlnP - grada, output).
224 : !! @param conv_vel The convection speed (cm/s).
225 : !! @param D The chemical diffusion coefficient (cm^2/s).
226 : !! @param mixing_type Set to semiconvective if convection operates (output).
227 : !! @param ierr Tracks errors (output).
228 0 : subroutine set_semiconvection(L, Lambda, m, T, P, Pr, beta, opacity, rho, alpha_semiconvection, &
229 : semiconvection_option, cgrav, Cp, gradr, grada, gradL, &
230 : gradL_composition_term, &
231 : gradT, Y_face, conv_vel, D, mixing_type, ierr)
232 65242 : use semiconvection
233 : type(auto_diff_real_star_order1), intent(in) :: L, Lambda, T, P, Pr, beta, opacity, rho
234 : type(auto_diff_real_star_order1), intent(in) :: Cp, gradr, grada, gradL
235 : character(len=*), intent(in) :: semiconvection_option
236 : real(dp), intent(in) :: alpha_semiconvection, cgrav, gradL_composition_term, m
237 : type(auto_diff_real_star_order1), intent(out) :: gradT, Y_face, conv_vel, D
238 : integer, intent(out) :: mixing_type, ierr
239 :
240 : call calc_semiconvection(L, Lambda, m, T, P, Pr, beta, opacity, rho, alpha_semiconvection, &
241 : semiconvection_option, cgrav, Cp, gradr, grada, gradL, &
242 : gradL_composition_term, &
243 0 : gradT, Y_face, conv_vel, D, mixing_type, ierr)
244 0 : end subroutine set_semiconvection
245 :
246 : !> Calculates the outputs of convective mixing length theory.
247 : !!
248 : !! @param MLT_option A string specifying which MLT option to use. Currently supported are Cox, Henyey, ML1, ML2, Mihalas.
249 : !! Note that 'TDC' is also a valid input and will return the Cox result.
250 : !! This is for use when falling back from TDC -> MLT, as Cox is the most-similar prescription to TDC.
251 : !! @param mixing_length_alpha The mixing length parameter.
252 : !! @param Henyey_MLT_nu_param The nu parameter in Henyey's MLT prescription.
253 : !! @param Henyey_MLT_y_param The y parameter in Henyey's MLT prescription.
254 : !! @param chiT dlnP/dlnT|rho
255 : !! @param chiRho dlnP/dlnRho|T
256 : !! @param Cp Specific heat at constant pressure (erg/g/K).
257 : !! @param grav The acceleration due to gravity (cm/s^2).
258 : !! @param Lambda The mixing length (cm).
259 : !! @param rho density (g/cm^3).
260 : !! @param T temperature (K).
261 : !! @param opacity opacity (cm^2/g)
262 : !! @param gradr The radiative temperature gradient dlnT/dlnP_{rad}
263 : !! @param grada The adiabatic temperature gradient dlnT/dlnP|s
264 : !! @param gradL The Ledoux temperature gradient dlnT/dlnP
265 : !! @param Gamma The convective efficiency parameter (output).
266 : !! @param gradT The temperature gradient dlnT/dlnP (output).
267 : !! @param Y_face The superadiabaticity (dlnT/dlnP - grada, output).
268 : !! @param conv_vel The convection speed (cm/s).
269 : !! @param D The chemical diffusion coefficient (cm^2/s).
270 : !! @param mixing_type Set to convective if convection operates (output).
271 : !! @param ierr Tracks errors (output).
272 67974 : subroutine set_MLT(MLT_option, mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param, &
273 : chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, &
274 : gradr, grada, gradL, &
275 : Gamma, gradT, Y_face, conv_vel, D, mixing_type, max_conv_vel, ierr)
276 : use mlt
277 : type(auto_diff_real_star_order1), intent(in) :: chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, gradr, grada, gradL
278 : character(len=*), intent(in) :: MLT_option
279 : real(dp), intent(in) :: mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param, max_conv_vel
280 :
281 : type(auto_diff_real_star_order1), intent(out) :: Gamma, gradT, Y_face, conv_vel, D
282 : integer, intent(out) :: mixing_type, ierr
283 :
284 : call calc_MLT(MLT_option, mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param, &
285 : chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, &
286 : gradr, grada, gradL, &
287 67974 : Gamma, gradT, Y_face, conv_vel, D, mixing_type, max_conv_vel, ierr)
288 :
289 2732 : end subroutine set_MLT
290 :
291 : end module turb
|
