Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2025 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 reconstructed_face_support
21 :
22 : use star_private_def
23 : use const_def, only: dp, ln10, pi4, clight, crad
24 : use auto_diff
25 : use kap_support, only: get_kap
26 :
27 : implicit none
28 :
29 : private
30 : public :: get_reconstructed_face_state_ad
31 : public :: get_reconstructed_face_eos_kap_ad
32 : public :: get_reconstructed_scale_height_ad
33 :
34 : contains
35 :
36 : ! Returns the MLT/TDC face thermodynamic state as
37 : ! auto_diff_real_star_order1 quantities, either from recomputed face
38 : ! EOS/opacity data or from the stored face quantities.
39 79994 : subroutine get_reconstructed_face_state_ad( &
40 : s, k, T_face, rho_face, P_face, energy_face, Cp_face, ChiRho_face, ChiT_face, grada_face, &
41 : opacity_face, scale_height_face, gradr_face, ierr)
42 : use star_utils, only: get_T_face, get_Peos_face, get_e_face, get_kap_face, get_rho_face, &
43 : get_ChiRho_face, get_ChiT_face, get_Cp_face, get_grada_face, get_scale_height_face, get_gradr_face
44 :
45 : type(star_info), pointer :: s
46 : integer, intent(in) :: k
47 : type(auto_diff_real_star_order1), intent(out) :: &
48 : T_face, rho_face, P_face, energy_face, Cp_face, ChiRho_face, ChiT_face, grada_face, &
49 : opacity_face, scale_height_face, gradr_face
50 : integer, intent(out) :: ierr
51 :
52 79994 : ierr = 0
53 79994 : if (s%use_face_reconstruction) then
54 0 : call ensure_reconstructed_face_state_ad(s, k, ierr)
55 0 : if (ierr /= 0) return
56 0 : T_face = s%reconstructed_T_face_ad(k)
57 0 : rho_face = s%reconstructed_rho_face_ad(k)
58 0 : P_face = s%reconstructed_P_face_ad(k)
59 0 : energy_face = s%reconstructed_energy_face_ad(k)
60 0 : Cp_face = s%reconstructed_Cp_face_ad(k)
61 0 : ChiRho_face = s%reconstructed_ChiRho_face_ad(k)
62 0 : ChiT_face = s%reconstructed_ChiT_face_ad(k)
63 0 : grada_face = s%reconstructed_grada_face_ad(k)
64 0 : opacity_face = s%reconstructed_opacity_face_ad(k)
65 0 : scale_height_face = s%reconstructed_scale_height_face_ad(k)
66 0 : gradr_face = s%reconstructed_gradr_face_ad(k)
67 : else
68 79994 : T_face = get_T_face(s, k)
69 79994 : P_face = get_Peos_face(s, k)
70 79994 : energy_face = get_e_face(s, k)
71 79994 : opacity_face = get_kap_face(s, k)
72 79994 : rho_face = get_rho_face(s, k)
73 79994 : ChiRho_face = get_ChiRho_face(s, k)
74 79994 : ChiT_face = get_ChiT_face(s, k)
75 79994 : Cp_face = get_Cp_face(s, k)
76 79994 : grada_face = get_grada_face(s, k)
77 79994 : scale_height_face = get_scale_height_face(s, k)
78 79994 : gradr_face = get_gradr_face(s, k)
79 : end if
80 : end subroutine get_reconstructed_face_state_ad
81 :
82 :
83 : ! Ensures that the recomputed MLT and TDC face thermodynamic quantities have
84 : ! been assembled and cached for face k.
85 0 : subroutine ensure_reconstructed_face_state_ad(s, k, ierr)
86 : type(star_info), pointer :: s
87 : integer, intent(in) :: k
88 : integer, intent(out) :: ierr
89 :
90 : type(auto_diff_real_star_order1) :: T_face, rho_face, P_face, energy_face, Cp_face
91 : type(auto_diff_real_star_order1) :: ChiRho_face, ChiT_face, grada_face, opacity_face
92 : type(auto_diff_real_star_order1) :: scale_height_face, gradr_face
93 : real(dp) :: csound_face
94 :
95 0 : ierr = 0
96 0 : if (s%reconstructed_face_state_valid(k)) return
97 :
98 : call build_reconstructed_face_state_ad( &
99 : s, k, T_face, rho_face, P_face, energy_face, Cp_face, ChiRho_face, ChiT_face, grada_face, &
100 0 : opacity_face, scale_height_face, gradr_face, csound_face, ierr)
101 0 : if (ierr /= 0) return
102 :
103 0 : s%reconstructed_T_face_ad(k) = T_face
104 0 : s%reconstructed_rho_face_ad(k) = rho_face
105 0 : s%reconstructed_P_face_ad(k) = P_face
106 0 : s%reconstructed_energy_face_ad(k) = energy_face
107 0 : s%reconstructed_Cp_face_ad(k) = Cp_face
108 0 : s%reconstructed_ChiRho_face_ad(k) = ChiRho_face
109 0 : s%reconstructed_ChiT_face_ad(k) = ChiT_face
110 0 : s%reconstructed_grada_face_ad(k) = grada_face
111 0 : s%reconstructed_opacity_face_ad(k) = opacity_face
112 0 : s%reconstructed_scale_height_face_ad(k) = scale_height_face
113 0 : s%reconstructed_gradr_face_ad(k) = gradr_face
114 0 : s%reconstructed_csound_face(k) = csound_face
115 0 : s%reconstructed_face_state_valid(k) = .true.
116 : end subroutine ensure_reconstructed_face_state_ad
117 :
118 : ! Reconstructs the face composition from either the current or the
119 : ! start-of-step composition, renormalizes xa_face, and derives zbar_face.
120 0 : subroutine get_face_composition(s, k, use_starting_comp, zbar_face, xa_face, ierr)
121 : use chem_lib, only: basic_composition_info
122 : use star_utils, only: get_face_weights
123 :
124 : type(star_info), pointer :: s
125 : integer, intent(in) :: k
126 : logical, intent(in) :: use_starting_comp
127 : real(dp), intent(out) :: zbar_face
128 : real(dp), intent(out) :: xa_face(:)
129 : integer, intent(out) :: ierr
130 :
131 : real(dp) :: alfa, beta, sum_xa
132 : real(dp) :: xh, xhe, z, abar, z2bar, z53bar, ye, mass_correction, sumx
133 :
134 0 : ierr = 0
135 0 : if (k == 1) then
136 0 : alfa = 1d0
137 0 : beta = 0d0
138 : else
139 0 : call get_face_weights(s, k, alfa, beta)
140 : end if
141 :
142 0 : if (use_starting_comp) then
143 0 : xa_face(1:s%species) = alfa*s%xa_start(1:s%species, k)
144 0 : if (k > 1) xa_face(1:s%species) = xa_face(1:s%species) + beta*s%xa_start(1:s%species, k-1)
145 : else
146 0 : xa_face(1:s%species) = alfa*s%xa(1:s%species, k)
147 0 : if (k > 1) xa_face(1:s%species) = xa_face(1:s%species) + beta*s%xa(1:s%species, k-1)
148 : end if
149 :
150 0 : sum_xa = sum(xa_face)
151 0 : if (sum_xa <= 0d0) then
152 0 : ierr = -1
153 0 : if (s%report_ierr) then
154 0 : !$OMP critical (reconstructed_face_report_ierr)
155 0 : write(*,*) 'get_face_composition: sum_xa <= 0 for k', k
156 : !$OMP end critical (reconstructed_face_report_ierr)
157 : end if
158 0 : return
159 : end if
160 0 : xa_face = xa_face/sum_xa
161 : call basic_composition_info( &
162 : s%species, s%chem_id, xa_face, xh, xhe, z, abar, zbar_face, &
163 0 : z2bar, z53bar, ye, mass_correction, sumx)
164 : end subroutine get_face_composition
165 :
166 :
167 : ! Builds the recomputed face EOS input state by wrapping T and rho to the
168 : ! face, reconstructing the face composition, and evaluating the EOS there.
169 0 : subroutine get_face_eos_inputs( &
170 : s, k, T_face, rho_face, eos_res, d_dlnd, d_dlnT, ierr)
171 : use eos_support, only: get_eos
172 : use eos_def, only: num_eos_basic_results, num_eos_d_dxa_results
173 : use star_utils, only: get_rho_face, get_T_face
174 :
175 : type(star_info), pointer :: s
176 : integer, intent(in) :: k
177 : type(auto_diff_real_star_order1), intent(out) :: T_face, rho_face
178 : real(dp), intent(out) :: eos_res(num_eos_basic_results), d_dlnd(num_eos_basic_results), d_dlnT(num_eos_basic_results)
179 : integer, intent(out) :: ierr
180 :
181 : real(dp) :: log10_T, log10_rho, zbar_face
182 0 : real(dp) :: eos_xa_face(s%species)
183 0 : real(dp) :: d_dxa(num_eos_d_dxa_results, s%species)
184 :
185 0 : ierr = 0
186 0 : T_face = get_T_face(s, k)
187 0 : rho_face = get_rho_face(s, k)
188 0 : if (T_face%val <= 0d0 .or. rho_face%val <= 0d0) then
189 0 : ierr = -1
190 0 : if (s%report_ierr) then
191 0 : !$OMP critical (reconstructed_face_report_ierr)
192 0 : write(*,*) 'get_face_eos_inputs: bad face T or rho for k', k, T_face%val, rho_face%val
193 : !$OMP end critical (reconstructed_face_report_ierr)
194 : end if
195 0 : return
196 : end if
197 :
198 0 : call get_face_composition(s, k, .false., zbar_face, eos_xa_face, ierr)
199 0 : if (ierr /= 0) return
200 :
201 0 : log10_T = log10(T_face%val)
202 0 : log10_rho = log10(rho_face%val)
203 :
204 : call get_eos( &
205 : s, k, eos_xa_face, rho_face%val, log10_rho, T_face%val, log10_T, &
206 0 : eos_res, d_dlnd, d_dlnT, d_dxa, ierr)
207 0 : if (ierr /= 0) then
208 0 : if (s%report_ierr) call write_face_eos_call_info(s, k, T_face, rho_face, zbar_face, eos_xa_face)
209 0 : return
210 : end if
211 : end subroutine get_face_eos_inputs
212 :
213 :
214 : ! Interpolates extra_opacity_factor to the face and applies the existing
215 : ! logT taper used to turn that factor on and off.
216 0 : subroutine get_face_opacity_factor(s, k, log10_T, opacity_factor_face)
217 : use star_utils, only: get_face_weights
218 :
219 : type(star_info), pointer :: s
220 : integer, intent(in) :: k
221 : real(dp), intent(in) :: log10_T
222 : real(dp), intent(out) :: opacity_factor_face
223 :
224 : real(dp) :: alfa, beta
225 :
226 0 : if (k == 1) then
227 0 : alfa = 1d0
228 0 : beta = 0d0
229 : else
230 0 : call get_face_weights(s, k, alfa, beta)
231 : end if
232 :
233 0 : opacity_factor_face = alfa*s%extra_opacity_factor(k)
234 0 : if (k > 1) opacity_factor_face = opacity_factor_face + beta*s%extra_opacity_factor(k-1)
235 0 : if (s%min_logT_for_opacity_factor_off > 0) then
236 0 : if (log10_T >= s%max_logT_for_opacity_factor_off .or. &
237 : log10_T <= s%min_logT_for_opacity_factor_off) then
238 0 : opacity_factor_face = 1d0
239 0 : else if (log10_T > s%max_logT_for_opacity_factor_on) then
240 : opacity_factor_face = 1d0 + (opacity_factor_face - 1d0)* &
241 : (log10_T - s%max_logT_for_opacity_factor_off)/ &
242 0 : (s%max_logT_for_opacity_factor_on - s%max_logT_for_opacity_factor_off)
243 0 : else if (log10_T < s%min_logT_for_opacity_factor_on) then
244 : opacity_factor_face = 1d0 + (opacity_factor_face - 1d0)* &
245 : (log10_T - s%min_logT_for_opacity_factor_off)/ &
246 0 : (s%min_logT_for_opacity_factor_on - s%min_logT_for_opacity_factor_off)
247 : end if
248 : end if
249 0 : end subroutine get_face_opacity_factor
250 :
251 :
252 : ! Returns the cached or newly built face EOS and opacity state as
253 : ! auto_diff_real_star_order1 quantities for the MLT/TDC solve,
254 : ! instead of using the stored face quantities.
255 0 : subroutine get_reconstructed_face_eos_kap_ad( &
256 : s, k, T_face, rho_face, P_face, Cp_face, ChiRho_face, ChiT_face, grada_face, opacity_face, ierr)
257 : type(star_info), pointer :: s
258 : integer, intent(in) :: k
259 : type(auto_diff_real_star_order1), intent(out) :: T_face, rho_face, P_face, Cp_face, ChiRho_face, ChiT_face, grada_face, opacity_face
260 : integer, intent(out) :: ierr
261 :
262 : ierr = 0
263 0 : call ensure_reconstructed_face_state_ad(s, k, ierr)
264 0 : if (ierr /= 0) return
265 :
266 0 : T_face = s%reconstructed_T_face_ad(k)
267 0 : rho_face = s%reconstructed_rho_face_ad(k)
268 0 : P_face = s%reconstructed_P_face_ad(k)
269 0 : Cp_face = s%reconstructed_Cp_face_ad(k)
270 0 : ChiRho_face = s%reconstructed_ChiRho_face_ad(k)
271 0 : ChiT_face = s%reconstructed_ChiT_face_ad(k)
272 0 : grada_face = s%reconstructed_grada_face_ad(k)
273 0 : opacity_face = s%reconstructed_opacity_face_ad(k)
274 : end subroutine get_reconstructed_face_eos_kap_ad
275 :
276 :
277 : ! Builds the full set of recomputed face thermodynamic quantities for the
278 : ! MLT and TDC solve from one EOS call and one opacity call at face k.
279 0 : subroutine build_reconstructed_face_state_ad( &
280 : s, k, T_face, rho_face, P_face, energy_face, Cp_face, ChiRho_face, ChiT_face, grada_face, &
281 : opacity_face, scale_height_face, gradr_face, csound_face, ierr)
282 : use eos_def, only: num_eos_basic_results, i_lnPgas, i_lnE, i_grad_ad, i_gamma1, i_Cp, i_chiRho, i_chiT, i_eta, i_lnfree_e
283 : use kap_def, only: num_kap_fracs
284 :
285 : type(star_info), pointer :: s
286 : integer, intent(in) :: k
287 : type(auto_diff_real_star_order1), intent(out) :: &
288 : T_face, rho_face, P_face, energy_face, Cp_face, ChiRho_face, ChiT_face, grada_face, &
289 : opacity_face, scale_height_face, gradr_face
290 : real(dp), intent(out) :: csound_face
291 : integer, intent(out) :: ierr
292 :
293 : real(dp) :: log10_T, log10_rho, kap_zbar_face, opacity_factor_face, csound2
294 : real(dp) :: eos_res(num_eos_basic_results), d_dlnd(num_eos_basic_results), d_dlnT(num_eos_basic_results)
295 : real(dp) :: dlnT_face(auto_diff_star_num_vars), dlnd_face(auto_diff_star_num_vars)
296 : real(dp) :: kap, dlnkap_dlnd, dlnkap_dlnT
297 0 : real(dp) :: kap_fracs(num_kap_fracs), kap_xa_face(s%species)
298 : type(auto_diff_real_star_order1) :: Pgas_face, gamma1_face, mlt_Pturb_ad, alpha
299 :
300 : ierr = 0
301 : call get_face_eos_inputs( &
302 0 : s, k, T_face, rho_face, eos_res, d_dlnd, d_dlnT, ierr)
303 0 : if (ierr /= 0) return
304 0 : log10_T = log10(T_face%val)
305 0 : log10_rho = log10(rho_face%val)
306 0 : call set_face_log_partials(T_face, rho_face, dlnT_face, dlnd_face)
307 :
308 0 : call set_face_ad_from_log(eos_res(i_lnPgas), d_dlnd(i_lnPgas), d_dlnT(i_lnPgas), dlnd_face, dlnT_face, Pgas_face)
309 0 : P_face = Pgas_face + crad*pow4(T_face)/3d0
310 0 : call set_face_ad_from_log(eos_res(i_lnE), d_dlnd(i_lnE), d_dlnT(i_lnE), dlnd_face, dlnT_face, energy_face)
311 0 : call set_face_ad_from_value(eos_res(i_Cp), d_dlnd(i_Cp), d_dlnT(i_Cp), dlnd_face, dlnT_face, Cp_face)
312 0 : call set_face_ad_from_value(eos_res(i_chiRho), d_dlnd(i_chiRho), d_dlnT(i_chiRho), dlnd_face, dlnT_face, ChiRho_face)
313 0 : call set_face_ad_from_value(eos_res(i_chiT), d_dlnd(i_chiT), d_dlnT(i_chiT), dlnd_face, dlnT_face, ChiT_face)
314 0 : call set_face_ad_from_value(eos_res(i_grad_ad), d_dlnd(i_grad_ad), d_dlnT(i_grad_ad), dlnd_face, dlnT_face, grada_face)
315 0 : call set_face_ad_from_value(eos_res(i_gamma1), d_dlnd(i_gamma1), d_dlnT(i_gamma1), dlnd_face, dlnT_face, gamma1_face)
316 0 : csound2 = gamma1_face%val*P_face%val/rho_face%val
317 0 : if (is_bad_num(csound2) .or. csound2 <= 0d0) then
318 0 : ierr = -1
319 0 : if (s%report_ierr) then
320 0 : !$OMP critical (reconstructed_face_report_ierr)
321 0 : write(*,*) 'build_reconstructed_face_state_ad: bad face csound squared for k', k, csound2
322 : !$OMP end critical (reconstructed_face_report_ierr)
323 : end if
324 0 : return
325 : end if
326 0 : csound_face = sqrt(csound2)
327 :
328 0 : call get_face_composition(s, k, s% use_starting_composition_for_kap, kap_zbar_face, kap_xa_face, ierr)
329 0 : if (ierr /= 0) return
330 0 : call get_face_opacity_factor(s, k, log10_T, opacity_factor_face)
331 :
332 : call get_kap( &
333 : s, k, kap_zbar_face, kap_xa_face, log10_rho, log10_T, &
334 : eos_res(i_lnfree_e), d_dlnd(i_lnfree_e), d_dlnT(i_lnfree_e), &
335 : eos_res(i_eta), d_dlnd(i_eta), d_dlnT(i_eta), &
336 0 : kap_fracs, kap, dlnkap_dlnd, dlnkap_dlnT, ierr)
337 0 : if (ierr /= 0) then
338 0 : if (s%report_ierr) call write_face_kap_call_info( &
339 0 : s, k, T_face, rho_face, kap_zbar_face, kap_xa_face, opacity_factor_face)
340 0 : return
341 : end if
342 0 : if (is_bad_num(kap) .or. kap <= 0d0) then
343 0 : ierr = -1
344 0 : if (s%report_ierr) then
345 0 : !$OMP critical (reconstructed_face_report_ierr)
346 0 : write(*,*) 'get_reconstructed_face_eos_kap_ad: bad face opacity for k', k, kap
347 : !$OMP end critical (reconstructed_face_report_ierr)
348 0 : call write_face_kap_call_info(s, k, T_face, rho_face, kap_zbar_face, kap_xa_face, opacity_factor_face)
349 : end if
350 0 : return
351 : end if
352 :
353 0 : kap = kap*opacity_factor_face
354 0 : if (s%opacity_max > 0d0 .and. kap > s%opacity_max) then
355 0 : kap = s%opacity_max
356 0 : dlnkap_dlnd = 0d0
357 0 : dlnkap_dlnT = 0d0
358 : end if
359 0 : if (s%opacity_min > 0d0 .and. kap < s%opacity_min) then
360 0 : kap = s%opacity_min
361 0 : dlnkap_dlnd = 0d0
362 0 : dlnkap_dlnT = 0d0
363 : end if
364 0 : call set_face_ad_from_value(kap, kap*dlnkap_dlnd, kap*dlnkap_dlnT, dlnd_face, dlnT_face, opacity_face)
365 :
366 : if (s% have_mlt_vc .and. s% okay_to_set_mlt_vc .and. s% include_mlt_Pturb_in_thermodynamic_gradients &
367 0 : .and. s% mlt_Pturb_factor > 0d0 .and. k > 1) then
368 0 : mlt_Pturb_ad = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*rho_face/3d0
369 0 : alpha = mlt_Pturb_ad/(P_face*gamma1_face)
370 0 : grada_face = grada_face*(P_face + mlt_Pturb_ad)/(P_face*(1d0 + alpha))
371 : end if
372 :
373 0 : call set_scale_height_from_face_state(s, k, P_face, rho_face, scale_height_face)
374 0 : call set_gradr_from_face_state(s, k, P_face, opacity_face, T_face, gradr_face)
375 : end subroutine build_reconstructed_face_state_ad
376 :
377 :
378 : ! Returns the cached or newly built face pressure scale height as an
379 : ! auto_diff_real_star_order1 quantity from the face EOS state.
380 0 : subroutine get_reconstructed_scale_height_ad(s, k, scale_height_face, ierr)
381 : type(star_info), pointer :: s
382 : integer, intent(in) :: k
383 : type(auto_diff_real_star_order1), intent(out) :: scale_height_face
384 : integer, intent(out) :: ierr
385 :
386 : ierr = 0
387 0 : call ensure_reconstructed_face_state_ad(s, k, ierr)
388 0 : if (ierr /= 0) return
389 0 : scale_height_face = s%reconstructed_scale_height_face_ad(k)
390 : end subroutine get_reconstructed_scale_height_ad
391 :
392 :
393 0 : subroutine set_scale_height_from_face_state(s, k, P_face, rho_face, scale_height_face)
394 : use auto_diff_support, only: wrap_r_00
395 :
396 : type(star_info), pointer :: s
397 : integer, intent(in) :: k
398 : type(auto_diff_real_star_order1), intent(in) :: P_face, rho_face
399 : type(auto_diff_real_star_order1), intent(out) :: scale_height_face
400 :
401 : real(dp) :: G
402 : type(auto_diff_real_star_order1) :: grav, scale_height2
403 :
404 0 : G = s%cgrav(k)
405 0 : grav = G*s%m_grav(k)/pow2(wrap_r_00(s,k))
406 0 : scale_height_face = P_face/(grav*rho_face)
407 0 : if (s%alt_scale_height_flag) then
408 0 : scale_height2 = sqrt(P_face/G)/rho_face
409 0 : if (scale_height2 < scale_height_face) scale_height_face = scale_height2
410 : end if
411 0 : end subroutine set_scale_height_from_face_state
412 :
413 :
414 0 : subroutine set_gradr_from_face_state(s, k, P_face, opacity_face, T_face, gradr_face)
415 : use auto_diff_support, only: wrap_L_00
416 :
417 : type(star_info), pointer :: s
418 : integer, intent(in) :: k
419 : type(auto_diff_real_star_order1), intent(in) :: P_face, opacity_face, T_face
420 : type(auto_diff_real_star_order1), intent(out) :: gradr_face
421 :
422 : real(dp) :: L_theta
423 : type(auto_diff_real_star_order1) :: L_face, Pr_face
424 :
425 0 : if (s%include_mlt_in_velocity_time_centering) then
426 : if (s%using_velocity_time_centering .and. &
427 0 : s%include_L_in_velocity_time_centering .and. &
428 : s%lnT(k) <= s%max_logT_for_include_P_and_L_in_velocity_time_centering*ln10) then
429 0 : L_theta = s%L_theta_for_velocity_time_centering
430 : else
431 0 : L_theta = 1d0
432 : end if
433 0 : L_face = L_theta*wrap_L_00(s, k) + (1d0 - L_theta)*s%L_start(k)
434 : else
435 0 : L_face = wrap_L_00(s, k)
436 : end if
437 :
438 0 : Pr_face = crad*pow4(T_face)/3d0
439 0 : gradr_face = P_face*opacity_face*L_face/(4d0*pi4*clight*s%m_grav(k)*s%cgrav(k)*Pr_face)
440 0 : end subroutine set_gradr_from_face_state
441 :
442 : ! Precomputes dlnT_face and dlnd_face for converting scalar d/dlnT and
443 : ! d/dlnd microphysics partials into star-order1 autodiff derivatives.
444 0 : subroutine set_face_log_partials(T_face, rho_face, dlnT_face, dlnd_face)
445 : type(auto_diff_real_star_order1), intent(in) :: T_face, rho_face
446 : real(dp), intent(out) :: dlnT_face(auto_diff_star_num_vars), dlnd_face(auto_diff_star_num_vars)
447 :
448 0 : dlnT_face = T_face%d1Array/T_face%val
449 0 : dlnd_face = rho_face%d1Array/rho_face%val
450 0 : end subroutine set_face_log_partials
451 :
452 :
453 : ! Converts a scalar value with d/dlnd and d/dlnT partials into an
454 : ! auto_diff_real_star_order1 quantity using the face chain rule.
455 0 : subroutine set_face_ad_from_value(value, dvalue_dlnd, dvalue_dlnT, dlnd_face, dlnT_face, quantity_ad)
456 : real(dp), intent(in) :: value, dvalue_dlnd, dvalue_dlnT
457 : real(dp), intent(in) :: dlnd_face(auto_diff_star_num_vars), dlnT_face(auto_diff_star_num_vars)
458 : type(auto_diff_real_star_order1), intent(out) :: quantity_ad
459 :
460 0 : quantity_ad = 0d0
461 0 : quantity_ad%val = value
462 0 : quantity_ad%d1Array = dvalue_dlnd*dlnd_face + dvalue_dlnT*dlnT_face
463 0 : end subroutine set_face_ad_from_value
464 :
465 :
466 : ! Same as set_face_ad_from_value, but for a quantity returned in
467 : ! logarithmic form by the microphysics routine.
468 0 : subroutine set_face_ad_from_log(log_value, dlog_dlnd, dlog_dlnT, dlnd_face, dlnT_face, quantity_ad)
469 : real(dp), intent(in) :: log_value, dlog_dlnd, dlog_dlnT
470 : real(dp), intent(in) :: dlnd_face(auto_diff_star_num_vars), dlnT_face(auto_diff_star_num_vars)
471 : type(auto_diff_real_star_order1), intent(out) :: quantity_ad
472 : real(dp) :: value
473 :
474 0 : value = exp(log_value)
475 0 : call set_face_ad_from_value(value, value*dlog_dlnd, value*dlog_dlnT, dlnd_face, dlnT_face, quantity_ad)
476 0 : end subroutine set_face_ad_from_log
477 :
478 :
479 : ! Writes the recomputed face EOS inputs that were passed to get_eos.
480 0 : subroutine write_face_eos_call_info(s, k, T_face, rho_face, zbar_face, xa_face)
481 : type(star_info), pointer :: s
482 : integer, intent(in) :: k
483 : type(auto_diff_real_star_order1), intent(in) :: T_face, rho_face
484 : real(dp), intent(in) :: zbar_face
485 : real(dp), intent(in) :: xa_face(:)
486 :
487 : integer :: j
488 : include 'formats'
489 :
490 0 : !$OMP critical (reconstructed_face_eos_call_info)
491 0 : write(*,'(A)')
492 0 : write(*,*) 'face EOS input info for k', k
493 0 : write(*,1) 'T_face', T_face%val
494 0 : write(*,1) 'rho_face', rho_face%val
495 0 : write(*,1) 'log10_T_face', log10(T_face%val)
496 0 : write(*,1) 'log10_rho_face', log10(rho_face%val)
497 0 : write(*,1) 'zbar_face', zbar_face
498 0 : write(*,1) 'sum(xa_face)', sum(xa_face)
499 0 : do j = 1, s%species
500 0 : write(*,2) 'xa_face ' // trim(s%nameofequ(j+s%nvar_hydro)), j, xa_face(j)
501 : end do
502 : !$OMP end critical (reconstructed_face_eos_call_info)
503 0 : end subroutine write_face_eos_call_info
504 :
505 :
506 : ! Writes the recomputed face opacity inputs that were passed to get_kap.
507 0 : subroutine write_face_kap_call_info(s, k, T_face, rho_face, zbar_face, xa_face, opacity_factor_face)
508 : type(star_info), pointer :: s
509 : integer, intent(in) :: k
510 : type(auto_diff_real_star_order1), intent(in) :: T_face, rho_face
511 : real(dp), intent(in) :: zbar_face
512 : real(dp), intent(in) :: xa_face(:)
513 : real(dp), intent(in) :: opacity_factor_face
514 :
515 : integer :: j
516 : include 'formats'
517 :
518 0 : !$OMP critical (reconstructed_face_kap_call_info)
519 0 : write(*,'(A)')
520 0 : write(*,*) 'face opacity input info for k', k
521 0 : write(*,1) 'T_face', T_face%val
522 0 : write(*,1) 'rho_face', rho_face%val
523 0 : write(*,1) 'log10_T_face', log10(T_face%val)
524 0 : write(*,1) 'log10_rho_face', log10(rho_face%val)
525 0 : write(*,1) 'zbar_face', zbar_face
526 0 : write(*,1) 'opacity_factor_face', opacity_factor_face
527 0 : write(*,1) 'sum(xa_face)', sum(xa_face)
528 0 : do j = 1, s%species
529 0 : write(*,2) 'xa_face ' // trim(s%nameofequ(j+s%nvar_hydro)), j, xa_face(j)
530 : end do
531 : !$OMP end critical (reconstructed_face_kap_call_info)
532 0 : end subroutine write_face_kap_call_info
533 :
534 : end module reconstructed_face_support
|
