Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2020 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_rsp2
21 :
22 : use star_private_def
23 : use const_def, only: dp, boltz_sigma, pi, clight, crad, ln10
24 : use utils_lib, only: is_bad
25 : use auto_diff
26 : use auto_diff_support
27 : use accurate_sum_auto_diff_star_order1
28 : use star_utils
29 :
30 : implicit none
31 :
32 : private
33 : public :: do1_rsp2_L_eqn
34 : public :: do1_turbulent_energy_eqn
35 : public :: do1_rsp2_Hp_eqn
36 : public :: compute_Eq_cell
37 : public :: compute_Uq_face
38 : public :: set_RSP2_vars
39 : public :: Hp_face_for_rsp2_val
40 : public :: Hp_face_for_rsp2_eqn, set_etrb_start_vars
41 : public :: RSP2_adjust_vars_before_call_solver
42 : public :: get_RSP2_alfa_beta_face_weights
43 :
44 : real(dp), parameter :: &
45 : x_ALFAP = 2.d0/3.d0, & ! Ptrb
46 : x_ALFAS = (1.d0/2.d0)*sqrt_2_div_3, & ! PII_face and Lc
47 : x_ALFAC = (1.d0/2.d0)*sqrt_2_div_3, & ! Lc
48 : x_CEDE = (8.d0/3.d0)*sqrt_2_div_3, & ! DAMP
49 : x_GAMMAR = 2.d0*sqrt(3.d0) ! DAMPR
50 :
51 : contains
52 :
53 0 : subroutine set_RSP2_vars(s,ierr)
54 : type (star_info), pointer :: s
55 : integer, intent(out) :: ierr
56 : type(auto_diff_real_star_order1) :: x
57 : integer :: k, op_err
58 : include 'formats'
59 0 : ierr = 0
60 0 : op_err = 0
61 0 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
62 : do k=1,s%nz
63 : ! Hp_face(k) <= 0 means it needs to be set. e.g., after read file
64 : if (s% Hp_face(k) <= 0) then
65 : s% Hp_face(k) = get_scale_height_face_val(s,k)
66 : s% xh(s% i_Hp,k) = s% Hp_face(k)
67 : end if
68 : x = compute_Y_face(s, k, op_err)
69 : if (op_err /= 0) ierr = op_err
70 : x = compute_PII_face(s, k, op_err)
71 : if (op_err /= 0) ierr = op_err
72 : !Pvsc skip?
73 : end do
74 : !$OMP END PARALLEL DO
75 0 : if (ierr /= 0) then
76 0 : if (s% report_ierr) write(*,2) 'failed in set_RSP2_vars loop 1', s% model_number
77 0 : return
78 : end if
79 0 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
80 : do k=1,s% nz
81 : x = compute_Chi_cell(s, k, op_err)
82 : if (op_err /= 0) ierr = op_err
83 : x = compute_Eq_cell(s, k, op_err)
84 : if (op_err /= 0) ierr = op_err
85 : x = compute_C(s, k, op_err) ! COUPL
86 : if (op_err /= 0) ierr = op_err
87 : x = compute_L_face(s, k, op_err) ! Lr, Lt, Lc
88 : if (op_err /= 0) ierr = op_err
89 : end do
90 : !$OMP END PARALLEL DO
91 0 : if (ierr /= 0) then
92 0 : if (s% report_ierr) write(*,2) 'failed in set_RSP2_vars loop 2', s% model_number
93 0 : return
94 : end if
95 0 : do k = 1, s% RSP2_num_outermost_cells_forced_nonturbulent
96 0 : s% Eq(k) = 0d0; s% Eq_ad(k) = 0d0
97 0 : s% Chi(k) = 0d0; s% Chi_ad(k) = 0d0
98 0 : s% COUPL(k) = 0d0; s% COUPL_ad(k) = 0d0
99 : !s% Ptrb(k) = 0d0;
100 0 : s% Lc(k) = 0d0; s% Lc_ad(k) = 0d0
101 0 : s% Lt(k) = 0d0; s% Lt_ad(k) = 0d0
102 : end do
103 0 : do k = s% nz + 1 - int(s% nz/s% RSP2_nz_div_IBOTOM) , s% nz
104 0 : s% Eq(k) = 0d0; s% Eq_ad(k) = 0d0
105 0 : s% Chi(k) = 0d0; s% Chi_ad(k) = 0d0
106 0 : s% COUPL(k) = 0d0; s% COUPL_ad(k) = 0d0
107 : !s% Ptrb(k) = 0d0;
108 0 : s% Lc(k) = 0d0; s% Lc_ad(k) = 0d0
109 0 : s% Lt(k) = 0d0; s% Lt_ad(k) = 0d0
110 : end do
111 : end subroutine set_RSP2_vars
112 :
113 :
114 0 : subroutine do1_rsp2_L_eqn(s, k, nvar, ierr)
115 : use star_utils, only: save_eqn_residual_info
116 : type (star_info), pointer :: s
117 : integer, intent(in) :: k, nvar
118 : integer, intent(out) :: ierr
119 : type(auto_diff_real_star_order1) :: &
120 : L_expected, L_actual,resid
121 : type(accurate_auto_diff_real_star_order1) :: L_sum
122 : real(dp) :: scale, residual, L_start_max
123 : logical :: test_partials
124 : include 'formats'
125 :
126 : !test_partials = (k == s% solver_test_partials_k)
127 0 : test_partials = .false.
128 0 : if (.not. s% RSP2_flag) then
129 0 : ierr = -1
130 0 : return
131 : end if
132 :
133 0 : ierr = 0
134 : !L_expected = compute_L_face(s, k, ierr)
135 : !if (ierr /= 0) return
136 0 : L_sum = s% Lr_ad(k)
137 0 : L_sum = L_sum + s% Lc_ad(k)
138 0 : L_sum = L_sum + s% Lt_ad(k)
139 0 : L_expected = L_sum
140 0 : L_actual = wrap_L_00(s, k)
141 0 : L_start_max = maxval(s% L_start(1:s% nz))
142 0 : scale = 1d0/L_start_max
143 0 : if (is_bad(scale)) then
144 0 : write(*,2) 'do1_rsp2_L_eqn scale', k, scale
145 0 : call mesa_error(__FILE__,__LINE__,'do1_rsp2_L_eqn')
146 : end if
147 0 : resid = (L_expected - L_actual)*scale
148 :
149 0 : residual = resid%val
150 0 : s% equ(s% i_equL, k) = residual
151 : if (test_partials) then
152 : s% solver_test_partials_val = residual
153 : end if
154 :
155 0 : call save_eqn_residual_info(s, k, nvar, s% i_equL, resid, 'do1_rsp2_L_eqn', ierr)
156 0 : if (ierr /= 0) return
157 :
158 : if (test_partials) then
159 : s% solver_test_partials_var = s% i_lnR
160 : s% solver_test_partials_dval_dx = resid%d1Array(i_lnR_00)
161 : write(*,4) 'do1_rsp2_L_eqn', s% solver_test_partials_var
162 : end if
163 : end subroutine do1_rsp2_L_eqn
164 :
165 :
166 0 : subroutine do1_rsp2_Hp_eqn(s, k, nvar, ierr)
167 : use star_utils, only: save_eqn_residual_info
168 : type (star_info), pointer :: s
169 : integer, intent(in) :: k, nvar
170 : integer, intent(out) :: ierr
171 : type(auto_diff_real_star_order1) :: &
172 : Hp_expected, Hp_actual,resid
173 : real(dp) :: residual, Hp_start
174 : logical :: test_partials
175 : include 'formats'
176 : !test_partials = (k == s% solver_test_partials_k)
177 0 : test_partials = .false.
178 :
179 0 : if (.not. s% RSP2_flag) then
180 0 : ierr = -1
181 0 : return
182 : end if
183 :
184 : ierr = 0
185 0 : Hp_expected = Hp_face_for_rsp2_eqn(s, k, ierr)
186 0 : if (ierr /= 0) return
187 0 : Hp_actual = wrap_Hp_00(s, k)
188 0 : Hp_start = s% Hp_face_start(k)
189 0 : resid = (Hp_expected - Hp_actual)/max(Hp_expected,Hp_actual)
190 :
191 0 : residual = resid%val
192 0 : s% equ(s% i_equ_Hp, k) = residual
193 : if (test_partials) then
194 : s% solver_test_partials_val = residual
195 : end if
196 :
197 0 : if (residual > 1d3) then
198 0 : !$omp critical (hydro_rsp2_1)
199 0 : write(*,2) 'residual', k, residual
200 0 : write(*,2) 'Hp_expected', k, Hp_expected%val
201 0 : write(*,2) 'Hp_actual', k, Hp_actual%val
202 0 : call mesa_error(__FILE__,__LINE__,'do1_rsp2_Hp_eqn')
203 : !$omp end critical (hydro_rsp2_1)
204 : end if
205 :
206 0 : call save_eqn_residual_info(s, k, nvar, s% i_equ_Hp, resid, 'do1_rsp2_Hp_eqn', ierr)
207 0 : if (ierr /= 0) return
208 :
209 : if (test_partials) then
210 : s% solver_test_partials_var = s% i_lnR
211 : s% solver_test_partials_dval_dx = resid%d1Array(i_lnR_00)
212 : write(*,4) 'do1_rsp2_Hp_eqn', s% solver_test_partials_var
213 : end if
214 :
215 : end subroutine do1_rsp2_Hp_eqn
216 :
217 :
218 0 : real(dp) function Hp_face_for_rsp2_val(s, k, ierr) result(Hp_face) ! cm
219 : type (star_info), pointer :: s
220 : integer, intent(in) :: k
221 : integer, intent(out) :: ierr
222 : type(auto_diff_real_star_order1) :: Hp_face_ad
223 : ierr = 0
224 0 : Hp_face_ad = Hp_face_for_rsp2_eqn(s, k, ierr)
225 0 : if (ierr /= 0) return
226 0 : Hp_face = Hp_face_ad%val
227 0 : end function Hp_face_for_rsp2_val
228 :
229 :
230 0 : function Hp_face_for_rsp2_eqn(s, k, ierr) result(Hp_face) ! cm
231 : type (star_info), pointer :: s
232 : integer, intent(in) :: k
233 : integer, intent(out) :: ierr
234 : type(auto_diff_real_star_order1) :: Hp_face
235 : type(auto_diff_real_star_order1) :: &
236 : rho_face, area, dlnPeos, &
237 : r_00, Peos_00, d_00, Peos_m1, d_m1, Peos_div_rho, &
238 : d_face, Peos_face, alt_Hp_face, A
239 : real(dp) :: alfa, beta
240 : include 'formats'
241 0 : ierr = 0
242 0 : if (k > s% nz) then
243 0 : Hp_face = 1d0 ! not used
244 0 : return
245 : end if
246 0 : if (k > 1 .and. .not. s% RSP2_assume_HSE) then
247 0 : call get_RSP2_alfa_beta_face_weights(s, k, alfa, beta)
248 0 : rho_face = alfa*wrap_d_00(s,k) + beta*wrap_d_m1(s,k)
249 0 : area = 4d0*pi*pow2(wrap_r_00(s,k))
250 0 : dlnPeos = wrap_lnPeos_m1(s,k) - wrap_lnPeos_00(s,k)
251 0 : Hp_face = -s% dm_bar(k)/(area*rho_face*dlnPeos)
252 : else
253 0 : r_00 = wrap_r_00(s, k) ! not time-centered in RSP
254 0 : d_00 = wrap_d_00(s, k)
255 0 : Peos_00 = wrap_Peos_00(s, k)
256 0 : if (k == 1) then
257 0 : Peos_div_rho = Peos_00/d_00
258 0 : Hp_face = pow2(r_00)*Peos_div_rho/(s% cgrav(k)*s% m(k))
259 : else
260 0 : d_m1 = wrap_d_m1(s, k)
261 0 : Peos_m1 = wrap_Peos_m1(s, k)
262 0 : call get_RSP2_alfa_beta_face_weights(s, k, alfa, beta)
263 0 : Peos_div_rho = alfa*Peos_00/d_00 + beta*Peos_m1/d_m1
264 0 : Hp_face = pow2(r_00)*Peos_div_rho/(s% cgrav(k)*s% m(k))
265 0 : if (k==-104) then
266 0 : write(*,3) 'RSP2 Hp P_div_rho Pdrho_00 Pdrho_m1', k, s% solver_iter, &
267 0 : Hp_face%val, Peos_div_rho%val, Peos_00%val/d_00%val, Peos_m1%val/d_m1%val
268 : !write(*,3) 'RSP2 Hp r2_div_Gm r_start r', k, s% solver_iter, &
269 : ! Hp_face%val, pow2(r_00%val)/(s% cgrav(k)*s% m(k)), &
270 : ! s% r_start(k), r_00%val
271 : end if
272 0 : if (s% alt_scale_height_flag) then
273 0 : call mesa_error(__FILE__,__LINE__,'Hp_face_for_rsp2_eqn: cannot use alt_scale_height_flag')
274 : ! consider sound speed*hydro time scale as an alternative scale height
275 : d_face = alfa*d_00 + beta*d_m1
276 : Peos_face = alfa*Peos_00 + beta*Peos_m1
277 : alt_Hp_face = sqrt(Peos_face/s% cgrav(k))/d_face
278 : if (alt_Hp_face%val < Hp_face%val) then ! blend
279 : A = pow2(alt_Hp_face/Hp_face) ! 0 <= A%val < 1
280 : Hp_face = A*Hp_face + (1d0 - A)*alt_Hp_face
281 : end if
282 : end if
283 : end if
284 : end if
285 0 : end function Hp_face_for_rsp2_eqn
286 :
287 :
288 0 : subroutine do1_turbulent_energy_eqn(s, k, nvar, ierr)
289 : use star_utils, only: set_energy_eqn_scal, save_eqn_residual_info
290 : type (star_info), pointer :: s
291 : integer, intent(in) :: k, nvar
292 : integer, intent(out) :: ierr
293 : ! for OLD WAY
294 : type(auto_diff_real_star_order1) :: &
295 : d_turbulent_energy_ad, Ptrb_dV_ad, dt_C_ad, dt_Eq_ad
296 : type(auto_diff_real_star_order1) :: w_00
297 : type(auto_diff_real_star_order1) :: tst, resid_ad, dt_dLt_dm_ad
298 : type(accurate_auto_diff_real_star_order1) :: esum_ad
299 : logical :: non_turbulent_cell, test_partials
300 : real(dp) :: residual, scal
301 : include 'formats'
302 : !test_partials = (k == s% solver_test_partials_k)
303 0 : test_partials = .false.
304 :
305 0 : ierr = 0
306 0 : w_00 = wrap_w_00(s,k)
307 :
308 : non_turbulent_cell = &
309 : s% mixing_length_alpha == 0d0 .or. &
310 : k <= s% RSP2_num_outermost_cells_forced_nonturbulent .or. &
311 0 : k > s% nz - int(s% nz/s% RSP2_nz_div_IBOTOM)
312 0 : if (.not. s% RSP2_flag) then
313 0 : resid_ad = w_00 - s% w_start(k) ! just hold w constant when not using RSP2
314 0 : else if (non_turbulent_cell) then
315 0 : resid_ad = w_00/s% csound(k) ! make w = 0
316 : else
317 0 : call setup_d_turbulent_energy(ierr); if (ierr /= 0) return ! erg g^-1 = cm^2 s^-2
318 0 : call setup_Ptrb_dV_ad(ierr); if (ierr /= 0) return ! erg g^-1
319 0 : call setup_dt_dLt_dm_ad(ierr); if (ierr /= 0) return ! erg g^-1
320 0 : call setup_dt_C_ad(ierr); if (ierr /= 0) return ! erg g^-1
321 0 : call setup_dt_Eq_ad(ierr); if (ierr /= 0) return ! erg g^-1
322 0 : call set_energy_eqn_scal(s, k, scal, ierr); if (ierr /= 0) return ! 1/(erg g^-1 s^-1)
323 : ! sum terms in esum_ad using accurate_auto_diff_real_star_order1
324 0 : esum_ad = d_turbulent_energy_ad
325 0 : esum_ad = esum_ad + Ptrb_dV_ad
326 0 : esum_ad = esum_ad + dt_dLt_dm_ad
327 0 : esum_ad = esum_ad - dt_C_ad
328 0 : esum_ad = esum_ad - dt_Eq_ad ! erg g^-1
329 0 : resid_ad = esum_ad
330 :
331 0 : if (k==-35 .and. s% solver_iter == 1) then
332 0 : write(*,3) 'RSP2 w dEt PdV dtC dtEq', k, s% solver_iter, &
333 0 : w_00%val, d_turbulent_energy_ad%val, Ptrb_dV_ad%val, dt_C_ad%val, dt_Eq_ad%val
334 : end if
335 :
336 0 : resid_ad = resid_ad*scal/s%dt ! to make residual unitless, must cancel out the dt in scal
337 :
338 : end if
339 :
340 0 : residual = resid_ad%val
341 0 : s% equ(s% i_detrb_dt, k) = residual
342 :
343 : if (test_partials) then
344 : tst = residual
345 : s% solver_test_partials_val = tst%val
346 : if (s% solver_iter == 12) &
347 : write(*,*) 'do1_turbulent_energy_eqn', s% solver_test_partials_var, s% lnd(k), tst%val
348 : end if
349 :
350 0 : call save_eqn_residual_info(s, k, nvar, s% i_detrb_dt, resid_ad, 'do1_turbulent_energy_eqn', ierr)
351 0 : if (ierr /= 0) return
352 :
353 : if (test_partials) then
354 : s% solver_test_partials_var = s% i_lnd
355 : s% solver_test_partials_dval_dx = tst%d1Array(i_lnd_00) ! xi0 good , xi1 partial 0, xi2 good. Af horrible.'
356 : write(*,*) 'do1_turbulent_energy_eqn', s% solver_test_partials_var, s% lnd(k)/ln10, tst%val
357 : end if
358 :
359 : contains
360 :
361 0 : subroutine setup_d_turbulent_energy(ierr) ! erg g^-1
362 : integer, intent(out) :: ierr
363 0 : ierr = 0
364 0 : d_turbulent_energy_ad = wrap_etrb_00(s,k) - get_etrb_start(s,k)
365 0 : end subroutine setup_d_turbulent_energy
366 :
367 : ! Ptrb_dV_ad = Ptrb_ad*dV_ad
368 0 : subroutine setup_Ptrb_dV_ad(ierr) ! erg g^-1
369 : use star_utils, only: calc_Ptrb_ad_tw
370 : integer, intent(out) :: ierr
371 : type(auto_diff_real_star_order1) :: Ptrb_ad, PT0, dV_ad, d_00
372 0 : call calc_Ptrb_ad_tw(s, k, Ptrb_ad, PT0, ierr)
373 0 : if (ierr /= 0) return
374 0 : d_00 = wrap_d_00(s,k)
375 0 : dV_ad = 1d0/d_00 - 1d0/s% rho_start(k)
376 0 : Ptrb_dV_ad = Ptrb_ad*dV_ad ! erg cm^-3 cm^-3 g^-1 = erg g^-1
377 : end subroutine setup_Ptrb_dV_ad
378 :
379 0 : subroutine setup_dt_dLt_dm_ad(ierr) ! erg g^-1
380 : integer, intent(out) :: ierr
381 : type(auto_diff_real_star_order1) :: Lt_00, Lt_p1
382 : real(dp) :: L_theta
383 : include 'formats'
384 0 : ierr = 0
385 0 : if (s% using_velocity_time_centering .and. &
386 : s% include_L_in_velocity_time_centering) then
387 0 : L_theta = s% L_theta_for_velocity_time_centering
388 : else
389 0 : L_theta = 1d0
390 : end if
391 0 : Lt_00 = L_theta*s% Lt_ad(k) + (1d0 - L_theta)*s% Lt_start(k)
392 0 : if (k == s% nz) then
393 0 : Lt_p1 = 0d0
394 : else
395 0 : Lt_p1 = L_theta*shift_p1(s% Lt_ad(k+1)) + (1d0 - L_theta)*s% Lt_start(k+1)
396 0 : if (ierr /= 0) return
397 : end if
398 0 : dt_dLt_dm_ad = (Lt_00 - Lt_p1)*s%dt/s%dm(k)
399 : end subroutine setup_dt_dLt_dm_ad
400 :
401 0 : subroutine setup_dt_C_ad(ierr) ! erg g^-1
402 : integer, intent(out) :: ierr
403 : type(auto_diff_real_star_order1) :: C
404 0 : C = s% COUPL_ad(k) ! compute_C(s, k, ierr) ! erg g^-1 s^-1
405 0 : if (ierr /= 0) return
406 0 : dt_C_ad = s%dt*C
407 : end subroutine setup_dt_C_ad
408 :
409 0 : subroutine setup_dt_Eq_ad(ierr) ! erg g^-1
410 : integer, intent(out) :: ierr
411 : type(auto_diff_real_star_order1) :: Eq_cell
412 0 : Eq_cell = s% Eq_ad(k) ! compute_Eq_cell(s, k, ierr) ! erg g^-1 s^-1
413 0 : if (ierr /= 0) return
414 0 : dt_Eq_ad = s%dt*Eq_cell
415 : end subroutine setup_dt_Eq_ad
416 :
417 : end subroutine do1_turbulent_energy_eqn
418 :
419 :
420 0 : subroutine get_RSP2_alfa_beta_face_weights(s, k, alfa, beta)
421 : type (star_info), pointer :: s
422 : integer, intent(in) :: k
423 : real(dp), intent(out) :: alfa, beta
424 : ! face_value = alfa*cell_value(k) + beta*cell_value(k-1)
425 0 : if (k == 1) call mesa_error(__FILE__,__LINE__,'bad k==1 for get_RSP2_alfa_beta_face_weights')
426 0 : if (s% RSP2_use_mass_interp_face_values) then
427 0 : alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k))
428 0 : beta = 1d0 - alfa
429 : else
430 0 : alfa = 0.5d0
431 0 : beta = 0.5d0
432 : end if
433 0 : end subroutine get_RSP2_alfa_beta_face_weights
434 :
435 :
436 0 : function compute_Y_face(s, k, ierr) result(Y_face) ! superadiabatic gradient [unitless]
437 : type (star_info), pointer :: s
438 : integer, intent(in) :: k
439 : integer, intent(out) :: ierr
440 : type(auto_diff_real_star_order1) :: Y_face
441 : type(auto_diff_real_star_order1) :: Hp_face, Y1, Y2, QQ_div_Cp_face, &
442 : r_00, d_00, Peos_00, Cp_00, T_00, chiT_00, chiRho_00, QQ_00, lnT_00, &
443 : r_m1, d_m1, Peos_m1, Cp_m1, T_m1, chiT_m1, chiRho_m1, QQ_m1, lnT_m1, &
444 : dlnT_dlnP, grad_ad_00, grad_ad_m1, grad_ad_face, dlnT, dlnP, alt_Y_face
445 : real(dp) :: dm_bar, alfa, beta
446 : include 'formats'
447 0 : ierr = 0
448 :
449 0 : if (k > s% nz) then
450 0 : Y_face = 0d0
451 0 : return
452 : end if
453 :
454 0 : if (k == 1 .or. s% mixing_length_alpha == 0d0) then
455 0 : Y_face = 0d0
456 0 : s% Y_face(k) = 0d0
457 0 : s% Y_face_ad(k) = 0d0
458 0 : return
459 : end if
460 :
461 0 : call get_RSP2_alfa_beta_face_weights(s, k, alfa, beta)
462 :
463 0 : if (s% RSP2_use_RSP_eqn_for_Y_face) then
464 :
465 0 : dm_bar = s% dm_bar(k)
466 0 : Hp_face = wrap_Hp_00(s,k)
467 0 : r_00 = wrap_r_00(s, k)
468 0 : d_00 = wrap_d_00(s, k)
469 0 : Peos_00 = wrap_Peos_00(s, k)
470 0 : Cp_00 = wrap_Cp_00(s, k)
471 0 : T_00 = wrap_T_00(s, k)
472 0 : chiT_00 = wrap_chiT_00(s, k)
473 0 : chiRho_00 = wrap_chiRho_00(s, k)
474 0 : QQ_00 = chiT_00/(d_00*T_00*chiRho_00)
475 0 : lnT_00 = wrap_lnT_00(s,k)
476 :
477 0 : r_m1 = wrap_r_m1(s, k)
478 0 : d_m1 = wrap_d_m1(s, k)
479 0 : Peos_m1 = wrap_Peos_m1(s, k)
480 0 : Cp_m1 = wrap_Cp_m1(s, k)
481 0 : T_m1 = wrap_T_m1(s, k)
482 0 : chiT_m1 = wrap_chiT_m1(s, k)
483 0 : chiRho_m1 = wrap_chiRho_m1(s, k)
484 0 : QQ_m1 = chiT_m1/(d_m1*T_m1*chiRho_m1)
485 0 : lnT_m1 = wrap_lnT_m1(s,k)
486 0 : QQ_div_Cp_face = alfa*QQ_00/Cp_00 + beta*QQ_m1/Cp_m1
487 : ! QQ units (g cm^-3 K)^-1 = g^-1 cm^3 K^-1
488 : ! Cp units erg g^-1 K^-1 = g cm^2 s^-2 g^-1 K^-1 = cm^2 s^-2 K^-1
489 : ! QQ/Cp units = (g^-1 cm^3 K^-1)/(cm^2 s^-2 K^-1)
490 : ! = g^-1 cm^3 K^-1 cm^-2 s^2 K
491 : ! = g^-1 cm s^2
492 : ! P units = erg cm^-3 = g cm^2 s^-2 cm^-3 = g cm^-1 s^-2
493 : ! QQ/Cp*P is unitless.
494 :
495 0 : Y1 = QQ_div_Cp_face*(Peos_m1 - Peos_00) - (lnT_m1 - lnT_00)
496 : ! Y1 unitless
497 :
498 0 : Y2 = 4d0*pi*pow2(r_00)*Hp_face*2d0/(1d0/d_00 + 1d0/d_m1)/dm_bar
499 : ! units = cm^2 cm / (cm^3 g^-1) / g
500 : ! = cm^2 cm cm^-3 g g^-1 = unitless
501 :
502 0 : Y_face = Y1*Y2 ! unitless
503 :
504 0 : if (k==-35) then
505 0 : write(*,3) 'RSP2 Y_face Y1 Y2', k, s% solver_iter, s% Y_face(k), Y1%val, Y2%val
506 0 : write(*,3) 'Peos', k, s% solver_iter, Peos_00%val
507 0 : write(*,3) 'Peos', k-1, s% solver_iter, Peos_m1%val
508 0 : write(*,3) 'QQ', k, s% solver_iter, QQ_00%val
509 0 : write(*,3) 'QQ', k-1, s% solver_iter, QQ_m1%val
510 0 : write(*,3) 'Cp', k, s% solver_iter, Cp_00%val
511 0 : write(*,3) 'Cp', k-1, s% solver_iter, Cp_m1%val
512 0 : write(*,3) 'lgT', k, s% solver_iter, lnT_00%val/ln10
513 0 : write(*,3) 'lgT', k-1, s% solver_iter, lnT_m1%val/ln10
514 0 : write(*,3) 'lgd', k, s% solver_iter, s% lnd(k)/ln10
515 0 : write(*,3) 'lgd', k-1, s% solver_iter, s% lnd(k-1)/ln10
516 : !call mesa_error(__FILE__,__LINE__,'compute_Y_face')
517 : end if
518 :
519 : else
520 :
521 0 : grad_ad_00 = wrap_grad_ad_00(s,k)
522 0 : grad_ad_m1 = wrap_grad_ad_m1(s,k)
523 0 : grad_ad_face = alfa*grad_ad_00 + beta*grad_ad_m1
524 0 : dlnT = wrap_lnT_m1(s,k) - wrap_lnT_00(s,k)
525 0 : dlnP = wrap_lnPeos_m1(s,k) - wrap_lnPeos_00(s,k)
526 0 : dlnT_dlnP = dlnT/dlnP
527 0 : if (is_bad(dlnT_dlnP%val)) then
528 0 : alt_Y_face = 0d0
529 0 : else if (s% use_Ledoux_criterion .and. s% calculate_Brunt_B) then
530 : ! gradL = grada + gradL_composition_term
531 0 : alt_Y_face = dlnT_dlnP - (grad_ad_face + s% gradL_composition_term(k))
532 : else
533 0 : alt_Y_face = dlnT_dlnP - grad_ad_face
534 : end if
535 0 : if (is_bad(alt_Y_face%val)) alt_Y_face = 0
536 0 : Y_face = alt_Y_face
537 :
538 : end if
539 :
540 0 : s% Y_face_ad(k) = Y_face
541 0 : s% Y_face(k) = Y_face%val
542 :
543 0 : end function compute_Y_face
544 :
545 :
546 0 : function compute_PII_face(s, k, ierr) result(PII_face) ! ergs g^-1 K^-1 (like Cp)
547 : type (star_info), pointer :: s
548 : integer, intent(in) :: k
549 : type(auto_diff_real_star_order1) :: PII_face
550 : integer, intent(out) :: ierr
551 : type(auto_diff_real_star_order1) :: Cp_00, Cp_m1, Cp_face, Y_face, T_00, T_m1
552 : type(auto_diff_real_star_order1) :: X, FL, scale, T_face, e_face, Peos_face, rho_face, h_face
553 : real(dp) :: ALFAS_ALFA, alfa, beta
554 : include 'formats'
555 0 : ierr = 0
556 0 : if (k > s% nz) then
557 0 : PII_face = 0d0
558 0 : return
559 : end if
560 0 : if (k == 1 .or. s% mixing_length_alpha == 0d0 .or. &
561 : k == s% nz) then ! just skip k == nz to be like RSP
562 0 : PII_face = 0d0
563 0 : s% PII(k) = 0d0
564 0 : s% PII_ad(k) = 0d0
565 0 : return
566 : end if
567 0 : Y_face = s% Y_face_ad(k) ! compute_Y_face(s, k, ierr)
568 0 : if (ierr /= 0) return
569 0 : Cp_00 = wrap_Cp_00(s, k)
570 0 : Cp_m1 = wrap_Cp_m1(s, k)
571 0 : T_00= wrap_T_00(s, k)
572 0 : T_m1 = wrap_T_m1(s, k)
573 0 : call get_RSP2_alfa_beta_face_weights(s, k, alfa, beta)
574 0 : Cp_face = alfa*Cp_00 + beta*Cp_m1 ! ergs g^-1 K^-1
575 0 : T_face = alfa*Cp_00 + beta*Cp_m1
576 0 : rho_face = alfa*wrap_d_00(s,k) + beta*wrap_d_m1(s,k)
577 0 : Peos_face = alfa*wrap_Peos_00(s,k) + beta*wrap_Peos_m1(s,k)
578 0 : e_face = alfa*wrap_e_00(s,k) + beta*wrap_e_m1(s,k)
579 0 : h_face = e_face + Peos_face/rho_face
580 0 : ALFAS_ALFA = x_ALFAS*s% mixing_length_alpha
581 0 : PII_face = ALFAS_ALFA*Cp_face*Y_face
582 :
583 0 : scale = 1d0
584 0 : if (Y_face > 0d0 .and. s% use_TDC_enthalpy_flux_limiter) then
585 : ! X = G/F
586 0 : X = (Cp_face*T_face/h_face)*ALFAS_ALFA* Y_face / sqrt_2_div_3
587 0 : FL = flux_limiter_function(X)
588 : ! Avoid 0/0 or tiny/tiny; for X ≈ 0, FL ≈ X so scale ~ 1 anyway.
589 0 : if (abs(X%val) >= 0.95d0) then
590 0 : scale = FL / X
591 : else
592 0 : scale = 1d0
593 : end if
594 : end if
595 :
596 0 : s% PII(k) = PII_face%val*scale%val
597 0 : s% PII_ad(k) = PII_face*scale
598 0 : if (k == -2 .and. s% PII(k) < 0d0) then
599 0 : write(*,2) 's% PII(k)', k, s% PII(k)
600 0 : write(*,2) 'Cp_face', k, Cp_face%val
601 0 : write(*,2) 'Y_face', k, Y_face%val
602 : !write(*,2) 'PII_face%val', k, PII_face%val
603 : !write(*,2) 'T_rho_face%val', k, T_rho_face%val
604 : !write(*,2) '', k,
605 : !write(*,2) '', k,
606 0 : call mesa_error(__FILE__,__LINE__,'compute_PII_face')
607 : end if
608 0 : end function compute_PII_face
609 :
610 0 : type(auto_diff_real_star_order1) function flux_limiter_function(X) result(FL) ! should be c2 continuous
611 : type(auto_diff_real_star_order1), intent(in) :: X
612 : real(dp), parameter :: X0 = 0.95_dp ! start of transition
613 : real(dp), parameter :: delta = 0.05_dp ! width of transition
614 : real(dp), parameter :: X1 = 1d0 !X0 + delta ! end of transition
615 :
616 : type(auto_diff_real_star_order1) :: s, p
617 :
618 : ! Region 1: purely linear, FL = X
619 0 : if (X%val < X0) then ! should not be encountered
620 0 : FL = X
621 :
622 : ! Region 3: saturated, FL = 1
623 0 : else if (X%val >= X1) then
624 0 : FL = 1.0_dp
625 :
626 : ! Region 2: smooth C² transition between the two
627 : else
628 : ! Normalized coordinate in [0,1]
629 0 : s = (X - X0) / (X1 - X0)
630 :
631 : ! Quintic "smootherstep" polynomial:
632 : ! p(s) = 10 s^3 - 15 s^4 + 6 s^5
633 : ! p(0)=0, p(1)=1, p'(0)=p'(1)=0, p''(0)=p''(1)=0
634 0 : p = pow3(s) * (10.0_dp + s * (-15.0_dp + 6.0_dp * s))
635 :
636 : ! Blend between line FL=X and flat FL=1
637 : ! At s=0: FL = X
638 : ! At s=1: FL = 1
639 : ! Because p', p'' vanish at 0 and 1, FL, FL', FL'' all match.
640 0 : FL = X + (1.0_dp - X) * p
641 : end if
642 0 : end function flux_limiter_function
643 :
644 0 : function compute_d_v_div_r(s, k, ierr) result(d_v_div_r) ! s^-1
645 : type (star_info), pointer :: s
646 : integer, intent(in) :: k
647 : type(auto_diff_real_star_order1) :: d_v_div_r
648 : integer, intent(out) :: ierr
649 : type(auto_diff_real_star_order1) :: v_00, v_p1, r_00, r_p1
650 : include 'formats'
651 0 : ierr = 0
652 0 : v_00 = wrap_v_00(s,k)
653 0 : v_p1 = wrap_v_p1(s,k)
654 0 : r_00 = wrap_r_00(s,k)
655 0 : r_p1 = wrap_r_p1(s,k)
656 0 : if (r_p1%val == 0d0) r_p1 = 1d0
657 0 : d_v_div_r = v_00/r_00 - v_p1/r_p1 ! units s^-1
658 0 : end function compute_d_v_div_r
659 :
660 :
661 0 : function compute_d_v_div_r_opt_time_center(s, k, ierr) result(d_v_div_r) ! s^-1
662 : type (star_info), pointer :: s
663 : integer, intent(in) :: k
664 : type(auto_diff_real_star_order1) :: d_v_div_r
665 : integer, intent(out) :: ierr
666 : type(auto_diff_real_star_order1) :: v_00, v_p1, r_00, r_p1
667 : include 'formats'
668 0 : ierr = 0
669 0 : v_00 = wrap_opt_time_center_v_00(s,k)
670 0 : v_p1 = wrap_opt_time_center_v_p1(s,k)
671 0 : r_00 = wrap_opt_time_center_r_00(s,k)
672 0 : r_p1 = wrap_opt_time_center_r_p1(s,k)
673 0 : if (r_p1%val == 0d0) r_p1 = 1d0
674 0 : d_v_div_r = v_00/r_00 - v_p1/r_p1 ! units s^-1
675 0 : end function compute_d_v_div_r_opt_time_center
676 :
677 :
678 0 : function wrap_Hp_cell(s, k) result(Hp_cell) ! cm
679 : type (star_info), pointer :: s
680 : integer, intent(in) :: k
681 : type(auto_diff_real_star_order1) :: Hp_cell
682 0 : Hp_cell = 0.5d0*(wrap_Hp_00(s,k) + wrap_Hp_p1(s,k))
683 0 : end function wrap_Hp_cell
684 :
685 :
686 0 : function Hp_cell_for_Chi(s, k, ierr) result(Hp_cell) ! cm
687 : type (star_info), pointer :: s
688 : integer, intent(in) :: k
689 : integer, intent(out) :: ierr
690 : type(auto_diff_real_star_order1) :: Hp_cell
691 : type(auto_diff_real_star_order1) :: d_00, Peos_00, rmid
692 : real(dp) :: mmid, cgrav_mid
693 : include 'formats'
694 0 : ierr = 0
695 :
696 0 : Hp_cell = wrap_Hp_cell(s, k)
697 0 : return
698 :
699 : d_00 = wrap_d_00(s, k)
700 : Peos_00 = wrap_Peos_00(s, k)
701 : if (k < s% nz) then
702 : rmid = 0.5d0*(wrap_r_00(s,k) + wrap_r_p1(s,k))
703 : mmid = 0.5d0*(s% m(k) + s% m(k+1))
704 : cgrav_mid = 0.5d0*(s% cgrav(k) + s% cgrav(k+1))
705 : else
706 : rmid = 0.5d0*(wrap_r_00(s,k) + s% r_center)
707 : mmid = 0.5d0*(s% m(k) + s% m_center)
708 : cgrav_mid = s% cgrav(k)
709 : end if
710 : Hp_cell = pow2(rmid)*Peos_00/(d_00*cgrav_mid*mmid)
711 : if (s% alt_scale_height_flag) then
712 : call mesa_error(__FILE__,__LINE__,'Hp_cell_for_Chi: cannot use alt_scale_height_flag')
713 : end if
714 : end function Hp_cell_for_Chi
715 :
716 :
717 0 : function compute_Chi_cell(s, k, ierr) result(Chi_cell)
718 : ! eddy viscosity energy (Kuhfuss 1986) [erg]
719 : type (star_info), pointer :: s
720 : integer, intent(in) :: k
721 : type(auto_diff_real_star_order1) :: Chi_cell
722 : integer, intent(out) :: ierr
723 : type(auto_diff_real_star_order1) :: &
724 : rho2, r6_cell, d_v_div_r, Hp_cell, w_00, d_00, r_00, r_p1
725 : real(dp) :: f, ALFAM_ALFA
726 : include 'formats'
727 0 : ierr = 0
728 0 : ALFAM_ALFA = s% RSP2_alfam*s% mixing_length_alpha
729 : if (ALFAM_ALFA == 0d0 .or. &
730 0 : k <= s% RSP2_num_outermost_cells_forced_nonturbulent .or. &
731 : k > s% nz - int(s% nz/s% RSP2_nz_div_IBOTOM)) then
732 0 : Chi_cell = 0d0
733 0 : if (k >= 1 .and. k <= s% nz) then
734 0 : s% Chi(k) = 0d0
735 0 : s% Chi_ad(k) = 0d0
736 : end if
737 : else
738 0 : Hp_cell = Hp_cell_for_Chi(s, k, ierr)
739 0 : if (ierr /= 0) return
740 0 : d_v_div_r = compute_d_v_div_r(s, k, ierr)
741 0 : if (ierr /= 0) return
742 0 : w_00 = wrap_w_00(s,k)
743 0 : d_00 = wrap_d_00(s,k)
744 0 : f = (16d0/3d0)*pi*ALFAM_ALFA/s% dm(k)
745 0 : rho2 = pow2(d_00)
746 0 : r_00 = wrap_r_00(s,k)
747 0 : r_p1 = wrap_r_p1(s,k)
748 0 : r6_cell = 0.5d0*(pow6(r_00) + pow6(r_p1))
749 0 : Chi_cell = f*rho2*r6_cell*d_v_div_r*Hp_cell*w_00
750 : ! units = g^-1 cm s^-1 g^2 cm^-6 cm^6 s^-1 cm
751 : ! = g cm^2 s^-2
752 : ! = erg
753 : end if
754 0 : s% Chi(k) = Chi_cell%val
755 0 : s% Chi_ad(k) = Chi_cell
756 :
757 0 : end function compute_Chi_cell
758 :
759 :
760 0 : function compute_Eq_cell(s, k, ierr) result(Eq_cell) ! erg g^-1 s^-1
761 : type (star_info), pointer :: s
762 : integer, intent(in) :: k
763 : type(auto_diff_real_star_order1) :: Eq_cell
764 : integer, intent(out) :: ierr
765 : type(auto_diff_real_star_order1) :: d_v_div_r, Chi_cell
766 : include 'formats'
767 0 : ierr = 0
768 : if (s% mixing_length_alpha == 0d0 .or. &
769 0 : k <= s% RSP2_num_outermost_cells_forced_nonturbulent .or. &
770 : k > s% nz - int(s% nz/s% RSP2_nz_div_IBOTOM)) then
771 0 : Eq_cell = 0d0
772 0 : if (k >= 1 .and. k <= s% nz) s% Eq_ad(k) = 0d0
773 : else
774 0 : Chi_cell = s% Chi_ad(k) ! compute_Chi_cell(s,k,ierr)
775 0 : if (ierr /= 0) return
776 0 : d_v_div_r = compute_d_v_div_r_opt_time_center(s, k, ierr)
777 0 : if (ierr /= 0) return
778 0 : Eq_cell = 4d0*pi*Chi_cell*d_v_div_r/s% dm(k) ! erg s^-1 g^-1
779 : end if
780 0 : s% Eq(k) = Eq_cell%val
781 0 : s% Eq_ad(k) = Eq_cell
782 0 : end function compute_Eq_cell
783 :
784 :
785 0 : function compute_Uq_face(s, k, ierr) result(Uq_face) ! cm s^-2, acceleration
786 : type (star_info), pointer :: s
787 : integer, intent(in) :: k
788 : type(auto_diff_real_star_order1) :: Uq_face
789 : integer, intent(out) :: ierr
790 : type(auto_diff_real_star_order1) :: Chi_00, Chi_m1, r_00
791 : include 'formats'
792 0 : ierr = 0
793 : if (s% mixing_length_alpha == 0d0 .or. &
794 0 : k <= s% RSP2_num_outermost_cells_forced_nonturbulent .or. &
795 : k > s% nz - int(s% nz/s% RSP2_nz_div_IBOTOM)) then
796 0 : Uq_face = 0d0
797 : else
798 0 : r_00 = wrap_opt_time_center_r_00(s,k)
799 0 : Chi_00 = s% Chi_ad(k) ! compute_Chi_cell(s,k,ierr)
800 0 : if (k > 1) then
801 : !Chi_m1 = shift_m1(compute_Chi_cell(s,k-1,ierr))
802 0 : Chi_m1 = shift_m1(s% Chi_ad(k-1))
803 0 : if (ierr /= 0) return
804 : else
805 0 : Chi_m1 = 0d0
806 : end if
807 0 : Uq_face = 4d0*pi*(Chi_m1 - Chi_00)/(r_00*s% dm_bar(k))
808 :
809 0 : if (k==-56) then
810 0 : write(*,3) 'RSP2 Uq chi_m1 chi_00 r', k, s% solver_iter, &
811 0 : Uq_face%val, Chi_m1%val, Chi_00%val, r_00%val
812 : end if
813 :
814 : end if
815 : ! erg g^-1 cm^-1 = g cm^2 s^-2 g^-1 cm^-1 = cm s^-2, acceleration
816 0 : s% Uq(k) = Uq_face%val
817 0 : end function compute_Uq_face
818 :
819 :
820 0 : function compute_Source(s, k, ierr) result(Source) ! erg g^-1 s^-1
821 : type (star_info), pointer :: s
822 : integer, intent(in) :: k
823 : type(auto_diff_real_star_order1) :: Source
824 : ! source_div_w assumes RSP2_source_seed == 0
825 : integer, intent(out) :: ierr
826 : type(auto_diff_real_star_order1) :: &
827 : w_00, T_00, d_00, Peos_00, Cp_00, chiT_00, chiRho_00, QQ_00, &
828 : Hp_face_00, Hp_face_p1, PII_face_00, PII_face_p1, PII_div_Hp_cell, &
829 : P_QQ_div_Cp
830 : include 'formats'
831 0 : ierr = 0
832 0 : w_00 = wrap_w_00(s, k)
833 0 : T_00 = wrap_T_00(s, k)
834 0 : d_00 = wrap_d_00(s, k)
835 0 : Peos_00 = wrap_Peos_00(s, k)
836 0 : Cp_00 = wrap_Cp_00(s, k)
837 0 : chiT_00 = wrap_chiT_00(s, k)
838 0 : chiRho_00 = wrap_chiRho_00(s, k)
839 0 : QQ_00 = chiT_00/(d_00*T_00*chiRho_00)
840 :
841 0 : Hp_face_00 = wrap_Hp_00(s,k)
842 0 : PII_face_00 = s% PII_ad(k) ! compute_PII_face(s, k, ierr)
843 0 : if (ierr /= 0) return
844 :
845 0 : if (k == s% nz) then
846 0 : PII_div_Hp_cell = PII_face_00/Hp_face_00
847 : else
848 0 : Hp_face_p1 = wrap_Hp_p1(s,k)
849 0 : if (ierr /= 0) return
850 : !PII_face_p1 = shift_p1(compute_PII_face(s, k+1, ierr))
851 0 : PII_face_p1 = shift_p1(s% PII_ad(k+1))
852 0 : if (ierr /= 0) return
853 0 : PII_div_Hp_cell = 0.5d0*(PII_face_00/Hp_face_00 + PII_face_p1/Hp_face_p1)
854 : end if
855 :
856 : ! Peos_00*QQ_00/Cp_00 = grad_ad if all perfect.
857 : !grad_ad_00 = wrap_grad_ad_00(s, k)
858 0 : P_QQ_div_Cp = Peos_00*QQ_00/Cp_00 ! use this to be same as RSP
859 0 : Source = (w_00 + s% RSP2_source_seed)*PII_div_Hp_cell*T_00*P_QQ_div_Cp
860 :
861 : ! PII units same as Cp = erg g^-1 K^-1
862 : ! P*QQ/Cp is unitless (see Y_face)
863 : ! Source units = (erg g^-1 K^-1) cm^-1 cm s^-1 K
864 : ! = erg g^-1 s^-1
865 :
866 0 : if (k==-109) then
867 0 : write(*,3) 'RSP2 Source w PII_div_Hp T_P_QQ_div_Cp', k, s% solver_iter, &
868 0 : Source%val, w_00%val, PII_div_Hp_cell%val, T_00%val*P_QQ_div_Cp% val
869 : !write(*,3) 'RSP2 PII_00 PII_p1 Hp_00 Hp_p1', k, s% solver_iter, &
870 : ! PII_face_00%val, PII_face_p1%val, Hp_face_00%val, Hp_face_p1%val
871 : end if
872 0 : s% SOURCE(k) = Source%val
873 :
874 0 : end function compute_Source
875 :
876 :
877 0 : function compute_D(s, k, ierr) result(D) ! erg g^-1 s^-1
878 : type (star_info), pointer :: s
879 : integer, intent(in) :: k
880 : type(auto_diff_real_star_order1) :: D
881 : type(auto_diff_real_star_order1) :: dw3, w_00
882 : integer, intent(out) :: ierr
883 : type(auto_diff_real_star_order1) :: Hp_cell
884 : include 'formats'
885 0 : ierr = 0
886 0 : if (s% mixing_length_alpha == 0d0) then
887 0 : D = 0d0
888 : else
889 0 : Hp_cell = wrap_Hp_cell(s,k)
890 0 : w_00 = wrap_w_00(s,k)
891 0 : dw3 = pow3(w_00) - pow3(s% RSP2_w_min_for_damping)
892 0 : D = (s% RSP2_alfad*x_CEDE/s% mixing_length_alpha)*dw3/Hp_cell
893 : ! units cm^3 s^-3 cm^-1 = cm^2 s^-3 = erg g^-1 s^-1
894 : end if
895 0 : if (k==-50) then
896 0 : write(*,3) 'RSP2 DAMP w Hp_cell dw3', k, s% solver_iter, &
897 0 : D%val, w_00%val, Hp_cell%val, dw3% val
898 : end if
899 0 : s% DAMP(k) = D%val
900 0 : end function compute_D
901 :
902 :
903 0 : function compute_Dr(s, k, ierr) result(Dr) ! erg g^-1 s^-1 = cm^2 s^-3
904 : type (star_info), pointer :: s
905 : integer, intent(in) :: k
906 : type(auto_diff_real_star_order1) :: Dr
907 : integer, intent(out) :: ierr
908 : type(auto_diff_real_star_order1) :: &
909 : w_00, T_00, d_00, Cp_00, kap_00, Hp_cell, POM2
910 : real(dp) :: gammar, alpha, POM
911 : include 'formats'
912 0 : ierr = 0
913 0 : alpha = s% mixing_length_alpha
914 0 : gammar = s% RSP2_alfar*x_GAMMAR
915 0 : if (gammar == 0d0) then
916 0 : Dr = 0d0
917 0 : s% DAMPR(k) = 0d0
918 0 : return
919 : end if
920 0 : w_00 = wrap_w_00(s,k)
921 0 : T_00 = wrap_T_00(s,k)
922 0 : d_00 = wrap_d_00(s,k)
923 0 : Cp_00 = wrap_Cp_00(s,k)
924 0 : kap_00 = wrap_kap_00(s,k)
925 0 : Hp_cell = wrap_Hp_cell(s,k)
926 0 : POM = 4d0*boltz_sigma*pow2(gammar/alpha) ! erg cm^-2 K^-4 s^-1
927 0 : POM2 = pow3(T_00)/(pow2(d_00)*Cp_00*kap_00)
928 : ! K^3 / ((g cm^-3)^2 (erg g^-1 K^-1) (cm^2 g^-1))
929 : ! K^3 / (cm^-4 erg K^-1) = K^4 cm^4 erg^-1
930 0 : Dr = get_etrb(s,k)*POM*POM2/pow2(Hp_cell)
931 : ! (erg cm^-2 K^-4 s^-1) (K^4 cm^4 erg^-1) cm^2 s^-2 cm^-2
932 : ! cm^2 s^-3 = erg g^-1 s^-1
933 0 : s% DAMPR(k) = Dr%val
934 0 : end function compute_Dr
935 :
936 :
937 0 : function compute_C(s, k, ierr) result(C) ! erg g^-1 s^-1
938 : type (star_info), pointer :: s
939 : integer, intent(in) :: k
940 : type(auto_diff_real_star_order1) :: C
941 : integer, intent(out) :: ierr
942 : type(auto_diff_real_star_order1) :: Source, D, Dr
943 : if (s% mixing_length_alpha == 0d0 .or. &
944 0 : k <= s% RSP2_num_outermost_cells_forced_nonturbulent .or. &
945 : k > s% nz - int(s% nz/s% RSP2_nz_div_IBOTOM)) then
946 0 : if (k >= 1 .and. k <= s% nz) then
947 0 : s% SOURCE(k) = 0d0
948 0 : s% DAMP(k) = 0d0
949 0 : s% DAMPR(k) = 0d0
950 0 : s% COUPL(k) = 0d0
951 0 : s% COUPL_ad(k) = 0d0
952 : end if
953 0 : C = 0d0
954 0 : return
955 : end if
956 0 : Source = compute_Source(s, k, ierr)
957 0 : if (ierr /= 0) return
958 0 : D = compute_D(s, k, ierr)
959 0 : if (ierr /= 0) return
960 0 : Dr = compute_Dr(s, k, ierr)
961 0 : if (ierr /= 0) return
962 0 : C = Source - D - Dr
963 0 : s% COUPL(k) = C%val
964 0 : s% COUPL_ad(k) = C
965 0 : end function compute_C
966 :
967 :
968 0 : function compute_L_face(s, k, ierr) result(L_face) ! erg s^-1
969 : type (star_info), pointer :: s
970 : integer, intent(in) :: k
971 : type(auto_diff_real_star_order1) :: L_face
972 : integer, intent(out) :: ierr
973 : type(auto_diff_real_star_order1) :: Lr, Lc, Lt
974 0 : call compute_L_terms(s, k, L_face, Lr, Lc, Lt, ierr)
975 0 : end function compute_L_face
976 :
977 :
978 0 : subroutine compute_L_terms(s, k, L, Lr, Lc, Lt, ierr)
979 : type (star_info), pointer, intent(in) :: s
980 : integer, intent(in) :: k
981 : type(auto_diff_real_star_order1), intent(out) :: L, Lr, Lc, Lt
982 : type(accurate_auto_diff_real_star_order1) :: L_sum
983 : integer, intent(out) :: ierr
984 : include 'formats'
985 0 : ierr = 0
986 0 : if (k > s% nz) then
987 0 : L = 0d0
988 0 : L%val = s% L_center
989 0 : Lr = 0d0
990 0 : Lc = 0d0
991 0 : Lt = 0d0
992 0 : return
993 : end if
994 0 : Lr = compute_Lr(s, k, ierr)
995 0 : if (ierr /= 0) return
996 0 : if (k == 1) then
997 0 : Lc = 0d0
998 0 : Lt = 0d0
999 : else
1000 0 : Lc = compute_Lc(s, k, ierr)
1001 0 : if (ierr /= 0) return
1002 0 : Lt = compute_Lt(s, k, ierr)
1003 0 : if (ierr /= 0) return
1004 : end if
1005 0 : L_sum = Lr
1006 0 : L_sum = L_sum + Lc
1007 0 : L_sum = L_sum + Lt
1008 0 : L = L_sum
1009 0 : s% Lr_ad(k) = Lr
1010 0 : s% Lc_ad(k) = Lc
1011 0 : s% Lt_ad(k) = Lt
1012 : end subroutine compute_L_terms
1013 :
1014 :
1015 0 : function compute_Lr(s, k, ierr) result(Lr) ! erg s^-1
1016 : type (star_info), pointer :: s
1017 : integer, intent(in) :: k
1018 : type(auto_diff_real_star_order1) :: Lr
1019 : integer, intent(out) :: ierr
1020 : type(auto_diff_real_star_order1) :: &
1021 : r_00, area, T_00, T400, Erad, T_m1, T4m1, &
1022 : kap_00, kap_m1, kap_face, diff_T4_div_kap, BW, BK
1023 : real(dp) :: alfa
1024 : include 'formats'
1025 0 : ierr = 0
1026 0 : if (k > s% nz) then
1027 0 : Lr = s% L_center
1028 : else
1029 0 : r_00 = wrap_r_00(s,k) ! not time centered
1030 0 : area = 4d0*pi*pow2(r_00)
1031 0 : T_00 = wrap_T_00(s,k)
1032 0 : T400 = pow4(T_00)
1033 0 : if (k == 1) then ! Lr(1) proportional to Erad in cell(1)
1034 0 : Erad = crad * T400
1035 0 : Lr = s% RSP2_Lsurf_factor * area * clight * Erad
1036 0 : s% Lr(k) = Lr%val
1037 0 : return
1038 : end if
1039 0 : T_m1 = wrap_T_m1(s,k)
1040 0 : T4m1 = pow4(T_m1)
1041 0 : alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k))
1042 0 : kap_00 = wrap_kap_00(s,k)
1043 0 : kap_m1 = wrap_kap_m1(s,k)
1044 0 : kap_face = alfa*kap_00 + (1d0 - alfa)*kap_m1
1045 0 : diff_T4_div_kap = (T4m1 - T400)/kap_face
1046 :
1047 0 : if (s% RSP2_use_Stellingwerf_Lr) then ! RSP style
1048 0 : BW = log(T4m1/T400)
1049 0 : if (abs(BW%val) > 1d-20) then
1050 0 : BK = log(kap_m1/kap_00)
1051 0 : if (abs(1d0 - BK%val/BW%val) > 1d-15 .and. abs(BW%val - BK%val) > 1d-15) then
1052 0 : diff_T4_div_kap = (T4m1/kap_m1 - T400/kap_00)/(1d0 - BK/BW)
1053 : end if
1054 : end if
1055 : end if
1056 0 : Lr = -crad*clight/3d0*diff_T4_div_kap*pow2(area)/s% dm_bar(k)
1057 : ! units (erg cm^-3 K^-4) (cm s^-1) (K^4 cm^-2 g cm^4) g^-1 = erg s^-1
1058 :
1059 : !s% xtra1_array(k) = s% T_start(k)
1060 : !s% xtra2_array(k) = T4m1%val - T400%val
1061 : !s% xtra3_array(k) = kap_face%val
1062 : !s% xtra4_array(k) = diff_T4_div_kap%val
1063 : !s% xtra5_array(k) = Lr%val/Lsun
1064 : !s% xtra6_array(k) = 1
1065 :
1066 : end if
1067 0 : s% Lr(k) = Lr%val
1068 0 : end function compute_Lr
1069 :
1070 :
1071 0 : function compute_Lc(s, k, ierr) result(Lc) ! erg s^-1
1072 : type (star_info), pointer :: s
1073 : integer, intent(in) :: k
1074 : type(auto_diff_real_star_order1) :: Lc
1075 : integer, intent(out) :: ierr
1076 : type(auto_diff_real_star_order1) :: Lc_div_w_face
1077 0 : Lc = compute_Lc_terms(s, k, Lc_div_w_face, ierr)
1078 0 : s% Lc(k) = Lc%val
1079 0 : end function compute_Lc
1080 :
1081 :
1082 0 : function compute_Lc_terms(s, k, Lc_div_w_face, ierr) result(Lc)
1083 : type (star_info), pointer :: s
1084 : integer, intent(in) :: k
1085 : type(auto_diff_real_star_order1) :: Lc, Lc_div_w_face
1086 : integer, intent(out) :: ierr
1087 : type(auto_diff_real_star_order1) :: r_00, area, &
1088 : T_m1, T_00, d_m1, d_00, w_m1, w_00, T_rho_face, PII_face, w_face
1089 : real(dp) :: ALFAC, ALFAS, alfa, beta
1090 : include 'formats'
1091 0 : ierr = 0
1092 : if (s% mixing_length_alpha == 0d0 .or. &
1093 0 : k <= s% RSP2_num_outermost_cells_forced_nonturbulent .or. &
1094 : k > s% nz - int(s% nz/s% RSP2_nz_div_IBOTOM)) then
1095 0 : Lc = 0d0
1096 0 : Lc_div_w_face = 1
1097 0 : return
1098 : end if
1099 0 : r_00 = wrap_r_00(s, k)
1100 0 : area = 4d0*pi*pow2(r_00)
1101 0 : T_m1 = wrap_T_m1(s, k)
1102 0 : T_00 = wrap_T_00(s, k)
1103 0 : d_m1 = wrap_d_m1(s, k)
1104 0 : d_00 = wrap_d_00(s, k)
1105 0 : w_m1 = wrap_w_m1(s, k)
1106 0 : w_00 = wrap_w_00(s, k)
1107 0 : call get_RSP2_alfa_beta_face_weights(s, k, alfa, beta)
1108 0 : T_rho_face = alfa*T_00*d_00 + beta*T_m1*d_m1
1109 0 : PII_face = s% PII_ad(k) ! compute_PII_face(s, k, ierr)
1110 0 : w_face = alfa*w_00 + beta*w_m1
1111 0 : ALFAC = x_ALFAC
1112 0 : ALFAS = x_ALFAS
1113 0 : Lc_div_w_face = area*(ALFAC/ALFAS)*T_rho_face*PII_face
1114 : ! units = cm^2 K g cm^-3 ergs g^-1 K^-1 = ergs cm^-1
1115 0 : Lc = w_face*Lc_div_w_face
1116 : ! units = cm s^-1 ergs cm^-1 = ergs s^-1
1117 0 : if (k == -458) then
1118 0 : write(*,2) 'Lc%val', k, Lc%val
1119 0 : write(*,2) 'w_face%val', k, w_face%val
1120 0 : write(*,2) 'Lc_div_w_face', k, Lc_div_w_face%val
1121 0 : write(*,2) 'PII_face%val', k, PII_face%val
1122 0 : write(*,2) 'T_rho_face%val', k, T_rho_face%val
1123 : !write(*,2) '', k,
1124 : !write(*,2) '', k,
1125 0 : call mesa_error(__FILE__,__LINE__,'compute_Lc_terms')
1126 : end if
1127 0 : end function compute_Lc_terms
1128 :
1129 :
1130 0 : function compute_Lt(s, k, ierr) result(Lt) ! erg s^-1
1131 : type (star_info), pointer :: s
1132 : integer, intent(in) :: k
1133 : type(auto_diff_real_star_order1) :: Lt
1134 : integer, intent(out) :: ierr
1135 : type(auto_diff_real_star_order1) :: r_00, area2, d_m1, d_00, &
1136 : rho2_face, Hp_face, w_m1, w_00, w_face, etrb_m1, etrb_00
1137 : real(dp) :: alpha_alpha_t, alfa, beta
1138 : include 'formats'
1139 0 : ierr = 0
1140 0 : if (k > s% nz) then
1141 0 : Lt = 0d0
1142 0 : return
1143 : end if
1144 0 : alpha_alpha_t = s% mixing_length_alpha*s% RSP2_alfat
1145 : if (alpha_alpha_t == 0d0 .or. &
1146 0 : k <= s% RSP2_num_outermost_cells_forced_nonturbulent .or. &
1147 : k > s% nz - int(s% nz/s% RSP2_nz_div_IBOTOM)) then
1148 0 : Lt = 0d0
1149 0 : s% Lt(k) = 0d0
1150 0 : return
1151 : end if
1152 0 : r_00 = wrap_r_00(s,k)
1153 0 : area2 = pow2(4d0*pi*pow2(r_00))
1154 0 : d_m1 = wrap_d_m1(s,k)
1155 0 : d_00 = wrap_d_00(s,k)
1156 0 : call get_RSP2_alfa_beta_face_weights(s, k, alfa, beta)
1157 0 : rho2_face = alfa*pow2(d_00) + beta*pow2(d_m1)
1158 0 : w_m1 = wrap_w_m1(s,k)
1159 0 : w_00 = wrap_w_00(s,k)
1160 0 : w_face = alfa*w_00 + beta*w_m1
1161 0 : etrb_m1 = wrap_etrb_m1(s,k)
1162 0 : etrb_00 = wrap_etrb_00(s,k)
1163 0 : Hp_face = wrap_Hp_00(s,k)
1164 : ! Ft = - alpha_t * rho_face * alpha * Hp_face * w_face * detrb/dr (thesis eqn 2.44)
1165 : ! replace dr by dm_bar/(area*rho_face)
1166 : ! Ft = - alpha_alpha_t * rho_face * Hp_face * w_face * (area*rho_face) * detrb/dm_bar
1167 : ! Lt = area * Ft
1168 : ! Lt = -alpha_alpha_t * (area*rho_face)**2 * Hp_face * w_face * (etrb(k-1) - etrb(k))/dm_bar
1169 0 : Lt = - alpha_alpha_t * area2 * rho2_face * Hp_face * w_face * (etrb_m1 - etrb_00) / s% dm_bar(k)
1170 : ! units = (cm^4) (g^2 cm^-6) (cm) (cm s^-1) (ergs g^-1) g^-1 = erg s^-1
1171 0 : s% Lt(k) = Lt%val
1172 0 : end function compute_Lt
1173 :
1174 :
1175 0 : subroutine set_etrb_start_vars(s, ierr)
1176 : type (star_info), pointer :: s
1177 : integer, intent(out) :: ierr
1178 : integer :: k
1179 : type(auto_diff_real_star_order1) :: Y_face, Lt
1180 : include 'formats'
1181 0 : ierr = 0
1182 0 : do k=1,s%nz
1183 0 : Y_face = compute_Y_face(s, k, ierr)
1184 0 : if (ierr /= 0) return
1185 0 : s% Y_face_start(k) = Y_face%val
1186 0 : Lt = compute_Lt(s, k, ierr)
1187 0 : if (ierr /= 0) return
1188 0 : s% Lt_start(k) = Lt%val
1189 0 : s% w_start(k) = s% w(k)
1190 0 : s% Hp_face_start(k) = s% Hp_face(k)
1191 : end do
1192 : end subroutine set_etrb_start_vars
1193 :
1194 :
1195 0 : subroutine RSP2_adjust_vars_before_call_solver(s, ierr) ! replaces check_omega in RSP
1196 : ! JAK OKRESLIC OMEGA DLA PIERWSZEJ ITERACJI
1197 : use micro, only: do_eos_for_cell
1198 : type (star_info), pointer :: s
1199 : integer, intent(out) :: ierr
1200 : real(dp) :: PII_div_Hp, QQ, SOURCE, Hp_cell, DAMP, POM, POM2, DAMPR, del, soln
1201 : !type(auto_diff_real_star_order1) :: x
1202 : integer :: k
1203 : include 'formats'
1204 0 : ierr = 0
1205 0 : if (s% mixing_length_alpha == 0d0) return
1206 :
1207 0 : !$OMP PARALLEL DO PRIVATE(k,PII_div_Hp,QQ,SOURCE,Hp_cell,DAMP,POM,POM2,DAMPR,del,soln) SCHEDULE(dynamic,2)
1208 : do k=s% RSP2_num_outermost_cells_forced_nonturbulent+1, &
1209 : s% nz - max(1,int(s% nz/s% RSP_nz_div_IBOTOM))
1210 :
1211 : if (s% w(k) > s% RSP2_w_min_for_damping) cycle
1212 :
1213 : PII_div_Hp = 0.5d0*(s% PII(k)/s% Hp_face(k) + s% PII(k+1)/s% Hp_face(k+1))
1214 : QQ = s% chiT(k)/(s% rho(k)*s% T(k)*s% chiRho(k))
1215 : SOURCE = PII_div_Hp*s% T(k)*s% Peos(k)*QQ/s% Cp(k)
1216 :
1217 : Hp_cell = 0.5d0*(s% Hp_face(k) + s% Hp_face(k+1))
1218 : DAMP = (s% RSP2_alfad*x_CEDE/s% mixing_length_alpha)/Hp_cell
1219 :
1220 : POM = 4d0*boltz_sigma*pow2(s% RSP2_alfar*x_GAMMAR/s% mixing_length_alpha)
1221 : POM2 = pow3(s% T(k))/(pow2(s% rho(k))*s% Cp(k)*s% opacity(k))
1222 : DAMPR = POM*POM2/pow2(Hp_cell)
1223 :
1224 : del = pow2(DAMPR) + 4d0*DAMP*SOURCE
1225 :
1226 : if (k==-35) then
1227 : write(*,2) 'del', k, del
1228 : write(*,2) 'DAMPR', k, DAMPR
1229 : write(*,2) 'DAMP', k, DAMP
1230 : write(*,2) 'SOURCE', k, SOURCE
1231 : write(*,2) 'POM', k, PII_div_Hp
1232 : write(*,2) 'POM2', k, s% T(k)*s% Peos(k)*QQ/s% Cp(k)
1233 : write(*,2) 's% Hp_face(k)', k, s% Hp_face(k)
1234 : write(*,2) 's% Hp_face(k+1)', k+1, s% Hp_face(k+1)
1235 : write(*,2) 's% PII(k)', k, s% PII(k)
1236 : write(*,2) 's% PII(k+1)', k+1, s% PII(k+1)
1237 : write(*,2) 's% Y_face(k)', k, s% Y_face(k)
1238 : write(*,2) 's% Y_face(k+1)', k+1, s% Y_face(k+1)
1239 : end if
1240 :
1241 : if (del < 0d0) cycle
1242 : soln = (-DAMPR + sqrt(del))/(2d0*DAMP)
1243 : if (k==-35) write(*,2) 'soln', k, soln
1244 : if (soln > 0d0) then
1245 : ! i tried soln = sqrt(soln) here. helps solver convergence, but hurts the model results.
1246 : if (s% RSP2_report_adjust_w) &
1247 : write(*,3) 'RSP2_adjust_vars_before_call_solver w', k, s% model_number, s% w(k), soln
1248 : s% w(k) = soln
1249 : end if
1250 :
1251 : end do
1252 : !$OMP END PARALLEL DO
1253 : end subroutine RSP2_adjust_vars_before_call_solver
1254 :
1255 : end module hydro_rsp2
|
