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 (k ==1 .and. s% use_RSP_L_eqn_outer_BC) then
257 0 : call set_RSP_Lsurf_BC(s, nvar, ierr)
258 0 : return
259 : end if
260 :
261 52304 : if (s% use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn) then
262 0 : call do1_gradT_eqn(s, k, nvar, ierr)
263 0 : return
264 : end if
265 :
266 52304 : if (s% use_dPrad_dm_form_of_T_gradient_eqn) then
267 0 : call do1_alt_dlnT_dm_eqn(s, k, nvar, ierr)
268 0 : return
269 : end if
270 :
271 : ! dT/dm = dP/dm * T/P * grad_T, grad_T = dlnT/dlnP from MLT.
272 : ! but use hydrostatic value for dP/dm in this.
273 : ! this is because of limitations of MLT for calculating grad_T.
274 : ! (MLT assumes hydrostatic equilibrium)
275 : ! see comment in K&W chpt 9.1.
276 :
277 52304 : call eval_dlnPdm_qhse(s, k, dlnPdm, Ppoint, ierr)
278 52304 : if (ierr /= 0) return
279 :
280 52304 : gradT = s% gradT_ad(k)
281 52304 : dlnTdm = dlnPdm*gradT
282 :
283 52304 : Tm1 = wrap_T_m1(s,k)
284 52304 : T00 = wrap_T_00(s,k)
285 52304 : dT = Tm1 - T00
286 52304 : alfa = s% dm(k-1)/(s% dm(k-1) + s% dm(k))
287 52304 : Tpoint = alfa*T00 + (1d0 - alfa)*Tm1
288 52304 : lnTdiff = dT/Tpoint ! use this in place of lnT(k-1)-lnT(k)
289 52304 : delm = (s% dm(k) + s% dm(k-1))/2
290 :
291 52304 : resid = delm*dlnTdm - lnTdiff
292 52304 : s% equ(i_equL, k) = resid%val
293 :
294 52304 : if (is_bad(s% equ(i_equL, k))) then
295 0 : ierr = -1
296 0 : if (s% report_ierr) write(*,2) 'equ(i_equL, k)', k, s% equ(i_equL, k)
297 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'hydro eqns')
298 0 : return
299 : write(*,2) 'equ(i_equL, k)', k, s% equ(i_equL, k)
300 : write(*,2) 'lnTdiff', k, lnTdiff
301 : write(*,2) 'delm', k, delm
302 : write(*,2) 'dlnPdm', k, dlnPdm
303 : write(*,2) 'gradT', k, gradT
304 : call mesa_error(__FILE__,__LINE__,'i_equL')
305 : end if
306 :
307 : if (test_partials) then
308 : s% solver_test_partials_val = s% equ(i_equL,k)
309 : end if
310 :
311 : call save_eqn_residual_info( &
312 52304 : s, k, nvar, i_equL, resid, 'do1_dlnT_dm_eqn', ierr)
313 :
314 52304 : end subroutine do1_dlnT_dm_eqn
315 :
316 :
317 :
318 0 : subroutine set_RSP_Lsurf_BC(s, nvar, ierr)
319 52304 : use const_def, only: crad, clight, pi4
320 : use eos_def
321 : use star_utils, only: save_eqn_residual_info, get_area_info_opt_time_center
322 : use auto_diff_support
323 : implicit none
324 :
325 : type(star_info), pointer :: s
326 : integer, intent(out) :: ierr
327 : integer, intent(in) :: nvar
328 :
329 : type(auto_diff_real_star_order1) :: L1_ad, r1_ad, area_ad, rhs_ad, lhs_ad, resid_ad, inv_R2
330 : type(auto_diff_real_star_order1) :: T_surf, Erad_ad
331 : integer :: i_equL
332 0 : real(dp) :: factor, scale, L_theta
333 : logical :: debug
334 :
335 0 : ierr = 0
336 0 : debug = .false.
337 :
338 0 : i_equL = s% i_equL
339 :
340 0 : if (s%nz < 1) then
341 0 : write(*,*) 'ERROR: Insufficient zones (nz < 1)'
342 0 : ierr = -1
343 : return
344 : end if
345 :
346 : if (debug) write(*,*) 'RSP zone 1 surface BC being set'
347 :
348 0 : call get_area_info_opt_time_center(s, 1, area_ad, inv_R2, ierr)
349 : ! no time centering the surface equations.
350 0 : L1_ad = wrap_L_00(s, 1)
351 0 : T_surf = wrap_T_00(s,1)
352 :
353 : if (debug) then
354 : write(*,*) 'T_surf =', T_surf%val, ' r_surf =', r1_ad%val, ' area =', area_ad%val
355 : end if
356 :
357 : ! rsp equation, zone 1
358 0 : rhs_ad = s%RSP2_Lsurf_factor * area_ad * clight * (crad * pow4(T_surf)) ! missing Lc at the moment, so only radiative surface
359 :
360 : if (debug) then
361 : write(*,*) 'RSP_Lsurf_factor =', s%RSP2_Lsurf_factor
362 : write(*,*) 'rhs_ad (RSP BC) =', rhs_ad%val
363 : end if
364 :
365 : ! residual
366 0 : lhs_ad = L1_ad
367 0 : resid_ad = lhs_ad - rhs_ad
368 :
369 0 : scale =maxval(s% L_start(1:s% nz))
370 0 : resid_ad = resid_ad / scale
371 :
372 : if (debug) then
373 : write(*,*) 'lhs (L1) =', lhs_ad%val
374 : write(*,*) 'scaled residual =', resid_ad%val
375 : end if
376 :
377 0 : s%equ(i_equL,1) = resid_ad%val
378 :
379 0 : if (is_bad(resid_ad%val)) then
380 0 : write(*,*) 'ERROR: NaN or Inf residual:', resid_ad%val
381 : ierr = -1
382 : end if
383 :
384 : call save_eqn_residual_info( &
385 0 : s, 1, nvar, i_equL, resid_ad, 'do1_dlnT_dm_eqn', ierr)
386 :
387 :
388 0 : end subroutine set_RSP_Lsurf_BC
389 :
390 : ! only used for dlnT_dm equation
391 52304 : subroutine eval_dlnPdm_qhse(s, k, & ! calculate the expected dlnPdm for HSE
392 : dlnPdm_qhse, Ppoint, ierr)
393 0 : use hydro_momentum, only: expected_HSE_grav_term
394 : type (star_info), pointer :: s
395 : integer, intent(in) :: k
396 : type(auto_diff_real_star_order1), intent(out) :: dlnPdm_qhse, Ppoint
397 : integer, intent(out) :: ierr
398 :
399 52304 : real(dp) :: alfa
400 : type(auto_diff_real_star_order1) :: grav, area, P00, Pm1, inv_R2, mlt_Ptrb00, mlt_Ptrbm1
401 : include 'formats'
402 :
403 : ierr = 0
404 :
405 : ! basic eqn is dP/dm = -G m / (4 pi r^4)
406 : ! divide by Ppoint to make it unitless
407 :
408 : ! for rotation, multiply gravity by factor fp. MESA 2, eqn 22.
409 52304 : call expected_HSE_grav_term(s, k, grav, area, ierr) ! note that expected_HSE_grav_term is negative
410 :
411 : ! mlt_pturb in thermodynamic gradients does not currently support time centering because it is timelagged.
412 : ! replace mlt_vc check with s% mlt_vc_old(k) >0 check.
413 : if ((s% have_mlt_vc .and. s% okay_to_set_mlt_vc) .and. s% include_mlt_Pturb_in_thermodynamic_gradients &
414 52304 : .and. s% mlt_Pturb_factor > 0d0) then
415 0 : if (k ==1) then
416 0 : mlt_Ptrb00 = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*wrap_d_00(s,k)/3d0
417 0 : mlt_Ptrbm1 = 0d0
418 : else
419 0 : mlt_Ptrb00 = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*wrap_d_00(s,k)/3d0
420 0 : mlt_Ptrbm1 = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*wrap_d_m1(s,k)/3d0
421 : end if
422 : else ! no mlt_pturb
423 52304 : mlt_Ptrb00 = 0d0
424 52304 : mlt_Ptrbm1 = 0d0
425 : end if
426 :
427 52304 : P00 = wrap_Peos_00(s,k)
428 :
429 : ! mlt Pturb doesn't support time centering yet.
430 52304 : if (s% using_velocity_time_centering) P00 = 0.5d0*(P00 + s% Peos_start(k))
431 :
432 52304 : if (k == 1) then
433 0 : Pm1 = 0d0
434 0 : Ppoint = P00 + mlt_Ptrb00
435 : else
436 52304 : Pm1 = wrap_Peos_m1(s,k)
437 52304 : if (s% using_velocity_time_centering) Pm1 = 0.5d0*(Pm1 + s% Peos_start(k-1)) ! pm1 wasn't time centered until now
438 52304 : Pm1 = Pm1 + mlt_Ptrbm1 ! include mlt Ptrb in k-1
439 52304 : P00 = P00 + mlt_Ptrb00 ! include mlt Ptrb in k
440 52304 : alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k))
441 52304 : Ppoint = alfa*P00 + (1d0-alfa)*Pm1
442 : end if
443 :
444 52304 : dlnPdm_qhse = grav/(area*Ppoint) ! note that expected_HSE_grav_term is negative
445 :
446 52304 : if (is_bad(dlnPdm_qhse%val)) then
447 0 : ierr = -1
448 0 : s% retry_message = 'eval_dlnPdm_qhse: is_bad(dlnPdm_qhse)'
449 0 : if (s% report_ierr) then
450 0 : !$OMP critical (hydro_vars_crit1)
451 0 : write(*,*) 'eval_dlnPdm_qhse: is_bad(dlnPdm_qhse)'
452 0 : stop
453 : !$OMP end critical (hydro_vars_crit1)
454 : end if
455 0 : if (s% stop_for_bad_nums) then
456 0 : write(*,2) 'dlnPdm_qhse', k, dlnPdm_qhse
457 0 : call mesa_error(__FILE__,__LINE__,'eval_dlnPdm_qhse')
458 : end if
459 : return
460 : end if
461 :
462 52304 : end subroutine eval_dlnPdm_qhse
463 :
464 : end module hydro_temperature
|
