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