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