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, T_00, T_m1
544 : type(auto_diff_real_star_order1) :: X, FL, scale, T_face, e_face, Peos_face, rho_face, h_face
545 0 : real(dp) :: ALFAS_ALFA, alfa, beta
546 : include 'formats'
547 0 : ierr = 0
548 0 : if (k > s% nz) then
549 0 : PII_face = 0d0
550 0 : return
551 : end if
552 0 : if (k == 1 .or. s% mixing_length_alpha == 0d0 .or. &
553 : k == s% nz) then ! just skip k == nz to be like RSP
554 0 : PII_face = 0d0
555 0 : s% PII(k) = 0d0
556 0 : s% PII_ad(k) = 0d0
557 0 : return
558 : end if
559 0 : Y_face = s% Y_face_ad(k) ! compute_Y_face(s, k, ierr)
560 0 : if (ierr /= 0) return
561 0 : Cp_00 = wrap_Cp_00(s, k)
562 0 : Cp_m1 = wrap_Cp_m1(s, k)
563 0 : T_00= wrap_T_00(s, k)
564 0 : T_m1 = wrap_T_m1(s, k)
565 0 : call get_RSP2_alfa_beta_face_weights(s, k, alfa, beta)
566 0 : Cp_face = alfa*Cp_00 + beta*Cp_m1 ! ergs g^-1 K^-1
567 0 : T_face = alfa*Cp_00 + beta*Cp_m1
568 0 : rho_face = alfa*wrap_d_00(s,k) + beta*wrap_d_m1(s,k)
569 0 : Peos_face = alfa*wrap_Peos_00(s,k) + beta*wrap_Peos_m1(s,k)
570 0 : e_face = alfa*wrap_e_00(s,k) + beta*wrap_e_m1(s,k)
571 0 : h_face = e_face + Peos_face/rho_face
572 0 : ALFAS_ALFA = x_ALFAS*s% mixing_length_alpha
573 0 : PII_face = ALFAS_ALFA*Cp_face*Y_face
574 :
575 0 : scale = 1d0
576 0 : if (Y_face > 0d0 .and. s% use_TDC_enthalpy_flux_limiter) then
577 : ! X = G/F
578 0 : X = (Cp_face*T_face/h_face)*ALFAS_ALFA* Y_face / sqrt_2_div_3
579 0 : FL = flux_limiter_function(X)
580 : ! Avoid 0/0 or tiny/tiny; for X ≈ 0, FL ≈ X so scale ~ 1 anyway.
581 0 : if (abs(X%val) >= 0.95d0) then
582 0 : scale = FL / X
583 : else
584 0 : scale = 1d0
585 : end if
586 : end if
587 :
588 0 : s% PII(k) = PII_face%val*scale%val
589 0 : s% PII_ad(k) = PII_face*scale
590 0 : if (k == -2 .and. s% PII(k) < 0d0) then
591 0 : write(*,2) 's% PII(k)', k, s% PII(k)
592 0 : write(*,2) 'Cp_face', k, Cp_face%val
593 0 : write(*,2) 'Y_face', k, Y_face%val
594 : !write(*,2) 'PII_face%val', k, PII_face%val
595 : !write(*,2) 'T_rho_face%val', k, T_rho_face%val
596 : !write(*,2) '', k,
597 : !write(*,2) '', k,
598 0 : call mesa_error(__FILE__,__LINE__,'compute_PII_face')
599 : end if
600 0 : end function compute_PII_face
601 :
602 0 : type(auto_diff_real_star_order1) function flux_limiter_function(X) result(FL) ! should be c2 continuous
603 : type(auto_diff_real_star_order1), intent(in) :: X
604 : real(dp), parameter :: X0 = 0.95_dp ! start of transition
605 : real(dp), parameter :: delta = 0.05_dp ! width of transition
606 : real(dp), parameter :: X1 = 1d0 !X0 + delta ! end of transition
607 :
608 : type(auto_diff_real_star_order1) :: s, p
609 :
610 : ! Region 1: purely linear, FL = X
611 0 : if (X%val < X0) then ! should not be encountered
612 0 : FL = X
613 :
614 : ! Region 3: saturated, FL = 1
615 0 : else if (X%val >= X1) then
616 0 : FL = 1.0_dp
617 :
618 : ! Region 2: smooth C² transition between the two
619 : else
620 : ! Normalized coordinate in [0,1]
621 0 : s = (X - X0) / (X1 - X0)
622 :
623 : ! Quintic "smootherstep" polynomial:
624 : ! p(s) = 10 s^3 - 15 s^4 + 6 s^5
625 : ! p(0)=0, p(1)=1, p'(0)=p'(1)=0, p''(0)=p''(1)=0
626 0 : p = pow3(s) * (10.0_dp + s * (-15.0_dp + 6.0_dp * s))
627 :
628 : ! Blend between line FL=X and flat FL=1
629 : ! At s=0: FL = X
630 : ! At s=1: FL = 1
631 : ! Because p', p'' vanish at 0 and 1, FL, FL', FL'' all match.
632 0 : FL = X + (1.0_dp - X) * p
633 : end if
634 0 : end function flux_limiter_function
635 :
636 0 : function compute_d_v_div_r(s, k, ierr) result(d_v_div_r) ! s^-1
637 : type (star_info), pointer :: s
638 : integer, intent(in) :: k
639 : type(auto_diff_real_star_order1) :: d_v_div_r
640 : integer, intent(out) :: ierr
641 : type(auto_diff_real_star_order1) :: v_00, v_p1, r_00, r_p1
642 : include 'formats'
643 0 : ierr = 0
644 0 : v_00 = wrap_v_00(s,k)
645 0 : v_p1 = wrap_v_p1(s,k)
646 0 : r_00 = wrap_r_00(s,k)
647 0 : r_p1 = wrap_r_p1(s,k)
648 0 : if (r_p1%val == 0d0) r_p1 = 1d0
649 0 : d_v_div_r = v_00/r_00 - v_p1/r_p1 ! units s^-1
650 0 : end function compute_d_v_div_r
651 :
652 :
653 0 : function compute_d_v_div_r_opt_time_center(s, k, ierr) result(d_v_div_r) ! s^-1
654 : type (star_info), pointer :: s
655 : integer, intent(in) :: k
656 : type(auto_diff_real_star_order1) :: d_v_div_r
657 : integer, intent(out) :: ierr
658 : type(auto_diff_real_star_order1) :: v_00, v_p1, r_00, r_p1
659 : include 'formats'
660 0 : ierr = 0
661 0 : v_00 = wrap_opt_time_center_v_00(s,k)
662 0 : v_p1 = wrap_opt_time_center_v_p1(s,k)
663 0 : r_00 = wrap_opt_time_center_r_00(s,k)
664 0 : r_p1 = wrap_opt_time_center_r_p1(s,k)
665 0 : if (r_p1%val == 0d0) r_p1 = 1d0
666 0 : d_v_div_r = v_00/r_00 - v_p1/r_p1 ! units s^-1
667 0 : end function compute_d_v_div_r_opt_time_center
668 :
669 :
670 0 : function wrap_Hp_cell(s, k) result(Hp_cell) ! cm
671 : type (star_info), pointer :: s
672 : integer, intent(in) :: k
673 : type(auto_diff_real_star_order1) :: Hp_cell
674 0 : Hp_cell = 0.5d0*(wrap_Hp_00(s,k) + wrap_Hp_p1(s,k))
675 0 : end function wrap_Hp_cell
676 :
677 :
678 0 : function Hp_cell_for_Chi(s, k, ierr) result(Hp_cell) ! cm
679 : type (star_info), pointer :: s
680 : integer, intent(in) :: k
681 : integer, intent(out) :: ierr
682 : type(auto_diff_real_star_order1) :: Hp_cell
683 : type(auto_diff_real_star_order1) :: d_00, Peos_00, rmid
684 0 : real(dp) :: mmid, cgrav_mid
685 : include 'formats'
686 0 : ierr = 0
687 :
688 0 : Hp_cell = wrap_Hp_cell(s, k)
689 0 : return
690 :
691 : d_00 = wrap_d_00(s, k)
692 : Peos_00 = wrap_Peos_00(s, k)
693 : if (k < s% nz) then
694 : rmid = 0.5d0*(wrap_r_00(s,k) + wrap_r_p1(s,k))
695 : mmid = 0.5d0*(s% m(k) + s% m(k+1))
696 : cgrav_mid = 0.5d0*(s% cgrav(k) + s% cgrav(k+1))
697 : else
698 : rmid = 0.5d0*(wrap_r_00(s,k) + s% r_center)
699 : mmid = 0.5d0*(s% m(k) + s% m_center)
700 : cgrav_mid = s% cgrav(k)
701 : end if
702 : Hp_cell = pow2(rmid)*Peos_00/(d_00*cgrav_mid*mmid)
703 : if (s% alt_scale_height_flag) then
704 : call mesa_error(__FILE__,__LINE__,'Hp_cell_for_Chi: cannot use alt_scale_height_flag')
705 : end if
706 : end function Hp_cell_for_Chi
707 :
708 :
709 0 : function compute_Chi_cell(s, k, ierr) result(Chi_cell)
710 : ! eddy viscosity energy (Kuhfuss 1986) [erg]
711 : type (star_info), pointer :: s
712 : integer, intent(in) :: k
713 : type(auto_diff_real_star_order1) :: Chi_cell
714 : integer, intent(out) :: ierr
715 : type(auto_diff_real_star_order1) :: &
716 : rho2, r6_cell, d_v_div_r, Hp_cell, w_00, d_00, r_00, r_p1
717 0 : real(dp) :: f, ALFAM_ALFA
718 : include 'formats'
719 0 : ierr = 0
720 0 : ALFAM_ALFA = s% RSP2_alfam*s% mixing_length_alpha
721 : if (ALFAM_ALFA == 0d0 .or. &
722 0 : k <= s% RSP2_num_outermost_cells_forced_nonturbulent .or. &
723 : k > s% nz - int(s% nz/s% RSP2_nz_div_IBOTOM)) then
724 0 : Chi_cell = 0d0
725 0 : if (k >= 1 .and. k <= s% nz) then
726 0 : s% Chi(k) = 0d0
727 0 : s% Chi_ad(k) = 0d0
728 : end if
729 : else
730 0 : Hp_cell = Hp_cell_for_Chi(s, k, ierr)
731 0 : if (ierr /= 0) return
732 0 : d_v_div_r = compute_d_v_div_r(s, k, ierr)
733 0 : if (ierr /= 0) return
734 0 : w_00 = wrap_w_00(s,k)
735 0 : d_00 = wrap_d_00(s,k)
736 0 : f = (16d0/3d0)*pi*ALFAM_ALFA/s% dm(k)
737 0 : rho2 = pow2(d_00)
738 0 : r_00 = wrap_r_00(s,k)
739 0 : r_p1 = wrap_r_p1(s,k)
740 0 : r6_cell = 0.5d0*(pow6(r_00) + pow6(r_p1))
741 0 : Chi_cell = f*rho2*r6_cell*d_v_div_r*Hp_cell*w_00
742 : ! units = g^-1 cm s^-1 g^2 cm^-6 cm^6 s^-1 cm
743 : ! = g cm^2 s^-2
744 : ! = erg
745 : end if
746 0 : s% Chi(k) = Chi_cell%val
747 0 : s% Chi_ad(k) = Chi_cell
748 :
749 0 : end function compute_Chi_cell
750 :
751 :
752 0 : function compute_Eq_cell(s, k, ierr) result(Eq_cell) ! erg g^-1 s^-1
753 : type (star_info), pointer :: s
754 : integer, intent(in) :: k
755 : type(auto_diff_real_star_order1) :: Eq_cell
756 : integer, intent(out) :: ierr
757 : type(auto_diff_real_star_order1) :: d_v_div_r, Chi_cell
758 : include 'formats'
759 0 : ierr = 0
760 : if (s% mixing_length_alpha == 0d0 .or. &
761 0 : k <= s% RSP2_num_outermost_cells_forced_nonturbulent .or. &
762 : k > s% nz - int(s% nz/s% RSP2_nz_div_IBOTOM)) then
763 0 : Eq_cell = 0d0
764 0 : if (k >= 1 .and. k <= s% nz) s% Eq_ad(k) = 0d0
765 : else
766 0 : Chi_cell = s% Chi_ad(k) ! compute_Chi_cell(s,k,ierr)
767 0 : if (ierr /= 0) return
768 0 : d_v_div_r = compute_d_v_div_r_opt_time_center(s, k, ierr)
769 0 : if (ierr /= 0) return
770 0 : Eq_cell = 4d0*pi*Chi_cell*d_v_div_r/s% dm(k) ! erg s^-1 g^-1
771 : end if
772 0 : s% Eq(k) = Eq_cell%val
773 0 : s% Eq_ad(k) = Eq_cell
774 0 : end function compute_Eq_cell
775 :
776 :
777 0 : function compute_Uq_face(s, k, ierr) result(Uq_face) ! cm s^-2, acceleration
778 : type (star_info), pointer :: s
779 : integer, intent(in) :: k
780 : type(auto_diff_real_star_order1) :: Uq_face
781 : integer, intent(out) :: ierr
782 : type(auto_diff_real_star_order1) :: Chi_00, Chi_m1, r_00
783 : include 'formats'
784 0 : ierr = 0
785 : if (s% mixing_length_alpha == 0d0 .or. &
786 0 : k <= s% RSP2_num_outermost_cells_forced_nonturbulent .or. &
787 : k > s% nz - int(s% nz/s% RSP2_nz_div_IBOTOM)) then
788 0 : Uq_face = 0d0
789 : else
790 0 : r_00 = wrap_opt_time_center_r_00(s,k)
791 0 : Chi_00 = s% Chi_ad(k) ! compute_Chi_cell(s,k,ierr)
792 0 : if (k > 1) then
793 : !Chi_m1 = shift_m1(compute_Chi_cell(s,k-1,ierr))
794 0 : Chi_m1 = shift_m1(s% Chi_ad(k-1))
795 0 : if (ierr /= 0) return
796 : else
797 0 : Chi_m1 = 0d0
798 : end if
799 0 : Uq_face = 4d0*pi*(Chi_m1 - Chi_00)/(r_00*s% dm_bar(k))
800 :
801 0 : if (k==-56) then
802 0 : write(*,3) 'RSP2 Uq chi_m1 chi_00 r', k, s% solver_iter, &
803 0 : Uq_face%val, Chi_m1%val, Chi_00%val, r_00%val
804 : end if
805 :
806 : end if
807 : ! erg g^-1 cm^-1 = g cm^2 s^-2 g^-1 cm^-1 = cm s^-2, acceleration
808 0 : s% Uq(k) = Uq_face%val
809 0 : end function compute_Uq_face
810 :
811 :
812 0 : function compute_Source(s, k, ierr) result(Source) ! erg g^-1 s^-1
813 : type (star_info), pointer :: s
814 : integer, intent(in) :: k
815 : type(auto_diff_real_star_order1) :: Source
816 : ! source_div_w assumes RSP2_source_seed == 0
817 : integer, intent(out) :: ierr
818 : type(auto_diff_real_star_order1) :: &
819 : w_00, T_00, d_00, Peos_00, Cp_00, chiT_00, chiRho_00, QQ_00, &
820 : Hp_face_00, Hp_face_p1, PII_face_00, PII_face_p1, PII_div_Hp_cell, &
821 : P_QQ_div_Cp
822 : include 'formats'
823 0 : ierr = 0
824 0 : w_00 = wrap_w_00(s, k)
825 0 : T_00 = wrap_T_00(s, k)
826 0 : d_00 = wrap_d_00(s, k)
827 0 : Peos_00 = wrap_Peos_00(s, k)
828 0 : Cp_00 = wrap_Cp_00(s, k)
829 0 : chiT_00 = wrap_chiT_00(s, k)
830 0 : chiRho_00 = wrap_chiRho_00(s, k)
831 0 : QQ_00 = chiT_00/(d_00*T_00*chiRho_00)
832 :
833 0 : Hp_face_00 = wrap_Hp_00(s,k)
834 0 : PII_face_00 = s% PII_ad(k) ! compute_PII_face(s, k, ierr)
835 0 : if (ierr /= 0) return
836 :
837 0 : if (k == s% nz) then
838 0 : PII_div_Hp_cell = PII_face_00/Hp_face_00
839 : else
840 0 : Hp_face_p1 = wrap_Hp_p1(s,k)
841 0 : if (ierr /= 0) return
842 : !PII_face_p1 = shift_p1(compute_PII_face(s, k+1, ierr))
843 0 : PII_face_p1 = shift_p1(s% PII_ad(k+1))
844 0 : if (ierr /= 0) return
845 0 : PII_div_Hp_cell = 0.5d0*(PII_face_00/Hp_face_00 + PII_face_p1/Hp_face_p1)
846 : end if
847 :
848 : ! Peos_00*QQ_00/Cp_00 = grad_ad if all perfect.
849 : !grad_ad_00 = wrap_grad_ad_00(s, k)
850 0 : P_QQ_div_Cp = Peos_00*QQ_00/Cp_00 ! use this to be same as RSP
851 0 : Source = (w_00 + s% RSP2_source_seed)*PII_div_Hp_cell*T_00*P_QQ_div_Cp
852 :
853 : ! PII units same as Cp = erg g^-1 K^-1
854 : ! P*QQ/Cp is unitless (see Y_face)
855 : ! Source units = (erg g^-1 K^-1) cm^-1 cm s^-1 K
856 : ! = erg g^-1 s^-1
857 :
858 0 : if (k==-109) then
859 0 : write(*,3) 'RSP2 Source w PII_div_Hp T_P_QQ_div_Cp', k, s% solver_iter, &
860 0 : Source%val, w_00%val, PII_div_Hp_cell%val, T_00%val*P_QQ_div_Cp% val
861 : !write(*,3) 'RSP2 PII_00 PII_p1 Hp_00 Hp_p1', k, s% solver_iter, &
862 : ! PII_face_00%val, PII_face_p1%val, Hp_face_00%val, Hp_face_p1%val
863 : end if
864 0 : s% SOURCE(k) = Source%val
865 :
866 0 : end function compute_Source
867 :
868 :
869 0 : function compute_D(s, k, ierr) result(D) ! erg g^-1 s^-1
870 : type (star_info), pointer :: s
871 : integer, intent(in) :: k
872 : type(auto_diff_real_star_order1) :: D
873 : type(auto_diff_real_star_order1) :: dw3, w_00
874 : integer, intent(out) :: ierr
875 : type(auto_diff_real_star_order1) :: Hp_cell
876 : include 'formats'
877 0 : ierr = 0
878 0 : if (s% mixing_length_alpha == 0d0) then
879 0 : D = 0d0
880 : else
881 0 : Hp_cell = wrap_Hp_cell(s,k)
882 0 : w_00 = wrap_w_00(s,k)
883 0 : dw3 = pow3(w_00) - pow3(s% RSP2_w_min_for_damping)
884 0 : D = (s% RSP2_alfad*x_CEDE/s% mixing_length_alpha)*dw3/Hp_cell
885 : ! units cm^3 s^-3 cm^-1 = cm^2 s^-3 = erg g^-1 s^-1
886 : end if
887 0 : if (k==-50) then
888 0 : write(*,3) 'RSP2 DAMP w Hp_cell dw3', k, s% solver_iter, &
889 0 : D%val, w_00%val, Hp_cell%val, dw3% val
890 : end if
891 0 : s% DAMP(k) = D%val
892 0 : end function compute_D
893 :
894 :
895 0 : function compute_Dr(s, k, ierr) result(Dr) ! erg g^-1 s^-1 = cm^2 s^-3
896 : type (star_info), pointer :: s
897 : integer, intent(in) :: k
898 : type(auto_diff_real_star_order1) :: Dr
899 : integer, intent(out) :: ierr
900 : type(auto_diff_real_star_order1) :: &
901 : w_00, T_00, d_00, Cp_00, kap_00, Hp_cell, POM2
902 0 : real(dp) :: gammar, alpha, POM
903 : include 'formats'
904 0 : ierr = 0
905 0 : alpha = s% mixing_length_alpha
906 0 : gammar = s% RSP2_alfar*x_GAMMAR
907 0 : if (gammar == 0d0) then
908 0 : Dr = 0d0
909 0 : s% DAMPR(k) = 0d0
910 0 : return
911 : end if
912 0 : w_00 = wrap_w_00(s,k)
913 0 : T_00 = wrap_T_00(s,k)
914 0 : d_00 = wrap_d_00(s,k)
915 0 : Cp_00 = wrap_Cp_00(s,k)
916 0 : kap_00 = wrap_kap_00(s,k)
917 0 : Hp_cell = wrap_Hp_cell(s,k)
918 0 : POM = 4d0*boltz_sigma*pow2(gammar/alpha) ! erg cm^-2 K^-4 s^-1
919 0 : POM2 = pow3(T_00)/(pow2(d_00)*Cp_00*kap_00)
920 : ! K^3 / ((g cm^-3)^2 (erg g^-1 K^-1) (cm^2 g^-1))
921 : ! K^3 / (cm^-4 erg K^-1) = K^4 cm^4 erg^-1
922 0 : Dr = get_etrb(s,k)*POM*POM2/pow2(Hp_cell)
923 : ! (erg cm^-2 K^-4 s^-1) (K^4 cm^4 erg^-1) cm^2 s^-2 cm^-2
924 : ! cm^2 s^-3 = erg g^-1 s^-1
925 0 : s% DAMPR(k) = Dr%val
926 0 : end function compute_Dr
927 :
928 :
929 0 : function compute_C(s, k, ierr) result(C) ! erg g^-1 s^-1
930 : type (star_info), pointer :: s
931 : integer, intent(in) :: k
932 : type(auto_diff_real_star_order1) :: C
933 : integer, intent(out) :: ierr
934 : type(auto_diff_real_star_order1) :: Source, D, Dr
935 : if (s% mixing_length_alpha == 0d0 .or. &
936 0 : k <= s% RSP2_num_outermost_cells_forced_nonturbulent .or. &
937 : k > s% nz - int(s% nz/s% RSP2_nz_div_IBOTOM)) then
938 0 : if (k >= 1 .and. k <= s% nz) then
939 0 : s% SOURCE(k) = 0d0
940 0 : s% DAMP(k) = 0d0
941 0 : s% DAMPR(k) = 0d0
942 0 : s% COUPL(k) = 0d0
943 0 : s% COUPL_ad(k) = 0d0
944 : end if
945 0 : C = 0d0
946 0 : return
947 : end if
948 0 : Source = compute_Source(s, k, ierr)
949 0 : if (ierr /= 0) return
950 0 : D = compute_D(s, k, ierr)
951 0 : if (ierr /= 0) return
952 0 : Dr = compute_Dr(s, k, ierr)
953 0 : if (ierr /= 0) return
954 0 : C = Source - D - Dr
955 0 : s% COUPL(k) = C%val
956 0 : s% COUPL_ad(k) = C
957 0 : end function compute_C
958 :
959 :
960 0 : function compute_L_face(s, k, ierr) result(L_face) ! erg s^-1
961 : type (star_info), pointer :: s
962 : integer, intent(in) :: k
963 : type(auto_diff_real_star_order1) :: L_face
964 : integer, intent(out) :: ierr
965 : type(auto_diff_real_star_order1) :: Lr, Lc, Lt
966 0 : call compute_L_terms(s, k, L_face, Lr, Lc, Lt, ierr)
967 0 : end function compute_L_face
968 :
969 :
970 0 : subroutine compute_L_terms(s, k, L, Lr, Lc, Lt, ierr)
971 : type (star_info), pointer, intent(in) :: s
972 : integer, intent(in) :: k
973 : type(auto_diff_real_star_order1), intent(out) :: L, Lr, Lc, Lt
974 : integer, intent(out) :: ierr
975 : include 'formats'
976 0 : ierr = 0
977 0 : if (k > s% nz) then
978 0 : L = 0d0
979 0 : L%val = s% L_center
980 0 : Lr = 0d0
981 0 : Lc = 0d0
982 0 : Lt = 0d0
983 0 : return
984 : end if
985 0 : Lr = compute_Lr(s, k, ierr)
986 0 : if (ierr /= 0) return
987 0 : if (k == 1) then
988 0 : Lc = 0d0
989 0 : Lt = 0d0
990 : else
991 0 : Lc = compute_Lc(s, k, ierr)
992 0 : if (ierr /= 0) return
993 0 : Lt = compute_Lt(s, k, ierr)
994 0 : if (ierr /= 0) return
995 : end if
996 0 : L = Lr + Lc + Lt
997 0 : s% Lr_ad(k) = Lr
998 0 : s% Lc_ad(k) = Lc
999 0 : s% Lt_ad(k) = Lt
1000 : end subroutine compute_L_terms
1001 :
1002 :
1003 0 : function compute_Lr(s, k, ierr) result(Lr) ! erg s^-1
1004 : type (star_info), pointer :: s
1005 : integer, intent(in) :: k
1006 : type(auto_diff_real_star_order1) :: Lr
1007 : integer, intent(out) :: ierr
1008 : type(auto_diff_real_star_order1) :: &
1009 : r_00, area, T_00, T400, Erad, T_m1, T4m1, &
1010 : kap_00, kap_m1, kap_face, diff_T4_div_kap, BW, BK
1011 0 : real(dp) :: alfa
1012 : include 'formats'
1013 0 : ierr = 0
1014 0 : if (k > s% nz) then
1015 0 : Lr = s% L_center
1016 : else
1017 0 : r_00 = wrap_r_00(s,k) ! not time centered
1018 0 : area = 4d0*pi*pow2(r_00)
1019 0 : T_00 = wrap_T_00(s,k)
1020 0 : T400 = pow4(T_00)
1021 0 : if (k == 1) then ! Lr(1) proportional to Erad in cell(1)
1022 0 : Erad = crad * T400
1023 0 : Lr = s% RSP2_Lsurf_factor * area * clight * Erad
1024 0 : s% Lr(k) = Lr%val
1025 0 : return
1026 : end if
1027 0 : T_m1 = wrap_T_m1(s,k)
1028 0 : T4m1 = pow4(T_m1)
1029 0 : alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k))
1030 0 : kap_00 = wrap_kap_00(s,k)
1031 0 : kap_m1 = wrap_kap_m1(s,k)
1032 0 : kap_face = alfa*kap_00 + (1d0 - alfa)*kap_m1
1033 0 : diff_T4_div_kap = (T4m1 - T400)/kap_face
1034 :
1035 0 : if (s% RSP2_use_Stellingwerf_Lr) then ! RSP style
1036 0 : BW = log(T4m1/T400)
1037 0 : if (abs(BW%val) > 1d-20) then
1038 0 : BK = log(kap_m1/kap_00)
1039 0 : if (abs(1d0 - BK%val/BW%val) > 1d-15 .and. abs(BW%val - BK%val) > 1d-15) then
1040 0 : diff_T4_div_kap = (T4m1/kap_m1 - T400/kap_00)/(1d0 - BK/BW)
1041 : end if
1042 : end if
1043 : end if
1044 0 : Lr = -crad*clight/3d0*diff_T4_div_kap*pow2(area)/s% dm_bar(k)
1045 : ! units (erg cm^-3 K^-4) (cm s^-1) (K^4 cm^-2 g cm^4) g^-1 = erg s^-1
1046 :
1047 : !s% xtra1_array(k) = s% T_start(k)
1048 : !s% xtra2_array(k) = T4m1%val - T400%val
1049 : !s% xtra3_array(k) = kap_face%val
1050 : !s% xtra4_array(k) = diff_T4_div_kap%val
1051 : !s% xtra5_array(k) = Lr%val/Lsun
1052 : !s% xtra6_array(k) = 1
1053 :
1054 : end if
1055 0 : s% Lr(k) = Lr%val
1056 0 : end function compute_Lr
1057 :
1058 :
1059 0 : function compute_Lc(s, k, ierr) result(Lc) ! erg s^-1
1060 : type (star_info), pointer :: s
1061 : integer, intent(in) :: k
1062 : type(auto_diff_real_star_order1) :: Lc
1063 : integer, intent(out) :: ierr
1064 : type(auto_diff_real_star_order1) :: Lc_div_w_face
1065 0 : Lc = compute_Lc_terms(s, k, Lc_div_w_face, ierr)
1066 0 : s% Lc(k) = Lc%val
1067 0 : end function compute_Lc
1068 :
1069 :
1070 0 : function compute_Lc_terms(s, k, Lc_div_w_face, ierr) result(Lc)
1071 : type (star_info), pointer :: s
1072 : integer, intent(in) :: k
1073 : type(auto_diff_real_star_order1) :: Lc, Lc_div_w_face
1074 : integer, intent(out) :: ierr
1075 : type(auto_diff_real_star_order1) :: r_00, area, &
1076 : T_m1, T_00, d_m1, d_00, w_m1, w_00, T_rho_face, PII_face, w_face
1077 0 : real(dp) :: ALFAC, ALFAS, alfa, beta
1078 : include 'formats'
1079 0 : ierr = 0
1080 : if (s% mixing_length_alpha == 0d0 .or. &
1081 0 : k <= s% RSP2_num_outermost_cells_forced_nonturbulent .or. &
1082 : k > s% nz - int(s% nz/s% RSP2_nz_div_IBOTOM)) then
1083 0 : Lc = 0d0
1084 0 : Lc_div_w_face = 1
1085 0 : return
1086 : end if
1087 0 : r_00 = wrap_r_00(s, k)
1088 0 : area = 4d0*pi*pow2(r_00)
1089 0 : T_m1 = wrap_T_m1(s, k)
1090 0 : T_00 = wrap_T_00(s, k)
1091 0 : d_m1 = wrap_d_m1(s, k)
1092 0 : d_00 = wrap_d_00(s, k)
1093 0 : w_m1 = wrap_w_m1(s, k)
1094 0 : w_00 = wrap_w_00(s, k)
1095 0 : call get_RSP2_alfa_beta_face_weights(s, k, alfa, beta)
1096 0 : T_rho_face = alfa*T_00*d_00 + beta*T_m1*d_m1
1097 0 : PII_face = s% PII_ad(k) ! compute_PII_face(s, k, ierr)
1098 0 : w_face = alfa*w_00 + beta*w_m1
1099 0 : ALFAC = x_ALFAC
1100 0 : ALFAS = x_ALFAS
1101 0 : Lc_div_w_face = area*(ALFAC/ALFAS)*T_rho_face*PII_face
1102 : ! units = cm^2 K g cm^-3 ergs g^-1 K^-1 = ergs cm^-1
1103 0 : Lc = w_face*Lc_div_w_face
1104 : ! units = cm s^-1 ergs cm^-1 = ergs s^-1
1105 0 : if (k == -458) then
1106 0 : write(*,2) 'Lc%val', k, Lc%val
1107 0 : write(*,2) 'w_face%val', k, w_face%val
1108 0 : write(*,2) 'Lc_div_w_face', k, Lc_div_w_face%val
1109 0 : write(*,2) 'PII_face%val', k, PII_face%val
1110 0 : write(*,2) 'T_rho_face%val', k, T_rho_face%val
1111 : !write(*,2) '', k,
1112 : !write(*,2) '', k,
1113 0 : call mesa_error(__FILE__,__LINE__,'compute_Lc_terms')
1114 : end if
1115 0 : end function compute_Lc_terms
1116 :
1117 :
1118 0 : function compute_Lt(s, k, ierr) result(Lt) ! erg s^-1
1119 : type (star_info), pointer :: s
1120 : integer, intent(in) :: k
1121 : type(auto_diff_real_star_order1) :: Lt
1122 : integer, intent(out) :: ierr
1123 : type(auto_diff_real_star_order1) :: r_00, area2, d_m1, d_00, &
1124 : rho2_face, Hp_face, w_m1, w_00, w_face, etrb_m1, etrb_00
1125 0 : real(dp) :: alpha_alpha_t, alfa, beta
1126 : include 'formats'
1127 0 : ierr = 0
1128 0 : if (k > s% nz) then
1129 0 : Lt = 0d0
1130 0 : return
1131 : end if
1132 0 : alpha_alpha_t = s% mixing_length_alpha*s% RSP2_alfat
1133 : if (alpha_alpha_t == 0d0 .or. &
1134 0 : k <= s% RSP2_num_outermost_cells_forced_nonturbulent .or. &
1135 : k > s% nz - int(s% nz/s% RSP2_nz_div_IBOTOM)) then
1136 0 : Lt = 0d0
1137 0 : s% Lt(k) = 0d0
1138 0 : return
1139 : end if
1140 0 : r_00 = wrap_r_00(s,k)
1141 0 : area2 = pow2(4d0*pi*pow2(r_00))
1142 0 : d_m1 = wrap_d_m1(s,k)
1143 0 : d_00 = wrap_d_00(s,k)
1144 0 : call get_RSP2_alfa_beta_face_weights(s, k, alfa, beta)
1145 0 : rho2_face = alfa*pow2(d_00) + beta*pow2(d_m1)
1146 0 : w_m1 = wrap_w_m1(s,k)
1147 0 : w_00 = wrap_w_00(s,k)
1148 0 : w_face = alfa*w_00 + beta*w_m1
1149 0 : etrb_m1 = wrap_etrb_m1(s,k)
1150 0 : etrb_00 = wrap_etrb_00(s,k)
1151 0 : Hp_face = wrap_Hp_00(s,k)
1152 : ! Ft = - alpha_t * rho_face * alpha * Hp_face * w_face * detrb/dr (thesis eqn 2.44)
1153 : ! replace dr by dm_bar/(area*rho_face)
1154 : ! Ft = - alpha_alpha_t * rho_face * Hp_face * w_face * (area*rho_face) * detrb/dm_bar
1155 : ! Lt = area * Ft
1156 : ! Lt = -alpha_alpha_t * (area*rho_face)**2 * Hp_face * w_face * (etrb(k-1) - etrb(k))/dm_bar
1157 0 : Lt = - alpha_alpha_t * area2 * rho2_face * Hp_face * w_face * (etrb_m1 - etrb_00) / s% dm_bar(k)
1158 : ! units = (cm^4) (g^2 cm^-6) (cm) (cm s^-1) (ergs g^-1) g^-1 = erg s^-1
1159 0 : s% Lt(k) = Lt%val
1160 0 : end function compute_Lt
1161 :
1162 :
1163 0 : subroutine set_etrb_start_vars(s, ierr)
1164 : type (star_info), pointer :: s
1165 : integer, intent(out) :: ierr
1166 : integer :: k
1167 : type(auto_diff_real_star_order1) :: Y_face, Lt
1168 : include 'formats'
1169 0 : ierr = 0
1170 0 : do k=1,s%nz
1171 0 : Y_face = compute_Y_face(s, k, ierr)
1172 0 : if (ierr /= 0) return
1173 0 : s% Y_face_start(k) = Y_face%val
1174 0 : Lt = compute_Lt(s, k, ierr)
1175 0 : if (ierr /= 0) return
1176 0 : s% Lt_start(k) = Lt%val
1177 0 : s% w_start(k) = s% w(k)
1178 0 : s% Hp_face_start(k) = s% Hp_face(k)
1179 : end do
1180 : end subroutine set_etrb_start_vars
1181 :
1182 :
1183 0 : subroutine RSP2_adjust_vars_before_call_solver(s, ierr) ! replaces check_omega in RSP
1184 : ! JAK OKRESLIC OMEGA DLA PIERWSZEJ ITERACJI
1185 : use micro, only: do_eos_for_cell
1186 : type (star_info), pointer :: s
1187 : integer, intent(out) :: ierr
1188 0 : real(dp) :: PII_div_Hp, QQ, SOURCE, Hp_cell, DAMP, POM, POM2, DAMPR, del, soln
1189 : !type(auto_diff_real_star_order1) :: x
1190 : integer :: k
1191 : include 'formats'
1192 0 : ierr = 0
1193 0 : if (s% mixing_length_alpha == 0d0) return
1194 :
1195 0 : !$OMP PARALLEL DO PRIVATE(k,PII_div_Hp,QQ,SOURCE,Hp_cell,DAMP,POM,POM2,DAMPR,del,soln) SCHEDULE(dynamic,2)
1196 : do k=s% RSP2_num_outermost_cells_forced_nonturbulent+1, &
1197 : s% nz - max(1,int(s% nz/s% RSP_nz_div_IBOTOM))
1198 :
1199 : if (s% w(k) > s% RSP2_w_min_for_damping) cycle
1200 :
1201 : PII_div_Hp = 0.5d0*(s% PII(k)/s% Hp_face(k) + s% PII(k+1)/s% Hp_face(k+1))
1202 : QQ = s% chiT(k)/(s% rho(k)*s% T(k)*s% chiRho(k))
1203 : SOURCE = PII_div_Hp*s% T(k)*s% Peos(k)*QQ/s% Cp(k)
1204 :
1205 : Hp_cell = 0.5d0*(s% Hp_face(k) + s% Hp_face(k+1))
1206 : DAMP = (s% RSP2_alfad*x_CEDE/s% mixing_length_alpha)/Hp_cell
1207 :
1208 : POM = 4d0*boltz_sigma*pow2(s% RSP2_alfar*x_GAMMAR/s% mixing_length_alpha)
1209 : POM2 = pow3(s% T(k))/(pow2(s% rho(k))*s% Cp(k)*s% opacity(k))
1210 : DAMPR = POM*POM2/pow2(Hp_cell)
1211 :
1212 : del = pow2(DAMPR) + 4d0*DAMP*SOURCE
1213 :
1214 : if (k==-35) then
1215 : write(*,2) 'del', k, del
1216 : write(*,2) 'DAMPR', k, DAMPR
1217 : write(*,2) 'DAMP', k, DAMP
1218 : write(*,2) 'SOURCE', k, SOURCE
1219 : write(*,2) 'POM', k, PII_div_Hp
1220 : write(*,2) 'POM2', k, s% T(k)*s% Peos(k)*QQ/s% Cp(k)
1221 : write(*,2) 's% Hp_face(k)', k, s% Hp_face(k)
1222 : write(*,2) 's% Hp_face(k+1)', k+1, s% Hp_face(k+1)
1223 : write(*,2) 's% PII(k)', k, s% PII(k)
1224 : write(*,2) 's% PII(k+1)', k+1, s% PII(k+1)
1225 : write(*,2) 's% Y_face(k)', k, s% Y_face(k)
1226 : write(*,2) 's% Y_face(k+1)', k+1, s% Y_face(k+1)
1227 : end if
1228 :
1229 : if (del < 0d0) cycle
1230 : soln = (-DAMPR + sqrt(del))/(2d0*DAMP)
1231 : if (k==-35) write(*,2) 'soln', k, soln
1232 : if (soln > 0d0) then
1233 : ! i tried soln = sqrt(soln) here. helps solver convergence, but hurts the model results.
1234 : if (s% RSP2_report_adjust_w) &
1235 : write(*,3) 'RSP2_adjust_vars_before_call_solver w', k, s% model_number, s% w(k), soln
1236 : s% w(k) = soln
1237 : end if
1238 :
1239 : end do
1240 : !$OMP END PARALLEL DO
1241 0 : end subroutine RSP2_adjust_vars_before_call_solver
1242 :
1243 : end module hydro_rsp2
|
