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