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