Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2025 Ebraheem Farag & 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 tdc_hydro
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 :: &
34 : compute_tdc_Eq_cell, compute_tdc_Uq_face, compute_tdc_Eq_div_w_face, &
35 : get_RSP2_alfa_beta_face_weights, set_viscosity_vars_TDC
36 :
37 : contains
38 :
39 : ! This routine is called to initialize eq and uq for TDC.
40 0 : subroutine set_viscosity_vars_TDC(s, ierr)
41 : type(star_info), pointer :: s
42 : integer, intent(out) :: ierr
43 : type(auto_diff_real_star_order1) :: x
44 : integer :: k, op_err
45 : include 'formats'
46 0 : ierr = 0
47 0 : op_err = 0
48 :
49 0 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
50 : do k = 1, s%nz
51 : ! Hp_face(k) <= 0 means it needs to be set. e.g., after read file
52 : if (s%Hp_face(k) <= 0) then
53 : ! this scale height for face is already calculated in TDC
54 : s%Hp_face(k) = get_scale_height_face_val(s, k) ! because this is called before s% scale_height(k) is updated in mlt_vars.
55 : end if
56 : end do
57 : !$OMP END PARALLEL DO
58 0 : if (ierr /= 0) then
59 0 : if (s%report_ierr) write (*, 2) 'failed in set_viscosity_vars_TDC loop 1', s%model_number
60 0 : return
61 : end if
62 0 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
63 : do k = 1, s%nz
64 : x = compute_Chi_cell(s, k, op_err)
65 : if (op_err /= 0) ierr = op_err
66 : x = compute_tdc_Eq_cell(s, k, op_err)
67 : if (op_err /= 0) ierr = op_err
68 : x = compute_tdc_Uq_face(s, k, op_err)
69 : if (op_err /= 0) ierr = op_err
70 : end do
71 : !$OMP END PARALLEL DO
72 0 : if (ierr /= 0) then
73 0 : if (s%report_ierr) write (*, 2) 'failed in set_viscosity_vars_TDC loop 2', s%model_number
74 0 : return
75 : end if
76 0 : if (.not. (s%v_flag .or. s%u_flag)) then ! set values 0 if not using v_flag or u_flag.
77 0 : do k = 1, s%nz
78 0 : s%Eq(k) = 0d0; s%Eq_ad(k) = 0d0
79 0 : s%Chi(k) = 0d0; s%Chi_ad(k) = 0d0
80 0 : s%Uq(k) = 0d0
81 : end do
82 : end if
83 : end subroutine set_viscosity_vars_TDC
84 :
85 0 : subroutine get_RSP2_alfa_beta_face_weights(s, k, alfa, beta)
86 : type(star_info), pointer :: s
87 : integer, intent(in) :: k
88 : real(dp), intent(out) :: alfa, beta
89 : ! face_value = alfa*cell_value(k) + beta*cell_value(k-1)
90 0 : if (k == 1) call mesa_error(__FILE__, __LINE__, 'bad k==1 for get_RSP2_alfa_beta_face_weights')
91 0 : if (s%TDC_hydro_use_mass_interp_face_values) then
92 0 : alfa = s%dq(k - 1)/(s%dq(k - 1) + s%dq(k))
93 0 : beta = 1d0 - alfa
94 : else
95 0 : alfa = 0.5d0
96 0 : beta = 0.5d0
97 : end if
98 0 : end subroutine get_RSP2_alfa_beta_face_weights
99 :
100 0 : function compute_d_v_div_r(s, k, ierr) result(d_v_div_r) ! s^-1
101 : type(star_info), pointer :: s
102 : integer, intent(in) :: k
103 : type(auto_diff_real_star_order1) :: d_v_div_r
104 : integer, intent(out) :: ierr
105 : type(auto_diff_real_star_order1) :: v_00, v_p1, r_00, r_p1, term1, term2
106 : logical :: dbg
107 : include 'formats'
108 0 : ierr = 0
109 0 : dbg = .false.
110 0 : v_00 = wrap_v_00(s, k)
111 0 : v_p1 = wrap_v_p1(s, k)
112 0 : r_00 = wrap_r_00(s, k)
113 0 : r_p1 = wrap_r_p1(s, k)
114 0 : if (r_p1%val == 0d0) r_p1 = 1d0
115 0 : d_v_div_r = v_00/r_00 - v_p1/r_p1 ! units s^-1
116 :
117 : ! Debugging output to trace values
118 : if (dbg .and. k == -63) then
119 : write (*, *) 'test d_v_div_r, k:', k
120 : write (*, *) 'v_00:', v_00%val, 'v_p1:', v_p1%val
121 : write (*, *) 'r_00:', r_00%val, 'r_p1:', r_p1%val
122 : write (*, *) 'd_v_div_r:', d_v_div_r%val
123 : end if
124 0 : end function compute_d_v_div_r
125 :
126 0 : function compute_rho_form_of_d_v_div_r(s, k, ierr) result(d_v_div_r)
127 : type(star_info), pointer :: s
128 : integer, intent(in) :: k
129 : integer, intent(out) :: ierr
130 : type(auto_diff_real_star_order1) :: d_v_div_r
131 : type(auto_diff_real_star_order1) :: r_cell, rho_cell, v_cell, dlnrho_dt
132 : real(dp) :: dm_cell
133 0 : ierr = 0
134 :
135 0 : r_cell = 0.5d0*(wrap_r_00(s, k) + wrap_r_p1(s, k))
136 0 : rho_cell = wrap_d_00(s, k)
137 0 : v_cell = wrap_v_00(s, k) ! cell-centred velocity (u_flag)
138 0 : dlnrho_dt = wrap_dxh_lnd(s, k)/s%dt ! (∂/∂t)lnρ
139 0 : dm_cell = s%dm(k) ! cell mass
140 :
141 0 : d_v_div_r = -dm_cell/(4d0*pi*rho_cell)*(dlnrho_dt/pow3(r_cell) + 3d0*v_cell/pow4(r_cell))
142 0 : end function compute_rho_form_of_d_v_div_r
143 :
144 0 : function compute_rho_form_of_d_v_div_r_opt_time_center(s, k, ierr) result(d_v_div_r) ! s^-1
145 : type(star_info), pointer :: s
146 : integer, intent(in) :: k
147 : integer, intent(out) :: ierr
148 : type(auto_diff_real_star_order1) :: d_v_div_r
149 : type(auto_diff_real_star_order1) :: r_cell, rho_cell, v_cell, dlnrho_dt
150 : real(dp) :: dm_cell
151 0 : ierr = 0
152 :
153 0 : r_cell = 0.5d0*(wrap_opt_time_center_r_00(s, k) + wrap_opt_time_center_r_p1(s, k))
154 0 : rho_cell = wrap_d_00(s, k)
155 0 : v_cell = wrap_opt_time_center_v_00(s, k) ! cell-centred velocity (u_flag)
156 0 : dlnrho_dt = wrap_dxh_lnd(s, k)/s%dt ! (∂/∂t)lnρ
157 0 : dm_cell = s%dm(k) ! cell mass
158 :
159 : d_v_div_r = -dm_cell/(4d0*pi*rho_cell)* &
160 : (dlnrho_dt/pow3(r_cell) &
161 0 : + 3d0*v_cell/pow4(r_cell))
162 0 : end function compute_rho_form_of_d_v_div_r_opt_time_center
163 :
164 :
165 0 : function compute_rho_form_of_d_v_div_r_face(s, k, ierr) result(d_v_div_r)
166 : type(star_info), pointer :: s
167 : integer, intent(in) :: k
168 : integer, intent(out) :: ierr
169 : type(auto_diff_real_star_order1) :: d_v_div_r
170 : type(auto_diff_real_star_order1) :: r_face, rho_face, v_face, dlnrho_dt
171 0 : real(dp) :: dm_bar
172 0 : ierr = 0
173 :
174 0 : r_face = wrap_r_00(s, k)
175 0 : rho_face = get_rho_face(s, k)
176 0 : if (s% v_flag .and. .not. s% u_flag) then
177 0 : v_face = wrap_v_00(s, k) ! face-centred velocity (v_flag)
178 0 : else if (s% u_flag .and. .not. s% v_flag) then
179 0 : v_face = s% u_face_ad(k) ! reconstructed face velocity (u_flag)
180 : end if
181 0 : dlnrho_dt = 0.5d0*(wrap_dxh_lnd(s, k) + wrap_dxh_lnd(s, k+1))/s%dt ! (∂/∂t)lnρ
182 :
183 0 : if (k >= 2) then
184 0 : dm_bar = 0.5d0*(s% dm(k) + s% dm(k-1))
185 : else
186 0 : dm_bar = 0.5d0*s% dm(k)
187 : end if
188 :
189 0 : d_v_div_r = -dm_bar/(4d0*pi*rho_face)*(dlnrho_dt/pow3(r_face) + 3d0*v_face/pow4(r_face))
190 0 : end function compute_rho_form_of_d_v_div_r_face
191 :
192 0 : function compute_rho_form_of_d_v_div_r_face_opt_time_center(s, k, ierr) result(d_v_div_r) ! s^-1
193 : type(star_info), pointer :: s
194 : integer, intent(in) :: k
195 : integer, intent(out) :: ierr
196 : type(auto_diff_real_star_order1) :: d_v_div_r
197 : type(auto_diff_real_star_order1) :: r_face, rho_face, v_face, dlnrho_dt
198 0 : real(dp) :: dm_bar
199 0 : ierr = 0
200 :
201 0 : r_face = wrap_opt_time_center_r_00(s, k)
202 0 : rho_face = get_rho_face(s, k)
203 :
204 0 : if (s% v_flag .and. .not. s% u_flag) then
205 0 : v_face = wrap_opt_time_center_v_00(s, k) ! face-centred velocity (v_flag)
206 0 : else if (s% u_flag .and. .not. s% v_flag) then
207 0 : v_face = s% u_face_ad(k) ! reconstructed face velocity (u_flag)
208 : end if
209 0 : dlnrho_dt = 0.5d0*(wrap_dxh_lnd(s, k) + wrap_dxh_lnd(s, k+1))/s%dt ! (∂/∂t)lnρ
210 :
211 0 : if (k >= 2) then
212 0 : dm_bar = 0.5d0*(s% dm(k) + s% dm(k-1))
213 : else
214 0 : dm_bar = 0.5d0*s% dm(k)
215 : end if
216 :
217 0 : d_v_div_r = -dm_bar/(4d0*pi*rho_face)*(dlnrho_dt/pow3(r_face) + 3d0*v_face/pow4(r_face))
218 0 : end function compute_rho_form_of_d_v_div_r_face_opt_time_center
219 :
220 :
221 0 : function compute_d_v_div_r_opt_time_center(s, k, ierr) result(d_v_div_r) ! s^-1
222 : type(star_info), pointer :: s
223 : integer, intent(in) :: k
224 : type(auto_diff_real_star_order1) :: d_v_div_r
225 : integer, intent(out) :: ierr
226 : type(auto_diff_real_star_order1) :: v_00, v_p1, r_00, r_p1
227 : include 'formats'
228 0 : ierr = 0
229 0 : v_00 = wrap_opt_time_center_v_00(s, k)
230 0 : v_p1 = wrap_opt_time_center_v_p1(s, k)
231 0 : r_00 = wrap_opt_time_center_r_00(s, k)
232 0 : r_p1 = wrap_opt_time_center_r_p1(s, k)
233 0 : if (r_p1%val == 0d0) r_p1 = 1d0
234 0 : d_v_div_r = v_00/r_00 - v_p1/r_p1 ! units s^-1
235 0 : end function compute_d_v_div_r_opt_time_center
236 :
237 0 : function wrap_Hp_cell(s, k) result(Hp_cell) ! cm , different than rsp2
238 : type(star_info), pointer :: s
239 : integer, intent(in) :: k
240 : type(auto_diff_real_star_order1) :: Hp1, Hp0, Hp_cell
241 0 : Hp0 = get_scale_height_face(s,k)
242 0 : Hp1 = 0d0
243 0 : if (k+1 < s%nz) then
244 0 : Hp1 = get_scale_height_face(s,k+1)
245 : end if
246 0 : Hp_cell = 0.5d0*(Hp0 + Hp1)
247 : !0.5d0*(wrap_Hp_00(s, k) + wrap_Hp_p1(s, k))
248 0 : end function wrap_Hp_cell
249 :
250 0 : function Hp_cell_for_Chi(s, k, ierr) result(Hp_cell) ! cm
251 : type(star_info), pointer :: s
252 : integer, intent(in) :: k
253 : integer, intent(out) :: ierr
254 : type(auto_diff_real_star_order1) :: Hp_cell
255 : type(auto_diff_real_star_order1) :: d_00, Peos_00, rmid
256 0 : real(dp) :: mmid, cgrav_mid
257 : include 'formats'
258 0 : ierr = 0
259 :
260 0 : Hp_cell = wrap_Hp_cell(s, k)
261 0 : return ! below is skipped, for now.
262 :
263 : d_00 = wrap_d_00(s, k)
264 : Peos_00 = wrap_Peos_00(s, k)
265 : if (k < s%nz) then
266 : rmid = 0.5d0*(wrap_r_00(s, k) + wrap_r_p1(s, k))
267 : mmid = 0.5d0*(s%m(k) + s%m(k + 1))
268 : cgrav_mid = 0.5d0*(s%cgrav(k) + s%cgrav(k + 1))
269 : else
270 : rmid = 0.5d0*(wrap_r_00(s, k) + s%r_center)
271 : mmid = 0.5d0*(s%m(k) + s%m_center)
272 : cgrav_mid = s%cgrav(k)
273 : end if
274 : Hp_cell = pow2(rmid)*Peos_00/(d_00*cgrav_mid*mmid)
275 : if (s%alt_scale_height_flag) then
276 : call mesa_error(__FILE__, __LINE__, 'Hp_cell_for_Chi: cannot use alt_scale_height_flag')
277 : end if
278 : end function Hp_cell_for_Chi
279 :
280 0 : function compute_Chi_cell(s, k, ierr) result(Chi_cell)
281 : ! eddy viscosity energy (Kuhfuss 1986) [erg]
282 : type(star_info), pointer :: s
283 : integer, intent(in) :: k
284 : type(auto_diff_real_star_order1) :: Chi_cell
285 : integer, intent(out) :: ierr
286 : type(auto_diff_real_star_order1) :: &
287 : rho2, r6_cell, d_v_div_r, Hp_cell, w_00, d_00, r_00, r_p1
288 0 : real(dp) :: f, ALFAM_ALFA
289 : logical :: dbg
290 : include 'formats'
291 0 : ierr = 0
292 0 : dbg = .false.
293 :
294 : ! check where we are getting alfam from.
295 0 : if (s%MLT_option == 'TDC' .and. .not. s%RSP2_flag) then
296 0 : ALFAM_ALFA = s%alpha_TDC_DampM*s%mixing_length_alpha
297 : else ! this is for safety, but probably is never called.
298 : ALFAM_ALFA = 0d0
299 : end if
300 :
301 0 : if (ALFAM_ALFA == 0d0 .or. &
302 : k > s% nz - s% TDC_num_innermost_cells_forced_nonturbulent) then
303 0 : Chi_cell = 0d0
304 0 : if (k >= 1 .and. k <= s%nz) then
305 0 : s%Chi(k) = 0d0
306 0 : s%Chi_ad(k) = 0d0
307 : end if
308 : else
309 0 : Hp_cell = Hp_cell_for_Chi(s, k, ierr)
310 0 : if (ierr /= 0) return
311 0 : if (s%u_flag .or. s%TDC_use_density_form_for_eddy_viscosity) then
312 : ! new density derivative term
313 0 : d_v_div_r = compute_rho_form_of_d_v_div_r(s, k, ierr)
314 : else
315 0 : d_v_div_r = compute_d_v_div_r(s, k, ierr)
316 : end if
317 0 : if (ierr /= 0) return
318 :
319 : ! don't need to check if mlt_vc > 0 here.
320 0 : if (s%MLT_option == 'TDC' .and. .not. s%RSP2_flag) then
321 0 : w_00 = s% mlt_vc_ad(k)/sqrt_2_div_3! same as info%A0 from TDC
322 : else ! normal RSP2
323 0 : w_00 = wrap_w_00(s, k)
324 : end if
325 0 : d_00 = wrap_d_00(s, k)
326 0 : f = (16d0/3d0)*pi*ALFAM_ALFA/s%dm(k)
327 0 : rho2 = pow2(d_00)
328 0 : r_00 = wrap_r_00(s, k)
329 0 : r_p1 = wrap_r_p1(s, k)
330 0 : r6_cell = 0.5d0*(pow6(r_00) + pow6(r_p1))
331 0 : Chi_cell = f*rho2*r6_cell*d_v_div_r*Hp_cell*w_00
332 : ! units = g^-1 cm s^-1 g^2 cm^-6 cm^6 s^-1 cm
333 : ! = g cm^2 s^-2
334 : ! = erg
335 :
336 : end if
337 0 : s%Chi(k) = Chi_cell%val
338 0 : s%Chi_ad(k) = Chi_cell
339 :
340 : if (dbg .and. k == -100) then
341 : write (*, *) ' s% ALFAM_ALFA', ALFAM_ALFA
342 : write (*, *) 'Hp_cell', Hp_cell%val
343 : write (*, *) 'd_v_div_r', d_v_div_r%val
344 : write (*, *) ' f', f
345 : write (*, *) 'w_00', w_00%val
346 : write (*, *) 'd_00 ', d_00%val
347 : write (*, *) 'rho2 ', rho2%val
348 : write (*, *) 'r_00', r_00%val
349 : write (*, *) 'r_p1 ', r_p1%val
350 : write (*, *) 'r6_cell', r6_cell%val
351 : end if
352 0 : end function compute_Chi_cell
353 :
354 0 : function compute_tdc_Eq_cell(s, k, ierr) result(Eq_cell) ! erg g^-1 s^-1
355 : type(star_info), pointer :: s
356 : integer, intent(in) :: k
357 : type(auto_diff_real_star_order1) :: Eq_cell
358 : integer, intent(out) :: ierr
359 : type(auto_diff_real_star_order1) :: d_v_div_r, Chi_cell
360 : include 'formats'
361 0 : ierr = 0
362 0 : if (s%mixing_length_alpha == 0d0 .or. &
363 : k > s% nz - s% TDC_num_innermost_cells_forced_nonturbulent) then
364 0 : Eq_cell = 0d0
365 0 : if (k >= 1 .and. k <= s%nz) s%Eq_ad(k) = 0d0
366 : else
367 0 : Chi_cell = s%Chi_ad(k) ! compute_Chi_cell(s,k,ierr)
368 0 : if (ierr /= 0) return
369 :
370 0 : if (s%u_flag .or. s%TDC_use_density_form_for_eddy_viscosity) then
371 : ! new density derivative term
372 0 : d_v_div_r = compute_rho_form_of_d_v_div_r_opt_time_center(s, k, ierr)
373 : else
374 0 : d_v_div_r = compute_d_v_div_r_opt_time_center(s, k, ierr)
375 : end if
376 :
377 0 : if (ierr /= 0) return
378 0 : Eq_cell = 4d0*pi*Chi_cell*d_v_div_r/s%dm(k) ! erg s^-1 g^-1
379 : end if
380 0 : s%Eq(k) = Eq_cell%val
381 0 : s%Eq_ad(k) = Eq_cell
382 0 : end function compute_tdc_Eq_cell
383 :
384 0 : function compute_tdc_Uq_face(s, k, ierr) result(Uq_face) ! cm s^-2, acceleration
385 : type(star_info), pointer :: s
386 : integer, intent(in) :: k
387 : type(auto_diff_real_star_order1) :: Uq_face
388 : integer, intent(out) :: ierr
389 : type(auto_diff_real_star_order1) :: Chi_00, Chi_m1, r_00
390 : include 'formats'
391 0 : ierr = 0
392 : if (s%mixing_length_alpha == 0d0 .or. &
393 0 : k <= s%RSP2_num_outermost_cells_forced_nonturbulent .or. &
394 : k > s%nz - int(s%nz/s%RSP2_nz_div_IBOTOM)) then
395 0 : Uq_face = 0d0
396 : else
397 0 : r_00 = wrap_opt_time_center_r_00(s, k)
398 :
399 : ! which do we adopt?
400 0 : Chi_00 = compute_Chi_cell(s, k, ierr) ! s% Chi_ad(k) XXX
401 : !Chi_00 = s% Chi_ad(k) ! compute_Chi_cell(s,k,ierr)
402 :
403 0 : if (k > 1) then
404 0 : Chi_m1 = shift_m1(compute_Chi_cell(s, k - 1, ierr))
405 : !Chi_m1 = shift_m1(s% Chi_ad(k-1)) XXX
406 0 : if (ierr /= 0) return
407 : else
408 0 : Chi_m1 = 0d0
409 : end if
410 0 : Uq_face = 4d0*pi*(Chi_m1 - Chi_00)/(r_00*s%dm_bar(k))
411 :
412 0 : if (k == -56) then
413 0 : write (*, 3) 'RSP2 Uq chi_m1 chi_00 r', k, s%solver_iter, &
414 0 : Uq_face%val, Chi_m1%val, Chi_00%val, r_00%val
415 : end if
416 :
417 : end if
418 : ! erg g^-1 cm^-1 = g cm^2 s^-2 g^-1 cm^-1 = cm s^-2, acceleration
419 0 : s%Uq(k) = Uq_face%val
420 0 : end function compute_tdc_Uq_face
421 :
422 : ! face centered variables for tdc update below
423 0 : function compute_Chi_div_w_face(s, k, ierr) result(Chi_face)
424 : ! eddy viscosity energy (Kuhfuss 1986) [erg]
425 : type(star_info), pointer :: s
426 : integer, intent(in) :: k
427 : type(auto_diff_real_star_order1) :: Chi_face
428 : integer, intent(out) :: ierr
429 : type(auto_diff_real_star_order1) :: &
430 : rho2, r6_face, d_v_div_r, Hp_face, w_00, d_00, r_00, r_p1
431 0 : real(dp) :: f, ALFAM_ALFA, dmbar
432 : logical :: dbg
433 : include 'formats'
434 0 : ierr = 0
435 0 : dbg = .false.
436 :
437 : ! check where we are getting alfam from.
438 0 : if (s%MLT_option == 'TDC' .and. .not. s%RSP2_flag) then
439 0 : ALFAM_ALFA = s%alpha_TDC_DampM*s%mixing_length_alpha
440 : else ! this is for safety, but probably is never called.
441 : ALFAM_ALFA = 0d0
442 : end if
443 :
444 0 : if (ALFAM_ALFA == 0d0 .or. &
445 : k > s%nz - s% TDC_num_innermost_cells_forced_nonturbulent) then
446 0 : Chi_face = 0d0
447 0 : if (k >= 1 .and. k <= s%nz) then
448 0 : s%Chi(k) = 0d0
449 0 : s%Chi_ad(k) = 0d0
450 : end if
451 : else
452 0 : Hp_face = get_scale_height_face(s,k) !Hp_cell_for_Chi(s, k, ierr)
453 0 : if (ierr /= 0) return
454 0 : if (s%u_flag .or. s%TDC_use_density_form_for_eddy_viscosity) then
455 : ! new density derivative term
456 0 : d_v_div_r = compute_rho_form_of_d_v_div_r_face(s, k, ierr)
457 : else
458 0 : d_v_div_r = compute_d_v_div_r_face(s, k, ierr)
459 : end if
460 0 : if (ierr /= 0) return
461 :
462 0 : if (k >= 2) then
463 0 : dmbar = 0.5d0*(s% dm(k) + s% dm(k-1))
464 : else
465 0 : dmbar = 0.5d0*s% dm(k)
466 : end if
467 0 : d_00 = get_rho_face(s, k)
468 0 : f = (16d0/3d0)*pi*ALFAM_ALFA/dmbar
469 0 : rho2 = pow2(d_00)
470 0 : r_00 = wrap_r_00(s, k)
471 : !r_p1 = wrap_r_p1(s, k)
472 0 : r6_face = pow6(r_00) !0.5d0*(pow6(r_00) + pow6(r_p1))
473 0 : Chi_face = f*rho2*r6_face*d_v_div_r*Hp_face!*w_00
474 : ! units = g^-1 cm s^-1 g^2 cm^-6 cm^6 s^-1 cm
475 : ! = g cm^2 s^-2
476 : ! = erg ! * (s / cm) - > [erg] * 1/w00
477 :
478 : end if
479 : !s%Chi(k) = Chi_face%val
480 : !s%Chi_ad(k) = Chi_face
481 :
482 : if (dbg .and. k == -100) then
483 : write (*, *) ' s% ALFAM_ALFA', ALFAM_ALFA
484 : write (*, *) 'Hp_face', Hp_face%val
485 : write (*, *) 'd_v_div_r', d_v_div_r%val
486 : write (*, *) ' f', f
487 : write (*, *) 'w_00', w_00%val
488 : write (*, *) 'd_00 ', d_00%val
489 : write (*, *) 'rho2 ', rho2%val
490 : write (*, *) 'r_00', r_00%val
491 : write (*, *) 'r_p1 ', r_p1%val
492 : write (*, *) 'r6_cell', r6_face%val
493 : end if
494 0 : end function compute_Chi_div_w_face
495 :
496 :
497 :
498 0 : function compute_tdc_Eq_div_w_face(s, k, ierr) result(Eq_face) ! erg g^-1 s^-1 * (cm^-1 s^1)
499 : type(star_info), pointer :: s
500 : integer, intent(in) :: k
501 : type(auto_diff_real_star_order1) :: Eq_face
502 : integer, intent(out) :: ierr
503 : type(auto_diff_real_star_order1) :: d_v_div_r, Chi_face
504 0 : real(dp) :: dmbar
505 : include 'formats'
506 0 : ierr = 0
507 0 : if (s%mixing_length_alpha == 0d0 .or. &
508 : k > s%nz - s% TDC_num_innermost_cells_forced_nonturbulent) then
509 0 : Eq_face = 0d0
510 0 : if (k >= 1 .and. k <= s%nz) s%Eq_ad(k) = 0d0
511 : else
512 0 : Chi_face = compute_Chi_div_w_face(s,k,ierr)
513 0 : if (ierr /= 0) return
514 :
515 0 : if (s%u_flag .or. s%TDC_use_density_form_for_eddy_viscosity) then
516 : ! new density derivative term
517 0 : d_v_div_r = compute_rho_form_of_d_v_div_r_face_opt_time_center(s, k, ierr)
518 : else
519 0 : d_v_div_r = compute_d_v_div_r_opt_time_center_face(s, k, ierr)
520 : end if
521 :
522 0 : if (k >= 2) then
523 0 : dmbar = 0.5d0*(s% dm(k) + s% dm(k-1))
524 : else
525 0 : dmbar = 0.5d0*s% dm(k)
526 : end if
527 :
528 0 : if (ierr /= 0) return
529 0 : Eq_face = 4d0*pi*Chi_face*d_v_div_r/dmbar ! erg s^-1 g^-1 * (cm^-1 s^1)
530 : end if
531 : !s%Eq(k) = Eq_face%val
532 : !s%Eq_ad(k) = Eq_face
533 0 : end function compute_tdc_Eq_div_w_face
534 :
535 0 : function compute_d_v_div_r_face(s, k, ierr) result(d_v_div_r) ! s^-1
536 : type(star_info), pointer :: s
537 : integer, intent(in) :: k
538 : type(auto_diff_real_star_order1) :: d_v_div_r
539 : integer, intent(out) :: ierr
540 : type(auto_diff_real_star_order1) :: v_00, v_p1, r_00, r_p1, term1, term2
541 : logical :: dbg
542 : include 'formats'
543 0 : ierr = 0
544 0 : dbg = .false.
545 :
546 0 : if (k >1) then
547 0 : v_00 = 0.5d0 * (wrap_v_00(s, k) + wrap_v_m1(s, k))
548 : else
549 0 : v_00 = 0.5d0 * wrap_v_00(s, k)
550 : end if
551 :
552 0 : v_p1 = 0.5d0 * (wrap_v_00(s, k) + wrap_v_p1(s, k))
553 :
554 0 : if (k >1) then
555 0 : r_00 = 0.5d0 * (wrap_r_00(s, k) + wrap_r_m1(s, k))
556 : else
557 0 : r_00 = 0.5d0 * wrap_r_00(s, k)
558 : end if
559 0 : r_p1 = 0.5d0*(wrap_r_00(s, k) + wrap_r_p1(s, k))
560 0 : if (r_p1%val == 0d0) r_p1 = 1d0
561 0 : d_v_div_r = v_00/r_00 - v_p1/r_p1 ! units s^-1
562 :
563 : ! Debugging output to trace values
564 : if (dbg .and. k == -63) then
565 : write (*, *) 'test d_v_div_r, k:', k
566 : write (*, *) 'v_00:', v_00%val, 'v_p1:', v_p1%val
567 : write (*, *) 'r_00:', r_00%val, 'r_p1:', r_p1%val
568 : write (*, *) 'd_v_div_r:', d_v_div_r%val
569 : end if
570 0 : end function compute_d_v_div_r_face
571 :
572 :
573 :
574 0 : function compute_d_v_div_r_opt_time_center_face(s, k, ierr) result(d_v_div_r) ! s^-1
575 : type(star_info), pointer :: s
576 : integer, intent(in) :: k
577 : type(auto_diff_real_star_order1) :: d_v_div_r
578 : integer, intent(out) :: ierr
579 : type(auto_diff_real_star_order1) :: v_00, v_p1, r_00, r_p1
580 : include 'formats'
581 0 : ierr = 0
582 :
583 0 : if (k >1) then
584 0 : v_00 = 0.5d0 * (wrap_opt_time_center_v_00(s, k) + wrap_opt_time_center_v_m1(s, k))
585 : else
586 0 : v_00 = 0.5d0 * wrap_opt_time_center_v_00(s, k)
587 : end if
588 :
589 0 : v_p1 = 0.5d0 * (wrap_opt_time_center_v_00(s, k) + wrap_opt_time_center_v_p1(s, k))
590 :
591 0 : if (k >1) then
592 0 : r_00 = 0.5d0 * (wrap_opt_time_center_r_00(s, k) + wrap_opt_time_center_r_m1(s, k))
593 : else
594 0 : r_00 = 0.5d0 * wrap_opt_time_center_r_00(s, k)
595 : end if
596 0 : r_p1 = 0.5d0*(wrap_opt_time_center_r_00(s, k) + wrap_opt_time_center_r_p1(s, k))
597 0 : if (r_p1%val == 0d0) r_p1 = 1d0
598 0 : d_v_div_r = v_00/r_00 - v_p1/r_p1 ! units s^-1
599 :
600 0 : if (r_p1%val == 0d0) r_p1 = 1d0
601 0 : d_v_div_r = v_00/r_00 - v_p1/r_p1 ! units s^-1
602 0 : end function compute_d_v_div_r_opt_time_center_face
603 :
604 : end module tdc_hydro
|
