Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2018-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 hydro_temperature
21 :
22 : use star_private_def
23 : use const_def, only: dp, ln10, pi4, crad, clight, convective_mixing
24 : use utils_lib, only: mesa_error, is_bad
25 : use auto_diff
26 : use auto_diff_support
27 :
28 : implicit none
29 :
30 : private
31 : public :: do1_alt_dlnT_dm_eqn
32 : public :: do1_gradT_eqn
33 : public :: do1_dlnT_dm_eqn
34 :
35 : contains
36 :
37 : ! just relate L_rad to T gradient.
38 : ! d_P_rad/dm = -<opacity_face>*L_rad/(clight*area^2) -- see, e.g., K&W (5.12)
39 : ! P_rad = (1/3)*crad*T^4
40 : ! d_P_rad/dm = (crad/3)*(T(k-1)^4 - T(k)^4)/dm_bar
41 : ! L_rad = L - L_non_rad, L_non_rad = L_start - L_rad_start
42 : ! L_rad_start = (-d_P_rad/dm_bar*clight*area^2/<opacity_face>)_start
43 52304 : subroutine do1_alt_dlnT_dm_eqn(s, k, nvar, ierr)
44 : use eos_def
45 : use star_utils, only: save_eqn_residual_info, get_face_weights
46 : type (star_info), pointer :: s
47 : integer, intent(in) :: k, nvar
48 : integer, intent(out) :: ierr
49 :
50 0 : real(dp) :: alfa, beta, scale, dm_bar
51 : type(auto_diff_real_star_order1) :: L_ad, r_00, area, area2, Lrad_ad, &
52 : kap_00, kap_m1, kap_face, d_P_rad_expected_ad, T_m1, T4_m1, T_00, T4_00, &
53 : P_rad_m1, P_rad_00, d_P_rad_actual_ad, resid
54 : type(auto_diff_real_star_order1) :: flxR, flxLambda
55 :
56 : integer :: i_equL
57 : logical :: dbg
58 : logical :: test_partials
59 :
60 : include 'formats'
61 0 : ierr = 0
62 0 : i_equL = s% i_equL
63 0 : if (i_equL == 0) return
64 :
65 0 : if (.not. s% use_dPrad_dm_form_of_T_gradient_eqn) then
66 0 : ierr = -1
67 0 : return
68 : end if
69 :
70 : !test_partials = (k == s% solver_test_partials_k)
71 0 : test_partials = .false.
72 :
73 0 : dbg = .false.
74 :
75 0 : call get_face_weights(s, k, alfa, beta)
76 :
77 0 : scale = s% energy_start(k)*s% rho_start(k)
78 0 : dm_bar = s% dm_bar(k)
79 0 : L_ad = wrap_L_00(s,k)
80 0 : r_00 = wrap_r_00(s,k)
81 0 : area = pi4*pow2(r_00); area2 = pow2(area)
82 :
83 : if (s% lnT(k)/ln10 <= s% max_logT_for_mlt &
84 : .and. s% mixing_type(k) == convective_mixing .and. s% gradr(k) > 0d0 &
85 0 : .and. abs(s% gradr(k) - s% gradT(k)) > abs(s% gradr(k))*1d-5) then
86 0 : Lrad_ad = L_ad*s% gradT_ad(k)/s% gradr_ad(k) ! C&G 14.109
87 : else
88 0 : Lrad_ad = L_ad
89 : end if
90 :
91 0 : kap_00 = wrap_kap_00(s,k)
92 0 : kap_m1 = wrap_kap_m1(s,k)
93 0 : kap_face = alfa*kap_00 + beta*kap_m1
94 0 : if (kap_face%val < s% min_kap_for_dPrad_dm_eqn) &
95 0 : kap_face = s% min_kap_for_dPrad_dm_eqn
96 :
97 : ! calculate expected d_P_rad from current L_rad
98 0 : d_P_rad_expected_ad = -dm_bar*kap_face*Lrad_ad/(clight*area2)
99 :
100 : ! calculate actual d_P_rad in current model
101 0 : T_m1 = wrap_T_m1(s,k); T4_m1 = pow4(T_m1)
102 0 : T_00 = wrap_T_00(s,k); T4_00 = pow4(T_00)
103 :
104 : !d_P_rad_expected = d_P_rad_expected*s% gradr_factor(k) !TODO(Pablo): check this
105 :
106 0 : P_rad_m1 = (crad/3._dp)*T4_m1
107 0 : P_rad_00 = (crad/3._dp)*T4_00
108 0 : d_P_rad_actual_ad = P_rad_m1 - P_rad_00
109 :
110 : ! enable flux-limited radiation transport derived by Levermore & Pomraning 1981
111 0 : s% flux_limit_R(k) = 0._dp
112 0 : s% flux_limit_lambda(k) =0._dp
113 0 : if (s% use_flux_limiting_with_dPrad_dm_form) then
114 : ! calculate the flux ratio R
115 : flxR = area * abs(T4_m1 - T4_00) / dm_bar / &
116 0 : (kap_face * 0.5_dp * (T4_m1 + T4_00))
117 :
118 0 : s% flux_limit_R(k) = flxR%val
119 :
120 : ! calculate the flux limiter lambda
121 0 : flxLambda = (6._dp + 3._dp*flxR) / (6._dp + (3._dp + flxR)*flxR)
122 :
123 0 : s% flux_limit_lambda(k) = flxLambda%val
124 :
125 : ! calculate d_P_rad given the flux limiter
126 0 : d_P_rad_expected_ad = d_P_rad_expected_ad / flxLambda
127 : end if
128 :
129 : ! residual
130 0 : resid = (d_P_rad_expected_ad - d_P_rad_actual_ad)/scale
131 0 : s% equ(i_equL, k) = resid%val
132 :
133 0 : if (is_bad(resid%val)) then
134 0 : !$OMP critical (star_alt_dlntdm_bad_num)
135 0 : write(*,2) 'resid%val', k, resid%val
136 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'do1_alt_dlnT_dm_eqn')
137 : !$OMP end critical (star_alt_dlntdm_bad_num)
138 : end if
139 :
140 : if (test_partials) then
141 : s% solver_test_partials_val = s% gradT(k)
142 : end if
143 :
144 : call save_eqn_residual_info( &
145 0 : s, k, nvar, i_equL, resid, 'do1_alt_dlnT_dm_eqn', ierr)
146 :
147 0 : if (test_partials) then
148 : s% solver_test_partials_var = 0
149 : s% solver_test_partials_dval_dx = 0
150 : write(*,*) 'do1_alt_dlnT_dm_eqn', s% solver_test_partials_var
151 : end if
152 :
153 : contains
154 :
155 : end subroutine do1_alt_dlnT_dm_eqn
156 :
157 :
158 0 : subroutine do1_gradT_eqn(s, k, nvar, ierr)
159 0 : use eos_def
160 : use star_utils, only: save_eqn_residual_info
161 : type (star_info), pointer :: s
162 : integer, intent(in) :: k, nvar
163 : integer, intent(out) :: ierr
164 :
165 : type(auto_diff_real_star_order1) :: &
166 : resid, gradT, dlnT, dlnP
167 : integer :: i_equL
168 : logical :: test_partials
169 :
170 : include 'formats'
171 0 : ierr = 0
172 :
173 : !test_partials = (k == s% solver_test_partials_k)
174 0 : test_partials = .false.
175 :
176 0 : i_equL = s% i_equL
177 0 : if (i_equL == 0) return
178 :
179 0 : gradT = s% gradT_ad(k)
180 0 : dlnT = wrap_lnT_m1(s,k) - wrap_lnT_00(s,k)
181 0 : dlnP = wrap_lnPeos_m1(s,k) - wrap_lnPeos_00(s,k)
182 :
183 0 : resid = gradT*dlnP - dlnT
184 0 : s% equ(i_equL, k) = resid%val
185 :
186 0 : if (is_bad(s% equ(i_equL, k))) then
187 0 : ierr = -1
188 0 : if (s% report_ierr) write(*,2) 'equ(i_equL, k)', k, s% equ(i_equL, k)
189 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'do1_gradT_eqn')
190 0 : return
191 : write(*,2) 'equ(i_equL, k)', k, s% equ(i_equL, k)
192 : write(*,2) 'gradT', k, gradT
193 : write(*,2) 'dlnT', k, dlnT
194 : write(*,2) 'dlnP', k, dlnP
195 : call mesa_error(__FILE__,__LINE__,'do1_gradT_eqn')
196 : end if
197 :
198 : if (test_partials) then
199 : s% solver_test_partials_val = s% equ(i_equL,k)
200 : end if
201 :
202 : call save_eqn_residual_info( &
203 0 : s, k, nvar, i_equL, resid, 'do1_gradT_eqn', ierr)
204 :
205 : !call set_xtras
206 :
207 : contains
208 :
209 : subroutine set_xtras
210 0 : use auto_diff_support
211 : use star_utils, only: get_Lrad
212 : type(auto_diff_real_star_order1) :: &
213 : T4m1, T400, kap_m1, kap_00, alfa, beta, kap_face, &
214 : diff_T4_div_kap
215 : T4m1 = pow4(wrap_T_m1(s,k))
216 : T400 = pow4(wrap_T_00(s,k))
217 : kap_m1 = wrap_kap_m1(s,k)
218 : kap_00 = wrap_kap_00(s,k)
219 : alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k))
220 : beta = 1d0 - alfa
221 : kap_face = alfa*kap_00 + beta*kap_m1
222 : diff_T4_div_kap = (T4m1 - T400)/kap_face
223 : s% xtra1_array(k) = s% T_start(k)
224 : s% xtra2_array(k) = T4m1%val - T400%val
225 : s% xtra3_array(k) = kap_face%val
226 : s% xtra4_array(k) = diff_T4_div_kap%val
227 : s% xtra5_array(k) = get_Lrad(s,k)
228 : s% xtra6_array(k) = 1
229 : end subroutine set_xtras
230 :
231 : end subroutine do1_gradT_eqn
232 :
233 :
234 52304 : subroutine do1_dlnT_dm_eqn(s, k, nvar, ierr)
235 : use eos_def
236 : use star_utils, only: save_eqn_residual_info
237 : type (star_info), pointer :: s
238 : integer, intent(in) :: k, nvar
239 : integer, intent(out) :: ierr
240 :
241 : type(auto_diff_real_star_order1) :: resid, &
242 : dlnPdm, Ppoint, gradT, dlnTdm, T00, Tm1, dT, Tpoint, lnTdiff
243 52304 : real(dp) :: delm, alfa
244 : integer :: i_equL
245 : logical :: test_partials
246 :
247 : include 'formats'
248 52304 : ierr = 0
249 :
250 : !test_partials = (k == s% solver_test_partials_k)
251 52304 : test_partials = .false.
252 :
253 52304 : i_equL = s% i_equL
254 52304 : if (i_equL == 0) return
255 :
256 52304 : if (s% use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn) then
257 0 : call do1_gradT_eqn(s, k, nvar, ierr)
258 0 : return
259 : end if
260 :
261 52304 : if (s% use_dPrad_dm_form_of_T_gradient_eqn) then
262 0 : call do1_alt_dlnT_dm_eqn(s, k, nvar, ierr)
263 0 : return
264 : end if
265 :
266 : ! dT/dm = dP/dm * T/P * grad_T, grad_T = dlnT/dlnP from MLT.
267 : ! but use hydrostatic value for dP/dm in this.
268 : ! this is because of limitations of MLT for calculating grad_T.
269 : ! (MLT assumes hydrostatic equilibrium)
270 : ! see comment in K&W chpt 9.1.
271 :
272 52304 : call eval_dlnPdm_qhse(s, k, dlnPdm, Ppoint, ierr)
273 52304 : if (ierr /= 0) return
274 :
275 52304 : gradT = s% gradT_ad(k)
276 52304 : dlnTdm = dlnPdm*gradT
277 :
278 52304 : Tm1 = wrap_T_m1(s,k)
279 52304 : T00 = wrap_T_00(s,k)
280 52304 : dT = Tm1 - T00
281 52304 : alfa = s% dm(k-1)/(s% dm(k-1) + s% dm(k))
282 52304 : Tpoint = alfa*T00 + (1d0 - alfa)*Tm1
283 52304 : lnTdiff = dT/Tpoint ! use this in place of lnT(k-1)-lnT(k)
284 52304 : delm = (s% dm(k) + s% dm(k-1))/2
285 :
286 52304 : resid = delm*dlnTdm - lnTdiff
287 52304 : s% equ(i_equL, k) = resid%val
288 :
289 52304 : if (is_bad(s% equ(i_equL, k))) then
290 0 : ierr = -1
291 0 : if (s% report_ierr) write(*,2) 'equ(i_equL, k)', k, s% equ(i_equL, k)
292 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'hydro eqns')
293 0 : return
294 : write(*,2) 'equ(i_equL, k)', k, s% equ(i_equL, k)
295 : write(*,2) 'lnTdiff', k, lnTdiff
296 : write(*,2) 'delm', k, delm
297 : write(*,2) 'dlnPdm', k, dlnPdm
298 : write(*,2) 'gradT', k, gradT
299 : call mesa_error(__FILE__,__LINE__,'i_equL')
300 : end if
301 :
302 : if (test_partials) then
303 : s% solver_test_partials_val = s% equ(i_equL,k)
304 : end if
305 :
306 : call save_eqn_residual_info( &
307 52304 : s, k, nvar, i_equL, resid, 'do1_dlnT_dm_eqn', ierr)
308 :
309 52304 : end subroutine do1_dlnT_dm_eqn
310 :
311 :
312 :
313 : ! only used for dlnT_dm equation
314 52304 : subroutine eval_dlnPdm_qhse(s, k, & ! calculate the expected dlnPdm for HSE
315 : dlnPdm_qhse, Ppoint, ierr)
316 52304 : use hydro_momentum, only: expected_HSE_grav_term
317 : type (star_info), pointer :: s
318 : integer, intent(in) :: k
319 : type(auto_diff_real_star_order1), intent(out) :: dlnPdm_qhse, Ppoint
320 : integer, intent(out) :: ierr
321 :
322 52304 : real(dp) :: alfa
323 : type(auto_diff_real_star_order1) :: grav, area, P00, Pm1
324 : include 'formats'
325 :
326 : ierr = 0
327 :
328 : ! basic eqn is dP/dm = -G m / (4 pi r^4)
329 : ! divide by Ppoint to make it unitless
330 :
331 : ! for rotation, multiply gravity by factor fp. MESA 2, eqn 22.
332 :
333 52304 : call expected_HSE_grav_term(s, k, grav, area, ierr)
334 52304 : if (ierr /= 0) return
335 :
336 52304 : P00 = wrap_Peos_00(s,k)
337 52304 : if (s% using_velocity_time_centering) P00 = 0.5d0*(P00 + s% Peos_start(k))
338 :
339 52304 : if (k == 1) then
340 0 : Pm1 = 0d0
341 0 : Ppoint = P00
342 : else
343 52304 : Pm1 = wrap_Peos_m1(s,k)
344 52304 : alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k))
345 52304 : Ppoint = alfa*P00 + (1d0-alfa)*Pm1
346 : end if
347 :
348 52304 : dlnPdm_qhse = grav/(area*Ppoint) ! note that expected_HSE_grav_term is negative
349 :
350 52304 : if (is_bad(dlnPdm_qhse%val)) then
351 0 : ierr = -1
352 0 : s% retry_message = 'eval_dlnPdm_qhse: is_bad(dlnPdm_qhse)'
353 0 : if (s% report_ierr) then
354 0 : !$OMP critical (hydro_vars_crit1)
355 0 : write(*,*) 'eval_dlnPdm_qhse: is_bad(dlnPdm_qhse)'
356 0 : stop
357 : !$OMP end critical (hydro_vars_crit1)
358 : end if
359 0 : if (s% stop_for_bad_nums) then
360 0 : write(*,2) 'dlnPdm_qhse', k, dlnPdm_qhse
361 0 : call mesa_error(__FILE__,__LINE__,'eval_dlnPdm_qhse')
362 : end if
363 0 : return
364 : end if
365 :
366 52304 : end subroutine eval_dlnPdm_qhse
367 :
368 : end module hydro_temperature
|
