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