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 TDC_alpha_D TDC turbulent damping parameter
87 : !! @param TDC_alpha_R TDC radiative damping parameter
88 : !! @param TDC_alpha_Pt 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 8677186 : subroutine set_TDC( &
115 : conv_vel_start, mixing_length_alpha, TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, 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, &
118 : max_conv_vel, Eq_div_w, grav, include_mlt_corr_to_TDC, TDC_alpha_C, &
119 : TDC_alpha_S, use_TDC_enthalpy_flux_limiter, energy, ierr)
120 0 : use tdc
121 : use tdc_support
122 : real(dp), intent(in) :: conv_vel_start, mixing_length_alpha, TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt
123 : real(dp), intent(in) :: dt, cgrav, m, scale, max_conv_vel, TDC_alpha_C, TDC_alpha_S
124 : type(auto_diff_real_star_order1), intent(in) :: &
125 : chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, scale_height, gradL, grada, Eq_div_w, grav, energy
126 : logical, intent(in) :: report, include_mlt_corr_to_TDC, use_TDC_enthalpy_flux_limiter
127 : type(auto_diff_real_star_order1),intent(out) :: conv_vel, Y_face, gradT, D
128 : integer, intent(out) :: tdc_num_iters, mixing_type, ierr
129 : type(tdc_info) :: info
130 : type(auto_diff_real_star_order1) :: L, Lambda, Gamma, h
131 : real(dp), parameter :: alpha_c = (1d0/2d0)*sqrt_2_div_3
132 : real(dp), parameter :: lower_bound_Z = -1d2
133 : real(dp), parameter :: upper_bound_Z = 1d2
134 : real(dp), parameter :: eps = 1d-2 ! Threshold in logY for separating multiple solutions.
135 : type(auto_diff_real_tdc) :: Zub, Zlb
136 : include 'formats'
137 :
138 : ! Do a call to MLT
139 : !grav = cgrav * m / pow2(r)
140 65242 : L = 64d0 * pi * boltz_sigma * pow4(T) * grav * pow2(r) * gradr / (3d0 * P * opacity)
141 65242 : Lambda = mixing_length_alpha * scale_height
142 : call set_MLT('Cox', mixing_length_alpha, 0d0, 0d0, &
143 : chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, &
144 : gradr, grada, gradL, &
145 65242 : Gamma, gradT, Y_face, conv_vel, D, mixing_type,1d99, ierr)
146 :
147 :
148 : ! Pack TDC info
149 65242 : info%report = report
150 65242 : info%include_mlt_corr_to_TDC = include_mlt_corr_to_TDC
151 65242 : info%use_TDC_enthalpy_flux_limiter = use_TDC_enthalpy_flux_limiter
152 65242 : info%mixing_length_alpha = mixing_length_alpha
153 65242 : info%TDC_alpha_D = TDC_alpha_D
154 65242 : info%TDC_alpha_R = TDC_alpha_R
155 65242 : info%TDC_alpha_Pt = TDC_alpha_Pt
156 65242 : info%dt = dt
157 65242 : info%L = convert(L)
158 65242 : info%gradL = convert(gradL)
159 65242 : info%grada = convert(grada)
160 65242 : info%c0 = convert(TDC_alpha_C * mixing_length_alpha * alpha_c * rho * T * Cp * 4d0 * pi * pow2(r))
161 65242 : info%L0 = convert((16d0*pi*crad*clight/3d0)*cgrav*m*pow4(T)/(P*opacity)) ! assumes QHSE for dP/dm
162 65242 : info%A0 = conv_vel_start/sqrt_2_div_3
163 65242 : info%h = energy + P/rho ! actual enthalpy
164 65242 : info%T = T
165 65242 : info%rho = rho
166 65242 : info%dV = dV
167 65242 : info%Cp = Cp
168 65242 : info%kap = opacity
169 65242 : info%Hp = scale_height
170 65242 : info%Gamma = Gamma
171 65242 : info%Eq_div_w = Eq_div_w
172 65242 : info%TDC_alpha_C = TDC_alpha_C
173 65242 : info%TDC_alpha_S = TDC_alpha_S
174 :
175 : ! Get solution
176 65242 : Zub = upper_bound_Z
177 65242 : Zlb = lower_bound_Z
178 65242 : call get_TDC_solution(info, scale, Zlb, Zub, conv_vel, Y_face, tdc_num_iters, ierr)
179 :
180 : ! Cap conv_vel at max_conv_vel_div_csound*cs
181 65242 : if (conv_vel%val > max_conv_vel) then
182 0 : conv_vel = max_conv_vel
183 : ! if max_conv_vel = csound,
184 : ! L = L0 * (gradL + Y) + c0 * Af * Y_env
185 : ! L = L0 * (gradL + Y) + c0 * sqrt_2_div_3 * csound * (Gamma / (1 + Gamma)) * Y
186 : ! L - L0 * gradL = Y * (L0 + c0 * sqrt_2_div_3 * csound * (Gamma / (1 + Gamma)))
187 0 : if (include_mlt_corr_to_TDC) then
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 * (info%Gamma / (1d0 + info%Gamma)))
190 : else
191 : Y_face = unconvert(info%L - info%L0 * info%gradL) / &
192 0 : (unconvert(info%L0) + unconvert(info%c0) * sqrt_2_div_3 * max_conv_vel)
193 : end if
194 : end if
195 :
196 : ! Unpack output
197 65242 : gradT = Y_face + gradL
198 65242 : D = conv_vel*scale_height*mixing_length_alpha/3d0 ! diffusion coefficient [cm^2/sec]
199 65242 : if (conv_vel > 0d0) then
200 8997 : mixing_type = convective_mixing
201 : else
202 56245 : mixing_type = no_mixing
203 : end if
204 47496176 : end subroutine set_TDC
205 :
206 : !> Calculates the outputs of semiconvective mixing theory.
207 : !!
208 : !! @param L Luminosity across a face (erg/s).
209 : !! @param Lambda The mixing length (cm).
210 : !! @param m Mass inside the face (g).
211 : !! @param T temperature (K).
212 : !! @param P pressure (erg/cm^3).
213 : !! @param Pr radiation pressure (erg/cm^3).
214 : !! @param beta ratio of gas pressure to radiation pressure.
215 : !! @param opacity opacity (cm^2/g).
216 : !! @param rho density (g/cm^3).
217 : !! @param alpha_semiconvection The semiconvective alpha parameter.
218 : !! @param semiconvection_option A string specifying which semiconvection theory to use. Currently supported are 'Langer_85 mixing; gradT = gradr' and 'Langer_85'.
219 : !! @param cgrav gravitational constant (erg*cm/g^2).
220 : !! @param Cp Specific heat at constant pressure (erg/g/K).
221 : !! @param gradr The radiative temperature gradient dlnT/dlnP_{rad}
222 : !! @param grada The adiabatic temperature gradient dlnT/dlnP|s
223 : !! @param gradL The Ledoux temperature gradient dlnT/dlnP
224 : !! @param gradL_composition_term The contribution of composition gradients to the Ledoux temperature gradient.
225 : !! @param gradT The temperature gradient dlnT/dlnP (output).
226 : !! @param Y_face The superadiabaticity (dlnT/dlnP - grada, output).
227 : !! @param conv_vel The convection speed (cm/s).
228 : !! @param D The chemical diffusion coefficient (cm^2/s).
229 : !! @param mixing_type Set to semiconvective if convection operates (output).
230 : !! @param ierr Tracks errors (output).
231 0 : subroutine set_semiconvection(L, Lambda, m, T, P, Pr, beta, opacity, rho, alpha_semiconvection, &
232 : semiconvection_option, cgrav, Cp, gradr, grada, gradL, &
233 : gradL_composition_term, &
234 : gradT, Y_face, conv_vel, D, mixing_type, ierr)
235 65242 : use semiconvection
236 : type(auto_diff_real_star_order1), intent(in) :: L, Lambda, T, P, Pr, beta, opacity, rho
237 : type(auto_diff_real_star_order1), intent(in) :: Cp, gradr, grada, gradL
238 : character(len=*), intent(in) :: semiconvection_option
239 : real(dp), intent(in) :: alpha_semiconvection, cgrav, gradL_composition_term, m
240 : type(auto_diff_real_star_order1), intent(out) :: gradT, Y_face, conv_vel, D
241 : integer, intent(out) :: mixing_type, ierr
242 :
243 : call calc_semiconvection(L, Lambda, m, T, P, Pr, beta, opacity, rho, alpha_semiconvection, &
244 : semiconvection_option, cgrav, Cp, gradr, grada, gradL, &
245 : gradL_composition_term, &
246 0 : gradT, Y_face, conv_vel, D, mixing_type, ierr)
247 0 : end subroutine set_semiconvection
248 :
249 : !> Calculates the outputs of convective mixing length theory.
250 : !!
251 : !! @param MLT_option A string specifying which MLT option to use. Currently supported are Cox, Henyey, ML1, ML2, Mihalas.
252 : !! Note that 'TDC' is also a valid input and will return the Cox result.
253 : !! This is for use when falling back from TDC -> MLT, as Cox is the most-similar prescription to TDC.
254 : !! @param mixing_length_alpha The mixing length parameter.
255 : !! @param Henyey_MLT_nu_param The nu parameter in Henyey's MLT prescription.
256 : !! @param Henyey_MLT_y_param The y parameter in Henyey's MLT prescription.
257 : !! @param chiT dlnP/dlnT|rho
258 : !! @param chiRho dlnP/dlnRho|T
259 : !! @param Cp Specific heat at constant pressure (erg/g/K).
260 : !! @param grav The acceleration due to gravity (cm/s^2).
261 : !! @param Lambda The mixing length (cm).
262 : !! @param rho density (g/cm^3).
263 : !! @param T temperature (K).
264 : !! @param opacity opacity (cm^2/g)
265 : !! @param gradr The radiative temperature gradient dlnT/dlnP_{rad}
266 : !! @param grada The adiabatic temperature gradient dlnT/dlnP|s
267 : !! @param gradL The Ledoux temperature gradient dlnT/dlnP
268 : !! @param Gamma The convective efficiency parameter (output).
269 : !! @param gradT The temperature gradient dlnT/dlnP (output).
270 : !! @param Y_face The superadiabaticity (dlnT/dlnP - grada, output).
271 : !! @param conv_vel The convection speed (cm/s).
272 : !! @param D The chemical diffusion coefficient (cm^2/s).
273 : !! @param mixing_type Set to convective if convection operates (output).
274 : !! @param ierr Tracks errors (output).
275 11283684 : subroutine set_MLT(MLT_option, mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param, &
276 : chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, &
277 : gradr, grada, gradL, &
278 : Gamma, gradT, Y_face, conv_vel, D, mixing_type, max_conv_vel, ierr)
279 : use mlt
280 : type(auto_diff_real_star_order1), intent(in) :: chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, gradr, grada, gradL
281 : character(len=*), intent(in) :: MLT_option
282 : real(dp), intent(in) :: mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param, max_conv_vel
283 :
284 : type(auto_diff_real_star_order1), intent(out) :: Gamma, gradT, Y_face, conv_vel, D
285 : integer, intent(out) :: mixing_type, ierr
286 :
287 : call calc_MLT(MLT_option, mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param, &
288 : chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, &
289 : gradr, grada, gradL, &
290 67974 : Gamma, gradT, Y_face, conv_vel, D, mixing_type, max_conv_vel, ierr)
291 :
292 67974 : end subroutine set_MLT
293 :
294 : end module turb
|
